文章来自公众号：微生态与微进化
Anosim、Adonis、MRPP等基于群落的组间差异分析可以快速的对分组的有效性进行评估。然而，有时候我们还想进一步知道不同区组的微生物群落差异在哪里，也即那些物种是显著差异的。这时候我们能想到的最简单的办法就是对所有物种按照分组进行显著性检验，这时候我们对于一个数据集进行了多重检验，则需要p值校正来获得更准确的结果。
在不同区组中寻找差异物种常用的两个工具是Metastats和LEfSe。抛开这两个工具本身，从算法原理上来说，Metastats实际上是非参数多重检验和p值校正的整合，而LEfSe则是Metastats和LDA判别的整合。当然，由于Metastats采用的非参数t检验，只能分析两个分组；而LEfSe则因为使用的Kruskal-Wallis秩和检验可以分析两个以上的分组。当我们明白了他们的原理，实际上可以不用拘泥于两个工具本身，可以自己在R中选择合适的方法来进行分析。
p值校正
假设检验是一种概率判断，因为小概率事件发生了所以我们拒绝假设。然而同时多次做这种概率判断，也会出错。例如当我们进行多重独立比较相关性时，假如有k个变量，那么需要进行n=k(k-1)/2个相关性分析，每个相关性均检验一次。在显著水平0.05（置信水平0.95）的情况下做出显著性判断，其正确概率为0.95，而n个独立检验均正确的概率为0.95n。若要使所有检验结果正确的概率大于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)
结果如下所示：
由热图可以看出，这些物种确实存在明显的组间差异，根据研究需求可以导出这些物种名并进行后续的系统发育与功能分析。示例数据下载请移步公众号。
