Package 'BPCells'

Description

Convenience functions for adding or multiplying each row / column of a matrix by a number.

Usage

add_rows(mat, vec)

add_cols(mat, vec)

multiply_rows(mat, vec)

multiply_cols(mat, vec)

Arguments

Value

Matrix-like object

Description

A matrix object can either be an input (i.e. a file on disk or a raw matrix in memory), or it can represent a delayed operation on one or more matrices. The all_matrix_inputs() getter and setter functions allow accessing the base-level input matrices as a list, and changing them. This is useful if you want to re-locate data on disk without losing your transformed BPCells matrix. (Note: experimental API; potentially subject to revisions).

Usage

all_matrix_inputs(x)

all_matrix_inputs(x) <- value

Arguments

Value

List of IterableMatrix objects. If a matrix m is itself an input object, then all_matrix_inputs(m) will return list(m).

Description

Binarize compares the matrix element values to the threshold value and sets the output elements to either zero or one. By default, element values greater than the threshold are set to one; otherwise, set to zero. When strict_inequality is set to FALSE, element values greater than or equal to the threshold are set to one. As an alternative, the <, <=, >, and >= operators are also supported.

Usage

binarize(mat, threshold = 0, strict_inequality = TRUE)

Arguments

Value

binarized IterableMatrix object

Description

Calling peaks from a pre-set list of tiles can be much faster than using dedicated peak-calling software like macs3. The resulting peaks are less precise in terms of exact coordinates, but should be sufficient for most analyses.

Usage

call_peaks_tile(
  fragments,
  chromosome_sizes,
  cell_groups = rep.int("all", length(cellNames(fragments))),
  effective_genome_size = NULL,
  peak_width = 200,
  peak_tiling = 3,
  fdr_cutoff = 0.01,
  merge_peaks = c("all", "group", "none")
)

Arguments

Details

Peak calling steps:

1. Estimate the genome-wide expected insertions per tile based on peak_width, effective_genome_size, and per-group read counts

2. Tile the genome with nonoverlapping tiles of size peak_width

3. For each tile and group, calculate p_value based on a Poisson model

4. Compute adjusted p-values using BH method and using the total number of tiles as the number of hypotheses tested.

5. Repeat steps 2-4 peak_tiling times, with evenly spaced offsets

6. If merge_peaks is "all" or "group": use merge_peaks_iterative() within each group to keep only the most significant of the overlapping candidate peaks

7. If merge_peaks is "all", perform a final round of merge_peaks_iterative(), prioritizing each peak by its within-group significance rank

Value

tibble with peak calls and the following columns:

• chr, start, end: genome coordinates

• group: group ID that this peak was identified in

• p_val, q_val: Poission p-value and BH-corrected p-value

• enrichment: Enrichment of counts in this peak compared to a genome-wide background

Description

Calculate the MD5 checksum of an IterableMatrix and return the checksum in hexidecimal format.

Usage

checksum(matrix)

Arguments

Details

checksum() converts the non-zero elements of the sparse input matrix to double precision, concatenates each element value with the element row and column index words, and uses these 16-byte blocks along with the matrix dimensions and row and column names to calculate the checksum. The checksum value depends on the storage order so column- and row-order matrices with the same element values give different checksum values. checksum() uses element and index values in little-endian CPU storage order. It converts to little-endian order on big-endian architecture although this has not been tested.

Value

MD5 checksum string in hexidecimal format.

Examples

library(Matrix)
library(BPCells)
m1 <- matrix(seq(1,12), nrow=3)
m2 <- as(m1, 'dgCMatrix')
m3 <- as(m2, 'IterableMatrix')
checksum(m3)

Description

Cluster an adjacency matrix

Usage

cluster_graph_leiden(snn, resolution = 0.001, seed = 12531, ...)

cluster_graph_louvain(snn, resolution = 1, seed = 12531)

cluster_graph_seurat(snn, resolution = 0.8, ...)

Arguments

Details

cluster_graph_leiden: Leiden graph clustering algorithm igraph::cluster_leiden()

cluster_graph_louvain: Louvain graph clustering algorithm igraph::cluster_louvain()

cluster_graph_seurat: Seurat's clustering algorithm Seurat::FindClusters()

Value

Factor vector containing the cluster assignment for each cell.

Description

Converts a vector of membership IDs into a sparse matrix

Usage

cluster_membership_matrix(groups, group_order = NULL)

Arguments

Value

cell x group matrix where an entry is 1 when a cell is in a given group

Description

Helper function for data on features to plot from a diverse set of data sources.

Usage

collect_features(
  source,
  features = NULL,
  gene_mapping = human_gene_mapping,
  n = 1
)

Arguments

Details

If source is a data.frame, features will be drawn from the columns. If source is a matrix object (IterableMatrix, dgCMatrix, or matrix), features will be drawn from rows.

Value

Data frame with one column for each feature requested

Description

Convert the type of a matrix

Usage

convert_matrix_type(matrix, type = c("uint32_t", "double", "float"))

Arguments

Value

IterableMatrix object

Description

BPCells fragments can be interconverted with GRanges and data.frame R objects. The main conversion method is R's builtin as() function, though the convert_to_fragments() helper is also available. For all R objects except GRanges, BPCells assumes a 0-based, end-exclusive coordinate system. (See genomic-ranges-like reference for details)

Usage

# Convert from R to BPCells
convert_to_fragments(x, zero_based_coords = !is(x, "GRanges"))
as(x, "IterableFragments")

# Convert from BPCells to R
as.data.frame(bpcells_fragments)
as(bpcells_fragments, "data.frame")
as(bpcells_fragments, "GRanges")

Arguments

Value

convert_to_fragments(): IterableFragments object

Description

These color palettes are derived from the ArchR color palettes, and provide large sets of distinguishable colors

Usage

discrete_palette(name, n = 1)

continuous_palette(name)

Arguments

Details

If the requested number of colors is too large, a new palette will be constructed via interpolation from the requested palette

Value

Character vector of hex color codes

Description

Extend genome ranges in a strand-aware fashion.

Usage

extend_ranges(
  ranges,
  upstream = 0,
  downstream = 0,
  metadata_cols = c("strand"),
  chromosome_sizes = NULL,
  zero_based_coords = !is(ranges, "GRanges")
)

Arguments

Details

