Remove Peaks

[sgibb]

This is a work-in-progress PR.
After IRanges were replaced by IRangesList in utils.removePeaks in d519ec5 I recognized a large slow down of this function. So I started to rewrite this function without IRanges (in favour of #201 I removed IRanges dependency completely, avoiding IRanges is faster than the recently merged solution by @hpages but uses a for loop again, see below).

library("MSnbase")
library("IRanges")
library("rbenchmark")


utils.removePeaksIRangesList <- function(int, t) {
    peakRanges <- as(IRanges(int > 0L), "IRangesList")
    sapply(peakRanges, function(x) {
        idx <- start(x):end(x)
        ## we get the indices of every peak in int
        if(all(int[idx]<=t))
            int[idx] <<- 0
        })
    int
}

utils.removePeaksIRanges <- function(int, t) {
    peakRanges <- IRanges(int > 0L)
    s <- start(peakRanges)
    e <- end(peakRanges)

    for (i in seq_along(s)) {
        idx <- s[[i]]:e[[i]]
        if (all(int[idx] <= t)) {
            int[idx] <- 0L
        }
    }
    int
}

utils.removePeaksIRangesOverlap <- function(int, t) {
    peakRanges <- IRanges(int > 0L)
    aboveThresh <- IRanges(int > t)
    belowThresh <- peakRanges %outside% aboveThresh
    s <- start(peakRanges[belowThresh])
    e <- end(peakRanges[belowThresh])
    int[unlist(mapply(seq, s, e))] <- 0L
    int
}

utils.removePeaksIRangesExtractList <- function(int, t) {
    peakRanges <- as(int > 0, "IRanges")
    peak_is_too_low <- max(extractList(int, peakRanges)) <= t
    replaceROWS(int, peakRanges[peak_is_too_low], 0L)
}

utils.removePeaksDiff <- function(int, t) {
    d <- diff(int > 0L)
    s <- which(d == 1L)
    e <- which(d == -1L)
    for (i in seq_along(s)) {
        j <- s[[i]]:e[[i]]
        if (all(int[j] <= t)) {
            int[j] <- 0L
        }
    }
    int
}

data(itraqdata)

int <- intensity(itraqdata)
thr <- 2e5


all.equal(utils.removePeaksIRangesList(int[[1]], thr),
          utils.removePeaksIRanges(int[[1]], thr))
# [1] TRUE
all.equal(utils.removePeaksIRanges(int[[1]], thr),
          utils.removePeaksIRangesOverlap(int[[1]], thr))
# [1] TRUE
all.equal(utils.removePeaksIRangesOverlap(int[[1]], thr),
          utils.removePeaksIRangesExtractList(int[[1]], thr))
# [1] TRUE
all.equal(utils.removePeaksIRangesExtractList(int[[1]], thr),
          utils.removePeaksDiff(int[[1]], thr))
# [1] TRUE

b <- lapply(int[1:5], function(i) {
    benchmark(IRangesList=utils.removePeaksIRangesList(i, thr),
              IRanges=utils.removePeaksIRanges(i, thr),
              Overlap=utils.removePeaksIRangesOverlap(i, thr),
              ExtractList=utils.removePeaksIRangesExtractList(i, thr),
              diff=utils.removePeaksDiff(i, thr),
              columns=c("test", "elapsed", "relative"),
              replications=10,
              order="relative")
})
b
# $X1
#          test elapsed relative
# 5        diff   0.002      1.0
# 2     IRanges   0.004      2.0
# 3     Overlap   0.028     14.0
# 4 ExtractList   0.029     14.5
# 1 IRangesList   0.937    468.5
# 
# $X10
#          test elapsed relative
# 5        diff   0.002      1.0
# 2     IRanges   0.004      2.0
# 3     Overlap   0.026     13.0
# 4 ExtractList   0.028     14.0
# 1 IRangesList   0.833    416.5
# 
# $X11
#          test elapsed relative
# 5        diff   0.002      1.0
# 2     IRanges   0.004      2.0
# 3     Overlap   0.025     12.5
# 4 ExtractList   0.029     14.5
# 1 IRangesList   0.781    390.5
# 
# $X12
#          test elapsed relative
# 5        diff   0.002      1.0
# 2     IRanges   0.005      2.5
# 3     Overlap   0.027     13.5
# 4 ExtractList   0.028     14.0
# 1 IRangesList   1.162    581.0
# 
# $X13
#          test elapsed relative
# 5        diff   0.003    1.000
# 2     IRanges   0.004    1.333
# 3     Overlap   0.029    9.667
# 4 ExtractList   0.029    9.667
# 1 IRangesList   1.264  421.333

I also wanted to address the issue #252 with this PR but unfortunately I won't have time the next days (maybe weeks) to finish this PR. If anybody wants to take over please feel free to do so.

[lgatto]

That's great, thanks. I won't have time either, but no hurry. Let me know if you prefer this PR to be merged now wait for #252.

[jorainer]

Nice @sgibb 👍

[sgibb]

@lgatto Because I won't be able to finish this in the near future and the current PR is an (non-destructive) improvement I merge it. But I am afraid nobody will work on #252 then ...

[lgatto]

Thanks @sgibb - I won't be able to commit time to #252 at the moment. Hopefully later though.

[hpages]

Be aware that utils.removePeaksDiff() doesn't exactly preserve the semantic of the original utils.removePeaks():

## Original utils.removePeaks() (slightly modified to work with BioC 3.7):
utils.removePeaks <- function(int, t) {
  peakRanges <- IRanges(sapply(int,">",0))
  sapply(seq_along(peakRanges),function(i) {
    x <- IRanges:::unlist_as_integer(peakRanges[i])
    ## we get the indices of every peak in int
    if(all(int[x]<=t))
      int[x] <<- 0
  })
  return(int)
}
> int1 <- c(2.5, -0.3, 0, 3.1, 4.2, -1, 2.2, 0)
> utils.removePeaks(int1, 2.5)
[1]  0.0 -0.3  0.0  3.1  4.2 -1.0  0.0  0.0
> utils.removePeaksDiff(int1, 2.5)
[1]  0.0  0.0  0.0  3.1  4.2 -1.0  2.2  0.0

Also it sometimes fails (the original utils.removePeaks() would fail only when the input is a zero-length numeric vector):

> int2 <- c(0, 0, 2.5, -0.3, 3.1, 4.2, -1, 2.2)
> utils.removePeaks(int2, 0.1)
[1]  0.0  0.0  2.5 -0.3  3.1  4.2 -1.0  2.2
> utils.removePeaksDiff(int2, 0.1)
Error in e[[i]] : subscript out of bounds

Note that with utils.removePeaksIRangesExtractList() you pay a fixed overhead for using S4 objects. That overhead penalizes any S4-based solution vs a non-S4-based solution when the input is small to medium-size. However it becomes neglectable when the input is big:

> system.time(res3a <- utils.removePeaksDiff(int3, 0.75))
   user  system elapsed 
 10.148   1.036  11.184 
> system.time(res3b <- utils.removePeaksIRangesExtractList(int3, 0.75))
   user  system elapsed 
  1.998   0.564   2.563 

H.

[lgatto]

Thanks @hpages for following up!

[sgibb]

@hpages thanks for this great comment! Indeed the utils.removePeaksDiff implementation contains a bug. A fixed version would be

utils.removePeaks <- function(int, t) {
    int <- c(0L, int, 0L)
    d <- diff(int > 0L)
    s <- which(d == 1L) + 1L
    e <- which(d == -1L)
    for (i in seq_along(s)) {
        j <- s[[i]]:e[[i]]
        if (all(int[j] <= t)) {
            int[j] <- 0L
        }
    }
    int[seq_len(length(int) - 2L) + 1L]
}

But you are right that the diff/for loop concept is just faster because it doesn't need to do S4-method-dispatch. I revoke all my chances (because of all the merge from bioc/github etc. it is easier to add another commit: ) and reapply your IRanges based implementation (but I removed the snake_case_variable_name 😉). Thanks for enlighten me.

@lgatto and @jotsetung sorry for this stupid PR and all its inconvenience (and I said it would be a non-destructive PR ...)