Optimising a calculation on every cumulative subset of a vector in R

user_1084470 0
Viewed 56

I have a collection of DNA sequencing reads of various lengths, sorted from longest to shortest. I would like to know the largest number of reads I can include in a set such that the N50 of that set is above some threshold t

For any given set of reads, the total amount of data is just the cumulative sum of the lengths of the reads. The N50 is defined as the length of the read such that half of the data are contained in reads at least that long.

I have a solution below, but it is slow for very large read sets. I tried vectorising it, but this was slower (probably because my threshold is usually relatively large, such that my solution below stops calculating fairly early on).

Here's a worked example:

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)

This works fine on small datasets, but my actual data are more like 5m reads that vary from ~200,000 to 1 in length (longer reads are rarer), and I'm interested in an N50 of 100,000, then it gets pretty slow.

This example is closer to something that's realistic. It takes ~15s on my desktop.

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)

So, I'm interested in any ideas, tips, or tricks to noticeably optimise this. It seems like this should be possible, but I'm out of ideas.

vectorization r optimization loops
2 Answers

As n is decreasing with i, you should use a 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))
  }
}

Then, just call

binSearch(1, nrow(df))
user_61030400
Related