当前位置：首页 > news >正文

MATLAB | 绘图复刻(三) | 分层聚类分析图：树状图+热图

news 来源：原创 2024/5/4 5:42:58

好久不见啊，今天时绘图复刻第三期，这期画的比较难应该文章也不会太短。。。
前段时间发现公众号SCIPainter发布了一期名为《如何对基因和蛋白质的表达丰度进行相关性分析》，其中有一幅图很好看：

于是我也复刻了一下。由于原文章没有提供数据，我这里随便捏造了点数据进行的绘图，以下是绘制效果：

大家光看图就能看出能展示哪些信息，明显行列相似度大的被调整的挨在一起，且能看出明显颜色分界处也正好是树状图主要分支的分界。

对了，重点需要安装Statistics and Machine Learning Toolbox即统计与机器学习工具箱！！！

0 注

本来发现Bioinformatics Toolbox工具箱就有分层聚类分析的绘制函数clustergram，但是画出来属实太丑了而且配色啊树状图啊统统很难调整：

嗯。。。自己写呗，直接开肝。

1 数据

这里假设X,Y数据为相同数据，计算出的相关性矩阵就是对称矩阵：

X=randn(20,20)+[(linspace(-1,2.5,20)').*ones(1,8),(linspace(.5,-.7,20)').*ones(1,5),(linspace(.9,-.2,20)').*ones(1,7)];
Y=X;
% 计算数据相关性系数
Data=corr(X,Y);

% 输入行列名称
colName={'FREM2','ALDH9A1','RBL1','AP2A2','HNRNPK','ATP1A1','ARPC3','SMG5','RPS27A',...
          'RAB8A','SPARC','DDX3X','EEF1D','EEF1B2','RPS11','RPL13','RPL34','GCN1','FGG','CCT3'};
rowName=colName;

如果是对称矩阵最终效果：

当然也可以X,Y数据为不同数据：

rng(1)
Y=randn(20,20)+[(linspace(-1,2.5,20)').*ones(1,8),(linspace(.5,-.7,20)').*ones(1,5),(linspace(.9,-.2,20)').*ones(1,7)];
X=randn(20,15)+[(linspace(-1,2.5,20)').*ones(1,4),(linspace(.5,-.7,20)').*ones(1,5),(linspace(.9,-.2,20)').*ones(1,6)];
% 计算数据相关性系数
Data=corr(X,Y);

% 输入行列名称
colName={'FREM2','ALDH9A1','RBL1','AP2A2','HNRNPK','ATP1A1','ARPC3','SMG5','RPS27A',...
          'RAB8A','SPARC','DDX3X','EEF1D','EEF1B2','RPS11','RPL13','RPL34','GCN1','FGG','CCT3'};
rowName={'A1','A2','A3','A4','A5','A6','A7','A8','A9','A10','A11','A12','A13','A14','A15'};

如果非对称矩阵最终效果：

2 坐标区域构建

创建坐标区域并调整各个坐标区域的位置Position，中间的热图，两个树状图以及colorbar都分别创建一个坐标区域：

其中比较少见的属性YAxisLocation是用来设置Y轴是在0点处还是坐标区域左侧或者右侧，本期需要设置为右侧，uistack用来调整图形对象层级关系，需要调整一部分坐标区域谁压着谁。

% =========================================================================
% 获取数据矩阵大小
[rows,cols]=size(Data);

fig=figure('Position',[500,200,800,750],'Name','slandarer');


% 坐标区域创建 =============================================================
% 热图坐标区域
placeMat=zeros(7,8);placeMat(2:7,2:7)=1;
axMain=subplot(7,8,find(placeMat'));
axMain.XLim=[1,cols]+[-.5,.5];
axMain.YLim=[1,rows]+[-.5,.5];
axMain.YAxisLocation='right';
axMain.YDir='reverse';
axMain.XTick=1:cols;
axMain.YTick=1:rows; 
hold on
% 树状图坐标区域
axTree1=subplot(7,8,(1:6).*8+1);
axTree1.Position(3)=axTree1.Position(3)+axTree1.Position(1)*4/5;
axTree1.Position(1)=axTree1.Position(1)/5;
axTree1.Position(3)=(axMain.Position(1)-axTree1.Position(1)-axTree1.Position(3))*4/5+axTree1.Position(3);
hold on
axTree2=subplot(7,8,2:7);
axTree2.Position(2)=axMain.Position(2)+axMain.Position(4)+(axTree2.Position(2)-axMain.Position(2)-axMain.Position(4))/5;
axTree2.Position(4)=axTree2.Position(4)+(1-axTree2.Position(2)-axTree2.Position(4))*4/5;
hold on
% colorbar坐标区域
axBar=subplot(7,8,(1:7).*8);
axBar.Position(1)=1/2;
axBar.Position(3)=.92/2;
axBar.Position(2)=axMain.Position(2)+axMain.Position(4)/2;
axBar.Position(4)=axMain.Position(2)+axMain.Position(4)-axBar.Position(2);
axBar.Color='none';
axBar.XColor='none';
axBar.YColor='none';
uistack(axBar,'bottom')
CM=colorbar;
hold on

3 树状图绘制

这里使用linkage获取集聚分层聚类树，再使用dendrogram函数将其绘制，dendrogram函数的Orientation属性可以调整树的方向：

tree1=linkage(Data,'average');
[treeHdl1,~,order1]=dendrogram(tree1,0,'Orientation','left');

设置为黑色并加粗：

% 设置树状图颜色
set(treeHdl1,'Color',[0,0,0]);
set(treeHdl1,'LineWidth',1);

但是注意到此时两个图不在同一figure:

此时我们就想到之前写过一篇《MATLAB | 如何复制figure图窗任意axes的全部信息？》中给了一个复制axes的函数，将其进行微调：

function axbag=copyAxes(fig,k,newAx)
% @author : slandarer
% 公众号  : slandarer随笔
% 知乎    : slandarer
% 
% 此段代码解析详见公众号 slandarer随笔 文章：
%《MATLAB | 如何复制figure图窗任意axes的全部信息？》
% https://mp.weixin.qq.com/s/3i8C78pv6Ok1cmEZYPMyWg
classList(length(fig.Children))=true;
for n=1:length(fig.Children)
    classList(n)=isa(fig.Children(n),'matlab.graphics.axis.Axes');
end
isaaxes=find(classList);
oriAx=fig.Children(isaaxes(end-k+1));
if isaaxes(end-k+1)-1<1||isa(fig.Children(isaaxes(end-k+1)-1),'matlab.graphics.axis.Axes')
    oriLgd=[];
else
    oriLgd=fig.Children(isaaxes(end-k+1)-1);
end
axbag=copyobj([oriAx,oriLgd],newAx.Parent);
axbag(1).Position=newAx.Position;
delete(newAx)
end

调用该函数复制axes并修饰：

tempFig=treeHdl1(1).Parent.Parent;
% 坐标区域二次修饰
axTree1=copyAxes(tempFig,1,axTree1);
axTree1.XColor='none';
axTree1.YColor='none';
axTree1.YDir='reverse';
axTree1.XTick=[];
axTree1.YTick=[];
axTree1.YLim=[1,rows]+[-.5,.5];
delete(tempFig)

另一个树状图代码几乎一模一样:

% 绘制上方树状图
tree2=linkage(Data.','average');
[treeHdl2,~,order2]=dendrogram(tree2,0,'Orientation','top');
set(treeHdl2,'Color',[0,0,0]);
set(treeHdl2,'LineWidth',1);
tempFig=treeHdl2(1).Parent.Parent;
axTree2=copyAxes(tempFig,1,axTree2);
axTree2.XColor='none';
axTree2.YColor='none';
axTree2.XTick=[];
axTree2.YTick=[];
axTree2.XLim=[1,cols]+[-.5,.5];

4 热图绘制

需要根据linkage及dendrogram的聚类结果调整相关系数矩阵行列顺序以及X,Y轴标签顺序：

% 绘制中央热图
Data=Data(order1,:);
Data=Data(:,order2);
imagesc(axMain,Data)
axMain.XTickLabel=colName(order2);
axMain.YTickLabel=rowName(order1);

5 绘制白线

% 绘制白线
LineX=repmat([[1,cols]+[-.5,.5],nan],[rows+1,1]).';
LineY=repmat((.5:1:(rows+.5)).',[1,3]).';
plot(axMain,LineX(:),LineY(:),'Color','w','LineWidth',1);
LineY=repmat([[1,rows]+[-.5,.5],nan],[cols+1,1]).';
LineX=repmat((.5:1:(cols+.5)).',[1,3]).';
plot(axMain,LineX(:),LineY(:),'Color','w','LineWidth',1);

6 调整colorbar范围和配色

范围使用clim调整，配色的话MATLAB自带的配色均可：

colormap(pink)
clim([min(min(Data)),max(max(Data))])

当然自己找一些数据进行插值也可以，这里提供了六种配色，可自行添加，原理就是自己找一些nx3大小的颜色列表，然后插值到上百行，颜色就渐变了：

% 调整colorbar
baseCM={[189, 53, 70;255,255,255; 97, 97, 97]./255,...
    [13,135,169;255,255,255;239,139,14]./255,...
    [ 28,127,119;255,255,255;204,157, 80]./255,...
    [130,130,255;255,255,255;255,133,133]./255,...
    [209,58,78;253,203,121;254,254,189;198,230,156;63,150,181]./255,...
    [243,166, 72;255,255,255;133,121,176]./255};
Cmap=baseCM{2};
Ci=1:size(Cmap,1);Cq=linspace(1,size(Cmap,1),200);% 插值到200行
Cmap=[interp1(Ci,Cmap(:,1),Cq,'linear')',...
     interp1(Ci,Cmap(:,2),Cq,'linear')',...
     interp1(Ci,Cmap(:,3),Cq,'linear')'];
colormap(Cmap)
clim([min(min(Data)),max(max(Data))])

调整Cmap=baseCM{2};中数字可以调整配色：

C1

C2

C3

C4

C5

C6

完整代码

% @author : slandarer
% gzh  : slandarer随笔
% zh    : slandarer

% 随机生成数据
rng(1)
X=randn(20,20)+[(linspace(-1,2.5,20)').*ones(1,8),(linspace(.5,-.7,20)').*ones(1,5),(linspace(.9,-.2,20)').*ones(1,7)];
Y=X;
% X=randn(20,15)+[(linspace(-1,2.5,20)').*ones(1,4),(linspace(.5,-.7,20)').*ones(1,5),(linspace(.9,-.2,20)').*ones(1,6)];
% 计算数据相关性系数
Data=corr(X,Y);

% 输入行列名称
colName={'FREM2','ALDH9A1','RBL1','AP2A2','HNRNPK','ATP1A1','ARPC3','SMG5','RPS27A',...
          'RAB8A','SPARC','DDX3X','EEF1D','EEF1B2','RPS11','RPL13','RPL34','GCN1','FGG','CCT3'};
rowName=colName;
% rowName={'A1','A2','A3','A4','A5','A6','A7','A8','A9','A10','A11','A12','A13','A14','A15'};


% =========================================================================
% 获取数据矩阵大小
[rows,cols]=size(Data);
fig=figure('Position',[500,200,800,750],'Name','slandarer','Color',[1,1,1]);

% 坐标区域创建 =============================================================
% 热图坐标区域
placeMat=zeros(7,8);placeMat(2:7,2:7)=1;
axMain=subplot(7,8,find(placeMat'));
axMain.XLim=[1,cols]+[-.5,.5];
axMain.YLim=[1,rows]+[-.5,.5];
axMain.YAxisLocation='right';
axMain.YDir='reverse';
axMain.XTick=1:cols;
axMain.YTick=1:rows;
hold on
% 树状图坐标区域
axTree1=subplot(7,8,(1:6).*8+1);
axTree1.Position(3)=axTree1.Position(3)+axTree1.Position(1)*4/5;
axTree1.Position(1)=axTree1.Position(1)/5;
axTree1.Position(3)=(axMain.Position(1)-axTree1.Position(1)-axTree1.Position(3))*4/5+axTree1.Position(3);
hold on
axTree2=subplot(7,8,2:7);
axTree2.Position(2)=axMain.Position(2)+axMain.Position(4)+(axTree2.Position(2)-axMain.Position(2)-axMain.Position(4))/5;
axTree2.Position(4)=axTree2.Position(4)+(1-axTree2.Position(2)-axTree2.Position(4))*4/5;
hold on
% colorbar坐标区域
axBar=subplot(7,8,(1:7).*8);
axBar.Position(1)=1/2;
axBar.Position(3)=.92/2;
axBar.Position(2)=axMain.Position(2)+axMain.Position(4)/2;
axBar.Position(4)=axMain.Position(2)+axMain.Position(4)-axBar.Position(2);
axBar.Color='none';
axBar.XColor='none';
axBar.YColor='none';
uistack(axBar,'bottom')
CM=colorbar;
hold on
% 图像绘制 =================================================================
% 绘制左侧树状图
tree1=linkage(Data,'average');
[treeHdl1,~,order1]=dendrogram(tree1,0,'Orientation','left');
% 设置树状图颜色
set(treeHdl1,'Color',[0,0,0]);
set(treeHdl1,'LineWidth',1);
tempFig=treeHdl1(1).Parent.Parent;
% 坐标区域二次修饰
axTree1=copyAxes(tempFig,1,axTree1);
axTree1.XColor='none';
axTree1.YColor='none';
axTree1.YDir='reverse';
axTree1.XTick=[];
axTree1.YTick=[];
axTree1.YLim=[1,rows]+[-.5,.5];
delete(tempFig)

% 绘制上方树状图
tree2=linkage(Data.','average');
[treeHdl2,~,order2]=dendrogram(tree2,0,'Orientation','top');
set(treeHdl2,'Color',[0,0,0]);
set(treeHdl2,'LineWidth',1);
tempFig=treeHdl2(1).Parent.Parent;
axTree2=copyAxes(tempFig,1,axTree2);
axTree2.XColor='none';
axTree2.YColor='none';
axTree2.XTick=[];
axTree2.YTick=[];
axTree2.XLim=[1,cols]+[-.5,.5];
delete(tempFig)
% 绘制中央热图
Data=Data(order1,:);
Data=Data(:,order2);
imagesc(axMain,Data)
axMain.XTickLabel=colName(order2);
axMain.YTickLabel=rowName(order1);
% 绘制白线
LineX=repmat([[1,cols]+[-.5,.5],nan],[rows+1,1]).';
LineY=repmat((.5:1:(rows+.5)).',[1,3]).';
plot(axMain,LineX(:),LineY(:),'Color','w','LineWidth',1);
LineY=repmat([[1,rows]+[-.5,.5],nan],[cols+1,1]).';
LineX=repmat((.5:1:(cols+.5)).',[1,3]).';
plot(axMain,LineX(:),LineY(:),'Color','w','LineWidth',1);
% 调整colorbar
baseCM={[189, 53, 70;255,255,255; 97, 97, 97]./255,...
    [13,135,169;255,255,255;239,139,14]./255,...
    [ 28,127,119;255,255,255;204,157, 80]./255,...
    [130,130,255;255,255,255;255,133,133]./255,...
    [209,58,78;253,203,121;254,254,189;198,230,156;63,150,181]./255,...
    [243,166, 72;255,255,255;133,121,176]./255};
Cmap=baseCM{5};
Ci=1:size(Cmap,1);Cq=linspace(1,size(Cmap,1),200);% 插值到200行
Cmap=[interp1(Ci,Cmap(:,1),Cq,'linear')',...
     interp1(Ci,Cmap(:,2),Cq,'linear')',...
     interp1(Ci,Cmap(:,3),Cq,'linear')'];
colormap(Cmap)
clim([min(min(Data)),max(max(Data))])



% 工具函数 =================================================================
% 复制坐标区域全部图形对象
function axbag=copyAxes(fig,k,newAx)
% @author : slandarer
% gzh  : slandarer随笔
% zh    : slandarer
% 
% 此段代码解析详见公众号 slandarer随笔 文章：
%《MATLAB | 如何复制figure图窗任意axes的全部信息？》
% https://mp.weixin.qq.com/s/3i8C78pv6Ok1cmEZYPMyWg
classList(length(fig.Children))=true;
for n=1:length(fig.Children)
    classList(n)=isa(fig.Children(n),'matlab.graphics.axis.Axes');
end
isaaxes=find(classList);
oriAx=fig.Children(isaaxes(end-k+1));
if isaaxes(end-k+1)-1<1||isa(fig.Children(isaaxes(end-k+1)-1),'matlab.graphics.axis.Axes')
    oriLgd=[];
else
    oriLgd=fig.Children(isaaxes(end-k+1)-1);
end
axbag=copyobj([oriAx,oriLgd],newAx.Parent);
axbag(1).Position=newAx.Position;
delete(newAx)
end