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

Question

我收集了各种长度的 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

随着 n 随 i 的增加而减少，您应该使用 binary search algorithm .

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))