From c8503b501e615d785ceec92904dc82a10c7f2128 Mon Sep 17 00:00:00 2001 From: Lucille Delisle Date: Thu, 11 Jun 2026 17:24:49 +0800 Subject: [PATCH 1/8] Add ExportGroupBW to export per-group bigwig files (#1887) Export coverage tracks as bigwig files, one per group of cells, adapted from the ArchR approach. Ported to the v2 branch: - replace SetIfNull (removed in v2) with %||% - add an explicit seqlengths argument, since v2 classes no longer store chromosome lengths; falls back to the object seqinfo when available and errors clearly when lengths are unavailable - drop the ineffective is.null(chromLengths) guard in favour of an NA check Co-authored-by: Tim Stuart <4591688+timoast@users.noreply.github.com> Co-Authored-By: Claude Fable 5 --- DESCRIPTION | 3 +- NAMESPACE | 2 + R/bigwig.R | 304 +++++++++++++++++++++++++++++++++++ man/ExportGroupBW.Rd | 72 +++++++++ tests/testthat/test-bigwig.R | 174 ++++++++++++++++++++ 5 files changed, 554 insertions(+), 1 deletion(-) create mode 100644 R/bigwig.R create mode 100644 man/ExportGroupBW.Rd create mode 100644 tests/testthat/test-bigwig.R diff --git a/DESCRIPTION b/DESCRIPTION index 68dc46f3..0eeebe0d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -65,10 +65,11 @@ Imports: Collate: 'RcppExports.R' 'Signac-package.R' + 'generics.R' + 'bigwig.R' 'classes.R' 'data.R' 'differential_accessibility.R' - 'generics.R' 'dimension_reduction.R' 'footprinting.R' 'fragments.R' diff --git a/NAMESPACE b/NAMESPACE index 0518296f..72f151ac 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -172,6 +172,7 @@ export(DensityScatter) export(DepthCor) export(DownsampleFeatures) export(EnrichedTerms) +export(ExportGroupBW) export(ExpressionPlot) export(Extend) export(FRiP) @@ -296,6 +297,7 @@ importFrom(GenomicRanges,makeGRangesFromDataFrame) importFrom(GenomicRanges,reduce) importFrom(GenomicRanges,resize) importFrom(GenomicRanges,seqnames) +importFrom(GenomicRanges,slidingWindows) importFrom(GenomicRanges,start) importFrom(GenomicRanges,strand) importFrom(GenomicRanges,tileGenome) diff --git a/R/bigwig.R b/R/bigwig.R new file mode 100644 index 00000000..3119e306 --- /dev/null +++ b/R/bigwig.R @@ -0,0 +1,304 @@ +#' @include generics.R +#' +NULL + +#' Export bigwig files for groups of cells +#' +#' Export coverage tracks as bigwig files, one per group of cells. Fragments +#' are split by group, the genome is tiled, and the number of fragments in each +#' tile is counted and (optionally) normalized before writing the bigwig files. +#' +#' @param object A Seurat object +#' @param assay Name of assay to use +#' @param group.by The metadata variable used to group the cells +#' @param idents Identities to include (defined by group.by parameter) +#' @param normMethod Normalization method for the bigwig files. Default 'RC'. +#' 'RC' will divide the number of fragments in a tile by the total number of +#' fragments in the group. A scaling factor of 10^4 will be applied. +#' 'ncells' will divide the number of fragments in a tile by the number of cells +#' in the group. 'none' will apply no normalization method. +#' A vector of values for each cell can also be passed as a metadata column +#' name. A scaling factor of 10^4 will be applied. +#' @param tileSize The size of the tiles in the bigwig file +#' @param minCells The minimum number of cells in a group to be exported +#' @param cutoff The maximum number of fragments in a given genomic tile +#' @param chromosome A vector of chromosomes to export. If \code{NULL}, use all +#' chromosomes present in \code{seqlengths}. +#' @param seqlengths Chromosome lengths used to define the genomic tiles. Can be +#' a named numeric vector of chromosome lengths, or any object with a +#' \code{seqlengths} method such as a \code{BSgenome} or +#' \code{\link[Seqinfo:Seqinfo]{Seqinfo}} object. If \code{NULL}, the chromosome +#' lengths stored in the object are used; note that these are frequently unset, +#' in which case \code{seqlengths} must be supplied. +#' @param outdir Directory to write output files (split bed files and bigwigs) +#' @param verbose Display messages +#' +#' @importFrom GenomicRanges GRanges slidingWindows +#' @importFrom IRanges IRanges +#' @importFrom future nbrOfWorkers +#' @importFrom future.apply future_lapply +#' @importFrom pbapply pblapply +#' @importFrom SeuratObject DefaultAssay +#' +#' @export +#' @concept bigwig +#' +#' @return Returns a list of paths to the bigwig files that were created +#' +#' @examples +#' \dontrun{ +#' ExportGroupBW(object, assay = "peaks") +#' } +ExportGroupBW <- function( + object, + assay = NULL, + group.by = NULL, + idents = NULL, + normMethod = "RC", + tileSize = 100, + minCells = 5, + cutoff = NULL, + chromosome = NULL, + seqlengths = NULL, + outdir = getwd(), + verbose = TRUE +) { + # Check if output directory exists + if (!dir.exists(paths = outdir)) { + dir.create(path = outdir) + } + if (length(x = Fragments(object = object)) == 0) { + stop("This object does not have Fragments, cannot generate bigwig.") + } + assay <- assay %||% DefaultAssay(object = object) + + # Resolve chromosome lengths, falling back to the object's seqinfo + chromLengths <- seqlengths %||% seqlengths(x = object) + if (!is.numeric(x = chromLengths)) { + # extract a named length vector from a Seqinfo, BSgenome, etc. + chromLengths <- Seqinfo::seqlengths(x = chromLengths) + } + if (is.null(x = chromLengths) || all(is.na(x = chromLengths))) { + stop( + "No chromosome lengths available. Supply them via the `seqlengths` ", + "argument (a named numeric vector, or a Seqinfo or BSgenome object)." + ) + } + # drop chromosomes with unknown length + chromLengths <- chromLengths[!is.na(x = chromLengths)] + + # Subset to requested chromosomes + if (!is.null(x = chromosome)) { + missing.chr <- setdiff(x = chromosome, y = names(x = chromLengths)) + if (length(x = missing.chr) > 0) { + stop( + "Requested chromosome(s) not found in seqlengths: ", + paste(missing.chr, collapse = ", ") + ) + } + chromLengths <- chromLengths[chromosome] + } + availableChr <- names(x = chromLengths) + chromSizes <- GRanges( + seqnames = availableChr, + ranges = IRanges( + start = rep(x = 1, times = length(x = availableChr)), + end = as.numeric(x = chromLengths) + ) + ) + obj.groups <- GetGroups( + object = object, + group.by = group.by, + idents = idents + ) + GroupsNames <- names( + x = table(obj.groups)[table(obj.groups) > minCells] + ) + # Check if output files already exist + lapply(X = GroupsNames, FUN = function(x) { + fn <- paste0(outdir, .Platform$file.sep, x, ".bed") + if (file.exists(fn)) { + message( + sprintf( + paste0( + "The group \"%s\" is already present in the destination folder ", + "and will be overwritten !" + ), + x + ) + ) + file.remove(fn) + } + }) + # Splitting fragments file for each ident in group.by + SplitFragments( + object = object, + assay = assay, + group.by = group.by, + idents = idents, + outdir = outdir, + file.suffix = "", + append = TRUE, + buffer_length = 256L, + verbose = verbose + ) + # Column to normalize by + if (!is.null(x = normMethod)) { + if (tolower(x = normMethod) %in% c("rc", "ncells", "none")) { + normBy <- normMethod + } else { + normBy <- object[[normMethod, drop = FALSE]] + } + } + + if (verbose) { + message("Creating tiles") + } + # Create tiles for each chromosome, from GenomicRanges + tiles <- unlist( + x = slidingWindows(x = chromSizes, width = tileSize, step = tileSize) + ) + if (verbose) { + message("Creating bigwig files at ", outdir) + } + # Run the creation of bigwig for each group of cells + if (nbrOfWorkers() > 1) { + mylapply <- future_lapply + } else { + mylapply <- ifelse(test = verbose, yes = pblapply, no = lapply) + } + covFiles <- mylapply( + GroupsNames, + FUN = CreateBWGroup, + availableChr, + chromLengths, + tiles, + normBy, + tileSize, + normMethod, + cutoff, + outdir + ) + return(covFiles) +} + +# Helper function for ExportGroupBW +# +# @param groupNamei The group to be exported +# @param availableChr Chromosomes to be processed +# @param chromLengths Chromosome lengths +# @param tiles The tiles object +# @param normBy A vector of values to normalize the cells by +# @param tileSize The size of the tiles in the bigwig file +# @param normMethod Normalization method for the bigwig files +# 'RC' will divide the number of fragments in a tile by the number of fragments +# in the group. A scaling factor of 10^4 will be applied +# 'ncells' will divide the number of fragments in a tile by the number of cells +# in the group. 'none' will apply no normalization method. A vector of values +# for each cell can also be passed as a meta.data column name. A scaling factor +# of 10^4 will be applied +# @param cutoff The maximum number of fragments in a given tile +# @param outdir The output directory for bigwig file +# +#' @importFrom GenomicRanges seqnames GRanges +#' @importFrom IRanges coverage +#' @importFrom S4Vectors mcols +#' @importFrom BiocGenerics start end +#' @importFrom Matrix sparseMatrix rowSums +CreateBWGroup <- function( + groupNamei, + availableChr, + chromLengths, + tiles, + normBy, + tileSize, + normMethod, + cutoff, + outdir +) { + if (!requireNamespace("rtracklayer", quietly = TRUE)) { + message( + "Please install rtracklayer. ", + "http://www.bioconductor.org/packages/rtracklayer/" + ) + return(NULL) + } + normMethod <- tolower(x = normMethod) + # Read the fragments file associated with the group + fragi <- rtracklayer::import( + paste0(outdir, .Platform$file.sep, groupNamei, ".bed"), + format = "bed" + ) + cellGroupi <- unique(x = fragi$name) + # Open the writing bigwig file + covFile <- file.path( + outdir, + paste0( + groupNamei, "-TileSize-", tileSize, "-normMethod-", normMethod, ".bw" + ) + ) + + covList <- lapply(X = seq_along(availableChr), FUN = function(k) { + fragik <- fragi[seqnames(x = fragi) == availableChr[k], ] + tilesk <- tiles[ + BiocGenerics::which( + S4Vectors::match(seqnames(x = tiles), availableChr[k], nomatch = 0) > 0 + ) + ] + if (length(x = fragik) == 0) { + tilesk$reads <- 0 + # If fragments + } else { + # N Tiles + nTiles <- chromLengths[availableChr[k]] / tileSize + # Add one tile if there are extra bases + if (nTiles %% 1 != 0) { + nTiles <- trunc(x = nTiles) + 1 + } + # Create Sparse Matrix + matchID <- S4Vectors::match(mcols(x = fragik)$name, cellGroupi) + + # For each tile of this chromosome, create start tile and end tile row, + # set the associated counts matching with the fragments + # This changes compared to ArchR version 1.0.2 + # See https://github.com/GreenleafLab/ArchR/issues/2214 + mat <- sparseMatrix( + i = c( + trunc(x = (start(x = fragik) - 1) / tileSize), + trunc(x = (end(x = fragik) - 1) / tileSize) + ) + 1, + j = as.vector(x = c(matchID, matchID)), + x = rep(x = 1, times = 2 * length(x = fragik)), + dims = c(nTiles, length(x = cellGroupi)) + ) + + # Max count for a cell in a tile is set to cutoff + if (!is.null(x = cutoff)) { + mat@x[mat@x > cutoff] <- cutoff + } + # Sum the cells + mat <- rowSums(x = mat) + tilesk$reads <- mat + # Normalization + if (!is.null(x = normMethod)) { + if (normMethod == "rc") { + tilesk$reads <- tilesk$reads * 10^4 / length(x = fragi$name) + } else if (normMethod == "ncells") { + tilesk$reads <- tilesk$reads / length(x = cellGroupi) + } else if (normMethod == "none") { + } else { + if (!is.null(x = normBy)) { + tilesk$reads <- tilesk$reads * 10^4 / sum(normBy[cellGroupi, 1]) + } + } + } + } + tilesk <- coverage(x = tilesk, weight = tilesk$reads)[[availableChr[k]]] + tilesk + }) + + names(x = covList) <- availableChr + covList <- as(object = covList, Class = "RleList") + rtracklayer::export.bw(object = covList, con = covFile) + return(covFile) +} diff --git a/man/ExportGroupBW.Rd b/man/ExportGroupBW.Rd new file mode 100644 index 00000000..9bb39ac1 --- /dev/null +++ b/man/ExportGroupBW.Rd @@ -0,0 +1,72 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bigwig.R +\name{ExportGroupBW} +\alias{ExportGroupBW} +\title{Export bigwig files for groups of cells} +\usage{ +ExportGroupBW( + object, + assay = NULL, + group.by = NULL, + idents = NULL, + normMethod = "RC", + tileSize = 100, + minCells = 5, + cutoff = NULL, + chromosome = NULL, + seqlengths = NULL, + outdir = getwd(), + verbose = TRUE +) +} +\arguments{ +\item{object}{A Seurat object} + +\item{assay}{Name of assay to use} + +\item{group.by}{The metadata variable used to group the cells} + +\item{idents}{Identities to include (defined by group.by parameter)} + +\item{normMethod}{Normalization method for the bigwig files. Default 'RC'. +'RC' will divide the number of fragments in a tile by the total number of +fragments in the group. A scaling factor of 10^4 will be applied. +'ncells' will divide the number of fragments in a tile by the number of cells +in the group. 'none' will apply no normalization method. +A vector of values for each cell can also be passed as a metadata column +name. A scaling factor of 10^4 will be applied.} + +\item{tileSize}{The size of the tiles in the bigwig file} + +\item{minCells}{The minimum number of cells in a group to be exported} + +\item{cutoff}{The maximum number of fragments in a given genomic tile} + +\item{chromosome}{A vector of chromosomes to export. If \code{NULL}, use all +chromosomes present in \code{seqlengths}.} + +\item{seqlengths}{Chromosome lengths used to define the genomic tiles. Can be +a named numeric vector of chromosome lengths, or any object with a +\code{seqlengths} method such as a \code{BSgenome} or +\code{\link[Seqinfo:Seqinfo]{Seqinfo}} object. If \code{NULL}, the chromosome +lengths stored in the object are used; note that these are frequently unset, +in which case \code{seqlengths} must be supplied.} + +\item{outdir}{Directory to write output files (split bed files and bigwigs)} + +\item{verbose}{Display messages} +} +\value{ +Returns a list of paths to the bigwig files that were created +} +\description{ +Export coverage tracks as bigwig files, one per group of cells. Fragments +are split by group, the genome is tiled, and the number of fragments in each +tile is counted and (optionally) normalized before writing the bigwig files. +} +\examples{ +\dontrun{ +ExportGroupBW(object, assay = "peaks") +} +} +\concept{bigwig} diff --git a/tests/testthat/test-bigwig.R b/tests/testthat/test-bigwig.R new file mode 100644 index 00000000..9e4b0429 --- /dev/null +++ b/tests/testthat/test-bigwig.R @@ -0,0 +1,174 @@ +test_that("CreateBWGroup works with single tile", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "createBW") + dir.create(outdir, showWarnings = FALSE) + fake.bed.data <- data.frame( + seqnames = rep("chr1", 5), + start = c(0, 10, 100, 110, 300), + end = c(100, 150, 200, 250, 500), + cell_name = rep("fake_cell", 5), + nb = 1:5 + ) + write.table( + fake.bed.data, file.path(outdir, "0.bed"), + col.names = FALSE, quote = FALSE, sep = "\t", + row.names = FALSE + ) + CreateBWGroup( + groupNamei = "0", + availableChr = "chr1", + chromLengths = c("chr1" = 249250621), + tiles = GRanges( + seqnames = "chr1", ranges = IRanges(start = 1, end = 249250621) + ), + normBy = NULL, + tileSize = 249250621, + normMethod = "RC", + cutoff = NULL, + outdir = outdir + ) + expect_equal(object = length(list.files(outdir)), expected = 2) + expect( + file.exists(file.path(outdir, "0-TileSize-249250621-normMethod-rc.bw")), + "File does not exist." + ) + bw <- rtracklayer::import.bw( + file.path(outdir, "0-TileSize-249250621-normMethod-rc.bw") + ) + expect_equal(object = bw$score, 20000) +}) + +test_that("CreateBWGroup works with 100bp tile", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "createBW2") + dir.create(outdir, showWarnings = FALSE) + fake.bed.data <- data.frame( + seqnames = rep("chr1", 5), + start = c(0, 10, 100, 110, 300), + end = c(100, 150, 200, 250, 500), + cell_name = rep("fake_cell", 5), + nb = 1:5 + ) + write.table( + fake.bed.data, file.path(outdir, "0.bed"), + col.names = FALSE, quote = FALSE, sep = "\t", + row.names = FALSE + ) + CreateBWGroup( + groupNamei = "0", + availableChr = "chr1", + chromLengths = c("chr1" = 249250621), + tiles = GRanges( + seqnames = "chr1", + ranges = IRanges( + start = seq(1, 249250621, 100), + end = c(seq(100, 249250621, 100), 249250621) + ) + ), + normBy = NULL, + tileSize = 100, + normMethod = "RC", + cutoff = NULL, + outdir = outdir + ) + expect_equal(object = length(list.files(outdir)), expected = 2) + expect( + file.exists(file.path(outdir, "0-TileSize-100-normMethod-rc.bw")), + "File does not exist." + ) + bw <- rtracklayer::import.bw( + file.path(outdir, "0-TileSize-100-normMethod-rc.bw") + ) + expect_equal(object = bw$score, c(6000, 8000, 2000, 0)) +}) + +test_that("CreateBWGroup works with seqlength equal to final pos", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "createBW3") + dir.create(outdir, showWarnings = FALSE) + fake.bed.data <- data.frame( + seqnames = rep("chr1", 5), + start = c(0, 10, 100, 110, 300), + end = c(100, 150, 200, 250, 500), + cell_name = rep("fake_cell", 5), + nb = 1:5 + ) + write.table( + fake.bed.data, file.path(outdir, "0.bed"), + col.names = FALSE, quote = FALSE, sep = "\t", + row.names = FALSE + ) + CreateBWGroup( + groupNamei = "0", + availableChr = "chr1", + chromLengths = c("chr1" = 500), + tiles = GRanges( + seqnames = "chr1", + ranges = IRanges(start = seq(1, 499, 100), end = c(seq(100, 500, 100))) + ), + normBy = NULL, + tileSize = 100, + normMethod = "RC", + cutoff = NULL, + outdir = outdir + ) + expect_equal(object = length(list.files(outdir)), expected = 2) + expect( + file.exists(file.path(outdir, "0-TileSize-100-normMethod-rc.bw")), + "File does not exist." + ) + bw <- rtracklayer::import.bw( + file.path(outdir, "0-TileSize-100-normMethod-rc.bw") + ) + expect_equal(object = bw$score, c(6000, 8000, 2000)) +}) + +test_that("ExportGroupBW works", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "ExportGroupBW") + if (dir.exists(paths = outdir)) { + unlink(x = outdir, recursive = TRUE) + } + dir.create(outdir, showWarnings = FALSE) + fpath <- system.file("extdata", "fragments.tsv.gz", package = "Signac") + cells <- colnames(x = atac_small) + names(x = cells) <- cells + frags <- CreateFragmentObject( + path = fpath, + cells = cells, + verbose = FALSE, + validate.fragments = FALSE + ) + Fragments(atac_small) <- frags + # split into two groups, each above minCells + groups <- rep(x = c("g1", "g2"), length.out = ncol(x = atac_small)) + atac_small <- AddMetaData( + object = atac_small, metadata = groups, col.name = "bw_group" + ) + # v2 objects do not store chromosome lengths, so supply them explicitly + covfiles <- ExportGroupBW( + object = atac_small, + group.by = "bw_group", + normMethod = "RC", + tileSize = 100, + minCells = 5, + seqlengths = c("chr1" = 1e6), + outdir = outdir, + verbose = FALSE + ) + # one bed file + one bigwig file per group + expect_length(object = covfiles, n = 2) + expect( + file.exists(file.path(outdir, "g1-TileSize-100-normMethod-rc.bw")), + "File does not exist." + ) + expect( + file.exists(file.path(outdir, "g2-TileSize-100-normMethod-rc.bw")), + "File does not exist." + ) + bw <- rtracklayer::import.bw( + file.path(outdir, "g1-TileSize-100-normMethod-rc.bw") + ) + expect_equal(object = length(seqlengths(bw)), expected = 1) + expect_equal(object = unname(seqlengths(bw)), expected = 1e6) +}) From 78b81b5c949792c61f5113077c9e95d87cab6351 Mon Sep 17 00:00:00 2001 From: timoast <4591688+timoast@users.noreply.github.com> Date: Thu, 11 Jun 2026 18:03:57 +0800 Subject: [PATCH 2/8] Rename ExportGroupBW to ExportBigwig and address review Rename: - ExportGroupBW -> ExportBigwig (function, docs, tests, NAMESPACE) Correctness: - treat normMethod = NULL as no normalization (previously crashed) - compute metadata-column normalization as the per-group sum of the column from object metadata rather than indexing by the bed file's cell barcodes (which differ after RenameCells, silently giving NaN); validate the column - clamp tile indices so fragment ends past the chromosome boundary land in the last tile instead of an out-of-bounds sparseMatrix index - sanitize group names (spaces, file separators) to match the bed files written by SplitFragments - check for rtracklayer up front instead of after splitting fragments Robustness / API: - add cleanup argument (default TRUE) to remove intermediate bed files - default outdir to tempdir() instead of getwd() - dir.create() with recursive = TRUE - minCells is now an inclusive minimum (>=) Docs: - clarify that tracks are Tn5 insertion (cut-site) coverage, not fragment pileup, and that each fragment contributes two insertion events Tests cover NULL/metadata-column normalization, boundary clamping, group names with spaces, and bed-file cleanup. Co-Authored-By: Claude Fable 5 --- NAMESPACE | 2 +- R/bigwig.R | 126 +++++++++++++------- man/{ExportGroupBW.Rd => ExportBigwig.Rd} | 43 ++++--- tests/testthat/test-bigwig.R | 133 +++++++++++++++++++++- 4 files changed, 246 insertions(+), 58 deletions(-) rename man/{ExportGroupBW.Rd => ExportBigwig.Rd} (51%) diff --git a/NAMESPACE b/NAMESPACE index 72f151ac..c4959ab1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -172,7 +172,7 @@ export(DensityScatter) export(DepthCor) export(DownsampleFeatures) export(EnrichedTerms) -export(ExportGroupBW) +export(ExportBigwig) export(ExpressionPlot) export(Extend) export(FRiP) diff --git a/R/bigwig.R b/R/bigwig.R index 3119e306..d74d86c3 100644 --- a/R/bigwig.R +++ b/R/bigwig.R @@ -4,24 +4,30 @@ NULL #' Export bigwig files for groups of cells #' -#' Export coverage tracks as bigwig files, one per group of cells. Fragments -#' are split by group, the genome is tiled, and the number of fragments in each -#' tile is counted and (optionally) normalized before writing the bigwig files. +#' Export Tn5 insertion coverage tracks as bigwig files, one per group of cells. +#' Fragments are split by group, the genome is tiled, and the number of fragment +#' ends (Tn5 insertion sites) falling in each tile is counted and (optionally) +#' normalized before writing the bigwig files. Note that each fragment +#' contributes two insertion events (one at each end), so the tracks represent +#' insertion-site coverage rather than fragment pileup. #' #' @param object A Seurat object #' @param assay Name of assay to use #' @param group.by The metadata variable used to group the cells #' @param idents Identities to include (defined by group.by parameter) #' @param normMethod Normalization method for the bigwig files. Default 'RC'. -#' 'RC' will divide the number of fragments in a tile by the total number of +#' 'RC' will divide the number of insertions in a tile by the total number of #' fragments in the group. A scaling factor of 10^4 will be applied. -#' 'ncells' will divide the number of fragments in a tile by the number of cells -#' in the group. 'none' will apply no normalization method. -#' A vector of values for each cell can also be passed as a metadata column -#' name. A scaling factor of 10^4 will be applied. +#' 'ncells' will divide the number of insertions in a tile by the number of +#' cells in the group. 'none' (or \code{NULL}) will apply no normalization. +#' The name of a metadata column can also be passed, in which case insertions +#' will be divided by the sum of that column over the cells in the group, with a +#' scaling factor of 10^4 applied. #' @param tileSize The size of the tiles in the bigwig file -#' @param minCells The minimum number of cells in a group to be exported -#' @param cutoff The maximum number of fragments in a given genomic tile +#' @param minCells The minimum number of cells in a group for it to be exported. +#' Groups with fewer than \code{minCells} cells are skipped. +#' @param cutoff The maximum number of insertions for a single cell in a given +#' genomic tile. Counts above this value are capped before summing across cells. #' @param chromosome A vector of chromosomes to export. If \code{NULL}, use all #' chromosomes present in \code{seqlengths}. #' @param seqlengths Chromosome lengths used to define the genomic tiles. Can be @@ -30,7 +36,10 @@ NULL #' \code{\link[Seqinfo:Seqinfo]{Seqinfo}} object. If \code{NULL}, the chromosome #' lengths stored in the object are used; note that these are frequently unset, #' in which case \code{seqlengths} must be supplied. -#' @param outdir Directory to write output files (split bed files and bigwigs) +#' @param outdir Directory to write output files (split bed files and bigwigs). +#' Defaults to a temporary directory. +#' @param cleanup Remove the intermediate per-group bed files after writing the +#' bigwig files. Default TRUE. #' @param verbose Display messages #' #' @importFrom GenomicRanges GRanges slidingWindows @@ -47,9 +56,9 @@ NULL #' #' @examples #' \dontrun{ -#' ExportGroupBW(object, assay = "peaks") +#' ExportBigwig(object, assay = "peaks") #' } -ExportGroupBW <- function( +ExportBigwig <- function( object, assay = NULL, group.by = NULL, @@ -60,17 +69,26 @@ ExportGroupBW <- function( cutoff = NULL, chromosome = NULL, seqlengths = NULL, - outdir = getwd(), + outdir = tempdir(), + cleanup = TRUE, verbose = TRUE ) { # Check if output directory exists if (!dir.exists(paths = outdir)) { - dir.create(path = outdir) + dir.create(path = outdir, recursive = TRUE) + } + if (!requireNamespace("rtracklayer", quietly = TRUE)) { + stop( + "Please install rtracklayer: ", + "BiocManager::install('rtracklayer')" + ) } if (length(x = Fragments(object = object)) == 0) { stop("This object does not have Fragments, cannot generate bigwig.") } assay <- assay %||% DefaultAssay(object = object) + # a NULL normMethod is equivalent to applying no normalization + normMethod <- normMethod %||% "none" # Resolve chromosome lengths, falling back to the object's seqinfo chromLengths <- seqlengths %||% seqlengths(x = object) @@ -111,8 +129,16 @@ ExportGroupBW <- function( group.by = group.by, idents = idents ) + # match the group name sanitization done by SplitFragments so that the bed + # file names we read back line up with the files that were written + group.cells <- names(x = obj.groups) + obj.groups <- gsub(pattern = " ", replacement = "_", x = obj.groups) + obj.groups <- gsub( + pattern = .Platform$file.sep, replacement = "_", x = obj.groups + ) + names(x = obj.groups) <- group.cells GroupsNames <- names( - x = table(obj.groups)[table(obj.groups) > minCells] + x = table(obj.groups)[table(obj.groups) >= minCells] ) # Check if output files already exist lapply(X = GroupsNames, FUN = function(x) { @@ -142,13 +168,25 @@ ExportGroupBW <- function( buffer_length = 256L, verbose = verbose ) - # Column to normalize by - if (!is.null(x = normMethod)) { - if (tolower(x = normMethod) %in% c("rc", "ncells", "none")) { - normBy <- normMethod - } else { - normBy <- object[[normMethod, drop = FALSE]] + # Determine the per-group normalization factor. For a metadata column we sum + # the values over the cells in each group here (keyed by group name), since + # the object cell names do not necessarily match the cell barcodes written to + # the split bed files. + if (tolower(x = normMethod) %in% c("rc", "ncells", "none")) { + normBy <- NULL + } else { + if (!normMethod %in% colnames(x = object[[]])) { + stop( + "normMethod must be one of 'RC', 'ncells', 'none', or the name of a ", + "metadata column. '", normMethod, "' is not a metadata column." + ) } + md <- object[[normMethod]] + normBy <- tapply( + X = md[names(x = obj.groups), 1], + INDEX = obj.groups, + FUN = sum + ) } if (verbose) { @@ -179,10 +217,18 @@ ExportGroupBW <- function( cutoff, outdir ) + # remove the intermediate split bed files (written for every group, including + # those below minCells) + if (cleanup) { + bedfiles <- file.path( + outdir, paste0(unique(x = obj.groups), ".bed") + ) + file.remove(bedfiles[file.exists(bedfiles)]) + } return(covFiles) } -# Helper function for ExportGroupBW +# Helper function for ExportBigwig # # @param groupNamei The group to be exported # @param availableChr Chromosomes to be processed @@ -262,11 +308,17 @@ CreateBWGroup <- function( # set the associated counts matching with the fragments # This changes compared to ArchR version 1.0.2 # See https://github.com/GreenleafLab/ArchR/issues/2214 + # Clamp tile indices so fragment ends sitting beyond the chromosome + # boundary (e.g. from read extension) fall in the last tile rather than + # producing an out-of-bounds matrix index. + startTile <- pmin( + trunc(x = (start(x = fragik) - 1) / tileSize) + 1, nTiles + ) + endTile <- pmin( + trunc(x = (end(x = fragik) - 1) / tileSize) + 1, nTiles + ) mat <- sparseMatrix( - i = c( - trunc(x = (start(x = fragik) - 1) / tileSize), - trunc(x = (end(x = fragik) - 1) / tileSize) - ) + 1, + i = c(startTile, endTile), j = as.vector(x = c(matchID, matchID)), x = rep(x = 1, times = 2 * length(x = fragik)), dims = c(nTiles, length(x = cellGroupi)) @@ -280,17 +332,15 @@ CreateBWGroup <- function( mat <- rowSums(x = mat) tilesk$reads <- mat # Normalization - if (!is.null(x = normMethod)) { - if (normMethod == "rc") { - tilesk$reads <- tilesk$reads * 10^4 / length(x = fragi$name) - } else if (normMethod == "ncells") { - tilesk$reads <- tilesk$reads / length(x = cellGroupi) - } else if (normMethod == "none") { - } else { - if (!is.null(x = normBy)) { - tilesk$reads <- tilesk$reads * 10^4 / sum(normBy[cellGroupi, 1]) - } - } + if (normMethod == "rc") { + tilesk$reads <- tilesk$reads * 10^4 / length(x = fragi$name) + } else if (normMethod == "ncells") { + tilesk$reads <- tilesk$reads / length(x = cellGroupi) + } else if (normMethod == "none") { + # no normalization + } else { + # normBy holds the per-group sum of the requested metadata column + tilesk$reads <- tilesk$reads * 10^4 / normBy[[groupNamei]] } } tilesk <- coverage(x = tilesk, weight = tilesk$reads)[[availableChr[k]]] diff --git a/man/ExportGroupBW.Rd b/man/ExportBigwig.Rd similarity index 51% rename from man/ExportGroupBW.Rd rename to man/ExportBigwig.Rd index 9bb39ac1..882c26a6 100644 --- a/man/ExportGroupBW.Rd +++ b/man/ExportBigwig.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/bigwig.R -\name{ExportGroupBW} -\alias{ExportGroupBW} +\name{ExportBigwig} +\alias{ExportBigwig} \title{Export bigwig files for groups of cells} \usage{ -ExportGroupBW( +ExportBigwig( object, assay = NULL, group.by = NULL, @@ -15,7 +15,8 @@ ExportGroupBW( cutoff = NULL, chromosome = NULL, seqlengths = NULL, - outdir = getwd(), + outdir = tempdir(), + cleanup = TRUE, verbose = TRUE ) } @@ -29,18 +30,21 @@ ExportGroupBW( \item{idents}{Identities to include (defined by group.by parameter)} \item{normMethod}{Normalization method for the bigwig files. Default 'RC'. -'RC' will divide the number of fragments in a tile by the total number of +'RC' will divide the number of insertions in a tile by the total number of fragments in the group. A scaling factor of 10^4 will be applied. -'ncells' will divide the number of fragments in a tile by the number of cells -in the group. 'none' will apply no normalization method. -A vector of values for each cell can also be passed as a metadata column -name. A scaling factor of 10^4 will be applied.} +'ncells' will divide the number of insertions in a tile by the number of +cells in the group. 'none' (or \code{NULL}) will apply no normalization. +The name of a metadata column can also be passed, in which case insertions +will be divided by the sum of that column over the cells in the group, with a +scaling factor of 10^4 applied.} \item{tileSize}{The size of the tiles in the bigwig file} -\item{minCells}{The minimum number of cells in a group to be exported} +\item{minCells}{The minimum number of cells in a group for it to be exported. +Groups with fewer than \code{minCells} cells are skipped.} -\item{cutoff}{The maximum number of fragments in a given genomic tile} +\item{cutoff}{The maximum number of insertions for a single cell in a given +genomic tile. Counts above this value are capped before summing across cells.} \item{chromosome}{A vector of chromosomes to export. If \code{NULL}, use all chromosomes present in \code{seqlengths}.} @@ -52,7 +56,11 @@ a named numeric vector of chromosome lengths, or any object with a lengths stored in the object are used; note that these are frequently unset, in which case \code{seqlengths} must be supplied.} -\item{outdir}{Directory to write output files (split bed files and bigwigs)} +\item{outdir}{Directory to write output files (split bed files and bigwigs). +Defaults to a temporary directory.} + +\item{cleanup}{Remove the intermediate per-group bed files after writing the +bigwig files. Default TRUE.} \item{verbose}{Display messages} } @@ -60,13 +68,16 @@ in which case \code{seqlengths} must be supplied.} Returns a list of paths to the bigwig files that were created } \description{ -Export coverage tracks as bigwig files, one per group of cells. Fragments -are split by group, the genome is tiled, and the number of fragments in each -tile is counted and (optionally) normalized before writing the bigwig files. +Export Tn5 insertion coverage tracks as bigwig files, one per group of cells. +Fragments are split by group, the genome is tiled, and the number of fragment +ends (Tn5 insertion sites) falling in each tile is counted and (optionally) +normalized before writing the bigwig files. Note that each fragment +contributes two insertion events (one at each end), so the tracks represent +insertion-site coverage rather than fragment pileup. } \examples{ \dontrun{ -ExportGroupBW(object, assay = "peaks") +ExportBigwig(object, assay = "peaks") } } \concept{bigwig} diff --git a/tests/testthat/test-bigwig.R b/tests/testthat/test-bigwig.R index 9e4b0429..19a0ab1b 100644 --- a/tests/testthat/test-bigwig.R +++ b/tests/testthat/test-bigwig.R @@ -123,9 +123,9 @@ test_that("CreateBWGroup works with seqlength equal to final pos", { expect_equal(object = bw$score, c(6000, 8000, 2000)) }) -test_that("ExportGroupBW works", { +test_that("ExportBigwig works", { skip_if_not_installed("rtracklayer") - outdir <- file.path(tempdir(), "ExportGroupBW") + outdir <- file.path(tempdir(), "ExportBigwig") if (dir.exists(paths = outdir)) { unlink(x = outdir, recursive = TRUE) } @@ -146,7 +146,7 @@ test_that("ExportGroupBW works", { object = atac_small, metadata = groups, col.name = "bw_group" ) # v2 objects do not store chromosome lengths, so supply them explicitly - covfiles <- ExportGroupBW( + covfiles <- ExportBigwig( object = atac_small, group.by = "bw_group", normMethod = "RC", @@ -171,4 +171,131 @@ test_that("ExportGroupBW works", { ) expect_equal(object = length(seqlengths(bw)), expected = 1) expect_equal(object = unname(seqlengths(bw)), expected = 1e6) + # intermediate bed files are cleaned up by default + expect_false(object = file.exists(file.path(outdir, "g1.bed"))) + expect_false(object = file.exists(file.path(outdir, "g2.bed"))) +}) + +test_that("ExportBigwig handles group names containing spaces", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "ExportBigwig_space") + if (dir.exists(paths = outdir)) { + unlink(x = outdir, recursive = TRUE) + } + dir.create(outdir, showWarnings = FALSE) + fpath <- system.file("extdata", "fragments.tsv.gz", package = "Signac") + cells <- colnames(x = atac_small) + names(x = cells) <- cells + frags <- CreateFragmentObject( + path = fpath, cells = cells, verbose = FALSE, validate.fragments = FALSE + ) + Fragments(atac_small) <- frags + # group labels with spaces must be sanitized consistently so the bed files + # written by SplitFragments are found again + groups <- rep(x = c("group a", "group b"), length.out = ncol(x = atac_small)) + atac_small <- AddMetaData( + object = atac_small, metadata = groups, col.name = "bw_group" + ) + expect_no_error( + ExportBigwig( + object = atac_small, + group.by = "bw_group", + seqlengths = c("chr1" = 1e6), + outdir = outdir, + verbose = FALSE + ) + ) + expect( + file.exists(file.path(outdir, "group_a-TileSize-100-normMethod-rc.bw")), + "File does not exist." + ) +}) + +test_that("CreateBWGroup clamps fragments beyond the chromosome end", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "createBW_clamp") + dir.create(outdir, showWarnings = FALSE) + # a fragment whose end (250) extends past the chromosome length (200) + fake.bed.data <- data.frame( + seqnames = "chr1", + start = 0, + end = 250, + cell_name = "fake_cell", + nb = 1 + ) + write.table( + fake.bed.data, file.path(outdir, "0.bed"), + col.names = FALSE, quote = FALSE, sep = "\t", + row.names = FALSE + ) + # should not error despite the out-of-bounds fragment end + expect_no_error( + CreateBWGroup( + groupNamei = "0", + availableChr = "chr1", + chromLengths = c("chr1" = 200), + tiles = GRanges( + seqnames = "chr1", ranges = IRanges(start = c(1, 101), end = c(100, 200)) + ), + normBy = NULL, + tileSize = 100, + normMethod = "RC", + cutoff = NULL, + outdir = outdir + ) + ) + bw <- rtracklayer::import.bw( + file.path(outdir, "0-TileSize-100-normMethod-rc.bw") + ) + # both ends counted, clamped within [1, 200]: 1 event in tile1, 1 in tile2, + # RC-normalized by the single fragment -> 10000 across the whole chromosome + expect_equal(object = bw$score, 10000) +}) + +test_that("ExportBigwig supports NULL and metadata-column normalization", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "ExportBigwig_norm") + if (dir.exists(paths = outdir)) { + unlink(x = outdir, recursive = TRUE) + } + dir.create(outdir, showWarnings = FALSE) + fpath <- system.file("extdata", "fragments.tsv.gz", package = "Signac") + cells <- colnames(x = atac_small) + names(x = cells) <- cells + frags <- CreateFragmentObject( + path = fpath, cells = cells, verbose = FALSE, validate.fragments = FALSE + ) + Fragments(atac_small) <- frags + atac_small$depth <- seq_len(length.out = ncol(x = atac_small)) + + # NULL normMethod must not error and is treated as "none" + expect_no_error( + ExportBigwig( + object = atac_small, + normMethod = NULL, + seqlengths = c("chr1" = 1e6), + outdir = outdir, + verbose = FALSE + ) + ) + expect( + file.exists( + file.path(outdir, "SeuratProject-TileSize-100-normMethod-none.bw") + ), + "File does not exist." + ) + + # metadata-column normalization must produce finite (non-NaN) scores + ExportBigwig( + object = atac_small, + normMethod = "depth", + seqlengths = c("chr1" = 1e6), + outdir = outdir, + verbose = FALSE + ) + bw <- rtracklayer::import.bw( + file.path(outdir, "SeuratProject-TileSize-100-normMethod-depth.bw") + ) + expect_true(all(is.finite(bw$score))) + expect_true(any(bw$score > 0)) }) From ba2de7a22be327317cdd4e9da9e7704eab1eedda Mon Sep 17 00:00:00 2001 From: timoast <4591688+timoast@users.noreply.github.com> Date: Thu, 11 Jun 2026 18:08:52 +0800 Subject: [PATCH 3/8] Validate ExportBigwig normalization column and empty group set - error clearly if a metadata column used for normalization is non-numeric - warn and return early if no group meets minCells, instead of silently producing no output after splitting fragments Co-Authored-By: Claude Fable 5 --- R/bigwig.R | 13 +++++++++++++ tests/testthat/test-bigwig.R | 37 ++++++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+) diff --git a/R/bigwig.R b/R/bigwig.R index d74d86c3..9ea9767f 100644 --- a/R/bigwig.R +++ b/R/bigwig.R @@ -140,6 +140,13 @@ ExportBigwig <- function( GroupsNames <- names( x = table(obj.groups)[table(obj.groups) >= minCells] ) + if (length(x = GroupsNames) == 0) { + warning( + "No groups contain at least minCells (", minCells, ") cells; ", + "nothing to export." + ) + return(list()) + } # Check if output files already exist lapply(X = GroupsNames, FUN = function(x) { fn <- paste0(outdir, .Platform$file.sep, x, ".bed") @@ -182,6 +189,12 @@ ExportBigwig <- function( ) } md <- object[[normMethod]] + if (!is.numeric(x = md[[1]])) { + stop( + "The '", normMethod, "' metadata column must be numeric to be used ", + "for normalization." + ) + } normBy <- tapply( X = md[names(x = obj.groups), 1], INDEX = obj.groups, diff --git a/tests/testthat/test-bigwig.R b/tests/testthat/test-bigwig.R index 19a0ab1b..84a5e73b 100644 --- a/tests/testthat/test-bigwig.R +++ b/tests/testthat/test-bigwig.R @@ -298,4 +298,41 @@ test_that("ExportBigwig supports NULL and metadata-column normalization", { ) expect_true(all(is.finite(bw$score))) expect_true(any(bw$score > 0)) + + # a non-numeric normalization column gives a clear error + atac_small$celltype <- rep(x = c("a", "b"), length.out = ncol(x = atac_small)) + expect_error( + object = ExportBigwig( + object = atac_small, + normMethod = "celltype", + seqlengths = c("chr1" = 1e6), + outdir = outdir, + verbose = FALSE + ), + regexp = "must be numeric" + ) +}) + +test_that("ExportBigwig warns when no group meets minCells", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "ExportBigwig_mincells") + dir.create(outdir, showWarnings = FALSE) + fpath <- system.file("extdata", "fragments.tsv.gz", package = "Signac") + cells <- colnames(x = atac_small) + names(x = cells) <- cells + frags <- CreateFragmentObject( + path = fpath, cells = cells, verbose = FALSE, validate.fragments = FALSE + ) + Fragments(atac_small) <- frags + expect_warning( + object = res <- ExportBigwig( + object = atac_small, + minCells = 1e6, + seqlengths = c("chr1" = 1e6), + outdir = outdir, + verbose = FALSE + ), + regexp = "minCells" + ) + expect_length(object = res, n = 0) }) From 1c15445c9ca629d3ab5b7e6e044f7e1ce2a24fdf Mon Sep 17 00:00:00 2001 From: timoast <4591688+timoast@users.noreply.github.com> Date: Thu, 11 Jun 2026 18:16:23 +0800 Subject: [PATCH 4/8] Use true per-group cell count for ncells normalization For assays with multiple fragment files, cells are identified in the split bed files by their raw fragment-file barcode, which can collide across files (e.g. reused 10x barcodes). Deriving the cell count from unique bed barcodes therefore undercounts cells and over-normalizes 'ncells' output. Pass the true per-group cell count from ExportBigwig down to CreateBWGroup and use it for 'ncells' normalization, falling back to the unique-barcode count when not supplied (direct CreateBWGroup calls). Document that the per-cell `cutoff` cap is still shared between colliding barcodes. Co-Authored-By: Claude Fable 5 --- R/bigwig.R | 39 ++++++++++++++++++++++++++---------- man/ExportBigwig.Rd | 5 ++++- tests/testthat/test-bigwig.R | 39 ++++++++++++++++++++++++++++++++++++ 3 files changed, 71 insertions(+), 12 deletions(-) diff --git a/R/bigwig.R b/R/bigwig.R index 9ea9767f..38885eaa 100644 --- a/R/bigwig.R +++ b/R/bigwig.R @@ -28,6 +28,9 @@ NULL #' Groups with fewer than \code{minCells} cells are skipped. #' @param cutoff The maximum number of insertions for a single cell in a given #' genomic tile. Counts above this value are capped before summing across cells. +#' Note that cells are identified by their fragment-file barcode, so when an +#' assay contains multiple fragment files this cap is shared between any cells +#' that have the same barcode in different files. #' @param chromosome A vector of chromosomes to export. If \code{NULL}, use all #' chromosomes present in \code{seqlengths}. #' @param seqlengths Chromosome lengths used to define the genomic tiles. Can be @@ -137,9 +140,11 @@ ExportBigwig <- function( pattern = .Platform$file.sep, replacement = "_", x = obj.groups ) names(x = obj.groups) <- group.cells - GroupsNames <- names( - x = table(obj.groups)[table(obj.groups) >= minCells] - ) + # true number of cells per group (used for 'ncells' normalization; deriving + # this from the bed files would undercount cells when fragment-file barcodes + # collide across multiple fragment files) + group.counts <- table(obj.groups) + GroupsNames <- names(x = group.counts[group.counts >= minCells]) if (length(x = GroupsNames) == 0) { warning( "No groups contain at least minCells (", minCells, ") cells; ", @@ -225,6 +230,7 @@ ExportBigwig <- function( chromLengths, tiles, normBy, + group.counts, tileSize, normMethod, cutoff, @@ -247,16 +253,19 @@ ExportBigwig <- function( # @param availableChr Chromosomes to be processed # @param chromLengths Chromosome lengths # @param tiles The tiles object -# @param normBy A vector of values to normalize the cells by +# @param normBy Per-group normalization factor (a named vector keyed by group) +# used when normMethod is a metadata column +# @param nCells Per-group cell counts (a named vector keyed by group) used for +# 'ncells' normalization. If NULL, the number of unique cell barcodes in the +# group's bed file is used. # @param tileSize The size of the tiles in the bigwig file # @param normMethod Normalization method for the bigwig files -# 'RC' will divide the number of fragments in a tile by the number of fragments +# 'RC' will divide the number of insertions in a tile by the number of fragments # in the group. A scaling factor of 10^4 will be applied -# 'ncells' will divide the number of fragments in a tile by the number of cells -# in the group. 'none' will apply no normalization method. A vector of values -# for each cell can also be passed as a meta.data column name. A scaling factor -# of 10^4 will be applied -# @param cutoff The maximum number of fragments in a given tile +# 'ncells' will divide the number of insertions in a tile by the number of cells +# in the group. 'none' will apply no normalization method. A meta.data column +# name can also be passed. A scaling factor of 10^4 will be applied +# @param cutoff The maximum number of insertions for a cell in a given tile # @param outdir The output directory for bigwig file # #' @importFrom GenomicRanges seqnames GRanges @@ -270,6 +279,7 @@ CreateBWGroup <- function( chromLengths, tiles, normBy, + nCells = NULL, tileSize, normMethod, cutoff, @@ -348,7 +358,14 @@ CreateBWGroup <- function( if (normMethod == "rc") { tilesk$reads <- tilesk$reads * 10^4 / length(x = fragi$name) } else if (normMethod == "ncells") { - tilesk$reads <- tilesk$reads / length(x = cellGroupi) + # use the true per-group cell count when supplied, otherwise fall back + # to the number of unique barcodes in the bed file + n.cells <- if (is.null(x = nCells)) { + length(x = cellGroupi) + } else { + nCells[[groupNamei]] + } + tilesk$reads <- tilesk$reads / n.cells } else if (normMethod == "none") { # no normalization } else { diff --git a/man/ExportBigwig.Rd b/man/ExportBigwig.Rd index 882c26a6..4fb0bc2b 100644 --- a/man/ExportBigwig.Rd +++ b/man/ExportBigwig.Rd @@ -44,7 +44,10 @@ scaling factor of 10^4 applied.} Groups with fewer than \code{minCells} cells are skipped.} \item{cutoff}{The maximum number of insertions for a single cell in a given -genomic tile. Counts above this value are capped before summing across cells.} +genomic tile. Counts above this value are capped before summing across cells. +Note that cells are identified by their fragment-file barcode, so when an +assay contains multiple fragment files this cap is shared between any cells +that have the same barcode in different files.} \item{chromosome}{A vector of chromosomes to export. If \code{NULL}, use all chromosomes present in \code{seqlengths}.} diff --git a/tests/testthat/test-bigwig.R b/tests/testthat/test-bigwig.R index 84a5e73b..2f220d12 100644 --- a/tests/testthat/test-bigwig.R +++ b/tests/testthat/test-bigwig.R @@ -38,6 +38,45 @@ test_that("CreateBWGroup works with single tile", { expect_equal(object = bw$score, 20000) }) +test_that("CreateBWGroup ncells uses the supplied per-group cell count", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "createBW_ncells") + dir.create(outdir, showWarnings = FALSE) + fake.bed.data <- data.frame( + seqnames = rep("chr1", 5), + start = c(0, 10, 100, 110, 300), + end = c(100, 150, 200, 250, 500), + cell_name = rep("fake_cell", 5), + nb = 1:5 + ) + write.table( + fake.bed.data, file.path(outdir, "0.bed"), + col.names = FALSE, quote = FALSE, sep = "\t", + row.names = FALSE + ) + # the bed has a single unique barcode, but nCells says the group has 5 cells; + # ncells normalization must divide by 5, not by 1 + CreateBWGroup( + groupNamei = "0", + availableChr = "chr1", + chromLengths = c("chr1" = 249250621), + tiles = GRanges( + seqnames = "chr1", ranges = IRanges(start = 1, end = 249250621) + ), + normBy = NULL, + nCells = c("0" = 5), + tileSize = 249250621, + normMethod = "ncells", + cutoff = NULL, + outdir = outdir + ) + bw <- rtracklayer::import.bw( + file.path(outdir, "0-TileSize-249250621-normMethod-ncells.bw") + ) + # 5 fragments -> 10 insertion events in the single tile, divided by 5 cells + expect_equal(object = bw$score, 2) +}) + test_that("CreateBWGroup works with 100bp tile", { skip_if_not_installed("rtracklayer") outdir <- file.path(tempdir(), "createBW2") From 9cfc33aeb3883548b862c575e11d9d69dbaa8137 Mon Sep 17 00:00:00 2001 From: timoast <4591688+timoast@users.noreply.github.com> Date: Thu, 11 Jun 2026 18:21:08 +0800 Subject: [PATCH 5/8] Use markdown formatting in ExportBigwig documentation Convert the roxygen docs from Rd-style \code{}/\link{} macros to markdown (backticks and [pkg::fn()] links), consistent with the rest of the package (Roxygen: list(markdown = TRUE)). Co-Authored-By: Claude Fable 5 --- R/bigwig.R | 23 +++++++++++------------ man/ExportBigwig.Rd | 17 ++++++++--------- 2 files changed, 19 insertions(+), 21 deletions(-) diff --git a/R/bigwig.R b/R/bigwig.R index 38885eaa..a1682e40 100644 --- a/R/bigwig.R +++ b/R/bigwig.R @@ -15,34 +15,33 @@ NULL #' @param assay Name of assay to use #' @param group.by The metadata variable used to group the cells #' @param idents Identities to include (defined by group.by parameter) -#' @param normMethod Normalization method for the bigwig files. Default 'RC'. -#' 'RC' will divide the number of insertions in a tile by the total number of +#' @param normMethod Normalization method for the bigwig files. Default `RC`. +#' `RC` will divide the number of insertions in a tile by the total number of #' fragments in the group. A scaling factor of 10^4 will be applied. -#' 'ncells' will divide the number of insertions in a tile by the number of -#' cells in the group. 'none' (or \code{NULL}) will apply no normalization. +#' `ncells` will divide the number of insertions in a tile by the number of +#' cells in the group. `none` (or `NULL`) will apply no normalization. #' The name of a metadata column can also be passed, in which case insertions #' will be divided by the sum of that column over the cells in the group, with a #' scaling factor of 10^4 applied. #' @param tileSize The size of the tiles in the bigwig file #' @param minCells The minimum number of cells in a group for it to be exported. -#' Groups with fewer than \code{minCells} cells are skipped. +#' Groups with fewer than `minCells` cells are skipped. #' @param cutoff The maximum number of insertions for a single cell in a given #' genomic tile. Counts above this value are capped before summing across cells. #' Note that cells are identified by their fragment-file barcode, so when an #' assay contains multiple fragment files this cap is shared between any cells #' that have the same barcode in different files. -#' @param chromosome A vector of chromosomes to export. If \code{NULL}, use all -#' chromosomes present in \code{seqlengths}. +#' @param chromosome A vector of chromosomes to export. If `NULL`, use all +#' chromosomes present in `seqlengths`. #' @param seqlengths Chromosome lengths used to define the genomic tiles. Can be #' a named numeric vector of chromosome lengths, or any object with a -#' \code{seqlengths} method such as a \code{BSgenome} or -#' \code{\link[Seqinfo:Seqinfo]{Seqinfo}} object. If \code{NULL}, the chromosome -#' lengths stored in the object are used; note that these are frequently unset, -#' in which case \code{seqlengths} must be supplied. +#' `seqlengths` method such as a `BSgenome` or [Seqinfo::Seqinfo()] object. If +#' `NULL`, the chromosome lengths stored in the object are used; note that these +#' are frequently unset, in which case `seqlengths` must be supplied. #' @param outdir Directory to write output files (split bed files and bigwigs). #' Defaults to a temporary directory. #' @param cleanup Remove the intermediate per-group bed files after writing the -#' bigwig files. Default TRUE. +#' bigwig files. Default `TRUE`. #' @param verbose Display messages #' #' @importFrom GenomicRanges GRanges slidingWindows diff --git a/man/ExportBigwig.Rd b/man/ExportBigwig.Rd index 4fb0bc2b..b837bb95 100644 --- a/man/ExportBigwig.Rd +++ b/man/ExportBigwig.Rd @@ -29,11 +29,11 @@ ExportBigwig( \item{idents}{Identities to include (defined by group.by parameter)} -\item{normMethod}{Normalization method for the bigwig files. Default 'RC'. -'RC' will divide the number of insertions in a tile by the total number of +\item{normMethod}{Normalization method for the bigwig files. Default \code{RC}. +\code{RC} will divide the number of insertions in a tile by the total number of fragments in the group. A scaling factor of 10^4 will be applied. -'ncells' will divide the number of insertions in a tile by the number of -cells in the group. 'none' (or \code{NULL}) will apply no normalization. +\code{ncells} will divide the number of insertions in a tile by the number of +cells in the group. \code{none} (or \code{NULL}) will apply no normalization. The name of a metadata column can also be passed, in which case insertions will be divided by the sum of that column over the cells in the group, with a scaling factor of 10^4 applied.} @@ -54,16 +54,15 @@ chromosomes present in \code{seqlengths}.} \item{seqlengths}{Chromosome lengths used to define the genomic tiles. Can be a named numeric vector of chromosome lengths, or any object with a -\code{seqlengths} method such as a \code{BSgenome} or -\code{\link[Seqinfo:Seqinfo]{Seqinfo}} object. If \code{NULL}, the chromosome -lengths stored in the object are used; note that these are frequently unset, -in which case \code{seqlengths} must be supplied.} +\code{seqlengths} method such as a \code{BSgenome} or \code{\link[Seqinfo:Seqinfo]{Seqinfo::Seqinfo()}} object. If +\code{NULL}, the chromosome lengths stored in the object are used; note that these +are frequently unset, in which case \code{seqlengths} must be supplied.} \item{outdir}{Directory to write output files (split bed files and bigwigs). Defaults to a temporary directory.} \item{cleanup}{Remove the intermediate per-group bed files after writing the -bigwig files. Default TRUE.} +bigwig files. Default \code{TRUE}.} \item{verbose}{Display messages} } From 3988ef11f7a420c1451c390c300151654652e15a Mon Sep 17 00:00:00 2001 From: timoast <4591688+timoast@users.noreply.github.com> Date: Thu, 11 Jun 2026 18:25:52 +0800 Subject: [PATCH 6/8] Default ExportBigwig outdir to the current working directory Co-Authored-By: Claude Fable 5 --- R/bigwig.R | 4 ++-- man/ExportBigwig.Rd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/bigwig.R b/R/bigwig.R index a1682e40..af12ea90 100644 --- a/R/bigwig.R +++ b/R/bigwig.R @@ -39,7 +39,7 @@ NULL #' `NULL`, the chromosome lengths stored in the object are used; note that these #' are frequently unset, in which case `seqlengths` must be supplied. #' @param outdir Directory to write output files (split bed files and bigwigs). -#' Defaults to a temporary directory. +#' Defaults to the current working directory. #' @param cleanup Remove the intermediate per-group bed files after writing the #' bigwig files. Default `TRUE`. #' @param verbose Display messages @@ -71,7 +71,7 @@ ExportBigwig <- function( cutoff = NULL, chromosome = NULL, seqlengths = NULL, - outdir = tempdir(), + outdir = getwd(), cleanup = TRUE, verbose = TRUE ) { diff --git a/man/ExportBigwig.Rd b/man/ExportBigwig.Rd index b837bb95..bb81d01d 100644 --- a/man/ExportBigwig.Rd +++ b/man/ExportBigwig.Rd @@ -15,7 +15,7 @@ ExportBigwig( cutoff = NULL, chromosome = NULL, seqlengths = NULL, - outdir = tempdir(), + outdir = getwd(), cleanup = TRUE, verbose = TRUE ) @@ -59,7 +59,7 @@ a named numeric vector of chromosome lengths, or any object with a are frequently unset, in which case \code{seqlengths} must be supplied.} \item{outdir}{Directory to write output files (split bed files and bigwigs). -Defaults to a temporary directory.} +Defaults to the current working directory.} \item{cleanup}{Remove the intermediate per-group bed files after writing the bigwig files. Default \code{TRUE}.} From d0f1fec2c1827b68d19d8d0f48524dc67543ed40 Mon Sep 17 00:00:00 2001 From: timoast <4591688+timoast@users.noreply.github.com> Date: Thu, 11 Jun 2026 18:30:02 +0800 Subject: [PATCH 7/8] Write ExportBigwig intermediate bed files to a temporary directory Separate the output location (bigwig files, outdir, defaults to the working directory) from the intermediate split bed files. Add a temp.dir argument (default tempdir()) controlling where the per-group bed files are written and cleaned up, so the working directory is not used as scratch space. Co-Authored-By: Claude Fable 5 --- R/bigwig.R | 37 +++++++++++++++++++++++------------- man/ExportBigwig.Rd | 8 ++++++-- tests/testthat/test-bigwig.R | 16 ++++++++++++---- 3 files changed, 42 insertions(+), 19 deletions(-) diff --git a/R/bigwig.R b/R/bigwig.R index af12ea90..bf38e45c 100644 --- a/R/bigwig.R +++ b/R/bigwig.R @@ -38,8 +38,10 @@ NULL #' `seqlengths` method such as a `BSgenome` or [Seqinfo::Seqinfo()] object. If #' `NULL`, the chromosome lengths stored in the object are used; note that these #' are frequently unset, in which case `seqlengths` must be supplied. -#' @param outdir Directory to write output files (split bed files and bigwigs). -#' Defaults to the current working directory. +#' @param outdir Directory to write the output bigwig files. Defaults to the +#' current working directory. +#' @param temp.dir Directory to write the intermediate per-group bed files. +#' Defaults to a temporary directory ([base::tempdir()]). #' @param cleanup Remove the intermediate per-group bed files after writing the #' bigwig files. Default `TRUE`. #' @param verbose Display messages @@ -72,13 +74,17 @@ ExportBigwig <- function( chromosome = NULL, seqlengths = NULL, outdir = getwd(), + temp.dir = tempdir(), cleanup = TRUE, verbose = TRUE ) { - # Check if output directory exists + # Check that the output and temporary directories exist if (!dir.exists(paths = outdir)) { dir.create(path = outdir, recursive = TRUE) } + if (!dir.exists(paths = temp.dir)) { + dir.create(path = temp.dir, recursive = TRUE) + } if (!requireNamespace("rtracklayer", quietly = TRUE)) { stop( "Please install rtracklayer: ", @@ -151,14 +157,14 @@ ExportBigwig <- function( ) return(list()) } - # Check if output files already exist + # Remove any pre-existing split bed files so we do not append to stale data lapply(X = GroupsNames, FUN = function(x) { - fn <- paste0(outdir, .Platform$file.sep, x, ".bed") + fn <- paste0(temp.dir, .Platform$file.sep, x, ".bed") if (file.exists(fn)) { message( sprintf( paste0( - "The group \"%s\" is already present in the destination folder ", + "The group \"%s\" is already present in the temporary folder ", "and will be overwritten !" ), x @@ -167,13 +173,13 @@ ExportBigwig <- function( file.remove(fn) } }) - # Splitting fragments file for each ident in group.by + # Split the fragment file for each group into the temporary directory SplitFragments( object = object, assay = assay, group.by = group.by, idents = idents, - outdir = outdir, + outdir = temp.dir, file.suffix = "", append = TRUE, buffer_length = 256L, @@ -233,13 +239,14 @@ ExportBigwig <- function( tileSize, normMethod, cutoff, - outdir + outdir, + temp.dir ) # remove the intermediate split bed files (written for every group, including # those below minCells) if (cleanup) { bedfiles <- file.path( - outdir, paste0(unique(x = obj.groups), ".bed") + temp.dir, paste0(unique(x = obj.groups), ".bed") ) file.remove(bedfiles[file.exists(bedfiles)]) } @@ -265,7 +272,9 @@ ExportBigwig <- function( # in the group. 'none' will apply no normalization method. A meta.data column # name can also be passed. A scaling factor of 10^4 will be applied # @param cutoff The maximum number of insertions for a cell in a given tile -# @param outdir The output directory for bigwig file +# @param outdir The output directory for the bigwig file +# @param bed.dir Directory containing the per-group bed file to read. If NULL, +# uses outdir. # #' @importFrom GenomicRanges seqnames GRanges #' @importFrom IRanges coverage @@ -282,7 +291,8 @@ CreateBWGroup <- function( tileSize, normMethod, cutoff, - outdir + outdir, + bed.dir = NULL ) { if (!requireNamespace("rtracklayer", quietly = TRUE)) { message( @@ -291,10 +301,11 @@ CreateBWGroup <- function( ) return(NULL) } + bed.dir <- bed.dir %||% outdir normMethod <- tolower(x = normMethod) # Read the fragments file associated with the group fragi <- rtracklayer::import( - paste0(outdir, .Platform$file.sep, groupNamei, ".bed"), + paste0(bed.dir, .Platform$file.sep, groupNamei, ".bed"), format = "bed" ) cellGroupi <- unique(x = fragi$name) diff --git a/man/ExportBigwig.Rd b/man/ExportBigwig.Rd index bb81d01d..132b8c88 100644 --- a/man/ExportBigwig.Rd +++ b/man/ExportBigwig.Rd @@ -16,6 +16,7 @@ ExportBigwig( chromosome = NULL, seqlengths = NULL, outdir = getwd(), + temp.dir = tempdir(), cleanup = TRUE, verbose = TRUE ) @@ -58,8 +59,11 @@ a named numeric vector of chromosome lengths, or any object with a \code{NULL}, the chromosome lengths stored in the object are used; note that these are frequently unset, in which case \code{seqlengths} must be supplied.} -\item{outdir}{Directory to write output files (split bed files and bigwigs). -Defaults to the current working directory.} +\item{outdir}{Directory to write the output bigwig files. Defaults to the +current working directory.} + +\item{temp.dir}{Directory to write the intermediate per-group bed files. +Defaults to a temporary directory (\code{\link[base:tempdir]{base::tempdir()}}).} \item{cleanup}{Remove the intermediate per-group bed files after writing the bigwig files. Default \code{TRUE}.} diff --git a/tests/testthat/test-bigwig.R b/tests/testthat/test-bigwig.R index 2f220d12..218ff3f9 100644 --- a/tests/testthat/test-bigwig.R +++ b/tests/testthat/test-bigwig.R @@ -169,6 +169,11 @@ test_that("ExportBigwig works", { unlink(x = outdir, recursive = TRUE) } dir.create(outdir, showWarnings = FALSE) + bdir <- file.path(tempdir(), "ExportBigwig_bed") + if (dir.exists(paths = bdir)) { + unlink(x = bdir, recursive = TRUE) + } + dir.create(bdir, showWarnings = FALSE) fpath <- system.file("extdata", "fragments.tsv.gz", package = "Signac") cells <- colnames(x = atac_small) names(x = cells) <- cells @@ -193,9 +198,10 @@ test_that("ExportBigwig works", { minCells = 5, seqlengths = c("chr1" = 1e6), outdir = outdir, + temp.dir = bdir, verbose = FALSE ) - # one bed file + one bigwig file per group + # one bigwig file per group written to outdir expect_length(object = covfiles, n = 2) expect( file.exists(file.path(outdir, "g1-TileSize-100-normMethod-rc.bw")), @@ -210,9 +216,11 @@ test_that("ExportBigwig works", { ) expect_equal(object = length(seqlengths(bw)), expected = 1) expect_equal(object = unname(seqlengths(bw)), expected = 1e6) - # intermediate bed files are cleaned up by default - expect_false(object = file.exists(file.path(outdir, "g1.bed"))) - expect_false(object = file.exists(file.path(outdir, "g2.bed"))) + # intermediate bed files go to temp.dir and are cleaned up by default + expect_false(object = file.exists(file.path(bdir, "g1.bed"))) + expect_false(object = file.exists(file.path(bdir, "g2.bed"))) + # no bed files left behind in the output directory + expect_length(object = list.files(outdir, pattern = "[.]bed$"), n = 0) }) test_that("ExportBigwig handles group names containing spaces", { From 462c509237c2f8096f559be7b17a617cd318e8e9 Mon Sep 17 00:00:00 2001 From: timoast <4591688+timoast@users.noreply.github.com> Date: Thu, 11 Jun 2026 18:54:20 +0800 Subject: [PATCH 8/8] Add ExportBigwig to NEWS and pkgdown, fix example - add a NEWS entry for ExportBigwig - use @concept visualization so it appears in the pkgdown reference index (which is built from has_concept()); 'bigwig' matched no section - make the example supply seqlengths, since v2 objects store no chromosome lengths and the previous example would error Co-Authored-By: Claude Fable 5 --- NEWS.md | 1 + R/bigwig.R | 10 ++++++++-- man/ExportBigwig.Rd | 10 ++++++++-- 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/NEWS.md b/NEWS.md index cbc68fbb..56babb2b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -60,6 +60,7 @@ regions side-by-side with shared formatting and axis labels * Added `ReadMQuad()` function to import output from [MQuad](https://github.com/single-cell-genetics/MQuad) for mitochondrial variant analysis * Added `ReadPWM()` function to read `.pwm` files from a directory into a `PWMatrixList` * Added `ReadJASPAR()` function to read motifs from a JASPAR-format file into a `PWMatrixList` +* Added `ExportBigwig()` function to export per-group Tn5 insertion coverage tracks as bigwig files Removed functions: diff --git a/R/bigwig.R b/R/bigwig.R index bf38e45c..35ea4944 100644 --- a/R/bigwig.R +++ b/R/bigwig.R @@ -54,13 +54,19 @@ NULL #' @importFrom SeuratObject DefaultAssay #' #' @export -#' @concept bigwig +#' @concept visualization #' #' @return Returns a list of paths to the bigwig files that were created #' #' @examples #' \dontrun{ -#' ExportBigwig(object, assay = "peaks") +#' # chromosome lengths can be supplied as a BSgenome object, a Seqinfo, or a +#' # named numeric vector +#' ExportBigwig( +#' object, +#' group.by = "celltype", +#' seqlengths = BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38 +#' ) #' } ExportBigwig <- function( object, diff --git a/man/ExportBigwig.Rd b/man/ExportBigwig.Rd index 132b8c88..84f22027 100644 --- a/man/ExportBigwig.Rd +++ b/man/ExportBigwig.Rd @@ -83,7 +83,13 @@ insertion-site coverage rather than fragment pileup. } \examples{ \dontrun{ -ExportBigwig(object, assay = "peaks") +# chromosome lengths can be supplied as a BSgenome object, a Seqinfo, or a +# named numeric vector +ExportBigwig( + object, + group.by = "celltype", + seqlengths = BSgenome.Hsapiens.UCSC.hg38::BSgenome.Hsapiens.UCSC.hg38 +) } } -\concept{bigwig} +\concept{visualization}