Utilizzo di multicore in R per analizzare i dati GWAS

Question

Sto utilizzando R per analizzare i dati di studio dell'associazione a livello dell'intero genoma. Ho circa 500.000 variabili predittive potenziali (polimorfismi a singolo nucleotide o SNP) e voglio testare l'associazione tra ciascuna di esse e un risultato continuo (in questo caso la concentrazione di lipoproteine ​​a bassa densità nel sangue).

Ho già scritto uno script che lo fa senza problemi. Per spiegare brevemente, ho un oggetto dati, chiamato "Dati". Ogni riga corrisponde a un particolare paziente nello studio. Ci sono colonne per età, sesso, indice di massa corporea (BMI) e concentrazione di LDL nel sangue. Ci sono anche mezzo milione di altre colonne con i dati SNP.

Attualmente sto usando un ciclo for per eseguire il modello lineare mezzo milione di volte, come mostrato:

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

Tutto questo funziona benissimo. Tuttavia, mi piacerebbe davvero accelerare la mia analisi. Ho quindi sperimentato i pacchetti multicore, DoMC e foreach.

La mia domanda è: qualcuno potrebbe aiutarmi ad adattare questo codice utilizzando lo schema foreach?

Sto eseguendo lo script su un server Linux che apparentemente dispone di 16 core disponibili. Ho provato a sperimentare con il pacchetto foreach e i risultati che ne ho ricavato sono stati comparativamente peggiori, il che significa che è necessario più lungo per eseguire l'analisi utilizzando foreach.

Per esempio, ho provato salvare gli oggetti modello lineare come mostrato:

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

questo richiede più di due volte più a lungo utilizzando solo un normale ciclo for. Qualunque consiglio su come farlo meglio o più rapidamente sarebbe apprezzato! Capisco che l'uso della versione parallela di lapply potrebbe essere un'opzione, ma non so come fare neanche questo.

Tutto il meglio,

Alex

Answer 1

per darvi un avvio: Se si utilizza Linux, si può fare l'approccio multicore contenuta all'interno del pacchetto parallel. Considerando che è necessario impostare il tutto quando si utilizza ad esempio il pacchetto foreach, non è più necessario con questo approccio. Il tuo codice sarebbe stato eseguito su 16 core semplicemente facendo:

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 

Qui si effettua una funzione che restituisce un vettore con le quantità volute, e applicare gli indici su questo. Non ho potuto controllare questo perché non ho accesso ai dati, ma dovrebbe funzionare.