title: Using multicore in R to analyse GWAS data
body: <p>I am using R to analyze genome-wide association study data. I have about 500,000 potential predictor variables (single-nucleotide polymorphisms, or SNPs) and want to test the association between each of them and a continuous outcome (in this case low-density lipoprotein concentration in the blood).</p>
  
  <p>I have already written a script that does this without problem. To briefly explain, I have a data object, called "Data". Each row corresponds to a particular patient in the study. There are columns for age, gender, body mass index (BMI), and blood LDL concentration. There are also half a million other columns with the SNP data.</p>
  
  <p>I am currently using a for loop to run the linear model half a million times, as shown:</p>
  
  <pre><code># Repeat loop half a million times
  for(i in 1:500000) {
  
  # Select the appropriate SNP
  SNP &lt;- Data[i]
  
  # For each iteration, perform linear regression adjusted for age, gender, and BMI and save the result in an object called "GenoMod"
  GenoMod  &lt;- 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] &lt;- summary(GenoMod)$coefficients["Geno","Pr(&gt;|t|)"]
  results[i,2] &lt;- summary(GenoMod)$coefficients["Geno","Estimate"]
  }
  </code></pre>
  
  <p>All of that works fine. However, I would really like to speed up my analysis. I've therefore been experimenting with the multicore, DoMC, and foreach packages.</p>
  
  <p>My question is, could someone please help me adapt this code using the foreach scheme?</p>
  
  <p>I am running the script on a Linux server that apparently has 16 cores available. I've tried experimenting with the foreach package, and my results using it have been comparatively worse, meaning that it takes <em>longer</em> to run the analysis using foreach.</p>
  
  <p>For example, I've tried saving the linear model objects as shown:</p>
  
  <pre><code>library(doMC)
  registerDoMC()
  results &lt;- foreach(i=1:500000) %dopar% { lm(bloodLDLlevel ~ SNP + Age + Gender + BMI, data = Data) }
  </code></pre>
  
  <p>This takes more than twice as long as using just a regular for loop. Any advice on how to do this better or more quickly would be appreciated! I understand that using the parallel version of lapply might be an option, but don't know how to do this either.</p>
  
  <p>All the best,</p>
  
  <p>Alex</p>