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

魁首哥

作者

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

今天就跟大家聊聊有关R语言中的Metastats分析如何理解，可能很多人都不太了解，为了让大家更加了解，小编给大家总结了以下内容，希望大家根据这篇文章可以有所收获。

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.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(iin2: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)

结果如下所示：