Note that ranges will be blocked from extending past the beginning of the chromosome (base 0), and if chromosome_sizes is given then they will also be blocked from extending past the end of the chromosome

Description

Get footprints around a set of genomic coordinates

Usage

footprint(
  fragments,
  ranges,
  zero_based_coords = !is(ranges, "GRanges"),
  cell_groups = rlang::rep_along(cellNames(fragments), "all"),
  cell_weights = rlang::rep_along(cell_groups, 1),
  flank = 125L,
  normalization_width = flank%/%10L
)

Arguments

Value

tibble::tibble() with columns group, position, and count, enrichment

Description

Check if two fragments objects are identical

Usage

fragments_identical(fragments1, fragments2)

Arguments

Value

boolean for whether the fragments objects are identical

Description

Conveniently look up the region of a gene by gene symbol. The value returned by this function can be used as the region argument for trackplot functions such as trackplot_coverage() or trackplot_gene()

Usage

gene_region(
  genes,
  gene_symbol,
  extend_bp = c(10000, 10000),
  gene_mapping = human_gene_mapping
)

Arguments

Value

List of chr, start, end positions for use with trackplot functions.

Description

ArchR-style gene activity scores are based on a weighted sum of each tile according to the signed distance from the tile to a gene body. This function calculates the signed distances according to ArchR's default parameters.

Usage

gene_score_tiles_archr(
  genes,
  chromosome_sizes = NULL,
  tile_width = 500,
  addArchRBug = FALSE
)

Arguments

Details

ArchR's tile distance algorithm works as follows

1. Genes are extended 5kb upstream

2. Genes are linked to any tiles 1kb-100kb upstream + downstream, but tiles beyond a neighboring gene are not considered

Value

Tibble with one range per tile, with additional metadata columns gene_idx (row index of the gene this tile corresponds to) and distance.

Distance is a signed distance calculated such that if the tile has a smaller start coordinate than the gene and the gene is on the + strand, distance will be negative. The distance of adjacent but non-overlapping regions is 1bp, counting up from there.

Description

Gene activity scores can be calculated as a distance-weighted sum of per-tile accessibility. The tile weights for each gene can be represented as a sparse matrix of dimension genes x tiles. If we multiply this weight matrix by a corresponding tile matrix (tiles x cells), then we can get a gene activity score matrix of genes x cells. gene_score_weights_archr() calculates the weight matrix (best if you have a pre-computed tile matrix), while gene_score_archr() provides a easy-to-use wrapper.

Usage

gene_score_weights_archr(
  genes,
  chromosome_sizes,
  blacklist = NULL,
  tile_width = 500,
  gene_name_column = "gene_id",
  addArchRBug = FALSE
)

gene_score_archr(
  fragments,
  genes,
  chromosome_sizes,
  blacklist = NULL,
  tile_width = 500,
  gene_name_column = "gene_id",
  addArchRBug = FALSE,
  tile_max_count = 4,
  scale_factor = 10000,
  tile_matrix_path = tempfile(pattern = "gene_score_tile_mat")
)

Arguments

Details

gene_score_weights_archr:

Given a set of tile coordinates and distances returned by gene_score_tiles_archr(), calculate a weight matrix of dimensions genes x tiles. This matrix can be multiplied with a tile matrix to obtain ArchR-compatible gene activity scores.

Value

gene_score_weights_archr

Weight matrix of dimension genes x tiles

gene_score_archr

Gene score matrix of dimension genes x cells.

Description

BPCells accepts a flexible set of genomic ranges-like objects as input, either GRanges, data.frame, lists, or character vectors. These objects must specify chromosome, start, and end coordinates along with optional metadata about each range. With the exception of GenomicRanges::GRanges objects, BPCells assumes all objects use a zero-based, end-exclusive coordinate system (see below for details).

Valid Range-like objects

BPCells can interpret the following types as ranges:

• list(), data.frame(), with columns:

  • chr: Character or factor of chromosome names

  • start: Start coordinates (0-based)

  • end: End coordinates (exclusive)

  • (optional) strand: "+"/"-" or TRUE/FALSE for pos/neg strand

  • (optional) Additional metadata as named list entries or data.frame columns

• GenomicRanges::GRanges

  • start(x) is interpreted as a 1-based start coordinate

  • end(x) is interpreted as an inclusive end coordinate

  • strand(x): "*" entries are interpeted as postive strand

  • (optional) mcols(x) holds additional metadata

• character

  • Given in format "chr1:1000-2000" or "chr1:1,000-2,000"

  • Uses 0-based, end-exclusive coordinate system

  • Cannot be used for ranges where additional metadata is required

Range coordinate systems

There are two main conventions for the coordinate systems:

One-based, end-inclusive ranges

• The first base of a chromosome is numbered 1

• The last base in a range is equal to the end coordinate

• e.g. 1-5 describes the first 5 bases of the chromosome

• Used in formats such as SAM, GTF

• In BPCells, used when reading or writing GenomicRanges::GRanges objects

Zero-based, end-exclusive ranges

• The first base of a chromosome is numbered 0

• The last base in a range is one less than the end coordinate

• e.g. 0-5 describes the first 5 bases of the chromosome

• Used in formats such as BAM, BED

• In BPCells, used for all other range objects

Description

Mapping of the canonical gene symbols corresponding to each unambiguous alias, previous symbol, ensembl ID, or entrez ID.

Usage

human_gene_mapping

mouse_gene_mapping

Format

human_gene_mapping

A named character vector. Names are aliases or IDs and values are the corresponding canonical gene symbol

mouse_gene_mapping

A named character vector. Names are aliases or IDs and values are the corresponding canonical gene symbol

Details

See the source code in data-raw/human_gene_mapping.R and data-raw/mouse_gene_mapping.R for exactly how these mappings were made.

Source

human_gene_mapping

http://ftp.ebi.ac.uk/pub/databases/genenames/hgnc/tsv/non_alt_loci_set.txt

mouse_gene_mapping

http://www.informatics.jax.org/downloads/reports/MGI_EntrezGene.rpt http://www.informatics.jax.org/downloads/reports/MRK_ENSEMBL.rpt

Description

Read a sparse matrix from a MatrixMarket file. This is a text-based format used by 10x, Parse, and others to store sparse matrices. Format details on the NIST website.

Usage

