【问题标题】:Using multicore in R to analyse GWAS data在 R 中使用多核分析 GWAS 数据
【发布时间】:2011-12-13 11:35:32
【问题描述】:

我正在使用 R 分析全基因组关联研究数据。我有大约 500,000 个潜在的预测变量(单核苷酸多态性,或 SNP),并想测试它们中的每一个与连续结果(在本例中为血液中的低密度脂蛋白浓度)之间的关联。

我已经编写了一个可以毫无问题地执行此操作的脚本。简单解释一下,我有一个数据对象,称为“数据”。每行对应于研究中的特定患者。有年龄、性别、体重指数 (BMI) 和血液 LDL 浓度列。还有 50 万列包含 SNP 数据。

我目前正在使用 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 方案调整此代码吗?

我在显然有 16 个内核可用的 Linux 服务器上运行脚本。我尝试过使用 foreach 包,但使用它的结果相对较差,这意味着使用 foreach 运行分析需要更长的时间

例如,我尝试保存线性模型对象,如下所示:

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

这比仅使用常规 for 循环花费的时间多一倍。任何有关如何更好或更快速地做到这一点的建议将不胜感激!我知道使用并行版本的 lapply 可能是一种选择,但也不知道该怎么做。

一切顺利,

亚历克斯

【问题讨论】:

  • 更新到 R 2.14 并使用 parallel 包。当我们这样做的时候,给我们一个可重复的例子来工作,肯定也会有所帮助。见this question
  • 乔里斯,谢谢你的建议。我在biomedcentral.com/content/pdf/1471-2105-9-390.pdf 找到了parallel 的文档,现在将阅读它。
  • 包裹 snowfall 是否已下架(不要打我,Dirk)?
  • @RomanLuštrik SMASH 我会做的:-)。 snowfall 应该与 parallel 包一起使用,但在 Linux 上你最好使用 multicore。比snowfall afaik 更容易工作。

标签: r foreach statistics multicore lapply


【解决方案1】:

给你一个启动:如果你使用 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

在这里,您创建了一个函数,该函数返回具有所需数量的向量,并在其上应用索引。我无法检查这个,因为我无权访问数据,但它应该可以工作。

【讨论】:

  • 乔里斯,感谢您的帮助!我已经实施了您的解决方案,并且似乎奏效了。我刚刚完成了一项之前需要 12 多个小时的工作,而它在 15 分钟内就从烤箱中出来了!现在我真希望我在三个月前问过你这个问题!
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2012-01-16
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多