使用并行 foreach 加速大型 for 循环

Question

我正在开展一个项目，在该项目中迭代小鼠基因组的 bin 并使用 GenomicRanges/rtracklayer 计算一些重叠。 bin 已经以 150 个基数/位置增量计算。对于那些不熟悉的人来说，这意味着对于恰好是最短的小鼠染色体之一的 Y 染色体，大约有 106,000 个 bin，而较大的染色体则包含大约 1,300,000 个 bin！ 这里的目标是遍历这些 bin，在两个方向上将 bin 的范围扩大 100,000 个位置，然后找出哪些基因与这些窗口重叠。我想计算从 bin 中心到扩展 bin 窗口中包含的基因开始的距离。

好的，这就是我写到这里的代码。它运行没有错误，并准确计算出我需要的东西。 这里的问题是它很慢，并且需要很长时间才能遍历 1M+ 个 bin。

progress <- function(w) {
  # just a function that prints out the window being processed
  cat(sprintf(paste0(as.character(Sys.time()), ": Window ", w, " completed!\n")))
}

extend <- function(x, upstream=0, downstream=0) {
# this will expand a `GenomicRanges` object range 
  if (any(strand(x) == "*"))
    warning("'*' ranges were treated as '+'")
  on_plus <- strand(x) == "+" | strand(x) == "*"
  new_start <- start(x) - ifelse(on_plus, upstream, downstream)
  new_end <- end(x) + ifelse(on_plus, downstream, upstream)
  ranges(x) <- IRanges(new_start, new_end)
  trim(x)
}


feature.overlap <- 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) <- window

  # compute signal window overlap with genes 
  overlaps <- subsetByOverlaps(genes, x.window)  

  if(length(overlaps) == 0){

    values <- 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 <- findOverlaps(x.window, genes)
    s.idx <- unique(subjectHits(hits))
    q.idx <- unique(queryHits(hits))    

    values <- 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 <- import.wig('~/file/linked/below.wig', format='wig', genome='mm9')
merged_wig <- keepSeqlevels(merged_wig, paste0('chr', c(seq(1,19), 'X', 'Y')))
chrY <- merged_wig[seqnames(merged_wig) == 'chrY'] 

# Generate gene info needed for computing overlap
library(TxDb.Mmusculus.UCSC.mm9.knownGene); library(Mus.musculus)
gene <- genes(TxDb.Mmusculus.UCSC.mm9.knownGene)
values(gene) <- merge(values(gene), as.data.frame(org.Mm.egSYMBOL), by='gene_id', all.x=T)
gene <- keepSeqlevels(gene, paste0('chr', c(seq(1,19), 'X', 'Y')))

# BEGIN LOOP GENOME WINDOWS *** TIME CONSUMING ***
window.overlaps <- list()
ptm <- 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]] <- result
  progress(i)
}
proc.time() - ptm
all.overlaps = do.call(rbind, window.overlaps)

上面的代码将在带有this 文件（88mb）的盒子中运行。

这是我尝试使用 foreach 和 doParallel 库来加速外观：

library(foreach)
library(doParallel)
cl<-makeCluster(8)
registerDoParallel(cl)
ptm <- proc.time()
ls<-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)

但是，此代码不起作用。返回的错误是Error: this S4 class is not subsettable，并且没有从progress() 产生输出。 已修复错误 - 请参阅编辑

同样，这里的目标是以更有效的方式编写此代码。一旦我有了values，我就可以轻松计算出我需要的指标。

任何帮助将不胜感激！谢谢！

编辑：我已经用 dopar 实现了一个有效的 foreach 循环，但它似乎比上面的实现还要慢。

library(foreach)
library(doParallel)
cl<-makeCluster(8)
registerDoParallel(cl)
ptm <- proc.time()
ls <- 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)

对于 100 个窗口，这大约需要 10 秒，而使用上述 for 循环处理的相同窗口需要 6 秒。

Answer 1

好的，我设法运行了你的代码。

• 对于 100 次运行，需要
  • 6 秒，ncores = 1
  • 7 秒，ncores = 2
  • 8 秒，ncores = 3
  • 10 秒，ncores = 4
• 对于 1000 次运行，需要
  • 56 秒顺序运行
  • 2 核 33 秒
  • 3 核 33 秒

所以，基本上，您看到的 7 秒是由于并行化的开销造成的。所以你做的越多，并行化就越有用。

然后，我们可以在这里验证的是，使用计算机上一半以上的内核通常是低效的。所以，如果你有 8 个，你应该只使用 4 个。

最后，如果你有一个很长的循环，不要使用.combine 参数，因为它会占用 100 个块并将它们组合起来，这很慢。如果可以，请在最后输入Reduce。