点在多边形内_玩儿转多组学的分析组合拳 | 泰森多边形制图全攻略

泰森多边形制图全攻略

亲爱的小伙伴们，在前两期的分享中，我们给大家介绍了一个功能强大的应用——

泰森多边形

——一款可以将多维分析一网打尽的可视化工具图形。

往期回顾

●泰森多边形 | 只要是多组学，您都需要这一款

●玩儿转扩增子测序 | 泰森多边形揭示为什么就你胖的奥秘

今天，凌恩生信团队就为大家带来这款分析作图的流程全攻略，一起看过来吧！

Voronoi图又叫泰森多边形或Dirichlet图,它是一组连续多边形组成，多边形的边界是由连接的垂直平分线组成。M个在平面上有区别的点，按照最近邻原则划分平面，每一个点与它最近邻的区域关联。

Voronoi树图通过使用加权重心Voronoi图递归划分凸多边形来可视化分层数据。多边形区域被节点的相对权重比例分割。 

Voronoi treemap

在之前的文献分享中，我们已经为大家介绍过Voronoi treemap在代谢通路和微生物多样性上的应用：

今天，我们就带着大家一起使用matlab来绘制如此高大上的voronoi treemap！

首先是输入数据格式，如下：

主程序

正五边形分割矩形

inf=importdata('data.txt'); %读取数据

data=inf.data;%获取数据中的数值信息

[t_row t_col]=size(data);%计算矩阵行数和列数

resolu=ceil(sqrt(t_row));   %nxn切割正方形

resolu_x=[];    resolu_y=[];    resolu_mat=[];

rc=1/(2*resolu);%正五边形半径

dy=2*rc;dx=rc*sqrt(3);

A=pi/3*[1:7];   rol=1;  cow=1;  num=0;  resolu_record=0;

for yk=-2:dy:2

    yfun=inline(['sqrt(3)*x/3+',num2str(yk)]);

    for xk=-2:dx:2

        xp=xk;      yp=yfun(xp);    resolu_record=resolu_record+1;

        resolu_x(resolu_record)=xp; resolu_mat(resolu_record,1)=xp;

        resolu_y(resolu_record)=yp; resolu_mat(resolu_record,2)=yp;

        if -rol

            plot([xp+1i*yp]+rc*exp(1i*A)*2/sqrt(3),'k','linewidth',2);

            hold on;

            a=real([xp+1i*yp]+rc*exp(1i*A)*2/sqrt(3));

            b=imag([xp+1i*yp]+rc*exp(1i*A)*2/sqrt(3));

            fill(a(1:6),b(1:6),'red');

            plot(xp,yp,'b.');

            %text(xp,yp,num2str(num))

        end

    end

end

生成kegg通路维诺图

site_rand=rand(500,layer1_num1,2);

layer1_x=site_rand(i,:,1);    layer1_y=site_rand(i,:,2);

Options.plot=1; %设置1表示画出维诺图

P = polytope(v); %生成边界

Options.pbound=P;

Pn=mpt_voronoi([layer1_x;layer1_y]',Options); %生成voronoi图信息

for j=1:length(layer1_x)

        [V,R] = extreme(Pn(j)); %这里的V就是第i个多边形的顶点序列

        sort_pos1 = convhull(V(:,1), V(:,2)); %获取多边形顶点

        [pX,pY] = poly2cw(V(sort_pos1,1),V(sort_pos1,2));  % make clockwise  顺时针转动

        area_info(i,j) = 100*polyarea(pX,pY);% 计算多边形面积

end

part_area_info=sort(area_info(i,:));

area_res(i)=sum(abs(part_area_info-value1'));   

对KEGG通路novonoi图填充正五边形基因信息

%绘图

for i=1:layer1_num1

    value_site=value1_site(i);

    [V,R] = extreme(Pn(final_area_site(i))); %这里的V就是第i大区域内多边形的顶点序列

    sort_pos1 = convhull(V(:,1), V(:,2));  

    [pX,pY] = poly2cw(V(sort_pos1,1),V(sort_pos1,2));    % 顺时针转动

    edge=[pX,pY];     patch(pX,pY,'white','edgecolor','black','FaceAlpha',.5,'linewidth',2); hold on;%voronoi图形填充颜色

    %第二层级基因信息

    layer2=tabulate(test(strcmp(test(:,2),layer1{value_site,1}),1));

    layer2_data=data(strcmp(test(:,2),layer1{value_site,1}));

    [layer2_num1 layer2_num2]=size(layer2);

    record=0;

    value2=sort(cell2mat(layer2(:,3)));

    [in,on]=inpolygon(resolu_x,resolu_y,pX,pY);%判断正五边形的中心点是否落在不规则多边形内

    site_select=resolu_mat(in&~on,:);      %选取落在多边形内部的中心点

    [aaa bbb]=size(site_select);

    layer2_x=site_select(:,1);

    layer2_y=site_select(:,2);

    Opart.plot=0; %设置1表示画出维诺图

    partP = polytope(edge); %生成边界

    Opart.pbound=partP;

part_Pn=mpt_voronoi([layer2_x,layer2_y],Opart); %落在多边形内部的中心点填充不规则多边形

获得多边形分割信息并填充颜色和绘图

最终绘制的图片

灰色表示表达没有变化基因，蓝色表示下调差异表达基因，红色表示上调差异表达基因，颜色深浅表示差异的倍数(Fold change)。 从图中可以明显的看出差异基因主要富集在哪些代谢通路上。

凌恩生物成立于2014年，专注组学技术在科研领域的应用与研究。公司成立以来，技术团队参与的项目成果成功发表在《Nature》《Cell》《PNAS》等国际顶端学术期刊。    

秉承“以客户需求为本，为客户创造价值”的服务宗旨；以高品质、高效率的技术服务，用心打造凌恩品牌，助力您的成功！