Usando multinúcleo en R para analizar los datos de GWAS

Question

Estoy usando R para analizar los datos del estudio de asociación del genoma completo. Tengo alrededor de 500,000 posibles variables predictoras (polimorfismos de un solo nucleótido o SNP) y quiero probar la asociación entre cada uno de ellos y un resultado continuo (en este caso, concentración de lipoproteínas de baja densidad en la sangre).

Ya he escrito una secuencia de comandos que hace esto sin ningún problema. Para explicar brevemente, tengo un objeto de datos, llamado "Datos". Cada fila corresponde a un paciente particular en el estudio. Hay columnas por edad, sexo, índice de masa corporal (IMC) y concentración de LDL en la sangre. También hay medio millón de otras columnas con los datos de SNP.

actualmente estoy usando un bucle for para ejecutar el modelo lineal de medio millón de veces, como se muestra:

# 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"] 
} 

Todo esto funciona bien. Sin embargo, realmente me gustaría acelerar mi análisis. Por lo tanto, he estado experimentando con los paquetes multinúcleo, DoMC y foreach.

Mi pregunta es, ¿podría alguien ayudarme a adaptar este código utilizando el esquema foreach?

Estoy ejecutando el script en un servidor Linux que aparentemente tiene 16 núcleos disponibles. Intenté experimentar con el paquete foreach, y mis resultados al usarlo han sido comparativamente peores, lo que significa que toma más tiempo para ejecutar el análisis utilizando foreach.

Por ejemplo, he tratado de guardar los objetos del modelo lineal como se muestra:

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

Esto toma más de dos veces más que utilizando sólo un habitual de bucle. ¡Cualquier consejo sobre cómo hacer esto mejor o más rápido sería apreciado! Entiendo que usar la versión paralela de lapply podría ser una opción, pero tampoco sé cómo hacer esto.

Todo lo mejor,

Alex

Answer 1

para darle un arranque: si utiliza Linux, que puede hacer el enfoque multicore contenida dentro del paquete parallel. Mientras que necesitabas configurar todo al usar, por ejemplo, el paquete foreach, ya no es necesario con este enfoque. Su código se ejecuta en 16 núcleos simplemente haciendo:

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 

Aquí haces una función que devuelve un vector con las cantidades deseadas, y aplicar los índices sobre esto. Sin embargo, no pude verificar esto ya que no tengo acceso a los datos, pero debería funcionar.