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..c4959ab1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -172,6 +172,7 @@ export(DensityScatter) export(DepthCor) export(DownsampleFeatures) export(EnrichedTerms) +export(ExportBigwig) 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/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 new file mode 100644 index 00000000..35ea4944 --- /dev/null +++ b/R/bigwig.R @@ -0,0 +1,400 @@ +#' @include generics.R +#' +NULL + +#' Export bigwig files for groups of cells +#' +#' 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 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 `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 `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 `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 +#' `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 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 +#' +#' @importFrom GenomicRanges GRanges slidingWindows +#' @importFrom IRanges IRanges +#' @importFrom future nbrOfWorkers +#' @importFrom future.apply future_lapply +#' @importFrom pbapply pblapply +#' @importFrom SeuratObject DefaultAssay +#' +#' @export +#' @concept visualization +#' +#' @return Returns a list of paths to the bigwig files that were created +#' +#' @examples +#' \dontrun{ +#' # 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, + assay = NULL, + group.by = NULL, + idents = NULL, + normMethod = "RC", + tileSize = 100, + minCells = 5, + cutoff = NULL, + chromosome = NULL, + seqlengths = NULL, + outdir = getwd(), + temp.dir = tempdir(), + cleanup = TRUE, + verbose = TRUE +) { + # 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: ", + "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) + 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 + ) + # 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 + # 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; ", + "nothing to export." + ) + return(list()) + } + # Remove any pre-existing split bed files so we do not append to stale data + lapply(X = GroupsNames, FUN = function(x) { + fn <- paste0(temp.dir, .Platform$file.sep, x, ".bed") + if (file.exists(fn)) { + message( + sprintf( + paste0( + "The group \"%s\" is already present in the temporary folder ", + "and will be overwritten !" + ), + x + ) + ) + file.remove(fn) + } + }) + # Split the fragment file for each group into the temporary directory + SplitFragments( + object = object, + assay = assay, + group.by = group.by, + idents = idents, + outdir = temp.dir, + file.suffix = "", + append = TRUE, + buffer_length = 256L, + verbose = verbose + ) + # 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]] + 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, + FUN = sum + ) + } + + 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, + group.counts, + tileSize, + normMethod, + cutoff, + outdir, + temp.dir + ) + # remove the intermediate split bed files (written for every group, including + # those below minCells) + if (cleanup) { + bedfiles <- file.path( + temp.dir, paste0(unique(x = obj.groups), ".bed") + ) + file.remove(bedfiles[file.exists(bedfiles)]) + } + return(covFiles) +} + +# Helper function for ExportBigwig +# +# @param groupNamei The group to be exported +# @param availableChr Chromosomes to be processed +# @param chromLengths Chromosome lengths +# @param tiles The tiles object +# @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 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 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 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 +#' @importFrom S4Vectors mcols +#' @importFrom BiocGenerics start end +#' @importFrom Matrix sparseMatrix rowSums +CreateBWGroup <- function( + groupNamei, + availableChr, + chromLengths, + tiles, + normBy, + nCells = NULL, + tileSize, + normMethod, + cutoff, + outdir, + bed.dir = NULL +) { + if (!requireNamespace("rtracklayer", quietly = TRUE)) { + message( + "Please install rtracklayer. ", + "http://www.bioconductor.org/packages/rtracklayer/" + ) + return(NULL) + } + bed.dir <- bed.dir %||% outdir + normMethod <- tolower(x = normMethod) + # Read the fragments file associated with the group + fragi <- rtracklayer::import( + paste0(bed.dir, .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 + # 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(startTile, endTile), + 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 (normMethod == "rc") { + tilesk$reads <- tilesk$reads * 10^4 / length(x = fragi$name) + } else if (normMethod == "ncells") { + # 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 { + # 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]]] + tilesk + }) + + names(x = covList) <- availableChr + covList <- as(object = covList, Class = "RleList") + rtracklayer::export.bw(object = covList, con = covFile) + return(covFile) +} diff --git a/man/ExportBigwig.Rd b/man/ExportBigwig.Rd new file mode 100644 index 00000000..84f22027 --- /dev/null +++ b/man/ExportBigwig.Rd @@ -0,0 +1,95 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bigwig.R +\name{ExportBigwig} +\alias{ExportBigwig} +\title{Export bigwig files for groups of cells} +\usage{ +ExportBigwig( + object, + assay = NULL, + group.by = NULL, + idents = NULL, + normMethod = "RC", + tileSize = 100, + minCells = 5, + cutoff = NULL, + chromosome = NULL, + seqlengths = NULL, + outdir = getwd(), + temp.dir = tempdir(), + cleanup = TRUE, + 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 \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. +\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.} + +\item{tileSize}{The size of the tiles in the bigwig file} + +\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 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.} + +\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::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 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}.} + +\item{verbose}{Display messages} +} +\value{ +Returns a list of paths to the bigwig files that were created +} +\description{ +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{ +# 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{visualization} diff --git a/tests/testthat/test-bigwig.R b/tests/testthat/test-bigwig.R new file mode 100644 index 00000000..218ff3f9 --- /dev/null +++ b/tests/testthat/test-bigwig.R @@ -0,0 +1,385 @@ +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 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") + 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("ExportBigwig works", { + skip_if_not_installed("rtracklayer") + outdir <- file.path(tempdir(), "ExportBigwig") + if (dir.exists(paths = outdir)) { + 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 + 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 <- ExportBigwig( + object = atac_small, + group.by = "bw_group", + normMethod = "RC", + tileSize = 100, + minCells = 5, + seqlengths = c("chr1" = 1e6), + outdir = outdir, + temp.dir = bdir, + verbose = FALSE + ) + # 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")), + "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) + # 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", { + 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)) + + # 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) +})