在R中优化对向量的每个累积子集进行计算

Question

在R中优化对向量的每个累积子集进行计算

rloopsoptimizationvectorization

3

我有一组各种长度的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)

这对小数据集没有问题，但我的实际数据更像是500万个读取结果，长度从约200,000到1不等（较长的读取结果较少），而我希望达到N50为100,000，那么速度就会变慢。

这个例子更接近于现实情况。在我的桌面上需要大约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)

因此，我对任何能够显著优化这个问题的想法、技巧或窍门都感兴趣。看起来这应该是可能的，但是我已经没有任何想法了。

- roblanf

2个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Damiano Fantini · Answer 1

由于您的数据按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.

- F. Privé · Answer 2

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