diff --git a/R/CoverageExperiment.R b/R/CoverageExperiment.R index b6a0268..ef30ef8 100644 --- a/R/CoverageExperiment.R +++ b/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( diff --git a/R/utils.R b/R/utils.R index 0fd6c3b..bafae64 100644 --- a/R/utils.R +++ b/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 } diff --git a/README.md b/README.md index cffcaf0..00ea877 100644 --- a/README.md +++ b/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)) + diff --git a/man/CoverageExperiment.Rd b/man/CoverageExperiment.Rd index 97699d8..75911e3 100644 --- a/man/CoverageExperiment.Rd +++ b/man/CoverageExperiment.Rd @@ -26,7 +26,7 @@ coarsen(x, window, ...) \S4method{CoverageExperiment}{BigWigFileList,GRangesList}( tracks, features, - width, + width = NULL, center = FALSE, scale = FALSE, ignore.strand = TRUE, @@ -47,7 +47,7 @@ coarsen(x, window, ...) \S4method{CoverageExperiment}{list,GRangesList}( tracks, features, - width, + width = NULL, center = FALSE, scale = FALSE, ignore.strand = TRUE, diff --git a/man/reexports.Rd b/man/reexports.Rd index 67be81f..3fa278b 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -19,9 +19,9 @@ These objects are imported from other packages. Follow the links below to see their documentation. \describe{ - \item{dplyr}{\code{\link[dplyr:reexports]{as_tibble}}} + \item{dplyr}{\code{\link[dplyr]{as_tibble}}} - \item{S4Vectors}{\code{\link[S4Vectors:aggregate-methods]{aggregate}}} + \item{S4Vectors}{\code{\link[S4Vectors]{aggregate}}} \item{tidyr}{\code{\link[tidyr]{expand}}} }} diff --git a/tests/testthat/test-compute-filter-import-plot.R b/tests/testthat/test-compute-filter-import-plot.R index e0f1457..a860d7b 100644 --- a/tests/testthat/test-compute-filter-import-plot.R +++ b/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)