Recount fragments over an explicit common ATAC peak universe

Use Signac::FeatureMatrix() to count local indexed fragments over one supplied set of genomic ranges. Rows and cells are explicitly aligned; no existing peak matrices are merged, padded, normalized or binarized.

Usage

RecountATAC(
  object,
  ranges,
  assay = NULL,
  cells = NULL,
  new.assay = NULL,
  process_n = 2000L,
  verbose = TRUE
)

Arguments

object: A Seurat object with a validated Signac ChromatinAssay and attached local fragment objects with explicit cell-to-barcode mappings.
ranges: A non-empty GRanges in 1-based closed coordinates. Its order is the output feature order. Intervals and optional names must be unique. The source assay must declare one verified genome label, and any labels on these ranges must agree. Unlabeled ranges explicitly inherit the source label; this does not verify their external provenance. Known lengths must agree, and source lengths for matching contigs are retained in the result. Chromosome spelling does not establish assembly.
assay: Source chromatin assay; defaults to the object's default assay.
cells: Unique source-assay cell IDs in requested matrix-return order. The default includes every source-assay cell. An empty selection is an error. Seurat writeback uses the original object cell order; the receipt records both requested and stored orders.
new.assay: Optional unused syntactically valid assay name. NULL returns a list; a name returns a Seurat object with a new assay, leaving the source assay and its analysis results unchanged. The constructor requires at least two peaks and two cells for this writeback route.
process_n: Number of regions per FeatureMatrix chunk.
verbose: Whether to print progress messages.

Value

With new.assay = NULL, a list containing sparse counts, aligned ranges named by the output peak IDs, and a receipt with original range names, cell/peak IDs, file and reference information, count semantics and dependency version. With new.assay specified, a Seurat object; its receipt is stored at Misc(object[[new.assay]], "scop_atac_recount").

Details

Counts are fragment-file record overlaps: a record contributes once to each overlapping feature, irrespective of its fifth-column read support. A record may overlap several peaks, including disjoint peaks. Counts are not split or reweighted, so matrix column sums are not necessarily unique fragments in peaks.

Requested raw barcodes and stored file/index hashes are checked with Signac's validators before counting. Different captures may reuse raw barcodes, but each final cell must map to exactly one file. No barcode suffixes, genome labels or chromosome styles are guessed. Only files needed by selected cells are read. Features on chromosomes present in each index are counted with keep_all_features = FALSE; ranges on absent chromosomes are explicitly mapped back as structural zero rows only when declared in the source assay's Seqinfo or observed in another selected file's index. Unknown contigs fail. A Seqinfo declaration establishes a reference name, not sequencing coverage; zero means no record in the supplied file, not biological inaccessibility. This avoids Signac 1.17.1's multi-chunk keep_all_features = TRUE zeroing bug without increasing process_n. A selected file with no matching requested chromosome fails clearly. Genuine zero-overlap features and cells are kept.

Writeback happens only after validation and counting succeed. The new assay contains the recounted matrix, supplied ranges, reference labels and selected fragment references. Assay initialization puts raw counts in its data layer; no normalized data, scaled data, variable features, gene annotation, motifs, links or old peak metadata are copied from the source. Existing reductions remain associated with their original assay; compute a new representation explicitly. The default assay is not changed. R dependencies follow thisutils::check_r() defaults. Genomic resources are not downloaded and external paths are not repaired.

Examples

(function() {
  directory <- tempfile("atac-recount-")
  dir.create(directory)
  on.exit(unlink(directory, recursive = TRUE), add = TRUE)
  path <- file.path(directory, "fragments.tsv")
  writeLines(c("chr1\t10\t20\tA\t5", "chr1\t110\t120\tB\t3"), path)
  compressed <- Rsamtools::bgzip(path, dest = paste0(path, ".gz"))
  Rsamtools::indexTabix(compressed, format = "bed")
  peaks <- GenomicRanges::GRanges(c("chr1:1-50", "chr1:101-150"))
  counts <- Matrix::sparseMatrix(
    i = 1:2, j = 1:2, x = 1, dims = c(2, 2),
    dimnames = list(c("p1", "p2"), c("cellA", "cellB"))
  )
  atac <- ReadATACData(counts,
    ranges = peaks, genome = "toy",
    fragments = compressed, fragment.cells = c(cellA = "A", cellB = "B"),
    verbose = FALSE
  )
  recounted <- RecountATAC(atac, ranges = peaks, verbose = FALSE)
  recounted$counts
})()
#> 2 x 2 sparse Matrix of class "dgCMatrix"
#>              cellA cellB
#> chr1-1-50        1     .
#> chr1-101-150     .     1