您好，登录后才能下订单哦！
密码登录
登录注册
点击 登录注册 即表示同意《亿速云用户服务条款》
今天就跟大家聊聊有关R语言中的Metastats分析如何理解，可能很多人都不太了解，为了让大家更加了解，小编给大家总结了以下内容，希望大家根据这篇文章可以有所收获。
Anosim、Adonis、MRPP等基于群落的组间差异分析可以快速的对分组的有效性进行评估。然而，有时候我们还想进一步知道不同区组的微生物群落差异在哪里，也即那些物种是显著差异的。这时候我们能想到的最简单的办法就是对所有物种按照分组进行显著性检验，这时候我们对于一个数据集进行了多重检验，则需要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)
结果如下所示：
看完上述内容，你们对R语言中的Metastats分析如何理解有进一步的了解吗？如果还想了解更多知识或者相关内容，请关注亿速云行业资讯频道，感谢大家的支持。
免责声明：本站发布的内容（图片、视频和文字）以原创、转载和分享为主，文章观点不代表本网站立场，如果涉及侵权请联系站长邮箱：is@yisu.com进行举报，并提供相关证据，一经查实，将立刻删除涉嫌侵权内容。