import_matrix_market(
  mtx_path,
  outdir = tempfile("matrix_market"),
  row_names = NULL,
  col_names = NULL,
  row_major = FALSE,
  tmpdir = tempdir(),
  load_bytes = 4194304L,
  sort_bytes = 1073741824L
)

import_matrix_market_10x(
  mtx_dir,
  outdir = tempfile("matrix_market"),
  feature_type = NULL,
  row_major = FALSE,
  tmpdir = tempdir(),
  load_bytes = 4194304L,
  sort_bytes = 1073741824L
)

Arguments

Details

Import MatrixMarket mtx files to the BPCells format. This implementation ensures fixed memory usage even for very large inputs by doing on-disk sorts. It will be much slower than hdf5 inputs, so only use MatrixMarket format when absolutely necessary.

As a rough speed estimate, importing the 17GB Parse 1M PBMC DGE_1M_PBMC.mtx file takes about 4 minutes and 1.3GB of RAM, producing a compressed output matrix of 1.5GB. mtx.gz files will be slower to import due to gzip decompression.

When importing from 10x mtx files, the row and column names can be read automatically using the import_matrix_market_10x() convenience function.

Value

MatrixDir object with the imported matrix

Description

Methods for IterableFragments objects

Usage

## S4 method for signature 'IterableFragments'
show(object)

cellNames(x)

cellNames(x, ...) <- value

chrNames(x)

chrNames(x, ...) <- value

Arguments

Details

• ⁠cellNames<-⁠ It is only possible to replace names, not add new names.

• ⁠chrNames<-⁠ It is only possible to replace names, not add new names.

Value

• cellNames() Character vector of cell names, or NULL if none are known

• chrNames(): Character vector of chromosome names, or NULL if none are known

Functions

• show(IterableFragments): Print IterableFragments

• cellNames(): Get cell names

• cellNames(x, ...) <- value: Set cell names

• chrNames(): Set chromosome names

• chrNames(x, ...) <- value: Set chromosome names

Description

Generic methods and built-in functions for IterableMatrix objects

Usage

matrix_type(x)

storage_order(x)

## S4 method for signature 'IterableMatrix'
show(object)

## S4 method for signature 'IterableMatrix'
t(x)

## S4 method for signature 'IterableMatrix,matrix'
x %*% y

## S4 method for signature 'IterableMatrix'
rowSums(x)

## S4 method for signature 'IterableMatrix'
colSums(x)

## S4 method for signature 'IterableMatrix'
rowMeans(x)

## S4 method for signature 'IterableMatrix'
colMeans(x)

## S4 method for signature 'IterableMatrix'
log1p(x)

log1p_slow(x)

## S4 method for signature 'IterableMatrix'
expm1(x)

expm1_slow(x)

## S4 method for signature 'IterableMatrix,numeric'
e1 ^ e2

pow_slow(x, exponent)

## S4 method for signature 'numeric,IterableMatrix'
e1 < e2

## S4 method for signature 'IterableMatrix,numeric'
e1 > e2

## S4 method for signature 'numeric,IterableMatrix'
e1 <= e2

## S4 method for signature 'IterableMatrix,numeric'
e1 >= e2

## S4 method for signature 'IterableMatrix'
round(x, digits = 0)

## S4 method for signature 'IterableMatrix,numeric'
e1 * e2

## S4 method for signature 'IterableMatrix,numeric'
e1 + e2

## S4 method for signature 'IterableMatrix,numeric'
e1 / e2

## S4 method for signature 'IterableMatrix,numeric'
e1 - e2

Arguments

Value

• t() Transposed object

• x %*% y: dense matrix result

• rowSums(): vector of row sums

• colSums(): vector of col sums

• rowMeans(): vector of row means

• colMeans(): vector of col means

Functions

• matrix_type(): Get the matrix data type (mat_uint32_t, mat_float, or mat_double for now)

• storage_order(): Get the matrix storage order ("row" or "col")

• show(IterableMatrix): Display an IterableMatrix

• t(IterableMatrix): Transpose an IterableMatrix

• x %*% y: Multiply by a dense matrix

• rowSums(IterableMatrix): Calculate rowSums

• colSums(IterableMatrix): Calculate colSums

• rowMeans(IterableMatrix): Calculate rowMeans

• colMeans(IterableMatrix): Calculate colMeans

• log1p(IterableMatrix): Calculate log(x + 1)

• log1p_slow(): Calculate log(x + 1) (non-SIMD version)

• expm1(IterableMatrix): Calculate exp(x) - 1

• expm1_slow(): Calculate exp(x) - 1 (non-SIMD version)

• e1^e2: Calculate x^y (elementwise)

• pow_slow(): Calculate x^y (elementwise, non-SIMD version)

• e1 < e2: Binarize matrix according to numeric < matrix comparison

• e1 > e2: Binarize matrix according to matrix > numeric comparison

• e1 <= e2: Binarize matrix according to numeric <= matrix comparison

• e1 >= e2: Binarize matrix according to matrix >= numeric comparison

• round(IterableMatrix): round to nearest integer (digits must be 0)

• e1 * e2: Multiply by a constant, or multiply rows by a vector length nrow(mat)

• e1 + e2: Add a constant, or row-wise addition with a vector length nrow(mat)

• e1 / e2: Divide by a constant, or divide rows by a vector length nrow(mat)

• e1 - e2: Subtract a constant, or row-wise subtraction with a vector length nrow(mat)

Description

Search for approximate nearest neighbors between cells in the reduced dimensions (e.g. PCA), and return the k nearest neighbors (knn) for each cell. Optionally, we can find neighbors between two separate sets of cells by utilizing both data and query.

Usage

knn_hnsw(
  data,
  query = NULL,
  k = 10,
  metric = c("euclidean", "cosine"),
  verbose = TRUE,
  threads = 1,
  ef = 100
)

knn_annoy(
  data,
  query = data,
  k = 10,
  metric = c("euclidean", "cosine", "manhattan", "hamming"),
  n_trees = 50,
  search_k = -1
)

Arguments

Details

knn_hnsw: Use RcppHNSW as knn engine

knn_annoy: Use RcppAnnoy as knn engine

Value

List of 2 matrices – idx for cell x K neighbor indices, dist for cell x K neighbor distances. If no query is given, nearest neighbors are found mapping the data matrix to itself, prohibiting self-neighbors

