R优化:在这种情况下如何避免for循环？

Question

我正在尝试在R中进行简单的基因组轨道交叉,并遇到主要的性能问题,可能与我使用for循环有关.

在这种情况下,我有100bp的预定义窗口,我试图计算mylist中的注释覆盖了多少窗口.从图形上看,它看起来像这样:

          0    100   200    300    400   500   600  
windows: |-----|-----|-----|-----|-----|-----|

mylist:    |-|   |-----------|

所以我写了一些代码来做到这一点,但它相当慢,并成为我的代码的瓶颈:

##window for each 100-bp segment    
windows <- numeric(6)

##second track
mylist = vector("list")
mylist[[1]] = c(1,20)
mylist[[2]] = c(120,320)


##do the intersection
for(i in 1:length(mylist)){
  st <- floor(mylist[[i]][1]/100)+1
  sp <- floor(mylist[[i]][2]/100)+1
  for(j in st:sp){       
    b <- max((j-1)*100, mylist[[i]][1])
    e <- min(j*100, mylist[[i]][2])
    windows[j] <- windows[j] + e - b + 1
  }
}

print(windows)
[1]  20  81 101  21   0   0

当然,这用于比我在此提供的示例大得多的数据集.通过一些剖析,我可以看到的瓶颈是在for循环,但我的笨拙尝试使用向量化它*应用功能,导致更多的运行速度慢一个数量级的代码.

我想我可以在C中写一些东西,但如果可能的话我想避免这样做.任何人都可以提出另一种方法来加速这种计算吗？

Answer 1

"正确"的做法是使用bioconductor IRanges包,它使用IntervalTree数据结构来表示这些范围.

将两个对象放在自己的IRanges对象中,然后使用该findOverlaps函数来获胜.

在这里得到它:

http://www.bioconductor.org/packages/release/bioc/html/IRanges.html

通过by,包的内部是用C语言编写的,所以它的速度非常快.

编辑

第二个想法,它并不像我建议的那样(一个单线程),但如果你在基因组间隔(或其他类型)工作,你绝对应该开始使用这个库...你可能需要做一些固定操作和东西.抱歉,没有时间提供确切的答案.

我只是认为将这个库指向你是很重要的.