Я использую R для анализа данных широкогеномного ассоциативного исследования. У меня есть около 500 000 потенциальных предикторных переменных (однонуклеотидных полиморфизмов, или SNPs), и я хочу проверить ассоциацию между каждой из них и непрерывным результатом (в данном случае концентрацией липопротеина низкой плотности в крови).
Я уже написал сценарий, который делает это без проблем. Чтобы кратко объяснить, у меня есть объект данных, называемый "Data". Каждая строка соответствует определенному пациенту в исследовании. Есть столбцы для возраста, пола, индекса массы тела (ИМТ) и концентрации ЛПНП в крови. Есть также полмиллиона других столбцов с данными о 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"]
}
Все это работает хорошо. Однако мне бы очень хотелось ускорить анализ. Поэтому я экспериментировал с пакетами multicore, 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 может быть вариантом, но не знаю, как это сделать.
Всего наилучшего,
Alex