title: Speeding up a large for loop using a parallel foreach
body: <p>I'm working on a project where I iterate over bins of the mouse genome and compute some overlaps using <code>GenomicRanges</code>/<code>rtracklayer</code>.  Bins are already computed in 150 base/position increments. For those who are not familiar, this means that for chromosomeY, which happens to be one of the shortest mouse chromosomes, there are approximately 106,000 bins and the larger chromosomes contain about 1,300,000 bins! <strong>The goal here is to iterate through the bins, expand the range of the bins 100,000 positions in both directions, then figure out which genes overlap these window. I want to compute the distance from the center of the bin to the start of the genes contained within the expanded bin window.</strong> </p>
  
  <p>Okay so here's the code I've written to this point. It runs without error and computes exactly what I need. <strong>The problem here is that it is slow and will take forever to iterate through 1M+ bins</strong>. </p>
  
  <pre><code>progress &lt;- function(w) {
    # just a function that prints out the window being processed
    cat(sprintf(paste0(as.character(Sys.time()), ": Window ", w, " completed!\n")))
  }
  
  extend &lt;- function(x, upstream=0, downstream=0) {
  # this will expand a `GenomicRanges` object range 
    if (any(strand(x) == "*"))
      warning("'*' ranges were treated as '+'")
    on_plus &lt;- strand(x) == "+" | strand(x) == "*"
    new_start &lt;- start(x) - ifelse(on_plus, upstream, downstream)
    new_end &lt;- end(x) + ifelse(on_plus, downstream, upstream)
    ranges(x) &lt;- IRanges(new_start, new_end)
    trim(x)
  }
  
  
  feature.overlap &lt;- function(x, window, genes, extend.upstream=100000, extend.downstream=100000) {
  
    # # test case
    # x = chrY; window = 2668 ; genes = gene; extend.upstream = 100000 ; extend.downstream = 100000
  
    # extend window of signal in both directions
    x.window = extend(x[window], extend.upstream, extend.downstream) 
    names(x.window) &lt;- window
  
    # compute signal window overlap with genes 
    overlaps &lt;- subsetByOverlaps(genes, x.window)  
  
    if(length(overlaps) == 0){
  
      values &lt;- data.frame(signal_window=names(x.window), 
                           signal_start=max(0, start(x.window)), 
                           signal_center=max(0, start(x.window)) + floor((width(x.window) - 1)/2), 
                           signal_end=end(x.window), 
                           signal_score=x.window$score,
                           symbol=NA, 
                           gene_id=NA,
                           gene_chr=NA,
                           gene_start=NA,
                           gene_end=NA,
                           gene_strand=NA)  
  
    } else {
  
      hits &lt;- findOverlaps(x.window, genes)
      s.idx &lt;- unique(subjectHits(hits))
      q.idx &lt;- unique(queryHits(hits))    
  
      values &lt;- data.frame(signal_window=names(x.window)[q.idx], 
                           signal_start=max(0, start(x.window)[q.idx]), 
                           signal_center=max(0, start(x.window)[q.idx]) + floor((width(x.window)[q.idx] - 1)/2), 
                           signal_end=end(x.window)[q.idx], 
                           signal_score=x.window$score[q.idx],
                           mcols(overlaps)[,c(2,1)], 
                           gene_chr=chrom(genes)[s.idx],
                           gene_start=ifelse(strand(genes)[s.idx] == '+', start(genes)[s.idx], end(genes)[s.idx]) ,
                           gene_end=end(genes)[s.idx],
                           gene_strand=strand(genes)[s.idx])
  
    }
  
    return(values)
  
  }
  
  # Import data
  library(rtracklayer)
  merged_wig &lt;- import.wig('~/file/linked/below.wig', format='wig', genome='mm9')
  merged_wig &lt;- keepSeqlevels(merged_wig, paste0('chr', c(seq(1,19), 'X', 'Y')))
  chrY &lt;- merged_wig[seqnames(merged_wig) == 'chrY'] 
  
  # Generate gene info needed for computing overlap
  library(TxDb.Mmusculus.UCSC.mm9.knownGene); library(Mus.musculus)
  gene &lt;- genes(TxDb.Mmusculus.UCSC.mm9.knownGene)
  values(gene) &lt;- merge(values(gene), as.data.frame(org.Mm.egSYMBOL), by='gene_id', all.x=T)
  gene &lt;- keepSeqlevels(gene, paste0('chr', c(seq(1,19), 'X', 'Y')))
  
  # BEGIN LOOP GENOME WINDOWS *** TIME CONSUMING ***
  window.overlaps &lt;- list()
  ptm &lt;- proc.time()
  for(i in 1:100) { # ideally 1:length(chrY) but this takes very long so I've only posted a few windows
    result = feature.overlap(chrY, i, gene, extend.upstream=100000, extend.downstream=100000)
    window.overlaps[[i]] &lt;- result
    progress(i)
  }
  proc.time() - ptm
  all.overlaps = do.call(rbind, window.overlaps)
  </code></pre>
  
  <p>The code above will run right our of the box with <a href="https://www.dropbox.com/s/xt4ay8rb77bq8ob/MEFES_K27AC.downsampled.sorted.wig?dl=0" rel="nofollow noreferrer"><strong>this</strong></a> file (88mb).   </p>
  
  <p>Here is my attempt to speed up the look using a <code>foreach</code> with the <code>doParallel</code> library:</p>
  
  <pre><code>library(foreach)
  library(doParallel)
  cl&lt;-makeCluster(8)
  registerDoParallel(cl)
  ptm &lt;- proc.time()
  ls&lt;-foreach(i = 1:100, chrY=chrY, gene=gene, .packages=c('rtracklayer', 'GenomicRanges')) %dopar% {
    result = feature.overlap(chrY, i, gene, extend.upstream=100000, extend.downstream=100000)
    progress(i)
    result
  }
  proc.time() - ptm
  stopCluster(cl)
  </code></pre>
  
  <p>But, this codes does not work. The error returned is <code>Error: this S4 class is not subsettable</code> and no output is produced from <code>progress()</code>. <strong>ERROR FIXED - SEE EDIT</strong> </p>
  
  <p>Again, the goal here is to write this in a more efficient manner. Once I have <code>values</code> I can easily compute the metrics I need. </p>
  
  <p>Any help would be much appreciated! Thank you! </p>
  
  <p><strong>EDIT</strong>: I've implement a working foreach loop with dopar, but it seems to be even slower than the above implementation. </p>
  
  <pre><code>library(foreach)
  library(doParallel)
  cl&lt;-makeCluster(8)
  registerDoParallel(cl)
  ptm &lt;- proc.time()
  ls &lt;- foreach(i = 1:100, .combine='rbind', .packages=c('rtracklayer', 'GenomicRanges')) %dopar% { 
    result = feature.overlap(chrY, i, gene, counts, extend.upstream=100000, extend.downstream=100000)
    progress(i)
    result
  }
  proc.time() - ptm
  stopCluster(cl)
  </code></pre>
  
  <p>For 100 windows this takes about 10 seconds whereas the same windows being processing with the for loop above takes 6 seconds.</p>