R语言中的Metastats分析如何理解

Anosim、Adonis、MRPP等基于群落的组间差异分析可以快速的对分组的有效性进行评估。然而，有时候我们还想进一步知道不同区组的微生物群落差异在哪里，也即那些物种是显著差异的。这时候我们能想到的最简单的办法就是对所有物种按照分组进行显著性检验，这时候我们对于一个数据集进行了多重检验，则需要p值校正来获得更准确的结果。

假设检验是一种概率判断，因为小概率事件发生了所以我们拒绝假设。然而同时多次做这种概率判断，也会出错。例如当我们进行多重独立比较相关性时，假如有k个变量，那么需要进行n=k(k-1)/2个相关性分析，每个相关性均检验一次。在显著水平0.05（置信水平0.95）的情况下做出显著性判断，其正确概率为0.95，而n个独立检验均正确的概率为0.95ⁿ。若要使所有检验结果正确的概率大于0.95，则需要调整显著水平或更常用的p值校正，一个常见的方法是Bonferroni校正，其原理为在同一数据集做n个独立的假设检验，那么每一个检验的显著水平应该为只有一个检验时的1/n。例如我们只做两个变量相关检验，那么显著水平0.05，假如同时做一个数据集5个变量相关检验，因为要检验10次，那么显著水平应为0.005，因此做Bonferroni校正后判断为显著的检验p值为原来p值的10倍。

在R中p值校正可以使用p.adjust()函数，其使用方法如下所示：

p.adjust(p, method=p.adjust.methods, n=length(p))

其中p为显著性检验的结果（为数值向量），n为独立检验次数，一般为length(p)，method为校正方法，常用的方法有"bonferroni"、"holm"、"hochberg"、"hommel"、"BH"、"fdr"、"BY"、"none"。其中刚刚提到的"bonferroni"最为保守，也即校正后p值变大最多，一般不是很常用，其余方法均为各种修正方法。

校正后的p值常称为q值，使用Benjamini-Hochberg（BH）方法校正的p值也称为错误发现率（false discovery rate，FDR）。

#读取抽平后的OTU_table和环境因子信息data=read.csv("otu_table.csv", header=TRUE, row.names=1)envir=read.table("environment.txt", header=TRUE)rownames(envir)=envir[,1]env=envir[,-1]#筛选高丰度物种并将物种数据标准化means=apply(data, 1, mean)otu=data[names(means[means>10]),]otu=t(otu)#根据地理距离聚类kms=kmeans(env, centers=3, nstart=22)Position=factor(kms$cluster)newotu=data.frame(Group=Position, otu)#进行多重Kruskal-Wallis秩和检验与p值校正pvalue=t(otu)[,1:2]colnames(pvalue)=c("p-value", "q-value")for (i in 2:ncol(newotu)) {  t=kruskal.test(newotu[,i]~newotu[,1])  pvalue[i-1,1]=t$p.value}pvalue[,2]=p.adjust(pvalue[,1], method="BH", n=nrow(pvalue))pvalue=pvalue[order(pvalue[,1]),]

#筛选q小于0.05的物种top=pvalue[pvalue[,2]<0.05,]topdata=cbind(Group=newotu[,1], newotu[,rownames(top)])#热图展示差异物种library(gplots)mycol=colorRampPalette(c("white", "blue","green", "red", "red"))(100)sidecol=c(rep("red",7), rep("green",12), rep("blue",3))heatmap.2(log(as.matrix(topdata[,-1])+1), Rowv=FALSE,          dendrogram="column",trace="none", col=mycol,          RowSideColors=sidecol, keysize = 1.2, key.title="",          cexRow = 1.2, cexCol = 0.5)