优化R中的每个向量的每个累积子集的计算

Question

i具有各种长度的DNA测序读数的集合，从最长到最短排序。我想知道我可以在一组中包含的最大读取，以使该集合的N50超过某个阈值t

对于任何给定的读取集，数据总数只是读取长度的累积总和。N50定义为读取的长度，使得一半的数据至少包含在读取中。

我有下面的解决方案，但是对于非常大的读取集来说，它很慢。我尝试对其进行矢量化，但这较慢（可能是因为我的阈值通常相对较大，因此下面的解决方案可以很早地停止计算）。

这是一个有效的例子：

df = data.frame(l = 100:1) # read lengths
df$cs = cumsum(df$l) # getting the cumulative sum is easy and quick
t = 95 # let's imagine that this is my threshold N50
for(i in 1:nrow(df)){
    N50 = df$l[min(which(df$cs>df$cs[i]/2))]
    if(N50 < t){ break }
}
# the loop will have gone one too far, so I subtract one
number.of.reads = as.integer(i-1)

这在小型数据集上正常工作，但是我的实际数据更像5m读取的读数从约200,000到1的长度（较长的读数更稀有），我对100,000的N50感兴趣，那么它变得很漂亮慢。

这个示例更接近现实的东西。我的桌面上需要约15秒。

l = ceiling(runif(100000, min = 0, max = 19999))
l = sort(l, decreasing = T)
df = data.frame(l = l)
df$cs = cumsum(df$l)
t = 18000
for(i in 1:nrow(df)){
    n = df$l[min(which(df$cs>df$cs[i]/2))]
    if(n < t){ break }
}
result = as.integer(i-1)

所以，我对任何想法，技巧或技巧感兴趣，以明显优化这一点。似乎这是可能的，但我没有想法。

Answer 1

as n随着 i的减少，您应该使用二进制搜索算法。

binSearch <- function(min, max) {
  print(mid <- floor(mean(c(min, max))))
  if (mid == min) {
    if (df$l[min(which(df$cs>df$cs[min]/2))] < t) {
      return(min - 1)
    } else {
      return(max - 1)
    }
  }
  n = df$l[min(which(df$cs>df$cs[mid]/2))]
  if (n >= t) {
    return(binSearch(mid, max))
  } else {
    return(binSearch(min, mid))
  }
}

然后，只需致电

binSearch(1, nrow(df))

Answer 2

由于您的数据按DNA/读取长度订购，也许您可以避免测试每一行。相反，您可以在每次迭代（例如使用while()）进行有限数量的行（相当间隔），因此可以逐渐接近解决方案。这应该使事情更快。只需确保一旦接近解决方案，就可以停止迭代。

这是您的解决方案

set.seed(111)
l = ceiling(runif(100000, min = 0, max = 19999))
l = sort(l, decreasing = T)
df = data.frame(l = l)
df$cs = cumsum(df$l)
t = 18000
for(i in 1:nrow(df)){
  n = df$l[min(which(df$cs>df$cs[i]/2))]
  if(n < t){ break }
}
result = as.integer(i-1)
result 
# 21216, in ~29 seconds

而不是测试每一行，让我们设置一个范围

i1 <- 1
i2 <- nrow(df)
i.range <- as.integer(seq(i1, i2, length.out = 10))

现在，仅测试这10行。通过重新定义范围来获取最接近的一个并"重点"。当您无法增加粒度时停止。

while(sum(duplicated(i.range))==0){
  for(i in 1:length(i.range)){
    N50 = df$l[min(which(df$cs>df$cs[i.range[i]]/2))]
    if(N50 < t){ break }
  }
  #update i1 and i2
  i1 <- i.range[(i-1)]
  i2 <- i.range[i]
  i.range <- as.integer(seq(i1, i2, length.out = 10))
}
i.range <- seq(i1, i2, by=1)
for(i in i.range){
  N50 = df$l[min(which(df$cs>df$cs[i]/2))]
  if(N50 < t){ break }
}
result <- as.integer(i-1)
result 
#21216, in ~ 0.06 seconds
Same result in a fraction of the time.