在R中使用多核來分析GWAS數據

我使用R來分析全基因組關聯研究數據。我有大約500,000個潛在的預測變量（單核苷酸多態性，或SNPs），並且想要測試它們各自與連續結果（在這種情況下爲血液中低密度脂蛋白濃度）之間的關聯。在R中使用多核來分析GWAS數據

我已經寫了一個沒有問題的腳本。簡要解釋一下，我有一個名爲「Data」的數據對象。每行對應於研究中的特定患者。有年齡，性別，體重指數（BMI）和血液LDL濃度的列。 SNP數據還有50萬列。

我目前使用for循環運行線性模型五十萬次，如圖所示：

# Repeat loop half a million times 
for(i in 1:500000) { 

# Select the appropriate SNP 
SNP <- Data[i] 

# For each iteration, perform linear regression adjusted for age, gender, and BMI and save the result in an object called "GenoMod" 
GenoMod <- lm(bloodLDLlevel ~ SNP + Age + Gender + BMI, data = Data) 

# For each model, save the p value and error for each SNP. I save these two data points in columns 1 and 2 of a matrix called "results" 
results[i,1] <- summary(GenoMod)$coefficients["Geno","Pr(>|t|)"] 
results[i,2] <- summary(GenoMod)$coefficients["Geno","Estimate"] 
}

所有這一切工作正常。不過，我真的很想加快我的分析。因此，我一直在試驗多核，DoMC和foreach軟件包。

我的問題是，有人可以幫助我使用foreach方案來調整此代碼嗎？

我在Linux服務器上運行腳本，顯然有16個可用的內核。我嘗試過使用foreach包進行試驗，並且使用它的結果一直比較糟糕，這意味着需要使用foreach來運行分析，需要更長的更長的。

舉例來說，我已經試過節約線性模型對象，如下所示：

library(doMC) 
registerDoMC() 
results <- foreach(i=1:500000) %dopar% { lm(bloodLDLlevel ~ SNP + Age + Gender + BMI, data = Data) }

這需要超過兩倍的時間多爲使用只是一個普通的for循環。任何建議如何更好或更快地做到這一點，將不勝感激！我知道使用lapply的平行版本可能是一種選擇，但不知道如何做到這一點。

一切順利，

亞歷

來源

2011-12-13 Alexander

更新到R 2.14並使用'parallel'軟件包。儘管我們處於這個階段，但給我們一個可重複使用的例子，肯定也會有所幫助。見[這個問題]（http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example） –

喬里斯，感謝您的意見。我已經在http://www.biomedcentral.com/content/pdf/1471-2105-9-390.pdf找到了'parallel'的文檔，現在將閱讀它。 – Alexander

桌子上有'snowfall'包裝嗎（不要打我，Dirk）？ –

爲了給你一個啓動：如果你使用Linux，你可以做包含parallel包內的multicore方法。而當你使用例如foreach包的時候你需要設置所有的東西，這對於這種方法不再需要。您的代碼將在16個內核通過簡單地做下運行：

require(parallel) 

mylm <- function(i){ 
    SNP <- Data[i] 
    GenoMod <- lm(bloodLDLlevel ~ SNP + Age + Gender + BMI, data = Data) 
    #return the vector 
    c(summary(GenoMod)$coefficients["Geno","Pr(>|t|)"], 
    summary(GenoMod)$coefficients["Geno","Estimate"]) 
} 

Out <- mclapply(1:500000, mylm,mc.cores=16) # returns list 
Result <- do.call(rbind,Out) # make list a matrix

在這裏，你作出這樣的返回與所需量的載體功能，並在此應用的指數。我無法檢查這個，因爲我沒有訪問數據，但它應該工作。

來源

2011-12-13 12:09:18

Joris，感謝您的幫助！我已經實施了你的解決方案，它看起來很有效。我剛剛完成了一項超過12小時的工作，並在15分鐘內從烤箱中出來！現在我只希望三個月前我問過你這個問題！ – Alexander

在R中使用多核來分析GWAS數據

回答

相關問題