feat: NULL width by default, with checks

R/CoverageExperiment.R

@@ -102,7 +102,7 @@ setMethod(
     "CoverageExperiment", 
     signature(tracks = "BigWigFileList", features = "GRangesList"), 
     function(
-        tracks, features, width, 
+        tracks, features, width = NULL, 
         center = FALSE, scale = FALSE, 
         ignore.strand = TRUE, 
         window = 1, 
@@ -115,7 +115,8 @@ setMethod(
         ## Extend and filter features
         tracks <- .set_seqinfo_bwfl(tracks)
         si <- seqinfo(tracks[[1]])
-        features <- lapply(features, .resize_granges, width = width, seqinfo = si) 
+        if (is.null(width)) .check_granges_widths(features)
+        features <- lapply(features, .resize_granges, width, si) 
 
         ## Prepare cData and rData
         cData <- data.frame(
@@ -274,7 +275,7 @@ setMethod(
     "CoverageExperiment", 
     signature(tracks = "list", features = "GRangesList"), 
     function(
-        tracks, features, width, 
+        tracks, features, width = NULL, 
         center = FALSE, scale = FALSE, 
         ignore.strand = TRUE, 
         window = 1, 
@@ -287,7 +288,8 @@ setMethod(
         ## Extend and filter features
         tracks <- .set_seqinfo(tracks)
         si <- seqinfo(tracks[[1]])[[1]]
-        features <- lapply(features, .resize_granges, width = width, seqinfo = si) 
+        if (is.null(width)) .check_granges_widths(features)
+        features <- lapply(features, .resize_granges, width, si) 
 
         ## Prepare cData and rData
         cData <- data.frame(

R/utils.R

@@ -37,14 +37,25 @@
     return(tracks)
 } 
 
-.resize_granges <- function(gr, width, seqinfo) {
+.resize_granges <- function(gr, w, seqinfo) {
     GenomeInfoDb::seqlevels(gr, pruning.mode = "coarse") <- GenomeInfoDb::seqlevels(seqinfo)
     GenomeInfoDb::seqinfo(gr) <- seqinfo
-    gr <- suppressWarnings(resize(gr, fix = 'center', width = width))
-    w <- width
+    if (is.null(w)) {
+        w <- width(gr) |> unlist() |> unique()
+    }
+    else {
+        gr <- suppressWarnings(resize(gr, fix = 'center', width = w))
+    }
     gr <- trim(gr)
     gr <- gr[width(gr) == w]
-    sort(gr)
+    gr
+}
+
+.check_granges_widths <- function(features) {
+    ws <- lapply(features, width) |> unlist() |> unique()
+    if (length(ws) > 1)
+        stop("Input GRanges do not all have the same width.
+  Use resize() to fix width prior to coverage data extraction.")
 }
 
 #' @importFrom IRanges NumericList
@@ -58,7 +69,9 @@
         m[which.flip, ] <- rev(m[which.flip, ])
         m <- as.matrix(m)
     }
-    m <- t(scale(t(m), center = center, scale = scale))
+    if (any(c(center, scale))) {
+        m <- t(scale(t(m), center = center, scale = scale))
+    }
     m
 }
 

README.md

@@ -52,7 +52,7 @@ CE |>
     filter(track %in% c('MNase', 'PolII')) |> 
     filter(features == 'TSSs') |> 
     aggregate() |> 
-    ggplot(aes(x = coord, y = mean)) + 
+    ggplot(aes(x = coord, y = mean, ymin = ci_low, ymax = ci_high, fill = track, col = track)) + 
     geom_ribbon(aes(ymin = ci_low, ymax = ci_high, fill = track), alpha = 0.2) + 
     geom_line(aes(col = track)) + 
     facet_grid(track ~ ., scales = "free") + 
@@ -66,7 +66,7 @@ CE |>
 ## Plot coverage over a single locus
 
 ```r
-CoverageExperiment(tracks, GRanges("II:450001-455000"), width = 5000) |> 
+CoverageExperiment(tracks, GRanges("II:450001-455000")) |> 
     expand() |> 
     ggplot(aes(x = coord, y = coverage)) + 
         geom_col(aes(fill = track, col = track)) + 

tests/testthat/test-compute-filter-import-plot.R

@@ -45,6 +45,7 @@ test_that("other CoverageExperiment methods work", {
     x <- BigWigFile(system.file("extdata", "RNA.fwd.bw", package = "tidyCoverage"))
     y <- import(system.file("extdata", "TSSs.bed", package = "tidyCoverage"))[1:200]
     expect_s4_class(CoverageExperiment(x, y, width = 100), "CoverageExperiment")
+    expect_s4_class(CoverageExperiment(x, y), "CoverageExperiment")
 
     # ~~~~~~~~ BigWigFile + GRangesList
     y <- GRangesList(tss = y, tss2 = y)
@@ -74,6 +75,7 @@ test_that("other CoverageExperiment methods work", {
     x <- import(system.file("extdata", "RNA.fwd.bw", package = "tidyCoverage"), as = 'Rle')
     y <- import(system.file("extdata", "TSSs.bed", package = "tidyCoverage"))
     expect_s4_class(CoverageExperiment(x, y, width = 100), "CoverageExperiment")
+    expect_s4_class(CoverageExperiment(x, y), "CoverageExperiment")
 
     # ~~~~~~~~ RleList + GRangesList
     y <- GRangesList(tss = y)