溫馨提示×

溫馨提示×

您好，登錄后才能下訂單哦！

密碼登錄×
登錄注冊×
其他方式登錄
點(diǎn)擊 登錄注冊 即表示同意《億速云用戶服務(wù)條款》

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

發(fā)布時(shí)間：2021-12-28 12:08:48 來源：億速云 閱讀：473 作者：柒染 欄目：大數(shù)據(jù)

今天就跟大家聊聊有關(guān)R語言中的Metastats分析如何理解，可能很多人都不太了解，為了讓大家更加了解，小編給大家總結(jié)了以下內(nèi)容，希望大家根據(jù)這篇文章可以有所收獲。

Anosim、Adonis、MRPP等基于群落的組間差異分析可以快速的對分組的有效性進(jìn)行評估。然而，有時(shí)候我們還想進(jìn)一步知道不同區(qū)組的微生物群落差異在哪里，也即那些物種是顯著差異的。這時(shí)候我們能想到的最簡單的辦法就是對所有物種按照分組進(jìn)行顯著性檢驗(yàn)，這時(shí)候我們對于一個(gè)數(shù)據(jù)集進(jìn)行了多重檢驗(yàn)，則需要p值校正來獲得更準(zhǔn)確的結(jié)果。

在不同區(qū)組中尋找差異物種常用的兩個(gè)工具是  Metastats  和  LEfSe  。拋開這兩個(gè)工具本身，從算法原理上來說，  Metastats  實(shí)際上是非參數(shù)多重檢驗(yàn)和  p  值校正的整合，而  LEfSe  則是  Metastats  和  LDA  判別的整合。當(dāng)然，由于  Metastats  采用的非參數(shù)  t  檢驗(yàn)，只能分析兩個(gè)分組；而  LEfSe  則因?yàn)槭褂玫?  Kruskal-Wallis  秩和檢驗(yàn)可以分析兩個(gè)以上的分組。當(dāng)我們明白了他們的原理，實(shí)際上可以不用拘泥于兩個(gè)工具本身，可以自己在  R  中選擇合適的方法來進(jìn)行分析。

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

p值校正

假設(shè)檢驗(yàn)是一種概率判斷，因?yàn)樾「怕适录l(fā)生了所以我們拒絕假設(shè)。然而同時(shí)多次做這種概率判斷，也會出錯(cuò)。例如當(dāng)我們進(jìn)行多重獨(dú)立比較相關(guān)性時(shí)，假如有k個(gè)變量，那么需要進(jìn)行n=k(k-1)/2個(gè)相關(guān)性分析，每個(gè)相關(guān)性均檢驗(yàn)一次。在顯著水平0.05（置信水平0.95）的情況下做出顯著性判斷，其正確概率為0.95，而n個(gè)獨(dú)立檢驗(yàn)均正確的概率為0.95n。若要使所有檢驗(yàn)結(jié)果正確的概率大于0.95，則需要調(diào)整顯著水平或更常用的p值校正，一個(gè)常見的方法是Bonferroni校正，其原理為在同一數(shù)據(jù)集做n個(gè)獨(dú)立的假設(shè)檢驗(yàn)，那么每一個(gè)檢驗(yàn)的顯著水平應(yīng)該為只有一個(gè)檢驗(yàn)時(shí)的1/n。例如我們只做兩個(gè)變量相關(guān)檢驗(yàn)，那么顯著水平0.05，假如同時(shí)做一個(gè)數(shù)據(jù)集5個(gè)變量相關(guān)檢驗(yàn)，因?yàn)橐獧z驗(yàn)10次，那么顯著水平應(yīng)為0.005，因此做Bonferroni校正后判斷為顯著的檢驗(yàn)p值為原來p值的10倍。

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

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

其中p為顯著性檢驗(yàn)的結(jié)果（為數(shù)值向量），n為獨(dú)立檢驗(yàn)次數(shù)，一般為length(p)，method為校正方法，常用的方法有"bonferroni"、"holm"、"hochberg"、"hommel"、"BH"、"fdr"、"BY"、"none"。其中剛剛提到的"bonferroni"最為保守，也即校正后p值變大最多，一般不是很常用，其余方法均為各種修正方法。

校正后的p值常稱為q值，使用Benjamini-Hochberg（BH）方法校正的p值也稱為錯(cuò)誤發(fā)現(xiàn)率（false discovery rate，F(xiàn)DR）。



接下來，我用相同數(shù)據(jù)為例，尋找不同分組間顯著差異的物種：
#讀取抽平后的OTU_table和環(huán)境因子信息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]#篩選高豐度物種并將物種數(shù)據(jù)標(biāo)準(zhǔn)化means=apply(data, 1, mean)otu=data[names(means[means>10]),]otu=t(otu)#根據(jù)地理距離聚類kms=kmeans(env, centers=3, nstart=22)Position=factor(kms$cluster)newotu=data.frame(Group=Position, otu)#進(jìn)行多重Kruskal-Wallis秩和檢驗(yàn)與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]),]

接下來我們可以篩選顯著差異的物種并進(jìn)行可視化：

#篩選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)

結(jié)果如下所示：

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

由熱圖可以看出，這些物種確實(shí)存在明顯的組間差異，根據(jù)研究需求可以導(dǎo)出這些物種名并進(jìn)行后續(xù)的系統(tǒng)發(fā)育與功能分析。

看完上述內(nèi)容，你們對R語言中的Metastats分析如何理解有進(jìn)一步的了解嗎？如果還想了解更多知識或者相關(guān)內(nèi)容，請關(guān)注億速云行業(yè)資訊頻道，感謝大家的支持。

向AI問一下細(xì)節(jié)

免責(zé)聲明：本站發(fā)布的內(nèi)容（圖片、視頻和文字）以原創(chuàng)、轉(zhuǎn)載和分享為主，文章觀點(diǎn)不代表本網(wǎng)站立場，如果涉及侵權(quán)請聯(lián)系站長郵箱：is@yisu.com進(jìn)行舉報(bào)，并提供相關(guān)證據(jù)，一經(jīng)查實(shí)，將立刻刪除涉嫌侵權(quán)內(nèi)容。

AI
助
手