coverageByTranscript weights ?

[Roleren]

For human data, using human NGS reads (1.4GB) as GRanges:

# Using reduced data that has a score column for duplicated reads,
#  every unique read is only once.
system.time(withWeights <- coverageByTranscript(a, cds, weight = "score"))
   user  system elapsed 
  3.987   0.683   4.667 
# Using raw data, without reducing the data-set by not creating 
# a $score column for number of duplicated reads
> system.time(withoutWeights <- coverageByTranscript(b, cds))
   user  system elapsed 
 37.956  21.586  59.547 
> identical(d, dd)
[1] TRUE

I just changed the coverageByTranscript function to allow weights, maybe a possible improvement?

[m-swirski]

The reason why GenomicFeatures developers didn't allow for weight in coverageByTranscript() is probably coverage method on GrangesList behavior. I elaborate on that in ORFik issues:
Roleren/ORFik#98

[Roleren]

I have this implemented as coverageByTranscriptW in ORFik for GRanges input.
Here

[lshep]

If you would like to make a pull request for the enhancement that would be fantastic.

[Roleren]

It should work to merge directly, I will add some more tests, and also add a note to remove the part:

 if ((is(x, "GAlignments") | is(x, "GAlignmentPairs"))
      & is.character(weight)) {

When coverage(), is fixed for GRangesList.

[lshep]

That would be great and when ready please put in the pull request

[Roleren]

As a note for you guys, I think coverage(GRangesList) should be changed to support meta columns asap, then much of my function is no longer needed.

[Roleren]

Pull request:
#24

[Roleren]

A new update, reading about other people also noting that IRanges::revElements is slow, I checked with profvis, and indeed this is the slow part of transcriptByCoverage:

I added a simple change in ORFik:::transcriptByCoverageW to make it much faster:

library(S4Vectors); library(IRanges)
#' Faster revElements for IRanges::RleList
#' @param x RleList
#' @return RleList (revered list elements)
revElementsF <- function(x) {
  b <- rev(x)
  b@unlistData <- rev(x@unlistData)
  return(rev(b))
}

# Make a large RleList of type integer-Rle
out <- sample(c(0L,1L,5L), size = 100000000, replace = T, prob = c(0.95, 0.4, 0.1))
out <- split(out, sort(sample(seq(length(out)/10), size = length(out), replace = T)))
out <- IRanges::RleList(out)
head(out)
system.time(a <- revElements(out))
system.time(b <- revElementsF(out))
identical(a, b)

# Adding longer names makes revElements even worse
names(out) <- paste0("ENST", sample(seq(100000, 100000 + length(out)), size = length(out), replace = T))
head(names(out))
system.time(a <- revElements(out))
system.time(b <- revElementsF(out))
identical(a, b)