在生物群落的研究中，经常会使用物种群落间的距离来评估样本间物种群落的差异程度，这种群落间的距离就是beta多样性。
用来评估物种群落距离的beta多样性指数有很多种，其中最常用的是Bray-Crutis距离和Unifrac距离。
Bray-Crutis距离在计算是同时考虑了物种在群落中是否存在以及物种在群落中的丰度，可以参考alpha多样性指数中的Shannon指数进行理解。
Bray-Crutis距离的计算结果在0-1之间，0代表两个样本的微生物群落没有任何差异，数值越接近1，两个样本微生物群落之间的差异越大。
R语言vegan包中的vegidst()函数可以很方便的计算多种beta多样性距离。
用于计算Bray-Crutis距离的输入文件为抽平后各样品OTU的丰度表格。
⚠️使用前要先处理一下，将第一行第一列位置的“#”删掉，并且删掉表格最后一列的分类学信息。
otu <- read.table('otu_table.txt',row.names = 1,header = T,sep = '\t')otu <- t(otu)
导入数据，并将数据进行转置，之后计算Bray-Crutis距离。
library(vegan)bray <- vegdist(otu,method = 'bray')
计算得到的距离文件为dist类型，需要将其转换为矩阵，才能将其保存到本地并进行下一步绘图。
bray <- as.matrix(bray)write.table(bray,'bray-crutis.txt',sep = '\t')
我们使用pheatmap包根据计算的结果绘制一个聚类热图。
library(pheatmap)library(RColorBrewer)pheatmap(bray,color = colorRampPalette(brewer.pal(7,'RdYlBu'))(100))
可以看出示例样本可以分为两类。
这里只提供了最简单的热图绘制代码，pheatmap绘图参数的详细调整方法可参考之前的推文。
Unifrac距离可以简单的理解为在普通beta多样性距离的基础上引入了OTU的进化距离。
Unifrac距离分为加权和不加权两种，Unweighted Unifrac在计算时只考虑物种在样本中是否存在，Weighted Unifrac在计算时同时考虑物种的存在及其丰度，类似于alpha多样性指数中的丰度度和多样性的区别。
可以使用phyloseq包计算Unifrac距离，phyloseq是集成在Bioconductor中的，第一次安装需要同时安装大量的依赖包。
if (!requireNamespace('BiocManager', quietly = TRUE))install.packages('BiocManager')BiocManager::install('phyloseq')
因为引入了OTU的进化距离，因此输入文件除了OTU的表格之外，还需要OTU代表序列的系统发育树。
这个可以使用QIIME分析自动生成的，也可以使用OTU的代表序列使用MEGA等软件自行绘制。
library(phyloseq)library(GUniFrac)tree <- read.tree('OTUs.tre')unifrac <- phyloseq(otu_table(otu,taxa_are_rows = F),phy_tree(tree))#加权UnifracW.unifrac <- distance(unifrac,method = 'wunifrac')#非加权UnifracU.unifrac <- distance(unifrac,method = 'unifrac')
W.unifrac <- as.matrix(W.unifrac)write.table(W.unifrac,'W.unifrac.txt',sep = '\t')U.unifrac <- as.matrix(U.unifrac)write.table(U.unifrac,'U.unifrac.txt',sep = '\t')结果的保存方式与Bray-Crutis距离一致。
因为Unifrac距离的结果同时包含加权和非加权两种，因而聚类热图如果还是按照Bray-Crutis距离的方式进行绘制就没什么新意了，这里的最终目的是在一幅图中同时展示加权和非加权的距离及其聚类模式。
经过思考之后打算使用corrplot包按照绘制相关性图的方式绘制Unifrac聚类热图.
首先须需要对计算得到的距离矩阵进行处理，使得矩阵的上三角为加权Unifrac，下三角为非加权Unifrac。
Unifrac <- W.unifracUnifrac[lower.tri(Unifrac)] <- U.unifrac[lower.tri(U.unifrac)]
接着我先绘制了一个简单的热图。
library(corrplot)corrplot(Unifrac,method = 'color',p.mat = Unifrac,insig = 'p-value',sig.level = -1,cl.lim = c(0,1),col = colorRampPalette(c('Red','white','navy'))(100),tl.col = 'black')
这幅图中只是分别在上三角和下三角展示了样品间的距离，并且标出了距离的具体数值，但是并没有进行聚类，虽然corrplot自身带有聚类功能，但是距离是按照整行和整列进行聚类，并无法对上下三角的数据进行分别聚类。
由于绘制一张热图横纵坐标对应样本的位置是只能有一种排列方式，虽然可以使用上下三角展示两种距离，但是无法在一张图中同时展示两种聚类方式。
想到了一种方法，就是分别画出上三角和下三角两幅图，之后在拼在一起。
可是遇到了不可抗力，暂时无法解决。
如果依然使用corrplot或者其它的热图绘制方法，是很容易得到两个带有不同聚类的上三角和下三角图像的，但是两个三角形图的要怎么能用R拼在一起呢？
用base plot的方式也能够得到两个三角形的图，也能够拼在一起，但是无法自动根据样品聚类模式调整各样本的位置，理论上倒是能解决这个问题，但是工程有点大啊，相当于写一个专门的R包来干这个事，性价比极低。
还有一种快捷的方式，就是输出两个png格式的带有聚类的上下三角图，之后用PS拼在一起，再手动修整样本名称和图例位置，我用这种方式画了一个，数据本身问题，聚类比较混乱，重在表达一个意思，大家看看就行了，但是这种方式没有可重复性，每次都要重新修一遍。
大家如果有好的想法和思路来解决上述问题，可以加我微信或者后台私信我，十分感谢！
联系客服
