#载入packages
library(vegan) #计算矩阵
library(pheatmap) #绘图, 虽然这个package也可以计算矩阵, 但方法不适于生物多样性分析
library(RColorBrewer) #这个package提供"渐变(seq)", "极端(div)", "离散(qual)" 三种好用的配色方案

#读入community 矩阵, 行为样本, 列为OTUs. 将输入的数据赋值给"fcommunity", 第一行row为变量名, 第一列为行名
fcommunity<-data.frame(read.table("xxxxx.csv", header=TRUE, sep="", row.names=1))

#
#(1)计算bray-curtis dissimilarity matrix, 并赋值给"fbray.dist", 算矩阵的时候OTU数量"0", 绘图时置换成"NA"
fbray.dist<-vegdist(fcommunity, method="bray", binary=FALSE, diag=FALSE, upper=FALSE, na.rm=FALSE)

#(2)已有矩阵: 由于计算多样性时, 样本测序深度(序列数)差异会造成α多样性差异(同理对β多样性也有影响?), 所以要按OTU最少的样本进行重抽样将数据均一化, MOTHUR可以实现这个功能, 所以我用MOTHUR计算的Bray-crutis dissimilarity distance. 在vegan中, 这个功能怎么实现?

fbray.dist<-as.dist(read.table("xxxxx.dist.csv", header=TRUE, row.names="groups", fill=TRUE))

#设置annotation:
ann_row=data.frame(SampleSites=factor(rep(c("xx1", "xx2", "xx3"), c(3, 3, 3))))
rownames(fcommunity)=paste("Groups", 1:9, sep="")
rownames(annotation_row)=paste("Groups", 1:9, sep="")

#设置annotation顔色
ann_colors=list(SampleSites=c(xx1="#FFA07A", xx2="#FFC125", xx3="#DDA0DD"))

#设置labels
labels_row=c("A11", "A12", "A13", "A21", "A22", "A23", "A31", "A32", "A33")

#绘图, 具体的参量设定请参考pheatmap的mannul
pheatmap (

fcommunityNA, color=colorRampPalette(brewer.pal(n=9, name="GnBu")[3:9])(3500),

breaks=0:3500, border_color="grey60", cellwidth=4, cellheight=30,

kmeans_k=NA, scale="none",

cluster_rows=TRUE, cluster_cols=FALSE, clustering_distance_rows=fbray.dist,

clustering_distance_cols="none", clustering_method="complete",

cutree_rows=2, treeheight_row=70,

legend=TRUE, legend_labels=0:3500,

annotation_row=ann_row, annotation_col=NA, annotation=NA, annotation_colors=ann_colors,

annotation_legend=TRUE,

drop_levels=TRUE, show_rownames=T, show_colnames=F, main=NA,

fontsize=12, fontsize_row=12, fontsize_col=4,

gaps_row=NULL, gaps_col=NULL,

labels_row=labels_row, labels_col=NULL,

filename="fheatmap.tiff", width=NA, height=NA, silent=FALSE

)