coveragePerTiling(..., as.data.table = TRUE, keep.names = TRUE) strips names

[m-swirski]

Hi, is it a desired behavior? Size of resulting data.table is very similar if you keep names as factor.

library(ORFik) attach(loadNamespace(ORFik)

only change is the if statement for keep.names within as.data.table if statement.

cpt_nm <- function (grl, reads, is.sorted = FALSE, keep.names = TRUE, as.data.table = FALSE, withFrames = FALSE, weight = "score") { if (!is.sorted) grl <- sortPerGroup(grl) if (is.numeric(weight) | (weight[1] %in% colnames(mcols(reads)))) { coverage <- coverageByTranscriptW(reads, grl, weight = weight) } else coverage <- coverageByTranscript(reads, grl) if (!keep.names) names(coverage) <- NULL if (as.data.table) { window_size <- unique(widthPerGroup(grl, FALSE)) count <- data.table(count = unlist(IntegerList(coverage), use.names = FALSE)) if (keep.names){ count[, :=(genes,factor(rep(names(grl), widthPerGroup(grl))))] } else count[, :=(genes, groupings(coverage))] if (length(window_size) != 1) { count[, :=(ones, rep.int(1L, length(genes)))] count[, :=(position, cumsum(ones)), by = genes] count$ones <- NULL } else { count[, :=(position, rep.int(seq.int(window_size), length(coverage)))] } if (withFrames) { count[, :=(frame, (position - 1)%%3)] } count[] return(count) } return(coverage) }

bam_file <- system.file("extdata", "ribo-seq.bam", package = "ORFik") footprints <- readBam(bam_file) gtf_file <- system.file("extdata", "annotations.gtf", package = "ORFik") txdb <- loadTxdb(gtf_file) tx <- exonsBy(txdb, by = "tx", use.names = TRUE)

result1 <- coveragePerTiling(tx, footprints, as.data.table = TRUE) result2 <- cpt_nm(tx, footprints, as.data.table = TRUE, keep.names = FALSE) result3 <- cpt_nm(tx, footprints, as.data.table = TRUE, keep.names = TRUE) identical(result1,result2) object.size(result1) object.size(result3)

[Roleren]

Hey, thanks for these tests, I just got married, so on honeymoon until 15th of August. Will check it out then.

[m-swirski]

Congratulations! Sure, I'll just keep opening new threads and resolving problems on my own, and you'll see to it then.

[Roleren]

I think any subsetting on names should be filtered before you tile the grl, And since you get the gene index, you can easily add it. Do you have a specific use case where it makes sense to keep the names in the tiled table ?

[m-swirski]

Hmm, my rationale was just that keeping names, instead of 'rank' numbers would make it easier to spot mistakes or the fact that genes got shuffled along the way - it can happen pretty easily while having GTF and TxDb based on it, because TxDb sorts by seqnames and strand (I don' really know why):

#pseudocode, just to show how it usually looks
strand(unlist(exons))
factor-Rle of length 5698 with 17 runs
  Lengths: 470 468 476 442 410 387 372 353 340 372 306 299 227 272 223 240  41
  Values :   +   -   +   -   +   -   +   -   +   -   +   -   +   -   +   -   +
Levels(3): + - *
strand(gtf_gr[gtf_gr$type == "exon"])
factor-Rle of length 5698 with 2975 runs
  Lengths:  4  2  2  3  2  1  1  1  1  2  1  2  1  1  1  1  1  1  1  1 ...  1  1  3  1  1  2  2  1  2  5  1  4  2  2  1  1  1  6  3 41
  Values :  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +  - ...  -  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +  -  +
Levels(3): + - *

And I think that in some cases subsetting by after using coveragePerTiling may be helpful, i.e. without as.data.table = TRUE, you just subset after tiling:

#pseudocode again
covs <- coveragePerTiling(grl, reads)
covs["mygene"]

and with data table it would be very similar if names were present:

covs <- coveragePerTiling(grl, reads, as.data.table = TRUE)
covs[genes == "mygene",]

It's a minor issue, of course, I just thought that having the names would be beneficial, given that a factor column is not significantly heavier than integer column.

[Roleren]

Will keep this open, but right now I think it's not worth implementing. Since it might lead to bugs downstream

[Roleren]

I've seen some more compelling reasons to have this as an option, will do some tests on how it affects size & speed

[m-swirski]

I observed no significant difference in performance when having gene names as factor (as in the function in the first post)

[Roleren]

We now allow names as long as drop.zero.dt is TRUE, and this is documented in the function. Still considering allowing it when drop.zero.dt is TRUE.