Description

Convert a KNN object (e.g. returned by knn_hnsw() or knn_annoy()) into a graph. The graph is represented as a sparse adjacency matrix.

Usage

knn_to_graph(knn, use_weights = FALSE, self_loops = TRUE)

knn_to_snn_graph(
  knn,
  min_val = 1/15,
  self_loops = FALSE,
  return_type = c("matrix", "list")
)

knn_to_geodesic_graph(knn, return_type = c("matrix", "list"), threads = 0L)

Arguments

Details

knn_to_graph Create a knn graph

knn_to_snn_graph Convert a knn object into a shared nearest neighbors adjacency matrix. This follows the algorithm that Seurat uses to compute SNN graphs

knn_to_geodesic_graph Convert a knn object into an undirected weighted graph, using the same geodesic distance estimation method as the UMAP package. This matches the output of umap._umap.fuzzy_simplicial_set from the umap-learn python package, used by default in scanpy.pp.neighbors. Because this only re-weights and symmetrizes the KNN graph, it will usually use less memory and return a sparser graph than knn_to_snn_graph which computes 2nd-order neighbors. Note: when cells don't have themselves listed as the nearest neighbor, results may differ slightly from umap._umap.fuzzy_simplicial_set, which assumes self is always successfully found in the approximate nearest neighbor search.

Value

knn_to_graph Sparse matrix (dgCMatrix) where mat[i,j] = distance from cell i to cell j, or 0 if cell j is not in the K nearest neighbors of i

knn_to_snn_graph

• return_type == "matrix": Sparse matrix (dgCMatrix) where mat[i,j] = jaccard index of the overlap in nearest neigbors between cell i and cell j, or 0 if the jaccard index is < min_val. Only the lower triangle is filled in, which is compatible with the BPCells clustering methods

• return_type == "list": List of 3 equal-length vectors i, j, and weight, along with an integer dim. These correspond to the rows, cols, and values of non-zero entries in the lower triangle adjacency matrix. dim is the total number of vertices (cells) in the graph

knn_to_geodesic_graph

• return_type == "matrix": Sparse matrix (dgCMatrix) where mat[i,j] = normalized similarity between cell i and cell j. Only the lower triangle is filled in, which is compatible with the BPCells clustering methods

• return_type == "list": List of 3 equal-length vectors i, j, and weight, along with an integer dim. These correspond to the rows, cols, and values of non-zero entries in the lower triangle adjacency matrix. dim is the total number of vertices (cells) in the graph

Description

Given a features x cells matrix, perform one-vs-all differential tests to find markers.

Usage

marker_features(mat, groups, method = "wilcoxon")

Arguments

Details

Tips for using the values from this function:

• Use dplyr::mutate() to add columns for e.g. adjusted p-value and log fold change.

• Use dplyr::filter() to get only differential genes above some given threshold

• To get adjusted p-values, use R p.adjust(), recommended method is "BH"

• To get log2 fold change: if your input matrix was already log-transformed, calculate (foreground_mean - background_mean)/log(2). If your input matrix was not log-transformed, calculate log2(forground_mean/background_mean)

Value

tibble with the following columns:

• foreground: Group ID used for the foreground

• background: Group ID used for the background (or NA if comparing to rest of cells)

• feature: ID of the feature

• p_val_raw: Unadjusted p-value for differential test

• foreground_mean: Average value in the foreground group

• background_mean: Average value in the background group

Description

Correct alias gene symbols, Ensembl IDs, and Entrez IDs to canonical gene symbols. This is useful for matching gene names between different datasets which might not always use the same gene naming conventions.

Usage

match_gene_symbol(query, subject, gene_mapping = human_gene_mapping)

canonical_gene_symbol(query, gene_mapping = human_gene_mapping)

Arguments

Value

match_gene_symbol

Integer vector of indices v such that subject[v] corresponds to the gene symbols in query

canonical_gene_symbol

Character vector of canonical gene symbols for each symbol in query

Description

BPCells matrices can be interconverted with Matrix package dgCMatrix sparse matrices, as well as base R dense matrices (though this may result in high memory usage for large matrices)

Usage

# Convert to R from BPCells
as(bpcells_mat, "dgCMatrix") # Sparse matrix conversion
as.matrix(bpcells_mat) # Dense matrix conversion

# Convert to BPCells from R
as(dgc_mat, "IterableMatrix")

Description

Calculate matrix stats

Usage

matrix_stats(
  matrix,
  row_stats = c("none", "nonzero", "mean", "variance"),
  col_stats = c("none", "nonzero", "mean", "variance"),
  threads = 0L
)

Arguments

Details

The statistics will be calculated in a single pass over the matrix, so this method is desirable to use for efficiency purposes compared to the more standard rowMeans or colMeans if multiple statistics are needed. The stats are ordered by complexity: nonzero, mean, then variance. All less complex stats are calculated in the process of calculating a more complicated stat. So to calculate mean and variance simultaneously, just ask for variance, which will compute mean and nonzero counts as a side-effect

Value

List of row_stats: matrix of n_stats x n_rows, col_stats: matrix of n_stats x n_cols

Description

Peak and tile matrix calculations can be sped up by reducing the number of cells. For cases where the outputs are going to be added together afterwards, this can provide a performance improvement

Usage

merge_cells(fragments, cell_groups)

Arguments

Description

Merge peaks according to ArchR's iterative merging algorithm. More details on the ArchR website

Usage

merge_peaks_iterative(peaks)

Arguments

Details

Properties of merged peaks:

• No peaks in the merged set overlap

• Peaks are prioritized according to their order in the original input

• The output peaks are a subset of the input peaks, with no peak boundaries changed

Value

tibble::tibble() with a nonoverlapping subset of the rows in peaks. All metadata columns are preserved

Description

min_scalar: Take minumum with a global constant

min_by_row: Take the minimum with a per-row constant

min_by_col: Take the minimum with a per-col constant

Usage

min_scalar(mat, val)

min_by_row(mat, vals)

min_by_col(mat, vals)

Arguments

Details

Take the minimum value of a matrix with a per-row, per-col, or global constant. This constant must be >0 to preserve sparsity of the matrix. This has the effect of capping the maximum value in the matrix.

Value

IterableMatrix

Description

Normalize an object representing genomic ranges

