用多核R分析GWAS数据

Question

我用R来分析全基因组的关联研究数据。我有大约500,000个潜在的预测变量(单核苷酸多态性，或SNPs)，并希望测试它们之间的联系和持续的结果(在这种情况下，低密度脂蛋白在血液中的浓度)。

我已经写了一个没有问题的脚本。简单地解释一下，我有一个数据对象，名为" data“。每一行对应于研究中的特定病人。列有年龄，性别，体重指数(BMI)和血LDL浓度。还有五十万其他列与SNP数据。

我目前正在使用for循环运行线性模型50万次，如下所示：

# 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的并行版本可能是一种选择，但也不知道如何做到这一点。

万事如意，

亚历克斯

Answer 1

要给您一个启动:如果您使用Linux，您可以执行包含在multicore包中的parallel方法。虽然您需要在使用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

在这里，您将创建一个函数，该函数返回一个具有所需数量的向量，并将索引应用于此。但我无法检查这一点，因为我没有访问数据的权限，但它应该能工作。