Использование многоядерности в R для анализа данных GWAS

Question

Я использую R для анализа данных исследования ассоциации всего генома.У меня около 500 000 потенциальных предикторов (однонуклеотидных полиморфизмов или SNP), и я хочу проверить связь между каждым из них и непрерывным исходом (в данном случае концентрация липопротеинов низкой плотности в крови).

Я уже написал скрипт, который делает это без проблем.Чтобы кратко объяснить, у меня есть объект данных, который называется «Данные».Каждый ряд соответствует конкретному пациенту в исследовании.Есть столбцы для возраста, пола, индекса массы тела (ИМТ) и концентрации ЛПНП в крови.Есть также полмиллиона других столбцов с данными 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?

Я выполняю сценарийна сервере Linux, который, очевидно, имеет 16 ядер.Я попытался поэкспериментировать с пакетом foreach, и мои результаты с его использованием были сравнительно хуже, то есть для выполнения анализа с использованием foreach требуется больше .

Например, япопытался сохранить объекты линейной модели, как показано:

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

Это занимает более чем вдвое больше времени, чем при использовании обычного цикла for.Любые советы о том, как сделать это лучше или быстрее, будут оценены!Я понимаю, что можно использовать параллельную версию lapply, но я тоже не знаю, как это сделать.

Всего наилучшего,

Alex

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

Здесь вы создаете функцию, которая возвращает вектор с требуемыми величинами, и применяете индексы над этим.Я не мог проверить это, хотя у меня нет доступа к данным, но это должно работать.