Usage

normalize_ranges(
  ranges,
  metadata_cols = character(0),
  zero_based_coords = !is(ranges, "GRanges"),
  n = 1
)

Arguments

Value

data frame with zero-based coordinates, and elements chr (factor), start (int), and end (int). If ranges does not have chr level information, chr levels are the sorted unique values of chr.

If strand is in metadata_cols, then the output strand element will be TRUE for positive strand, and FALSE for negative strand. (Converted from a character vector of "+"/"-" if necessary)

Description

Count fragments by nucleosomal size

Usage

nucleosome_counts(fragments, nucleosome_width = 147)

Arguments

Details

Shorter than nucleosome_width is subNucleosomal, nucleosome_width to 2*nucleosome_width-1 is monoNucleosomal, and anything longer is multiNucleosomal. The sum of all fragments is given as nFrags

Value

List with names subNucleosomal, monoNucleosomal, multiNucleosomal, and nFrags, containing the count vectors of fragments in each class per cell.

Description

10x fragment files come in a bed-like format, with columns chr, start, end, cell_id, and pcr_duplicates. Unlike a standard bed format, the format from cellranger has an inclusive end-coordinate, meaning the end coordinate itself is what should be counted as the tagmentation site, rather than offset by 1.

Usage

open_fragments_10x(path, comment = "#", end_inclusive = TRUE)

write_fragments_10x(
  fragments,
  path,
  end_inclusive = TRUE,
  append_5th_column = FALSE
)

Arguments

Details

open_fragments_10x

No disk operations will take place until the fragments are used in a function

write_fragments_10x

Fragments will be written to disk immediately, then returned in a readable object.

Value

10x fragments file object

Description

Read/write a 10x feature matrix

Usage

open_matrix_10x_hdf5(path, feature_type = NULL, buffer_size = 16384L)

write_matrix_10x_hdf5(
  mat,
  path,
  barcodes = colnames(mat),
  feature_ids = rownames(mat),
  feature_names = rownames(mat),
  feature_types = "Gene Expression",
  feature_metadata = list(),
  buffer_size = 16384L,
  chunk_size = 1024L,
  gzip_level = 0L,
  type = c("uint32_t", "double", "float", "auto")
)

Arguments

Details

The 10x format makes use of gzip compression for the matrix data, which can slow down read performance. Consider writing into another format if the read performance is important to you.

Input matrices must be in column-major storage order, and if the rownames and colnames are not set, names must be provided for the relevant metadata parameters. Some of the metadata parameters are not read by default in BPCells, but it is possible to export them for use with other tools.

Value

BPCells matrix object

Description

Read or write a sparse matrix from an anndata hdf5 file. These functions will automatically transpose matrices when converting to/from the AnnData format. This is because the AnnData convention stores cells as rows, whereas the R convention stores cells as columns. If this behavior is undesired, call t() manually on the matrix inputs and outputs of these functions.

Usage

open_matrix_anndata_hdf5(path, group = "X", buffer_size = 16384L)

write_matrix_anndata_hdf5(
  mat,
  path,
  group = "X",
  buffer_size = 16384L,
  chunk_size = 1024L,
  gzip_level = 0L
)

Arguments

Value

AnnDataMatrixH5 object, with cells as the columns.

Description

Use this function to order regioins prior to calling peak_matrix() or tile_matrix().

Usage

order_ranges(ranges, chr_levels, sort_by_end = TRUE)

Arguments

Value

Numeric vector analagous to the order function. Provides an index selection that will reorder the input ranges to be sorted by chr, end, start

Description

Calculate ranges x cells overlap matrix

Usage

peak_matrix(
  fragments,
  ranges,
  mode = c("insertions", "fragments", "overlaps"),
  zero_based_coords = !is(ranges, "GRanges"),
  explicit_peak_names = TRUE
)

Arguments

Value

Iterable matrix object with dimension ranges x cells. When saved, the column names of the output matrix will be in the format chr1:500-1000, where start and end coords are given in a 0-based coordinate system.

mode options

• "insertions": Start and end coordinates are separately overlapped with each peak

• "fragments": Like "insertions", but each fragment can contribute at most 1 count to each peak, even if both the start and end coordinates overlap

• "overlaps": Like "fragments", but an overlap is also counted if the fragment fully spans the peak even if neither the start or end falls within the peak

Note

When calculating the matrix directly from a fragments tsv, it's necessary to first call select_chromosomes() in order to provide the ordering of chromosomes to expect while reading the tsv.

Description

Plot feature levels per group or cluster as a grid of dots. Dots are colored by z-score normalized average expression, and sized by percent non-zero.

Usage

plot_dot(
  source,
  features,
  groups,
  group_order = NULL,
  gene_mapping = human_gene_mapping,
  colors = c("lightgrey", "#4682B4"),
  return_data = FALSE,
  apply_styling = TRUE
)

Arguments

Description

Plot one or more features by coloring cells in a UMAP plot.

Usage

plot_embedding(
  source,
  embedding,
  features = NULL,
  quantile_range = c(0.01, 0.99),
  randomize_order = TRUE,
  smooth = NULL,
  smooth_rounds = 3,
  gene_mapping = human_gene_mapping,
  size = NULL,
  rasterize = FALSE,
  raster_pixels = 512,
  legend_continuous = c("auto", "quantile", "value"),
  labels_quantile_range = TRUE,
  colors_continuous = c("lightgrey", "#4682B4"),
  legend_discrete = TRUE,
  labels_discrete = TRUE,
  colors_discrete = discrete_palette("stallion"),
  return_data = FALSE,
  return_plot_list = FALSE,
  apply_styling = TRUE
)

Arguments

Details

Smoothing

Smoothing is performed as follows: first, the smoothing matrix is normalized so the sum of incoming weights to every cell is 1. Then, the raw data values are repeatedly multiplied by the smoothing matrix and re-scaled so the average value stays the same.

Value

By default, returns a ggplot2 object with all the requested features plotted in a grid. If return_data or return_plot_list is called, the return value will match that argument.

Description

Plot the distribution of fragment lengths, with length in basepairs on the x-axis, and proportion of fragments on the y-axis. Typical plots will show 10-basepair periodicity, as well as humps spaced at multiples of a nucleosome width (about 150bp).

Usage

