|
生科云網(wǎng)址：https://www.
本文由微科盟胡天龍根據(jù)實踐經(jīng)驗而整理，希望對大家有幫助,。
微科盟原創(chuàng)微文,，歡迎轉(zhuǎn)發(fā)轉(zhuǎn)載，轉(zhuǎn)載須注明來源《微生態(tài)》公眾號,。 共現(xiàn)網(wǎng)絡(luò) (co-occurrence network) 和群落構(gòu)建 (assembly process) 分析日益成為微生物生態(tài)學分析中重要的組成部分,，成為目前文章發(fā)表的熱點技術(shù)。這些分析中往往會有大量的循環(huán)語句,，如果順序執(zhí)行的話往往非常耗時,。R語言的并行計算可以有效地解決耗時的問題，加快程序運行的速度,。本文主要介紹如何利用R語言并行計算spearman相關(guān)系數(shù),，加快共現(xiàn)網(wǎng)絡(luò)（co-occurrence network）構(gòu)建速度。后期的的文章將逐一介紹群落構(gòu)建分析的并行計算,，包括R語言并行計算beta-NTI (beta-nearest taxon index),、β多樣性零偏差 (null deviation of beta diversity) 和Raup-Crick距離 (Raup-Crick dissimilarity)。利用spearman相關(guān)性分析是構(gòu)建共現(xiàn)網(wǎng)絡(luò)的重要方法,，但由于OTU table往往有成千上萬行,，用R自帶的corr.test()函數(shù)計算較為費時，嚴重制約我們的分析速度,。對spearman相關(guān)性分析進行并行化運行可大大節(jié)省計算時間,，為此我們手寫了spearman相關(guān)性分析函數(shù)來實現(xiàn)并行化運行,。為方便講解，本文以我們常見的OTU table 數(shù)據(jù)為例（聯(lián)系您所添加的任一微科盟組學老師即可免費領(lǐng)取）,，對OTU進行兩兩spearman相關(guān)性分析,，獲得相關(guān)系數(shù)r和顯著性p值。我們將自己手寫的函數(shù)network_construct()與psych包中的corr.test()函數(shù)兩者運行時間和計算的結(jié)果進行了比較,，我們自己的函數(shù)network_construct()計算時間遠遠少于corr.test()函數(shù)且結(jié)果相同,，具體的R代碼見下文。如果對您有幫助,，請三連一波哦~
點贊,，在看，轉(zhuǎn)發(fā)?。,。?/span> 我們經(jīng)常使用corr.test () 函數(shù)計算OTU之間的相關(guān)性,，但該函數(shù)在面對較多的OTU時速度較慢,。 library(psych) system.time(corr.test(otu[,1:500],use="pairwise",method="spearman",adjust="fdr",alpha=.05)) |
運行時間如下：（elapsed欄為程序運行的時間） 
圖1 原函數(shù)運行時間 為了加快計算速度，我們根據(jù)相關(guān)系數(shù)矩陣的計算特點,，進行了并行化重寫,，代碼如下：#otu_table行名為樣點名，列名為OTU名 #threads為使用的CPU核數(shù) #本函數(shù)默認使用'fdr’進行P值矯正,，可參考corr.test()函數(shù)的R幫助文檔 network_construct <- function(otu_table,threads){ library(foreach) library(doParallel) library(abind) otu_table2 <- apply(otu_table,2,rank) r <- function(rx,ry){ n <- length(rx) lxy <- sum((rx-mean(rx))*(ry-mean(ry))) lxx <- sum((rx-mean(rx))^2) lyy <- sum((ry-mean(ry))^2) r <- lxy/sqrt(lxx*lyy) t <- (r * sqrt(n - 2))/sqrt(1 - r^2) p <- -2 * expm1(pt(abs(t), (n - 2), log.p = TRUE)) return(c(r,p)) } arraybind <- function(...){ abind(...,along = 3,force.array=TRUE) } nc <- ncol(otu_table) registerDoParallel(cores = threads) corr <- foreach (i = 1:nc,.combine = "arraybind") %dopar%{ corr1 <- matrix(rep(0,2*nc),nrow = 2,ncol=nc) for(j in 1:nc) { if(j > i) corr1[,j] <- r(otu_table2[,i],otu_table2[,j]) } corr <- corr1 } rr <- corr[1,,] rr <- rr+t(rr) diag(rr) <- 1 pp <- corr[2,,] lp <- lower.tri(pp) pa <- pp[lp] pa <- p.adjust(pa, "fdr") pp[lower.tri(pp, diag = FALSE)] <- pa pp <- pp+t(pp) rownames(pp) <- colnames(otu_table) colnames(pp) <- colnames(otu_table) rownames(rr) <- colnames(otu_table) colnames(rr) <- colnames(otu_table) return(list(r = rr,p = pp)) } |
使用我們自己構(gòu)造的函數(shù)后,，可利用10核進行計算，運行時間如下： 
圖2 自建函數(shù)十核運行時間 我們可以看到使用corr.test () 函數(shù)進行計算需要313秒,，而使用R代碼進行并行計算僅需0.99秒,。我們再比較一下不同線程情況下所需的時間。單核運行代碼如下： system.time(network_construct(otu[,1:500],1)) |
圖3 自建函數(shù)單核運行時間 即使使用單核進行計算速度,，也要比corr.test () 函數(shù)快上很多,，畢竟該函數(shù)里面有很多判斷語句， 計算更加費時,。當我們加大數(shù)據(jù)量時,，單核和多核的區(qū)別就更加明顯。下圖是單核與10核運行的比較： 
圖4 自建函數(shù)單核與十核運行時間的比較 我們可以看到計算量越大,，多核的性能就越優(yōu)越,。當我們有成千上萬個OTU時可以節(jié)省很多時間。那么我們運行的結(jié)果與corr.test () 函數(shù)的結(jié)果是否一致呢,？ res1<- corr.test(otu[,1:50],use="pairwise",method="spearman",adjust="fdr",alpha=.05) res2 <- network_construct(otu[,1:50],10) res1$r[1:5,1:5] res2$r[1:5,1:5] |
將我們構(gòu)建的函數(shù)network_construct () 與psych包中的corr.test () 的結(jié)果進行比較,，結(jié)果如下： 
圖5 自建函數(shù)與原函數(shù)結(jié)果比較
我們計算的結(jié)果與原函數(shù)的結(jié)果完全相同,，可以放心使用哈,，現(xiàn)在就可以把我們的函數(shù)復制過去用于自己的項目啦,。
|
