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

好久不见啊，今天时绘图复刻第三期，这期画的比较难应该文章也不会太短。。。
前段时间发现公众号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;
1
2
3
4
5
6
7
8
9

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

当然也可以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'};
1
2
3
4
5
6
7
8
9
10

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

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

3 树状图绘制

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

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

设置为黑色并加粗：

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

但是注意到此时两个图不在同一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 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

调用该函数复制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)
1
2
3
4
5
6
7
8
9
10

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

% 绘制上方树状图
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];
1
2
3
4
5
6
7
8
9
10
11
12

4 热图绘制

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

% 绘制中央热图
Data=Data(order1,:);
Data=Data(:,order2);
imagesc(axMain,Data)
axMain.XTickLabel=colName(order2);
axMain.YTickLabel=rowName(order1);
1
2
3
4
5
6

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);
1
2
3
4
5
6
7

6 调整colorbar范围和配色

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

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

当然自己找一些数据进行插值也可以，这里提供了六种配色，可自行添加，原理就是自己找一些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))])
1
2
3
4
5
6
7
8
9
10
11
12
13
14

调整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 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

---