plot_fragment_length(
  fragments,
  max_length = 500,
  return_data = FALSE,
  apply_styling = TRUE
)

Arguments

Value

Numeric vector where index i contans the number of length-i fragments

Description

Plots read count rank vs. number of reads on a log-log scale.

Usage

plot_read_count_knee(
  read_counts,
  cutoff = NULL,
  return_data = FALSE,
  apply_styling = TRUE
)

Arguments

Details

Performs logarithmic downsampling to reduce the number of points plotted

Value

ggplot2 plot object

Description

Plot the footprinting around TF motif sites

Usage

plot_tf_footprint(
  fragments,
  motif_positions,
  cell_groups = rlang::rep_along(cellNames(fragments), "all"),
  flank = 250L,
  smooth = 0L,
  zero_based_coords = !is(genes, "GRanges"),
  colors = discrete_palette("stallion"),
  return_data = FALSE,
  apply_styling = TRUE
)

Arguments

See Also

footprint(), plot_tss_profile()

Description

Plot the enrichmment of insertions relative to transcription start sites (TSS). Typically, this plot shows strong enrichment of insertions near a TSS, and a small bump downstream around 220bp downstream of the TSS for the +1 nucleosome.

Usage

plot_tss_profile(
  fragments,
  genes,
  cell_groups = rlang::rep_along(cellNames(fragments), "all"),
  flank = 2000L,
  smooth = 0L,
  zero_based_coords = !is(genes, "GRanges"),
  colors = discrete_palette("stallion"),
  return_data = FALSE,
  apply_styling = TRUE
)

Arguments

See Also

footprint(), plot_tf_footprint()

Description

Density scatter plot with log10(fragment_count) on the x-axis and TSS enrichment on the y-axis. This plot is most useful to select which cell barcodes in an experiment correspond to high-quality cells

Usage

plot_tss_scatter(
  atac_qc,
  min_frags = NULL,
  min_tss = NULL,
  bins = 100,
  apply_styling = TRUE
)

Arguments

Description

Rename cells by adding a prefix to the names. Most commonly this will be a sample name. All cells will recieve the exact text of prefix added to the beginning, so any separator characters like "_" must be included in the given prefix. Use this prior to merging fragments from different experiments with c() in order to help prevent cell name clashes.

Usage

prefix_cell_names(fragments, prefix)

Arguments

Value

Fragments object with prefixed names

Description

Calculate ArchR-compatible per-cell QC statistics

Usage

qc_scATAC(fragments, genes, blacklist)

Arguments

Details

This implementation mimics ArchR's default parameters. For uses requiring more flexibility to tweak default parameters, the best option is to re-implement this function with required changes. Output columns of data.frame:

• cellName: cell name for each cell

• nFrags: number of fragments per cell

• subNucleosomal, monoNucleosomal, multiNucleosomal: number of fragments of size 1-146bp, 147-254bp, and 255bp + respectively. equivalent to ArchR's nMonoFrags, nDiFrags, nMultiFrags respectively

• TSSEnrichment: AvgInsertInTSS / max(AvgInsertFlankingTSS, 0.1), where AvgInsertInTSS is ReadsInTSS / 101 (window size), and AvgInsertFlankingTSS is ReadsFlankingTSS / (100*2) (window size). The max(0.1) ensures that very low-read cells do not get assigned spuriously high TSSEnrichment.

• ReadsInPromoter: Number of reads from 2000bp upstream of TSS to 101bp downstream of TSS

• ReadsInBlacklist: Number of reads in the provided blacklist region

• ReadsInTSS: Number of reads overlapping the 101bp centered around each TSS

• ReadsFlankingTSS: Number of reads overlapping 1901-2000bp +/- each TSS

Differences from ArchR: Note that ArchR by default uses a different set of annotations to derive TSS sites and promoter sites. This function uses just one annotation for gene start+end sites, so must be called twice to exactly re-calculate the ArchR QC stats.

ArchR's PromoterRatio and BlacklistRatio are not included in the output, as they can be easily calculated from ReadsInPromoter / nFrags and ReadsInBlacklist / nFrags. Similarly, ArchR's NucleosomeRatio can be calculated as (monoNucleosomal + multiNucleosomal) / subNucleosomal.

Value

data.frame with QC data

Description

Given a set of genomic ranges, find the distance to the nearest neighbors both upstream and downstream.

Usage

range_distance_to_nearest(
  ranges,
  addArchRBug = FALSE,
  zero_based_coords = !is(ranges, "GRanges")
)

Arguments

Value

A 2-column data.frame with columns upstream and downstream, containing the distances to the nearest neighbor in the respective directions. For ranges on + or * strand, distance is calculated as:

• upstream = max(start(range) - end(upstreamNeighbor), 0)

• downstream = max(start(downstreamNeighbor) - end(range), 0)

For ranges on - strand, the definition of upstream and downstream is flipped. Note that this definition of distance is one off from GenomicRanges::distance(), as ranges that neighbor but don't overlap are given a distance of 1 rather than 0.

Description

Bed files can contain peak or blacklist annotations. These utilities help read thos annotations

Usage

read_bed(
  path,
  additional_columns = character(0),
  backup_url = NULL,
  timeout = 300
)

read_encode_blacklist(
  dir,
  genome = c("hg38", "mm10", "hg19", "dm6", "dm3", "ce11", "ce10"),
  timeout = 300
)

Arguments

Details

read_bed

Read a bed file from disk or a url.

read_encode_blacklist

Downloads the Boyle Lab blacklist, as described in https://doi.org/10.1038/s41598-019-45839-z

Value

Data frame with coordinates using the 0-based convention.

See Also

read_gtf(), read_gencode_genes()

Description

Read gene annotations from gtf format into a data frame. The source can be a URL, a gtf file on disk, or a gencode release version.

Usage

read_gtf(
  path,
  attributes = c("gene_id"),
  tags = character(0),
  features = c("gene"),
  keep_attribute_column = FALSE,
  backup_url = NULL,
  timeout = 300
)

read_gencode_genes(
  dir,
  release = "latest",
  annotation_set = c("basic", "comprehensive"),
  gene_type = "lncRNA|protein_coding|IG_.*_gene|TR_.*_gene",
  attributes = c("gene_id", "gene_type", "gene_name"),
  tags = character(0),
  features = c("gene"),
  timeout = 300
)

