5

我正在尝试在 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 循环中，但是我笨拙地尝试使用 *apply 函数对其进行矢量化导致代码运行速度要慢一个数量级。

我想我可以用 C 写一些东西，但如果可能的话，我想避免这种情况。任何人都可以提出另一种可以加快计算速度的方法吗？

4

5 回答 5

6

“正确”的做法是使用 bioconductorIRanges包，它使用 IntervalTree 数据结构来表示这些范围。

将两个对象都放在自己的IRanges对象中，然后您将使用该findOverlaps功能获胜。

在这里获取：

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

顺便说一句，包的内部是用 C 编写的，所以它超级快。

编辑

再想一想，它并不像我建议的那样扣篮（单线），但如果你正在使用基因组间隔（或其他类型），你绝对应该开始使用这个库......你可能需要做一些设置操作和东西。对不起，虽然没有时间提供确切的答案。

我只是认为向您指出这个库很重要。

于 2010-03-25T20:25:44.450 回答
4

所以我不完全确定为什么第三个和第四个窗口不是 100 和 20，因为这对我来说更有意义。这是该行为的一个衬里：

Reduce('+', lapply(mylist, function(x) hist(x[1]:x[2], breaks = (0:6) * 100, plot = F)$counts))

请注意，您需要在 中指定上限breaks，但如果您事先不知道，应该不难再通过一次来获得它。

于 2010-03-25T20:24:25.097 回答
4

好的，所以我在这方面浪费了太多时间，但仍然只得到了 3 倍的加速。谁能打败这个？

编码：

my <- do.call(rbind,mylist)
myFloor <- floor(my/100)
myRem <- my%%100
#Add intervals, over counting interval endpoints
counts <- table(do.call(c,apply(myFloor,1,function(r) r[1]:r[2])))
windows[as.numeric(names(counts))+1] <- counts*101

#subtract off lower and upper endpoints
lowerUncovered <- tapply(myRem[,1],myFloor[,1],sum)
windows[as.numeric(names(lowerUncovered))+1]  <-  windows[as.numeric(names(lowerUncovered))+1]  - lowerUncovered
upperUncovered <- tapply(myRem[,2],myFloor[,2],function(x) 100*length(x) - sum(x))
windows[as.numeric(names(upperUncovered))+1]  <-  windows[as.numeric(names(upperUncovered))+1] - upperUncovered

考试：

mylist = vector("list")
for(i in 1:20000){
    d <- round(runif(1,,500))
    mylist[[i]] <- c(d,d+round(runif(1,,700)))
}

windows <- numeric(200)


new_code <-function(){
    my <- do.call(rbind,mylist)
    myFloor <- floor(my/100)
    myRem <- my%%100
    counts <- table(do.call(c,apply(myFloor,1,function(r) r[1]:r[2])))
    windows[as.numeric(names(counts))+1] <- counts*101

    lowerUncovered <- tapply(myRem[,1],myFloor[,1],sum)
    windows[as.numeric(names(lowerUncovered))+1]  <-  windows[as.numeric(names(lowerUncovered))+1]  - lowerUncovered

    upperUncovered <- tapply(myRem[,2],myFloor[,2],function(x) 100*length(x) - sum(x))
    windows[as.numeric(names(upperUncovered))+1]  <-  windows[as.numeric(names(upperUncovered))+1] - upperUncovered

    #print(windows)
}


#old code
old_code <- function(){
    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)
}

system.time(old_code())
system.time(new_code())

结果：

> system.time(old_code())
   user  system elapsed 
  2.403   0.021   2.183 
> system.time(new_code())
   user  system elapsed 
  0.739   0.033   0.588

很郁闷，系统时间基本都是0，但是观察到的时间却是那么的大。我敢打赌，如果你真的降到 C，你会得到 50-100 倍的加速。

于 2010-03-26T08:01:45.030 回答
1

我想我让它变得更加复杂...... System.time 并没有帮助我在这么小的数据集中进行性能评估。 

windows <- numeric(6)

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


library(plyr)

l_ply(mylist, function(x) {
sapply((floor(x[1]/100)+1) : (floor(x[2]/100)+1), function(z){
    eval.parent(parse(text=paste("windows[",z,"] <- ", 
        min(z*100, x[2]) - max((z-1)*100, x[1]) + 1,sep="")),sys.nframe())
    })          
})

print(windows)

编辑

消除的修改eval

g <- llply(mylist, function(x) {
ldply((floor(x[1]/100)+1) : (floor(x[2]/100)+1), function(z){
        t(matrix(c(z,min(z*100, x[2]) - max((z-1)*100, x[1]) + 1),nrow=2))
    })          
})

for(i in 1:length(g)){
    windows[unlist(g[[i]][1])] <- unlist(g[[i]][2])
}
于 2010-03-25T18:52:17.050 回答
0

我没有一个好主意，但是您可以摆脱内部循环，并加快速度。请注意，如果一个窗口完全落在 mylist 间隔内，那么您只需将 100 添加到相应的windows元素。所以只有st-th 和sp-th 窗口需要特殊处理。

  windows <- numeric(100)
  for(i in 1:length(mylist)){ 
    win <- mylist[[i]]         # for cleaner code
    st <- floor(win[1]/100)+1 
    sp <- floor(win[2]/100)+1 
    # start and stop are within the same window
    if (sp == st){
      windows[st] <- windows[st] + (win[2]%%100) - (win[1]%%100) +1 
    }
    # start and stop are in separate windows - take care of edges
    if (sp > st){
      windows[st] <- windows[st] + 100 - (win[1]%%100) + 1
      windows[sp] <- windows[sp] + (win[2]%%100)
    }
    # windows completely inside win
    if (sp > st+1){
      windows[(st+1):(sp-1)] <- windows[(st+1):(sp-1)] + 100
    }       
  }

我生成了一个更大的列表：

  cuts <- sort(sample(1:10000, 70))  # random interval endpoints
  mylist <- split(cuts, gl(35,2))

此版本的 1000 次重复得到 1.08 秒，而原始版本的 1000 次重复得到 1.72 秒。对于真实数据，加速将取决于间隔是否mylist往往比 100 长得多。

顺便说一句，可以将内部循环重写为一个单独的函数，然后将lapply其 over mylist，但这并不能使它更快地工作。

于 2010-03-25T19:17:05.393 回答