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

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<xp && xp<rol && -cow<yp && yp<cow;

            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); %落在多边形内部的中心点填充不规则多边形

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