read_gencode_transcripts(
  dir,
  release = "latest",
  transcript_choice = c("MANE_Select", "Ensembl_Canonical", "all"),
  annotation_set = c("basic", "comprehensive"),
  gene_type = "lncRNA|protein_coding|IG_.*_gene|TR_.*_gene",
  attributes = c("gene_id", "gene_type", "gene_name", "transcript_id"),
  features = c("transcript", "exon"),
  timeout = 300
)

Arguments

Details

read_gtf

Read gtf from a file or URL

read_gencode_genes

Read gene annotations directly from GENCODE. The file name will vary depending on the release and annotation set requested, but will be of the format gencode.v42.annotation.gtf.gz. GENCODE currently recommends the basic set: https://www.gencodegenes.org/human/. In release 42, both the comprehensive and basic sets had identical gene-level annotations, but the comprehensive set had additional transcript variants annotated.

read_gencode_transcripts

Read transcript models from GENCODE, for use with trackplot_gene()

Value

Data frame with coordinates using the 0-based convention. Columns are:

• chr

• source

• feature

• start

• end

• score

• strand

• frame

• attributes (optional; named according to listed attributes)

• tags (named according to listed tags)

See Also

read_bed(), read_encode_blacklist()

Description

Read chromosome sizes from UCSC and return as a tibble with one row per chromosome. The underlying data is pulled from here: https://hgdownload.soe.ucsc.edu/downloads.html

Usage

read_ucsc_chrom_sizes(
  dir,
  genome = c("hg38", "mm39", "mm10", "mm9", "hg19"),
  keep_chromosomes = "chr[0-9]+|chrX|chrY",
  timeout = 300
)

Description

Rotate ggplot x axis labels

Usage

rotate_x_labels(degrees = 45)

Arguments

Description

Calculate pearson residuals of a negative binomial sctransform model. Normalized values are calculated as (X - mu) / sqrt(mu + mu^2/theta). mu is calculated as cell_read_counts * gene_beta.

Usage

sctransform_pearson(
  mat,
  gene_theta,
  gene_beta,
  cell_read_counts,
  min_var = -Inf,
  clip_range = c(-10, 10),
  columns_are_cells = TRUE,
  slow = FALSE
)

Arguments

Details

The parameterization used is somewhat simplified compared to the original SCTransform paper, in particular it uses a linear-scale rather than log-scale to represent the cell_read_counts and gene_beta variables. It also does not support the addition of arbitrary cell metadata (e.g. batch) to add to the negative binomial regression.

Value

IterableMatrix

Description

Subset, translate, or reorder cell IDs

Usage

select_cells(fragments, cell_selection)

Arguments

Description

Subset, translate, or reorder chromosome IDs

Usage

select_chromosomes(fragments, chromosome_selection)

Arguments

Description

Fragments can be subset based on overlapping (or not overlapping) a set of regions

Usage

select_regions(
  fragments,
  ranges,
  invert_selection = FALSE,
  zero_based_coords = !is(ranges, "GRanges")
)

Arguments

Value

Fragments object filtered according to the selected regions

Description

Adjust labels and heights on trackplots. Labels are set as facet labels in ggplot2, and heights are additional properties read by trackplot_combine() to determine relative height of input plots.

Usage

set_trackplot_label(plot, labels)

set_trackplot_height(plot, height)

get_trackplot_height(plot)

Arguments

Value

set_trackplot_label: ggplot object with adjusted facet labels

set_trackplot_height: ggplot object with adjusted trackplot height

get_trackplot_height: ggplot2::unit object with height setting

Description

Shifts start or end of fragments by a fixed amount, which can be useful to correct the Tn5 offset.

Usage

shift_fragments(fragments, shift_start = 0L, shift_end = 0L)

Arguments

Details

The correct Tn5 offset is +/- 4bp since the Tn5 cut sites on opposite strands are offset by 9bp. However, +4/-5 bp is often applied to bed-format files, since the end coordinate in bed files is 1 past the last basepair of the sequenced DNA fragment. This results in a bed-like format except with inclusive end coordinates.

Value

Shifted fragments object

Description

Subset fragments by length

Usage

subset_lengths(fragments, min_len = 0L, max_len = NA_integer_)

Arguments

Details

Fragment length is calculated as end-start

Value

Fragments object

Description

Use the C++ Spectra solver (same as RSpectra package), in order to compute the largest k values and corresponding singular vectors. Empirically, memory usage is much lower than using irlba::irlba(), likely due to avoiding R garbage creation while solving due to the pure-C++ solver. This documentation is a slightly-edited version of the RSpectra::svds() documentation.

Usage

svds(A, k, nu = k, nv = k, opts = list(), threads=0L, ...)

Arguments

Details

When RSpectra is installed, this function will just add a method to RSpectra::svds() for the IterableMatrix class.

The opts argument is a list that can supply any of the following parameters:

ncv

Number of Lanzcos basis vectors to use. More vectors will result in faster convergence, but with greater memory use. ncv must be satisfy $k < ncv \le p$ where p = min(m, n). Default is min(p, max(2*k+1, 20)).

tol

Precision parameter. Default is 1e-10.

maxitr

Maximum number of iterations. Default is 1000.

center

Either a logical value (TRUE/FALSE), or a numeric vector of length $n$. If a vector $c$ is supplied, then SVD is computed on the matrix $A - 1c'$, in an implicit way without actually forming this matrix. center = TRUE has the same effect as center = colMeans(A). Default is FALSE. Ignored in BPCells

scale

