绘制效果代码说明Mantel test 相关性热图的构造复刻自 R 语言代码:Houyun Huang(2021). linkET: Everything is Linkable. R package version 0.0.3.我为我开发的 MATLAB 绘图工具 SHeatmap 添加了新的辅助函数 SMantelLink 使得 MATLAB 也能具备绘制这样图像的能力。因此想要绘制首先需要下载工具函数包可在以下仓库获取【gitee】https://gitee.com/slandarer/matlab-special-heatmap【github】https://github.com/slandarer/MATLAB-special-heatmap【fileexchange】https://www.mathworks.com/matlabcentral/fileexchange/125520-special-heatmapMantel test在给出绘图教程前先给出 Mantel test 的原理和代码实现这样大家可以在其他地方使用 Mantel test。Mantel test 被用来检验两个距离矩阵的相关性基本实现是将两个矩阵下三角部分(不含对角线)拉成向量计算相关系数 r, 随机同时交换距离矩阵2的行列之后将矩阵2下三角部分再次拉成长条对两个列向量进行相关性检验记录r p r e m r_{prem}rprem​值看有多少次因为交换行列导致r p r e m r_{prem}rprem​值对比 r 值上升计算上升次数/总次数作为 p 值同时为了和主流计算保持一致(如 R 语言的vegan包)采用加 1 校正r 值越高说明相关性越强p 值越低说明越显著。以下是 Mantel test 实现代码mantel工具包还有一个拿存储空间换运行速度的mantelFast函数运行速度是下面这个函数的10倍可以工具包里瞅瞅。function[rho,pval,rperm]mantel(dist1,dist2,varargin)% MANTEL - Mantel test for matrix correlation% rho mantel(dist1, dist2) performs Mantel test with 999 permutations% using Pearson correlation.%% rho mantel(dist1, dist2, nperm) specifies number of permutations.%% rho mantel(dist1, dist2, method) specifies correlation method% (Pearson or Spearman, default Pearson).%% rho mantel(dist1, dist2, nperm, method) sets both.%% [rho, pval] mantel(___) also returns the p-value (one-sided).%% [rho, pval, rperm] mantel(___) returns the permuted r values.% References% [1] Mantel N. The detection of disease clustering and a generalized regression approach.% Cancer Res. 1967 Feb;27(2):209-20. PMID: 6018555.% [2] Borcard, D. Legendre, P. (2012) Is the Mantel correlogram powerful enough to be% useful in ecological analysis? A simulation study. Ecology 93: 1473-1481.% [3] Legendre, P. and Legendre, L. (2012) Numerical Ecology. 3rd English Edition. Elsevier.% Parse input argumentsswitchnargincase2nperm999;methodPearson;case3ifisnumeric(varargin{1})npermvarargin{1};methodPearson;elsenperm999;methodvarargin{1};endcase4npermvarargin{1};methodvarargin{2};end% Extract lower-triangular vectorsNsize(dist1,1);indtril(true(N),-1);V1dist1(ind);V2dist2(ind);% Observed correlationrhocorr(V1,V2,Type,method);% Permutation testrpermzeros(nperm,1);fori1:nperm permrandperm(N);dist2_permdist2(perm,perm);V2_permdist2_perm(ind);rperm(i)corr(V1,V2_perm,Type,method);end% p-value (one-sided, positive correlation)pval(sum(rpermrho)1)/(nperm1);end以下说明该代码计算结果与 R 语言相同MATLAB 计算rng(1)X1rand([100,20]);X2rand([100,5]);writematrix(X1,X1.csv);writematrix(X2,X2.csv);D1squareform(pdist(X1));D2squareform(pdist(X2));[rho,pval]mantel(D1,D2,9999,Pearson)计算得到rho 0.0438, pval 0.0827R 语言计算dataX1-read.csv(X1.csv,headerFALSE,row.namesNULL)dataX2-read.csv(X2.csv,headerFALSE,row.namesNULL)D1-dist(dataX1)D2-dist(dataX2)result-mantel(D1,D2,methodpearson,permutations9999)print(result)计算得到Mantel statistic r: 0.04375, Significance: 0.0828 , 计算结果相似。当然 Mantel test 具有随机性建议加大随机置换次数并设置随机数种子。mantelFast 函数就是把循环部分完全换成 MATLAB 矩阵运算但是在数据量超过 2e8 时容易会内存不足% Permutation test[~,perm]sort(rand(N,nperm),1);[indR,indC]find(tril(true(N),-1));V2dist2(sub2ind([N,N],perm(indR,:),perm(indC,:)));rpermcorr(V1,V2,Type,method);Mantel test 相关性热图数据准备数据来源Väre, H., Ohtonen, R. and Oksanen, J. (1995) Effects of reindeer grazing on understorey vegetation in dry Pinus sylvestris forests. Journal of Vegetation Science 6, 523–530.%% Load data (加载数据)load(lichenData.mat)% Load pre-saved data package (加载预存的数据包)Data1varechem.Variables;% Environmental matrix (环境因子矩阵)Data2varespec.Variables;% Species composition matrix (物种组成矩阵)labelsvarechem.Properties.VariableNames;% Environmental variable names (环境变量名称)% Define species groups: 44 columns into 4 groups (将44个物种列分为4组)groupName{Spec01,Spec02,Spec03,Spec04};groupzeros(1,size(Data2,2));group(1:7)1;% Group 1: columns 1-7 (第1组列1-7)group(8:18)2;% Group 2: columns 8-18 (第2组列8-18)group(19:37)3;% Group 3: columns 19-37 (第3组列19-37)group(38:44)4;% Group 4: columns 38-44 (第4组列38-44)热图绘制axgca;%% Draw heatmap[rho,pval]corr(Data1);objHMSHeatmap(ax,rho,Format,sq);objHM.draw();% Display significance stars: p 0.05 *, p 0.01 **, p 0.001 ***objHM.setText()objHM.showStars(pval,Levels,[0.05,0.01,0.001],CorrLabel,off)objHM.setVarName(labels)% Set row/column labels (设置行列标签)objHM.setType(tril0);% Show only lower triangle without diagonal (仅显示下三角矩阵不含对角线)% Adjust label positions (调整标签位置)objHM.setRowLabelLocation(left)objHM.setColLabelLocation(bottom)% Show all labels, including those that were previously hidden. (显示所有标签包括被隐藏的标签)objHM.setRowLabel(Visible,on)objHM.setColLabel(Visible,on)delete(objHM.Colorbar)Mantel 链接绘制% Create Mantel link object with env data, species data and groups (创建Mantel链接对象传入环境数据、物种数据及分组信息)objMLSMantelLink(ax,Data1,Data2,Group,group);objML.GroupNamegroupName;% Set group names (设置组名)objML.LegendLocationwest;% Place legend on the left (图例置于左侧)objML.Layouttriu;% Links placed in upper triangle (链接采用上三角布局)objML.draw()% Render the links (渲染链接图)在draw()之前很多参数都可以设置比如为了每次结果一致建议提前设置随机数种子且可以设置更大的置换次数rng(123456)objML.NumPerm9999;再比如计算距离矩阵的方式 (以下是默认计算方式)更多计算方式参见 MATLAB 官网pdist函数的文档这里不赘述:objML.Distance1euclidean;objML.Distance2(ZI,ZJ)sum(abs(ZI-ZJ),2)./sum(ZIZJ,2);同时Curvature可以设置链接线的弯曲程度和方向比如为负数就是反向弯曲为0就是直线LinkBendMode可以设置为simple就是所有线弯曲方向一致objML.Curvature 0;objML.Curvature-1/3;objML.LinkBendModesimple;左下角布局与字体、配色修改%% Figure and axesfigfigure(Units,normalized,Position,[.05,.15,.72,.72]);axaxes(Parent,fig,Position,[.06,.05,.88,.9]);%% Draw heatmap[rho,pval]corr(Data1);objHMSHeatmap(ax,rho,Format,sq);objHM.draw();objHM.setText()objHM.showStars(pval,Levels,[0.05,0.01,0.001],CorrLabel,off)objHM.setVarName(labels)objHM.setType(triu0);objHM.setRowLabelLocation(right)objHM.setColLabelLocation(top)objHM.setRowLabel(Visible,on)objHM.setColLabel(Visible,on)delete(objHM.Colorbar)% Apply a custom colormap with 25 colors (应用自定义 25 色 colormap)colormap(slanCM(102,25))% Adjust font properties for labels (调整标签字体)set([objHM.rowLabelHdl,objHM.colLabelHdl],FontSize,14,FontName,Helvetica)%% Draw mantel linksobjMLSMantelLink(ax,Data1,Data2,Group,group);objML.GroupNamegroupName;objML.Layouttril;% Customize colors (自定义颜色)objML.PColor[0,64,115;79,148,204;224,224,224]./255;objML.NodeColor1[184,207,248]./255;objML.NodeColor2[184,207,248]./255;objML.draw()% Adjust legend and group label fonts (调整图例和组标签字体)set(objML.legendTitleHdl,FontName,Helvetica)set(objML.legendTickLabelHdl,FontSize,13,FontName,Helvetica)set(objML.groupLabelHdl,FontSize,14,FontName,Helvetica)结以上绘图代码需要下载工具函数才能运行【gitee】https://gitee.com/slandarer/matlab-special-heatmap【github】https://github.com/slandarer/MATLAB-special-heatmap【fileexchange】https://www.mathworks.com/matlabcentral/fileexchange/125520-special-heatmap