title: Using %dopar% with a custom function
body: <p>So I've got this function meant to group measurements from multiple probes that fall into defined regions.  </p>
  
  <pre><code>HMkit.dmr&lt;-function(Mat,Classes,method.fdr=c("BH","bonferroni"),probe.features) {
  
     #Annotate first...
     require(plyr)
     require(dplyr)
  
  
     #Filter matrix for testing and stuff...
  
     message("Setting up merged table")
     Mat2&lt;-Mat[match(probe.features$probe,rownames(Mat)),]
  
     #Split by classes
  
     if(!is.factor(Classes)) {
         Classes&lt;-as.factor(Classes)
     }
     Class.1&lt;-levels(Classes)[[1]]
     Class.2&lt;-levels(Classes)[[2]]
  
     C1.Mat&lt;-Mat2[,Classes==Class.1]
     C2.Mat&lt;-Mat2[,Classes==Class.2]
  
     #Summarise and run wilcoxon's test for each dmr...
     num.regions&lt;-length(unique(as.character(probe.features$region.id)))
     pvals.vec&lt;-numeric(length=num.regions)
     unique.regions&lt;-unique(as.character(probe.features$region.id))
     message(num.regions)
     Meds.1&lt;-numeric(length=num.regions);Meds.2&lt;-numeric(length=num.regions)
  
     for (i in 1:num.regions) {
         region&lt;-probe.features%&gt;%filter(region.id %in% unique.regions[[i]])
         Set1.Mat&lt;-as.numeric(C1.Mat[rownames(C1.Mat) %in% region$probe,])
         Set2.Mat&lt;-as.numeric(C2.Mat[rownames(C2.Mat) %in% region$probe,])
         pvals.vec[[i]]&lt;-wilcox.test(Set1.Mat,Set2.Mat)$p.value
         Meds.1[[i]]&lt;-median(Set1.Mat)
         Meds.2[[i]]&lt;-median(Set2.Mat)
         message(i)
     }
  
      #Output frame
      dmrs.frame&lt;-data.frame(region=unique.regions,pval=pvals.vec,G1=Meds.1,G2=Meds.2,dB=Meds.1-Meds.2)
      dmrs.frame$q.val&lt;-p.adjust(dmrs.frame$pval,method=method.fdr)
      groups.ids&lt;-levels(Classes)
      return(list(dmrs=dmrs.frame,groups=groups.ids))
  }
  </code></pre>
  
  <p>The code basically splits a matrix into two groups by samples and then pulls in the values of all probes that are defined as being in a region, calls a wilcox.test and a median summarisation step and uses it to populate vectors created beforehand. </p>
  
  <p>I have tried to replace the for in the for loop with doparallel function in the foreach package but have not been able to get it to populate the vector with the correct outcomes. I want to know how to correctly use parallelisation with the function above - either by modifying the for loop, or by modifying the function call so regions are broken down into chunks that are processed in parallel. </p>
  
  <p>Example objects follow below...</p>
  
  <pre><code>Mat&lt;-matrix(runif(200,0,1), ncol=10,nrow=20)
  rownames(Mat)&lt;-paste0("p",1:20)
  colnames(Mat)&lt;-paste0("S",1:10)
  
  Classes&lt;-as.character(c(rep("G1",6),rep("G2",4)))
  probe.features&lt;-data.frame(probe=paste0("p",1:20),region.id=c(rep("R1",5),rep("R2",3),rep("R3",4),rep("R5",4),rep("R6",4))
  </code></pre>
  
  <p>and the function is run using </p>
  
  <pre><code>x&lt;-HMkit.dmr(Mat,Classes,method.fdr=c("BH"),probe.features=probe.features)
  </code></pre>
  
  <p>In practise, there are 30,000 regions I am looking at, and want to parallelise the function across multiple cores on windows because serial execution can take up to 40 minutes. How do I do this? </p>
  
  <p>Addendum - I tried to do this with </p>
  
  <pre><code> library(doParallel)
   ncores&lt;-2
   Cl&lt;-makeCluster(2)
   registerDoParallel(Cl)
  x&lt;-foreach(i=1:length(unique(probe.features$region.id)),packages=c("plyr","dplyr"))%dopar%HMkit.dmr(Mat,Classes,probe.features=probe.features,method.fdr="BH")
  </code></pre>
  
  <p>However, doing that just returned two copies of the same results as the serial function, what I want it to do is break down regions in probe.features$region.id into chunks that go to different cores. </p>