Either a logical value (TRUE/FALSE), or a numeric vector of length $n$. If a vector $s$ is supplied, then SVD is computed on the matrix $(A - 1c')S$, where $c$ is the centering vector and $S = diag(1/s)$. If scale = TRUE, then the vector $s$ is computed as the column norm of $A - 1c'$. Default is FALSE. Ignored in BPCells

Value

A list with the following components:

References

Qiu Y, Mei J (2022). RSpectra: Solvers for Large-Scale Eigenvalue and SVD Problems. R package version 0.16-1, https://CRAN.R-project.org/package=RSpectra.

Description

Calculate ranges x cells tile overlap matrix

Usage

tile_matrix(
  fragments,
  ranges,
  zero_based_coords = !is(ranges, "GRanges"),
  explicit_tile_names = FALSE
)

Arguments

Value

Iterable matrix object with dimension ranges x cells. When saved, the column names will be in the format chr1:500-1000, where start and end coords are given in a 0-based coordinate system.

Note

When calculating the matrix directly from a fragments tsv, it's necessary to first call select_chromosomes() in order to provide the ordering of chromosomes to expect while reading the tsv.

Description

Combines multiple track plots of the same region into a single grid. Uses the patchwork package to perform the alignment.

Usage

trackplot_combine(
  tracks,
  side_plot = NULL,
  title = NULL,
  side_plot_width = 0.3
)

Arguments

Value

A plot object with aligned genome plots. Each aligned row has the text label, y-axis, and plot body. The relative height of each row is given by heights. A shared title and x-axis are put at the top.

See Also

trackplot_coverage(), trackplot_gene(), trackplot_loop(), trackplot_scalebar()

Description

Plot a pseudobulk genome track, showing the number of fragment insertions across a region for each cell type or group.

Usage

trackplot_coverage(
  fragments,
  region,
  groups,
  cell_read_counts,
  group_order = NULL,
  bins = 500,
  clip_quantile = 0.999,
  colors = discrete_palette("stallion"),
  legend_label = "group",
  zero_based_coords = !is(region, "GRanges"),
  return_data = FALSE
)

Arguments

Value

Returns a combined plot of pseudobulk genome tracks. For compatability with draw_trackplot_grid(), the extra attribute ⁠$patches$labels⁠ will be added to specify the labels for each track. If return_data or return_plot_list is TRUE, the return value will be modified accordingly.

See Also

trackplot_combine(), trackplot_gene(), trackplot_loop(), trackplot_scalebar()

Description

Plot transcript models

Usage

trackplot_gene(
  transcripts,
  region,
  exon_size = 2.5,
  gene_size = 0.5,
  label_size = 11 * 0.8/ggplot2::.pt,
  track_label = "Genes",
  return_data = FALSE
)

Arguments

Value

Plot of gene locations

See Also

trackplot_combine(), trackplot_coverage(), trackplot_loop(), trackplot_scalebar()

Description

Plot loops

Usage

trackplot_loop(
  loops,
  region,
  color_by = NULL,
  colors = c("#bfd3e6", "#8c96c6", "#88419d", "#4d004b"),
  allow_truncated = TRUE,
  curvature = 0.75,
  track_label = "Links",
  return_data = FALSE
)

Arguments

Value

Plot of loops connecting genomic coordinates

See Also

trackplot_combine(), trackplot_coverage(), trackplot_gene(), trackplot_scalebar()

Description

Plots a human-readable scale bar and coordinates of the region being plotted

Usage

trackplot_scalebar(region, font_pt = 11)

Arguments

Value

Plot with coordinates and scalebar for plotted genomic region

See Also

trackplot_combine(), trackplot_coverage(), trackplot_gene(), trackplot_loop()

Description

Transpose the storage order for a matrix

Usage

transpose_storage_order(
  matrix,
  outdir = tempfile("transpose"),
  tmpdir = tempdir(),
  load_bytes = 4194304L,
  sort_bytes = 1073741824L
)

Arguments

Details

This re-sorts the entries of a matrix to change the storage order from row-major to col-major. For large matrices, this can be slow – around 2 minutes to transpose a 500k cell RNA-seq matrix The default load_bytes (4MiB) and sort_bytes (1GiB) parameters allow ~85GB of data to be sorted with two passes through the data, and ~7.3TB of data to be sorted in three passes through the data.

Value

MatrixDir object with a copy of the input matrix, but the storage order flipped

Description

BPCells fragments can be read/written in compressed (bitpacked) or uncompressed form in a variety of storage locations: in memory (as an R object), in an hdf5 file, or in a directory on disk (containing binary files).

Usage

write_fragments_memory(fragments, compress = TRUE)

write_fragments_dir(
  fragments,
  dir,
  compress = TRUE,
  buffer_size = 1024L,
  overwrite = FALSE
)

open_fragments_dir(dir, buffer_size = 1024L)

write_fragments_hdf5(
  fragments,
  path,
  group = "fragments",
  compress = TRUE,
  buffer_size = 8192L,
  chunk_size = 1024L,
  overwrite = FALSE,
  gzip_level = 0L
)

open_fragments_hdf5(path, group = "fragments", buffer_size = 16384L)

Arguments

Details

Saving in a directory on disk is a good default for local analysis, as it provides the best I/O performance and lowest memory usage. The HDF5 format allows saving within existing hdf5 files to group data together, and the in memory format provides the fastest performance in the event memory usage is unimportant.

Value

Fragment object

Description

BPCells matrices are stored in sparse format, meaning only the non-zero entries are stored. Matrices can store integer counts data or decimal numbers (float or double). See details for more information.

Usage

write_matrix_memory(mat, compress = TRUE)

write_matrix_dir(
  mat,
  dir,
  compress = TRUE,
  buffer_size = 8192L,
  overwrite = FALSE
)

open_matrix_dir(dir, buffer_size = 8192L)

write_matrix_hdf5(
  mat,
  path,
  group,
  compress = TRUE,
  buffer_size = 8192L,
  chunk_size = 1024L,
  overwrite = FALSE,
  gzip_level = 0L
)

open_matrix_hdf5(path, group, buffer_size = 16384L)

Arguments

Details

Storage locations

Matrices can be stored in a directory on disk, in memory, or in an HDF5 file. Saving in a directory on disk is a good default for local analysis, as it provides the best I/O performance and lowest memory usage. The HDF5 format allows saving within existing hdf5 files to group data together, and the in memory format provides the fastest performance in the event memory usage is unimportant.

Bitpacking Compression

For typical RNA counts matrices holding integer counts, this bitpacking compression will result in 6-8x less space than an R dgCMatrix, and 4-6x smaller than a scipy csc_matrix. The compression will be more effective when the count values in the matrix are small, and when the rows of the matrix are sorted by rowMeans. In tests on RNA-seq data optimal ordering could save up to 40% of storage space. On non-integer data only the row indices are compressed, not the values themselves so space savings will be smaller.

For non-integer data matrices, bitpacking compression is much less effective, as it can only be applied to the indexes of each entry but not the values. There will still be some space savings, but far less than for counts matrices.

Value

BPCells matrix object