Package 'bixverse'

Description

⁠%>>%⁠ chains pipeline steps. Either side can be a ScStep or a ScPipeline; the result is always a ScPipeline.

Usage

lhs %>>% rhs

Arguments

Value

A ScPipeline.

Examples

# either side may be a step, the result is always a pipeline
step_hvg_sc(hvg_no = 30L) %>>% step_pca_sc(no_pcs = 10L)

Description

This method allows you to add ADT counts to a SingleCellsMultiModal. Assumes the transcriptomics data has already been ingested and the cells to keep are known. The method subsets the ADT counts to the kept cells, applies the requested normalisation (CLR or DSB), populates the "var_adt" table in the DuckDB, and attaches an ADTCounts to the class.

Usage

add_adt_counts_sc(object, adt_counts, method = c("clr", "dsb"), ...)

Arguments

Value

Returns a SingleCellsMultiModal with the ADT data added.

Examples

# a CLR-normalised ADT layer on top of an ingested RNA modality
rna <- generate_single_cell_test_data()
adt <- generate_single_cell_test_data_adt()
dir <- tempfile("bixverse_mm")
dir.create(dir)
object <- load_r_data(
  SingleCellsMultiModal(dir_data = dir),
  counts = rna$counts,
  obs = rna$obs,
  var = rna$var,
  sc_qc_param = params_sc_min_quality(min_unique_genes = 5L),
  .verbose = FALSE
)
object <- add_adt_counts_sc(object, adt_counts = adt$counts, method = "clr")
get_adt_names(object)

unlink(dir, recursive = TRUE, force = TRUE)

Description

This function will replace the meta data within the given object

Usage

add_new_metadata(object, new_metadata, ...)

Arguments

Value

The object with updated metadata.

Examples

# swap in a metadata table that carries an extra batch column
set.seed(42)
counts <- matrix(rpois(60, 20), nrow = 10, ncol = 6)
rownames(counts) <- sprintf("gene_%i", 1:10)
colnames(counts) <- sprintf("sample_%i", 1:6)
meta <- data.table::data.table(
  sample_id = colnames(counts),
  case_control = rep(c("case", "control"), each = 3)
)
object <- BulkDge(raw_counts = counts, meta_data = meta)
new_meta <- data.table::copy(meta)[, batch := rep(c("b1", "b2"), 3)]
object <- add_new_metadata(object, new_metadata = new_meta)
head(S7::prop(object, "meta_data"))

Description

This function adds cell type composition information to the nhoods_info slot within the miloR object. For each neighbourhood, it calculates the proportion of the majority cell type and identifies which cell type is most abundant. This is useful for annotating differential abundance results with the cellular composition of each neighbourhood.

Usage

add_nhoods_info(x, cell_info)

## S3 method for class 'miloR'
add_nhoods_info(x, cell_info)

Arguments

Value

Modified miloR object with updated nhoods_info containing majority_celltype and majority_prop columns.

Examples

# tag each neighbourhood with its majority cell type
sc <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 500L,
    n_genes = 50L,
    n_samples = 6L,
    sample_bias = "even"
  )
)
milo <- get_miloR_abundances_sc(
  sc,
  sample_id_col = "sample_id",
  miloR_params = params_sc_miloR(k_refine = 10L),
  .verbose = FALSE
)
design_df <- data.frame(
  grp = rep(c("a", "b"), each = 3),
  row.names = sprintf("sample_%i", 1:6)
)
milo <- test_nhoods(milo, design = ~grp, design_df = design_df)
milo <- add_nhoods_info(milo, cell_info = get_sc_obs(sc)$cell_grp)
head(get_differential_abundance_res(milo))

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Add an obs table derived from a method to the SingleCells.

Usage

add_sc_new_obs(object, obs_data)

Arguments

Value

The class with updated obs table in the DuckDB

Examples

# cluster memberships carry their own cell_idx, so they join straight on
sc <- demo_single_cells()
clusters <- fast_cluster_sc(sc, resolutions = 1.0, .verbose = FALSE)
sc <- add_sc_new_obs(sc, get_data(clusters))
head(get_sc_obs(sc), 3)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

This function will add a modality as an affinity matrix for subsequent fusion to the object. To note: all of the modalities need to have the same number of rows, i.e., samples!

Usage

add_snf_data_modality(object, data, data_name, params = NULL)

Arguments

Value

The class with added adjacency matrix for this data.

Examples

# add a categorical clinical modality to a continuous one
set.seed(42)
continuous <- matrix(rnorm(120), nrow = 12, ncol = 10)
rownames(continuous) <- sprintf("sample_%02i", 1:12)
colnames(continuous) <- sprintf("feature_%i", 1:10)
clinical <- data.table::data.table(
  sample_id = rownames(continuous),
  sex = factor(sample(c("M", "F"), 12, replace = TRUE)),
  stage = factor(sample(c("I", "II", "III"), 12, replace = TRUE))
)
object <- SimilarityNetworkFusion(
  data = continuous,
  data_name = "continuous",
  snf_params = params_snf(k = 3L)
)
object <- add_snf_data_modality(object, clinical, data_name = "clinical")
dim(get_snf_adjcacency_mat(object, name = "clinical"))

Description

Reads one or more obs columns from a SingleCells object and stores them in the reference's labels slot for use by transfer_labels_symphony(). The provided sc_object must have cells_to_keep matching the cells used at reference-build time; this is enforced by a length check against nrow(z_corr).

Usage

add_symphony_labels(reference, sc_object, columns, overwrite = FALSE)

Arguments

Value

The reference with updated labels.

Examples

# attach obs labels to a reference built without them
ref <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 500L,
    n_genes = 50L,
    n_batches = 2L
  )
)
symphony_ref <- build_symphony_ref(
  ref,
  batch_column = "batch_index",
  hvg = get_hvg(ref) + 1L,
  harmony_params = params_sc_harmony(k = 10L),
  no_pcs = 10L,
  .verbose = FALSE
)
symphony_ref <- add_symphony_labels(
  symphony_ref,
  sc_object = ref,
  columns = "cell_grp"
)
head(get_symphony_labels(symphony_ref))

unlink(ref@dir_data, recursive = TRUE, force = TRUE)

Description

This class helps dealing with h5ad objects from Python AnnData. You have several options will allow for retrieval of the underlying data.

Methods

Public methods

• AnnDataParser$new()

• AnnDataParser$get_obs_table()

• AnnDataParser$get_var_info()

• AnnDataParser$get_raw_counts()

• AnnDataParser$get_key_data()

• AnnDataParser$clone()

AnnDataParser$new()

Initialises the Anndata Parser.

Usage

AnnDataParser$new(h5_path)

Arguments

Returns

Returns the initialised class.

AnnDataParser$get_obs_table()

Returns the observation table with all the data from the h5ad file.

Usage

AnnDataParser$get_obs_table()

Returns

data.table. The found observations are returned. The pandas index will be named sample_id. Remaining columns (if found) will be returned as factors due to the way the data is stored in h5.

AnnDataParser$get_var_info()

Returns the variable table with all the data from the h5ad file.

Usage

AnnDataParser$get_var_info()

Returns

data.table. The found observations are returned. The pandas index will be named var_id. Remaining columns (if found) will be returned as factors due to the way the data is stored in h5.

AnnDataParser$get_raw_counts()

Returns the counts that are stored in X slot of the anndata object.

Usage

AnnDataParser$get_raw_counts()

Returns

Returns the count matrix with samples = columns and rows = features.

AnnDataParser$get_key_data()

Wrapper function that returns a list of the stored count data and the metadata found in the h5ad file.

Usage

AnnDataParser$get_key_data()

Returns

List with following elements:

• metadata - metadata from the respective h5ad file

• var_info - metadata on the variables from the respective h5ad file.

• counts - counts that were found in the h5ad file.

AnnDataParser$clone()

The objects of this class are cloneable with this method.

Usage

AnnDataParser$clone(deep = FALSE)

Arguments

Examples

# parse counts and obs back out of a small h5ad file
set.seed(42)
counts <- matrix(rpois(50, 5), nrow = 10, ncol = 5)
obs <- data.table::data.table(sample_id = sprintf("cell_%i", 1:10))
var <- data.table::data.table(var_id = sprintf("gene_%i", 1:5))
h5_path <- tempfile(fileext = ".h5ad")
write_h5ad_sc_dense(h5_path, counts, obs, var, .verbose = FALSE)
parser <- AnnDataParser$new(h5_path)
dim(parser$get_raw_counts())
head(parser$get_obs_table())
unlink(h5_path)

Description

Runs each step in order. The first step receives object; subsequent steps receive the result of the previous step. The chain is validated against the class of object before anything runs, see validate_pipeline(). Errors propagate; nothing is caught.

Usage

apply_pipeline(pipeline, object)

Arguments

Value

The object after all steps have run.

Examples

# HVG then PCA, run in order on a freshly loaded object
sc <- demo_single_cells(prepped = FALSE)
p <- step_hvg_sc(hvg_no = 30L, .verbose = FALSE) %>>%
  step_pca_sc(no_pcs = 10L, .verbose = FALSE)
sc <- apply_pipeline(p, sc)
dim(get_pca_factors(sc))

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Splits object by group_col, applies pipeline to each subset, and returns a named list of processed SingleCellsSubsets. Useful for per-sample / per-cell-type re-analysis where the same chain (HVG, PCA, neighbours, clusters, ...) is run on each group, e.g. sample-pure metacell generation followed by an external merge.

Usage

apply_pipeline_per_group(
  pipeline,
  object,
  group_col,
  groups = NULL,
  progress = FALSE
)

Arguments

Value

Named list of processed objects, names being the group values. Usually SingleCellsSubset, or MetaCells if the pipeline ends on step_metacells_sc(), in which case merge_meta_cells() puts them back together.

Examples

# the same chain re-run inside each cell type
sc <- demo_single_cells(prepped = FALSE)
p <- sc_pipeline() %>>% step_hvg_sc(hvg_no = 20L, .verbose = FALSE)
res <- apply_pipeline_per_group(p, sc, group_col = "cell_grp")
names(res)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Compresses counts_cells.bin into counts.bxa with zstd. The gene-based file is not archived, as it is a transpose of the cell-based one. Normalised values are only kept where they cannot be recomputed from the raw counts. load_existing() restores both binaries automatically when it finds only the archive. Save the in-memory data first with save_sc_exp_to_disk().

Usage

archive_sc_exp(object, level = 3L, remove_bins = TRUE, .verbose = TRUE)

Arguments

Value

A list with n_cells, nnz, n_norm_stored and archive_bytes, invisibly.

Examples

sc <- demo_single_cells()
dir <- sc@dir_data
save_sc_exp_to_disk(sc, type = "rds")
archive_sc_exp(sc, .verbose = FALSE)
sc <- load_existing(SingleCells(dir_data = dir), .verbose = FALSE)

unlink(dir, recursive = TRUE, force = TRUE)

Description

Turns a regulatory potential matrix built elsewhere, for example with nichenetr's construct_ligand_target_matrix(), into a LigandTargetInfluence object, so it can be scored with ligand_activity_scores().

Usage

as_ligand_target_influence(x, ligands_as_cols = FALSE, params = NULL)

Arguments

Value

A LigandTargetInfluence object. Each ligand is its own seed.

Examples

# genes x ligands, as nichenetr returns it
ltm <- matrix(
  c(0.5, 0.0, 0.1, 0.4, 0.3, 0.1),
  nrow = 3,
  dimnames = list(c("G1", "G2", "G3"), c("L1", "L2"))
)
inf <- as_ligand_target_influence(ltm, ligands_as_cols = TRUE)
ligand_activity_scores(inf, gene_sets = list(set_A = "G1"))

Description

The hard tier. Called at the entry of anything that hands cached indices to Rust or writes a derived result back into the object, where a stale artefact means either a crash in the binding layer or silently mis-aligned biology.

Usage

assert_sc_state(
  x,
  artefacts,
  modality = "rna",
  .var.name = checkmate::vname(x),
  add = NULL
)

Arguments

Value

Invisibly x when the assertion holds, otherwise an error.

Description

The cluster-level ScType path (score_clusters()) stamps one label on every cell of a cluster, so a minority population sharing a Leiden community with a bigger one gets absorbed without anything flagging it. This runs the scoring per cell instead: the score matrix is smoothed over the sNN graph via label spreading (Zhou et al.), then each cell takes its own argmax. Cells whose best score falls below score_floor come back as NA.

Pass cluster_col to also get the per-cluster composition (purity, entropy, runner-up cell type) and the hybrid assignment, where clusters at or above purity_threshold keep the blanket cluster-level call and mixed clusters fall back to the per-cell calls.

Usage

assign_sc_type(
  object,
  sc_type_res,
  cluster_col = NULL,
  sctype_cell_params = params_sctype_cells(),
  .verbose = TRUE
)

Arguments

Value

An ScTypeCellResults results class.

References

Ianevski et al., Nat Comm, 2022. Zhou et al., NIPS, 2004.

Examples

# per cell calls from the ScType scores
sc <- demo_single_cells()
markers <- data.table::data.table(
  cell_type = rep(sprintf("cell_type_%i", 1:3), each = 10),
  gene_id = sprintf("gene_%02d", 1:30)
)
cell_markers <- prepare_cell_markers(obj = sc, marker_df = markers)
scores <- calc_sc_type_scores(
  sc,
  cell_marker_list = cell_markers,
  .verbose = FALSE
)
res <- assign_sc_type(sc, sc_type_res = scores, .verbose = FALSE)
table(res$assignments)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Calculates an AUC-type score akin to AUCell across the gene sets, see Aibar et al. Three statistics are on offer, all consuming the same within-cell ranking but weighting it differently. "recovery" (default) is the recovery-curve AUC under a rank cutoff, i.e. the AUCell statistic of Aibar, et al., and is top-heavy: only genes inside the top max_rank of the cell contribute. "wilcox" is the AUC derived from the Mann-Whitney U statistic over the full ranking; its null sits at 0.5 for any gene set size, which makes it a good fit for pathway activity. "ap" is average precision, the most top-heavy of the three, but its null tracks the gene set prevalence, so raw values are not comparable across gene sets of different size unless standardise is on.

Usage

aucell_sc(
  object,
  gs_list,
  aucell_params = params_sc_aucell(),
  streaming = NULL,
  .verbose = TRUE
)

Arguments

Value

AUCell results in form of a matrix that is cells x gene sets or as ScMatrixRes pending the input.

References

Aibar, et al., Nat Methods, 2017

Examples

# recovery curve AUC for two marker programmes
sc <- demo_single_cells()
gs_list <- list(
  type_1 = get_gene_names(sc)[1:10],
  type_2 = get_gene_names(sc)[11:20]
)
res <- aucell_sc(sc, gs_list = gs_list, .verbose = FALSE)
dim(res)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Runs a linear batch correction over the data regressing out batch effects and adding normalised_counts_corrected to the object. Should these counts be found by calculate_dge_hedges(), they will be used for calculations of effect sizes based on Hedge's G.

Usage

batch_correction_bulk_dge(
  object,
  contrast_column,
  batch_col,
  scale_genes = FALSE,
  no_hvg_genes = 2500L
)

Arguments

Value

Returns the class with additional data added to the outputs.

Examples

# regress out a batch column and keep the corrected counts
syn <- synthetic_bulk_cor_matrix()
meta <- data.table::data.table(
  sample_id = colnames(syn$counts),
  case_control = rep(c("case", "control"), each = 50),
  batch = rep(c("b1", "b2"), times = 50)
)
object <- BulkDge(raw_counts = syn$counts, meta_data = meta)
object <- qc_bulk_dge(object, group_col = "case_control", .verbose = FALSE)
object <- normalise_bulk_dge(
  object,
  group_col = "case_control",
  .verbose = FALSE
)
object <- calculate_pca_bulk_dge(object, no_hvg_genes = 500L)
object <- batch_correction_bulk_dge(
  object,
  contrast_column = "case_control",
  batch_col = "batch",
  no_hvg_genes = 500L
)
get_outputs(object)$normalised_counts_corrected[1:3, 1:3]

Description

This function implements the batch-balanced k-nearest neighbour algorithm from Polański, et al. Briefly, the algorithm generates a KNN index on a per batch basis and identifies the neighbours of cells for each individual index. Subsequently, it leverages UMAP connectivity calculations to reduce spurious connections. For more details, please refer to Polański, et al.

Usage

bbknn_sc(
  object,
  batch_column,
  no_neighbours_to_keep = 5L,
  embd_to_use = "pca",
  no_embd_to_use = NULL,
  bbknn_params = params_sc_bbknn(),
  seed = 42L,
  .verbose = TRUE
)

Arguments

Value

The object with added kNN matrix based on BBKNN and the graph based on the returned connectivities of the algorithm.

References

Polański, et al., Bioinformatics, 2020

Examples

# batch balanced kNN, replacing the graph the demo object carries
sc <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 600L, n_genes = 50L, n_batches = 3L
  )
)
sc <- bbknn_sc(
  remove_knn(sc),
  batch_column = "batch_index",
  no_neighbours_to_keep = 5L,
  .verbose = FALSE
)
dim(get_knn_mat(sc))

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

The last SCENIC step. Each regulon gets its own threshold derived from the shape of its AUC distribution across cells, and a cell counts as on when its score sits strictly above that threshold.

The thresholds come back alongside the calls, so you can inspect them, override any that look wrong and re-apply the comparison yourself. SCENIC does the same, it writes them to an editable file between scoring and assignment.

Regulons flagged as not bimodal fell back to mean + 2 * sd. If most of your regulons land there, the AUC distributions are too flat to separate, which usually points at the scoring statistic rather than the thresholding. Check you used auc_type = "recovery", see params_sc_aucell().

Usage

binarise_regulon_activity(
  auc_matrix,
  binarise_params = params_scenic_binarise(),
  .verbose = TRUE
)

Arguments

Value

A list with:

• binary - Logical matrix of cells x regulons, TRUE where the regulon is on.

• thresholds - data.table with one row per regulon, holding the threshold, whether it was bimodal and the number of cells called on.

References

Aibar, et al., Nat Methods, 2017

Examples

# on/off calls for a bimodal and a unimodal regulon
set.seed(7L)
auc <- cbind(
  regulon_a = c(rnorm(50, 0.1, 0.02), rnorm(50, 0.4, 0.02)),
  regulon_b = rnorm(100, 0.2, 0.05)
)
rownames(auc) <- sprintf("cell_%i", 1:100)
binarise_regulon_activity(auc, .verbose = FALSE)$thresholds

Description

Draws random gene sets across a log-spaced grid of set sizes, fits gamma tails to the resulting enrichment scores and smooths the fitted parameters across sizes, see Lachmann, et al. No gene set library enters here, so one calibration serves every library you score against that signature. This is where nearly all of the runtime sits; scoring afterwards costs one gamma tail evaluation per pathway.

The returned null is a plain list, so it survives saveRDS() and can be cached against a signature.

Usage

blitzgsea_calibrate(stats, blitz_params = params_blitzgsea())

Arguments

Value

An object of class BlitzGseaNull, a list with the following elements:

• anchor_sizes - Numeric vector. The anchor set sizes, ascending.

• shape_pos - Numeric vector. Smoothed positive-tail gamma shape.

• scale_pos - Numeric vector. Smoothed positive-tail gamma scale.

• shape_neg - Numeric vector. Smoothed negative-tail gamma shape.

• scale_neg - Numeric vector. Smoothed negative-tail gamma scale.

• pos_ratio - Numeric vector. Smoothed fraction of positive null scores at each anchor.

• ks_pos - Float. Mean goodness-of-fit p-value for the positive tail.

• ks_neg - Float. Mean goodness-of-fit p-value for the negative tail.

• centred - Boolean. Whether the signature was centred.

• n_genes - Integer. Size of the signature it was calibrated on.

References

Lachmann, et al., Bioinformatics, 2022

Examples

# calibrate a null once, to be reused across gene set libraries
set.seed(42L)
stats <- stats::setNames(rnorm(500), sprintf("gene_%03i", 1:500))
null_model <- blitzgsea_calibrate(
  stats,
  blitz_params = params_blitzgsea(permutations = 1000L, anchors = 10L)
)
null_model

Description

Bonsai reconstructs a tree over the cells in which every cell is a leaf and the internal nodes are inferred ancestral states, with branch lengths that carry the amount of change between them. Unlike a kNN graph it uses each measurement's error bar, which is what Sanity provides: the raw counts go through Sanity first, for posterior log fold changes with error bars, and Bonsai builds the tree on those. The tree is then laid out in 2D. For details, please refer to de Groot, et al. and Breda, et al.

No HVG selection needed. By default every gene goes in, streamed through Sanity in chunks, and only the genes with enough signal over their own noise (min_signal_to_noise in params_sc_bonsai()) are kept for the tree. That is the gene selection of the Bonsai paper, and it keeps memory at one chunk plus the survivors. Pass hvg to restrict the candidates.

Runtime grows a little faster than linearly with the number of cells. In bonsai-rs's own benchmarks the search took about 70 seconds at 10,000 cells and 210 seconds at 25,000 on ten cores. Sanity comes on top, linear in the number of genes it has to fit; on the CPU that is the larger share once all genes go in.

The object itself is not touched. Store the leaf coordinates with set_bonsai_embedding() if you want them next to the other embeddings.

On MetaCells every metacell is a leaf. Their raw counts are sums over their cells, and a sum of Poisson counts is Poisson, so Sanity treats a metacell exactly as it treats a cell: its total counts are the library size, and a bigger metacell gets tighter error bars, which Bonsai weights by. That is the way to large data: 100k cells in metacells of 50 is a 2,000-leaf tree. What you give up is any structure inside a metacell.

Usage

bonsai_sc(
  object,
  hvg = NULL,
  bonsai_params = params_sc_bonsai(),
  .verbose = TRUE
)

Arguments

Value

A BonsaiTree S3 object with:

• nodes - data.table with node, parent (NA for the root), branch, is_leaf, cell_id (NA for the inferred ancestors), n_cells (the cells behind a leaf: 1 for single cells, the originating cells for a metacell; NA for the ancestors), x and y. The leaves come first, in the order of the object's cells or metacells.

• loglik - The loglikelihood of the final tree.

• steps - data.table with the loglikelihood after each search step and the step's wall time in seconds.

• timings - data.table with the wall time in seconds of each stage: sanity, ingest, bonsai (the whole search), layout, and total, the whole call as R saw it.

• genes_used - The genes the tree was built on.

• genes_dropped - The candidate genes left out: no counts in the cells, ill-conditioned Sanity posteriors, or a signal-to-noise ratio below min_signal_to_noise.

• cell_idx - The cells or metacells the tree was built over (0-indexed).

• layout, hyperbolic - The current layout.

• params - The parameters of the run.

References

de Groot, et al., Nat. Biotechnol., 2026; Breda, et al., Nat. Biotechnol., 2021.

Examples

# a tree over the demo cells, genes selected by their signal-to-noise
sc <- demo_single_cells(
  prepped = FALSE,
  syn_data_params = params_sc_synthetic_data(n_cells = 200L, n_genes = 50L)
)
tree <- bonsai_sc(sc, .verbose = FALSE)
tree

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Turns the TF to gene table into the gene sets you hand to aucell_sc(). tf_to_genes_motif_enrichment() only flags which pairs survived the CisTarget leading edge, it does not remove anything, so this is the step that applies that filter.

Following SCENIC, the TF is added back to its own target list and regulons below min_genes are dropped. Note the CisTarget leading edge usually cuts hard: a module of a few hundred candidate targets typically ends up as a regulon of a few dozen.

Usage

build_regulons(
  x,
  use_leading_edge = TRUE,
  add_tf = TRUE,
  min_genes = 10L,
  mode = c("activating", "repressing", "both"),
  .verbose = TRUE
)

## S3 method for class 'ScenicGrn'
build_regulons(
  x,
  use_leading_edge = TRUE,
  add_tf = TRUE,
  min_genes = 10L,
  mode = c("activating", "repressing", "both"),
  .verbose = TRUE
)

Arguments

Value

A named list of character vectors, one per regulon.

References

Aibar, et al., Nat Methods, 2017

Examples

# turn the signed TF to gene table into regulon gene sets
sc <- demo_single_cells()
grn <- scenic_grn_sc(
  sc,
  tf_ids = sprintf("gene_%02d", 1:5),
  scenic_params = params_scenic(
    min_counts = 1L,
    learner_params = list(n_trees = 20L)
  ),
  .verbose = FALSE
)
grn <- identify_tf_to_genes(
  grn,
  method = "top_k",
  k_tfs = 3L,
  .verbose = FALSE
)
grn <- tf_to_genes_correlations(grn, object = sc, .verbose = FALSE)
lengths(build_regulons(
  grn,
  use_leading_edge = FALSE,
  min_genes = 3L,
  .verbose = FALSE
))

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Runs sparse PCA over the provided HVGs, runs Harmony for batch correction, and compresses the result into the cached terms used at query time. The Harmony version is auto-detected from the class of harmony_params. Optionally snapshots one or more obs columns (e.g. cell type annotations) into the reference's labels slot for downstream label transfer.

Usage

build_symphony_ref(
  object,
  batch_column,
  additional_batch_columns = NULL,
  hvg,
  harmony_params = params_sc_harmony(),
  pca_params = params_sc_pca(),
  no_pcs = 30L,
  slim = FALSE,
  label_columns = NULL,
  seed = 42L,
  .verbose = TRUE
)

Arguments

Value

A SymphonyReference object.

Examples

# two batch reference, cell type labels snapshotted along the way
ref <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 500L,
    n_genes = 50L,
    n_batches = 2L,
    batch_effect_strength = "medium"
  )
)
symphony_ref <- build_symphony_ref(
  ref,
  batch_column = "batch_index",
  hvg = get_hvg(ref) + 1L,
  harmony_params = params_sc_harmony(k = 10L),
  no_pcs = 10L,
  label_columns = "cell_grp",
  .verbose = FALSE
)
symphony_ref

unlink(ref@dir_data, recursive = TRUE, force = TRUE)

Description

This is a helper function that can be used to create a BulkDge object (see BulkDge()) directly from h5ad objects.

Usage

bulk_dge_from_h5ad(h5_path, .verbose = TRUE)

Arguments

Value

BulkDge object.

Examples

# round trip a synthetic count matrix through a temporary h5ad
syn <- synthetic_bulk_cor_matrix()
h5_path <- tempfile(fileext = ".h5ad")
write_h5ad_sc_dense(
  f_path = h5_path,
  counts = t(syn$counts),
  obs = data.table::data.table(sample_id = colnames(syn$counts)),
  var = data.table::data.table(var_id = rownames(syn$counts)),
  .verbose = FALSE
)
object <- bulk_dge_from_h5ad(h5_path, .verbose = FALSE)
object

unlink(h5_path)

Description

Class for applying various co-expression module detection methods on top of bulk RNAseq data.

Usage

BulkCoExp(raw_data, meta_data, variable_info = NULL)

Arguments

Value

Returns the BulkCoExp class for further operations.

Properties

raw_data

A numerical matrix of the provided raw data.

meta_data

A data.table with the meta-information about the samples.

variable_info

An optional data.table containing the variable info.

processed_data

A list in which various types of processed data will be stored.

outputs

A list in which key outputs will be stored.

params

A (nested) list that will store all the parameters of the applied function.

final_results

A data.table that will contain the final results.

Examples

# co-expression class over synthetic bulk counts (samples x genes)
syn <- synthetic_bulk_cor_matrix()
mat <- log1p(t(syn$counts))
meta <- data.table::data.table(sample_id = rownames(mat))
object <- BulkCoExp(raw_data = mat, meta_data = meta)
object

Description

Class for coordinating differential gene expression analyses with subsequent GSE in a structured format. The filtered counts, library sizes and normalisation factors are stored in the class; get_dge_list() turns them into an edgeR DGEList on demand.

Usage

BulkDge(
  raw_counts,
  meta_data,
  variable_info = NULL,
  alternative_gene_id = NULL
)

Arguments

Value

Returns the BulkDge class for further operations.

Properties

raw_counts

A numerical matrix of the provided raw data.

meta_data

A data.table with the meta-information about the samples.

variable_info

An optional data.table containing the variable info.

outputs

A list in which key outputs will be stored.

plots

A list with the plots that are generated during subsequent QC steps.

params

A (nested) list that will store all the parameters of the applied function.

final_results

A list in which final results will be stored.

Examples

# DGE class over synthetic bulk counts (genes x samples)
syn <- synthetic_bulk_cor_matrix()
meta <- data.table::data.table(
  sample_id = colnames(syn$counts),
  case_control = rep(c("case", "control"), each = 50)
)
object <- BulkDge(raw_counts = syn$counts, meta_data = meta)
object

Description

This function will plot various alphas to highlight the most interesting alpha parameters akin to the implementation of contrastive PCA in Python.

Usage

c_pca_plot_alphas(
  object,
  label_column = NULL,
  min_alpha = 0.1,
  max_alpha = 100,
  n_alphas = 10L,
  .verbose = TRUE
)

Arguments

Value

A ggplot showing the impact of various alpha parameters on the samples in form of 2D plots.

References

Abid, et al., Nature Communications, 2018

Examples

# how the group separation responds to alpha
cpca_data <- synthetic_c_pca_data()
target <- t(cpca_data$target)
background <- t(cpca_data$background)
meta <- data.table::data.table(
  sample_id = rownames(target),
  grp = cpca_data$target_labels
)
obj <- BulkCoExp(target, meta)
obj <- preprocess_bulk_coexp(obj, .verbose = FALSE)
obj <- contrastive_pca_processing(obj, background, .verbose = FALSE)
c_pca_plot_alphas(obj, label_column = "grp", n_alphas = 6L, .verbose = FALSE)

Description

Rust-based version of blitzGSEA, see Lachmann, et al. Instead of permuting per pathway, it calibrates a gamma null once for the signature and then reads each pathway's p-value straight off the fitted tail. That makes the cost independent of how many pathways you throw at it, which is where it pulls away from permutation methods on large libraries.

Hand back a null_model from blitzgsea_calibrate() to reuse a calibration across libraries. Left as NULL, one is calibrated on the fly.

Usage

calc_blitzgsea(
  stats,
  pathways,
  blitz_params = params_blitzgsea(),
  null_model = NULL
)

Arguments

Value

A data.table with the results from the blitzGSEA with the following columns:

• pathway_name - Character. The name of the pathway.

• es - Float. The enrichment score for this pathway.

• nes - Float. The normalised enrichment score, the signed normal quantile of the one-sided tail probability.

• pvals - Float. The two-sided p-value from the gamma approximation.

• fdr - Float. The Benjamini-Hochberg adjusted p-value.

• sidak - Float. The Sidak adjusted p-value.

• size - Integer. The size of the pathway after intersection.

• leading_edge - List of character vectors with the leading edge genes.

References

Lachmann, et al., Bioinformatics, 2022

Examples

# blitzGSEA reading p-values off a fitted gamma null
set.seed(42L)
stats <- stats::setNames(rnorm(500), sprintf("gene_%03i", 1:500))
stats[1:20] <- stats[1:20] + 2
pathways <- list(
  up_set = sprintf("gene_%03i", 1:20),
  random_set = sprintf("gene_%03i", 200:240)
)
res <- calc_blitzgsea(
  stats,
  pathways,
  blitz_params = params_blitzgsea(permutations = 1000L, anchors = 10L)
)
res[, c("pathway_name", "es", "nes", "pvals")]

Description

To leverage the quality metrics from Persad, et al., we need the diffusion coordinates to then calculate if a cell is a dense or sparse region of the manifold, its compactness and separation to other meta cells. To do so, generate a diffusion map on the original data based on the approach of SEACells and add the data to the object, see Persad, et al.

Usage

calc_diffusion_coordinates(
  object,
  knn_data,
  n_dcs = 10L,
  k_density = 150L,
  seed = 42L,
  .verbose = TRUE
)

Arguments

Value

The class with the diffusion map coordinates, density distance and region attached.

References

Persad, et al. Nat Biotechnol, 2023

Examples

# diffusion map off the source kNN graph, giving each meta cell a region
sc <- demo_single_cells()
mc <- generate_bt_meta_cells_sc(
  sc,
  sc_meta_cell_params = params_sc_bt_metacells(target_no_metacells = 50L),
  .verbose = FALSE
)
mc <- calc_diffusion_coordinates(
  mc,
  knn_data = get_knn_obj(sc),
  .verbose = FALSE
)
table(mc[["density_region"]]$density_region)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Rust-based version of the fgsea simple and multi-level algorithm. Initially, the simple method is run. For low p-values, the multi-level method is used to estimate lower p-values than possible just based on the permutations, see Korotkevich, et al.

Usage

calc_fgsea(
  stats,
  pathways,
  nperm = 1000L,
  gsea_params = params_gsea(),
  seed = 123L
)

Arguments

Value

A data.table with the results from the GSEA with the following columns:

• es - Float. The enrichment score for this pathway.

• nes - Float. The normalised enrichment score for this pathway.

• pvals - Float. The p-value for this pathway.

• n_more_extreme - Integer. Number of permutation that had more extreme enrichment scores than the actual.

• size - Integer. The size of the pathway.

• pathway_name - Character. The name of the pathway.

• leading_edge - List of character vectors with the leading edge genes.

• fdr - Float. The adjusted pval.

References

Korotkevich, et al., bioRxiv

Examples

# fgsea with the multi-level refinement of low p-values
set.seed(42L)
stats <- stats::setNames(rnorm(500), sprintf("gene_%03i", 1:500))
stats[1:20] <- stats[1:20] + 2
pathways <- list(
  up_set = sprintf("gene_%03i", 1:20),
  random_set = sprintf("gene_%03i", 200:240)
)
res <- calc_fgsea(stats, pathways, nperm = 1000L)
res[, c("pathway_name", "es", "pvals")]

Description

Rust-based version of the fgsea simple algorithm. This one is permutation- based and similar to the traditional implementation, but leverages some clear tricks to be way faster, see Korotkevich, et al.

Usage

calc_fgsea_simple(
  stats,
  pathways,
  nperm = 2000L,
  gsea_params = params_gsea(),
  seed = 123L
)

Arguments

Value

A data.table with the results from the GSEA with the following columns:

• es - Float. The enrichment score for this pathway.

• nes - Float. The normalised enrichment score for this pathway.

• pvals - Float. The p-value for this pathway.

• n_more_extreme - Integer. Number of permutation that had more extreme enrichment scores than the actual.

• size - Integer. The size of the pathway.

• pathway_name - Character. The name of the pathway.

• leading_edge - List of character vectors with the leading edge genes.

• fdr - Float. The adjusted pval.

References

Korotkevich, et al., bioRxiv

Examples

# fgsea simple against a signature with one enriched set
set.seed(42L)
stats <- stats::setNames(rnorm(500), sprintf("gene_%03i", 1:500))
stats[1:20] <- stats[1:20] + 2
pathways <- list(
  up_set = sprintf("gene_%03i", 1:20),
  random_set = sprintf("gene_%03i", 200:240)
)
res <- calc_fgsea_simple(stats, pathways, nperm = 1000L)
res[, c("pathway_name", "es", "pvals")]

Description

Rust-based version of the traditional permutation-based GSEA algorithm.

Usage

calc_gsea_traditional(
  stats,
  pathways,
  nperm = 2000L,
  gsea_params = params_gsea(),
  seed = 123L
)

Arguments

Value

A data.table with the results from the GSEA with the following columns:

• es - Float. The enrichment score for this pathway.

• nes - Float. The normalised enrichment score for this pathway.

• pvals - Float. The p-value for this pathway.

• n_more_extreme - Integer. Number of permutation that had more extreme enrichment scores than the actual.

• size - Integer. The size of the pathway.

• pathway_name - Character. The name of the pathway.

• leading_edge - List of character vectors with the leading edge genes.

• fdr - Float. The adjusted pval.

Examples

# permutation GSEA against a signature with one enriched set
set.seed(42L)
stats <- stats::setNames(rnorm(500), sprintf("gene_%03i", 1:500))
stats[1:20] <- stats[1:20] + 2
pathways <- list(
  up_set = sprintf("gene_%03i", 1:20),
  random_set = sprintf("gene_%03i", 200:240)
)
res <- calc_gsea_traditional(stats, pathways, nperm = 1000L)
res[, c("pathway_name", "es", "pvals")]

Description

Implementation of the bixverse version of the gene set variation analysis (GSVA), see Hänzelmann, et al.

Usage

calc_gsva(
  exp,
  pathways,
  kernel = c("gaussian", "poisson", "none"),
  gaussian = deprecated(),
  gsva_params = params_gsva(),
  .verbose = FALSE
)

Arguments

Value

A matrix of shape pathways (that passed the thresholds) x samples.

References

see Hänzelmann, et al. Bmc Bioinformatics, 2013

Examples

# per-sample GSVA scores for two gene sets
set.seed(123L)
exp_mat <- matrix(
  rnorm(200 * 10),
  nrow = 200,
  dimnames = list(sprintf("gene_%03i", 1:200), sprintf("sample_%i", 1:10))
)
pathways <- list(
  set_a = sprintf("gene_%03i", 1:20),
  set_b = sprintf("gene_%03i", 50:80)
)
dim(calc_gsva(exp_mat, pathways))

Description

Helper function to compare the results of two SingleCellNearestNeighbour against each other. The first one can serve as a reference (ground truth) and you can compare against the second one.

Usage

calc_knn_metrics(ref_knn, query_knn)

Arguments

Value

A list with:

• matches - The intersecting indices between the reference and query kNN for each sample. In an ideal match up should be equal to k.

• distance_ratio - The distance ratio. Calculates sum(dist_query) / sum(dist_ref) per sample. Indicates how much worse the reference is.

• final_recall - The final recall across all samples.

• final_ratio - The final distance ratio across all samples.

Examples

# recall of an annoy index against the default one
set.seed(42L)
embd <- matrix(rnorm(500 * 10), nrow = 500)
rownames(embd) <- sprintf("cell_%03d", 1:500)
ref <- generate_sc_knn(embd, .verbose = FALSE)
query <- generate_sc_knn(
  embd,
  neighbours_params = params_sc_neighbours(
    knn = list(knn_method = "annoy")
  ),
  .verbose = FALSE
)
calc_knn_metrics(ref, query)$final_recall

Description

This function will calculate the compactness and separation of your metacells on the manifold (defined by the diffusion map). You must have run calc_diffusion_coordinates() before calling this function. The idea is that compactness indicates how tight the metacell spans the manifold, whereas separation indicates how well the different metacells span the manifold.

Usage

calc_manifold_metrics(object)

Arguments

Value

The class with the compactness and separation scores added.

References

Persad, et al. Nat Biotechnol, 2023

Examples

# compactness and separation on top of the diffusion coordinates
sc <- demo_single_cells()
mc <- generate_bt_meta_cells_sc(
  sc,
  sc_meta_cell_params = params_sc_bt_metacells(target_no_metacells = 50L),
  .verbose = FALSE
)
mc <- calc_diffusion_coordinates(
  mc,
  knn_data = get_knn_obj(sc),
  .verbose = FALSE
)
mc <- calc_manifold_metrics(mc)
head(mc[[c("compactness", "separation")]])

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Scores a meta cell assignment for compactness, i.e. the mean variance across the diffusion components over the cells of each meta cell. Lower is better. Unlike calc_manifold_metrics(), this works on any assignment, so meta cells generated outside bixverse (or on a different object) can be scored against an existing diffusion map.

Usage

calc_meta_cell_compactness(dcs, assignment)

Arguments

Value

Named numeric vector with one compactness value per meta cell, named and ordered like split() orders assignment.

References

Persad, et al., Nat Biotechnol, 2023

Examples

set.seed(1L)
dcs <- matrix(rnorm(200L), ncol = 2L)
assignment <- rep(sprintf("mc_%i", 1:10), each = 10L)
calc_meta_cell_compactness(dcs, assignment)

Description

A potential metric to see how well the meta cells are aggregated is their cell type purity. This helper function helps to plot the meta-cell purity based on annotated cell types. These can be also just unsupervised memberships to graph-based clustering, etc.

Usage

calc_meta_cell_purity(
  object,
  original_cell_type,
  add_additional_info = c("none", "top_label", "top_two_labels"),
  add_entropy = FALSE
)

Arguments

Details

Ties for the top (or second) label resolve to whichever label sorts first, as the labels are factorised internally. The entropy is the Shannon entropy of the label distribution within a meta cell, divided by ⁠log(<number of distinct labels in original_cell_type>)⁠, so it sits in ⁠[0, 1]⁠ and stays comparable between meta cells. It is 0 if there is only a single label in the data.

Value

The MetaCells with added columns to the observation table:

• mc_purity - Fraction of the meta cell's cells that carry the most abundant label. Always added.

• mc_top_label - Name of the most abundant label. Added for add_additional_info %in% c("top_label", "top_two_labels").

• mc_second_label - Name of the second most abundant label, NA for a pure meta cell. Added for add_additional_info = "top_two_labels".

• mc_second_frac - Fraction of the meta cell's cells carrying that second label, 0 for a pure meta cell. Added for add_additional_info = "top_two_labels".

• mc_entropy - Normalised Shannon entropy of the label distribution, see details. Added for add_entropy = TRUE.

Examples

# purity of the meta cells against the cell types they were built from
sc <- demo_single_cells()
mc <- generate_bt_meta_cells_sc(
  sc,
  sc_meta_cell_params = params_sc_bt_metacells(target_no_metacells = 50L),
  .verbose = FALSE
)
cell_types <- sc[["cell_grp"]]$cell_grp
mc <- calc_meta_cell_purity(mc, original_cell_type = cell_types)
head(mc[[c("meta_cell_id", "mc_purity")]])

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Scores a meta cell assignment for separation, i.e. the Euclidean distance from each meta cell centroid in diffusion space to the nearest other meta cell centroid. Higher is better. Unlike calc_manifold_metrics(), this works on any assignment, so meta cells generated outside bixverse (or on a different object) can be scored against an existing diffusion map.

Usage

calc_meta_cell_separation(dcs, assignment)

Arguments

Value

Named numeric vector with one separation value per meta cell, named and ordered like split() orders assignment. A single meta cell yields Inf.

References

Persad, et al., Nat Biotechnol, 2023

Examples

set.seed(1L)
dcs <- matrix(rnorm(200L), ncol = 2L)
assignment <- rep(sprintf("mc_%i", 1:10), each = 10L)
calc_meta_cell_separation(dcs, assignment)

Description

Rust-based version of the mitch multi-contrast enrichment, see Kaspi and Ziemann. Takes in a matrix representing the contrasts you wish to test against.

Usage

calc_mitch(contrast_mat, gene_set_list, min_size = 5L)

Arguments

Value

A data.table with the Mitch enrichment results.

References

Kaspi and Ziemann, Bmc Genomics, 2020

Examples

# multi-contrast enrichment over two contrasts
set.seed(10L)
contrast_mat <- matrix(
  rnorm(300 * 2),
  ncol = 2,
  dimnames = list(sprintf("gene_%03i", 1:300), c("contrast_a", "contrast_b"))
)
contrast_mat[1:30, ] <- contrast_mat[1:30, ] + 1.5
gene_sets <- list(
  hit_set = sprintf("gene_%03i", 1:30),
  bg_set = sprintf("gene_%03i", 100:150)
)
res <- calc_mitch(contrast_mat, gene_sets)
res[, c("pathway_names", "manova_pval", "s_dist")]

Description

Implements the approach from

Usage

calc_sc_type_scores(
  object,
  cell_marker_list,
  sensitivity = TRUE,
  weight_floor = NULL,
  .verbose = TRUE
)

Arguments

Value

An ScTypeResults results class

Examples

# ScType scores against the planted marker genes
sc <- demo_single_cells()
markers <- data.table::data.table(
  cell_type = rep(sprintf("cell_type_%i", 1:3), each = 10),
  gene_id = sprintf("gene_%02d", 1:30)
)
cell_markers <- prepare_cell_markers(obj = sc, marker_df = markers)
res <- calc_sc_type_scores(
  sc,
  cell_marker_list = cell_markers,
  .verbose = FALSE
)
dim(get_scores(res))

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Scores all samples against one up-regulated gene set with an optional paired down-regulated set. When n_permutations > 0, additionally runs a permutation test and includes empirical p-values.

Usage

calc_singscore(
  ranks,
  up_set,
  down_set = NULL,
  center_score = TRUE,
  known_direction = TRUE,
  n_permutations = 0L,
  seed = 42L
)

Arguments

Value

A data.table with one row per sample. Columns: total_score, total_dispersion, optionally up_score, up_dispersion, down_score, down_dispersion, and (when permutations are run) pval. With permutations, the null distribution is attached as attr(., "null_distribution").

References

Foroutan et al., BMC Bioinformatics, 2018.

Examples

# score every sample against one up-regulated gene set
set.seed(123L)
exp_mat <- matrix(
  rnorm(200 * 10),
  nrow = 200,
  dimnames = list(sprintf("gene_%03i", 1:200), sprintf("sample_%i", 1:10))
)
ranks <- calc_singscore_rank(exp_mat)
head(calc_singscore(ranks, up_set = sprintf("gene_%03i", 1:20)))

Description

Scores all samples against many up-regulated gene sets with optional paired down-regulated sets. Down sets are paired with up sets by name; unmatched names are dropped.

Usage

calc_singscore_multi(
  ranks,
  up_pathways,
  down_pathways = NULL,
  center_score = TRUE,
  known_direction = TRUE,
  min_size = 1L,
  max_size = 500L
)

Arguments

Value

A named list with two matrices, scores and dispersions, each of shape gene_sets × samples.

References

Foroutan et al., BMC Bioinformatics, 2018.

Examples

# score every sample against several up-regulated gene sets
set.seed(123L)
exp_mat <- matrix(
  rnorm(200 * 10),
  nrow = 200,
  dimnames = list(sprintf("gene_%03i", 1:200), sprintf("sample_%i", 1:10))
)
ranks <- calc_singscore_rank(exp_mat)
pathways <- list(
  set_a = sprintf("gene_%03i", 1:20),
  set_b = sprintf("gene_%03i", 50:80)
)
res <- calc_singscore_multi(ranks, up_pathways = pathways)
dim(res$scores)

Description

Column-ranks an expression matrix. If stable_genes is provided, uses the stable-genes ranking method (unit-normalised ranks against a small set of stable genes); otherwise standard column ranks. Sets attr(., "stable") so downstream singscore functions know which bounds formula to use.

Usage

calc_singscore_rank(exp, stable_genes = NULL)

Arguments

Value

A rank matrix with the same shape as exp and attr(., "stable") set to TRUE or FALSE.

References

1.) Foroutan et al., BMC Bioinformatics, 2018.; 2.) Bhuva, et al., Nucleic Acids Res., 2020

Examples

# column ranks of an expression matrix, ready for singscore
set.seed(123L)
exp_mat <- matrix(
  rnorm(200 * 10),
  nrow = 200,
  dimnames = list(sprintf("gene_%03i", 1:200), sprintf("sample_%i", 1:10))
)
ranks <- calc_singscore_rank(exp_mat)
dim(ranks)

Description

Implementation of the bixverse version of the single sample gene set enrichment analysis (ssGSEA), see Barbie et al.

Usage

calc_ssgsea(exp, pathways, ssgsea_params = params_ssgsea(), .verbose = FALSE)

Arguments

Value

A matrix of shape pathways (that passed the thresholds) x samples.

References

Barbie et al., Nature, 2009

Examples

# per-sample ssGSEA scores for two gene sets
set.seed(123L)
exp_mat <- matrix(
  rnorm(200 * 10),
  nrow = 200,
  dimnames = list(sprintf("gene_%03i", 1:200), sprintf("sample_%i", 1:10))
)
pathways <- list(
  set_a = sprintf("gene_%03i", 1:20),
  set_b = sprintf("gene_%03i", 50:80)
)
round(calc_ssgsea(exp_mat, pathways)[, 1:3], 3)

Description

This is a deprecated method and will raise an error. Please use calculate_dge_limma() and calculate_dge_hedges().

Usage

calculate_all_dges(
  object,
  contrast_column,
  contrast_list = NULL,
  filter_column = NULL,
  co_variates = NULL,
  small_sample_correction = NULL,
  ...,
  .verbose = TRUE
)

Arguments

Value

Throws an error, as it is deprecated.

Description

Computes the average silhouette width (ASW) on batch labels using pairwise Euclidean distances in the embedding space. For each cell, the function estimates the mean distance to cells of the same batch (a) and the mean distance to cells of the nearest other batch (b), then computes s = (b - a) / max(a, b).

Values near 0 indicate good batch mixing, values near 1 indicate batch separation, and negative values suggest overcorrection. This metric is best suited for embedding-based correction methods (e.g. Harmony, fastMNN). For graph-based methods like BBKNN, consider using calculate_lisi_sc() instead.

Usage

calculate_batch_asw_sc(
  object,
  batch_column,
  embd_to_use = "pca",
  max_cells = 5000L,
  seed = 42L,
  .verbose = TRUE
)

Arguments

Value

A BatchSilhouetteScores object with the following elements

• per_cell - Per-cell silhouette scores in ⁠[-1, 1]⁠.

• mean_asw - Mean silhouette width across all cells.

• median_asw - Median silhouette width across all cells.

• n_batches - Number of batches in the data.

• embedding_used - Which embedding the ASW was computed on.

Examples

# batch silhouette width on the PCA embedding
sc <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 600L, n_genes = 50L, n_batches = 3L
  )
)
calculate_batch_asw_sc(
  sc,
  batch_column = "batch_index",
  .verbose = FALSE
)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Average silhouette width on cell type labels in the embedding, rescaled to ⁠[0, 1]⁠ via (s + 1) / 2 as in scIB. Higher values mean cell types stay separated after correction. Counterpart to calculate_batch_asw_sc() on the biology side.

Usage

calculate_cell_type_asw_sc(
  object,
  cell_type_column,
  embd_to_use = "pca",
  max_cells = 5000L,
  seed = 42L,
  .verbose = TRUE
)

Arguments

Value

A CellTypeAswScores object with the following elements

• per_cell - Per-cell rescaled silhouette scores in ⁠[0, 1]⁠.

• mean_asw - Mean rescaled silhouette width.

• median_asw - Median rescaled silhouette width.

• n_cell_types - Number of cell types.

• embedding_used - Which embedding the ASW was computed on.

References

Luecken, et al., Nat. Methods, 2022

Examples

# cell type silhouette width on the PCA embedding
sc <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 600L, n_genes = 50L, n_batches = 3L
  )
)
calculate_cell_type_asw_sc(
  sc,
  cell_type_column = "cell_grp",
  .verbose = FALSE
)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

This function calculates the critical value for a given ontology similarity matrix.

Usage

calculate_critical_value(x, alpha, permutations = 100000L, seed = 10101L)

Arguments

Value

The critical value.

Examples

# critical value of a Wang similarity matrix at alpha 0.1
onto <- data.table::data.table(
  parent = c("a", "b", "b", "b", "c"),
  child = c("b", "c", "d", "e", "f"),
  type = c("part_of", "part_of", "part_of", "is_a", "is_a")
)
weights <- c(part_of = 0.8, is_a = 0.6)
sim_mat <- calculate_wang_sim_mat(onto, weights = weights)
calculate_critical_value(sim_mat, alpha = 0.1)

Description

This function will calculate the Hedge's G effect size on the normalised counts. Should batch-corrected counts be found, these will be used. At a minimum you will need to provide contrast_column that can be found in the meta-data. If you do not provide a vector of contrasts that you wish to test for, every permutation of groups represented in that column will be tested against each other.

Usage

calculate_dge_hedges(
  object,
  contrast_column,
  contrast_list = NULL,
  filter_column = NULL,
  .verbose = TRUE
)

Arguments

Value

Returns the class with additional data added to the outputs.

Examples

# effect sizes over every contrast in the case_control column
syn <- synthetic_bulk_cor_matrix()
meta <- data.table::data.table(
  sample_id = colnames(syn$counts),
  case_control = rep(c("case", "control"), each = 50)
)
object <- BulkDge(raw_counts = syn$counts, meta_data = meta)
object <- qc_bulk_dge(object, group_col = "case_control", .verbose = FALSE)
object <- normalise_bulk_dge(
  object,
  group_col = "case_control",
  .verbose = FALSE
)
object <- calculate_dge_hedges(
  object,
  contrast_column = "case_control",
  .verbose = FALSE
)
head(get_dge_effect_sizes(object))

Description

This function will apply the Limma Voom DGE workflow. At a minimum you will need to provide contrast_column that can be found in the meta-data. If you do not provide a vector of contrasts that you wish to test for, every permutation of groups represented in that column will be tested against each other.

Usage

calculate_dge_limma(
  object,
  contrast_column,
  contrast_list = NULL,
  filter_column = NULL,
  co_variates = NULL,
  limma_params = params_limma_voom(),
  .verbose = TRUE
)

Arguments

Value

Returns the class with additional data added to the outputs.

Examples

# limma voom over every contrast in the case_control column
syn <- synthetic_bulk_cor_matrix()
meta <- data.table::data.table(
  sample_id = colnames(syn$counts),
  case_control = rep(c("case", "control"), each = 50)
)
object <- BulkDge(raw_counts = syn$counts, meta_data = meta)
object <- qc_bulk_dge(object, group_col = "case_control", .verbose = FALSE)
object <- normalise_bulk_dge(
  object,
  group_col = "case_control",
  .verbose = FALSE
)
object <- calculate_dge_limma(
  object,
  contrast_column = "case_control",
  .verbose = FALSE
)
head(get_dge_limma_voom(object))

Description

This functions can take a given NetworkDiffusions object and calculates an AUC and generates a Z-score based on random permutation of random_aucs for test for statistical significance if desired.

Usage

calculate_diffusion_auc(
  object,
  hit_nodes,
  auc_iters = 10000L,
  random_aucs = 1000L,
  permutation_test = FALSE,
  seed = 42L
)

Arguments

Value

List with AUC and Z-score as the two named elements if permutations test set to TRUE; otherwise just the AUC.

Examples

# AUROC of the diffusion score against two known hit nodes
set.seed(42)
g <- igraph::sample_pa(15, directed = FALSE)
edges <- data.table::setDT(igraph::as_data_frame(g))[, `:=`(
  from = sprintf("node_%i", from),
  to = sprintf("node_%i", to)
)]
object <- NetworkDiffusions(edges, weighted = FALSE, directed = FALSE)
object <- diffuse_seed_nodes(object, c(node_1 = 1, node_3 = 1), "max")
calculate_diffusion_auc(object, hit_nodes = c("node_2", "node_4"))

Description

This function takes two matrices in and calculate on a per column basis the Hedge's G effect size and the standard error. These results can be subsequently used for meta-analyses or other approaches.

Usage

calculate_effect_size(
  mat_a,
  mat_b,
  small_sample_correction = NULL,
  .verbose = TRUE
)

Arguments

Value

x, robustly scaled.

Examples

# Hedge's G between two groups of ten samples
set.seed(42)
mat_a <- matrix(rnorm(100), nrow = 10, ncol = 10)
mat_b <- matrix(rnorm(100, mean = 1), nrow = 10, ncol = 10)
colnames(mat_a) <- colnames(mat_b) <- sprintf("gene_%i", 1:10)
res <- calculate_effect_size(mat_a, mat_b, .verbose = FALSE)
head(res$effect_sizes)

Description

For each cell type, restricts the kNN graph to that cell type and takes the fraction of its cells in the largest connected component. A cell type split across batches after correction falls apart into several components and scores low. 1 means every cell type is one connected piece.

Usage

calculate_graph_connectivity_sc(object, cell_type_column, .verbose = TRUE)

Arguments

Value

A GraphConnectivityScores object with the following elements

• per_cell_type - Named numeric. Connectivity per cell type.

• mean_connectivity - Mean connectivity across cell types.

• median_connectivity - Median connectivity across cell types.

References

Luecken, et al., Nat. Methods, 2022

Examples

# connectivity of each cell type in the kNN graph
sc <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 600L, n_genes = 50L, n_batches = 3L
  )
)
calculate_graph_connectivity_sc(sc, cell_type_column = "cell_grp")

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

this function will calculate the information content of each provided term based on a list of ancestors, which is a named list of terms with ancestor identifiers as their values. Can be calculated using get_ontology_ancestry(). The information content is calculated as ⁠-log2(number descendant/total terms in the ontology)⁠. More information can be found here.

Usage

calculate_information_content(ancestor_list)

Arguments

Value

A named list of each term and their information content as values.

Examples

# information content from the descendant counts of a toy ontology
onto <- data.table::data.table(
  parent = c("a", "b", "b", "b", "c"),
  child = c("b", "c", "d", "e", "f")
)
ancestry <- get_ontology_ancestry(onto)
unlist(calculate_information_content(ancestry$descendants))

Description

Runs the batch mixing and (if cell type labels are given) the biological conservation metrics in one go and returns one row per call, so results across correction methods can be rbind-ed into one table. Every column is on ⁠[0, 1]⁠ (PCR comparison can go negative) and higher is better, the scIB convention:

Batch mixing:

• kbet_accept - ⁠1 - ⁠ kBET rejection rate.

• batch_asw - ⁠mean(1 - |s|)⁠ over the per-cell batch silhouettes.

• ilisi - Normalised iLISI.

• pcr_comparison - (pre - post) / pre of the batch PCR.

Biological conservation:

• clisi - Normalised cLISI.

• cell_type_asw - Rescaled cell type silhouette width.

• graph_connectivity - Mean graph connectivity over cell types.

The kNN metrics read the kNN graph currently stored in the object, so recompute the neighbours on the corrected embedding first. Embedding metrics are NA if embd_to_use = NULL (e.g. BBKNN, which only returns a graph).

Usage

calculate_integration_metrics_sc(
  object,
  batch_column,
  cell_type_column = NULL,
  embd_to_use = "pca",
  max_cells = 5000L,
  seed = 42L,
  .verbose = TRUE
)

Arguments

Value

A one-row data.table with the columns embedding, kbet_accept, batch_asw, ilisi, pcr_comparison, clisi, cell_type_asw and graph_connectivity.

References

Luecken, et al., Nat. Methods, 2022

Examples

# all metrics on the uncorrected PCA
sc <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 600L, n_genes = 50L, n_batches = 3L
  )
)
calculate_integration_metrics_sc(
  sc,
  batch_column = "batch_index",
  cell_type_column = "cell_grp",
  .verbose = FALSE
)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

This function calculates the k-nearest neighbour batch-effect test (kBET). Briefly, the function leverages a Chi-Square statistic to calculate the differences in batch proportions observed in the neighbourhood of a given cell with the overall batch proportions. If the test is significant for that cell it indicates poor mixing for that cell specifically. Large number of positive tests indicate bad mixing overall. For more details, please see Büttner et al.

Usage

calculate_kbet_sc(object, batch_column, threshold = 0.05, .verbose = TRUE)

Arguments

Value

A KbetScores object with the following elements

• kbet_score - Proportion of significant tests over all cells. 0 indicates perfect mixing, 1 indicates no mixing between batches.

• significant_tests - Logical vector indicating for which cells the test was below the threshold.

• p_values - The p-values from the Chi-Square test.

• chi_square_stats - Per-cell Chi-Square statistics.

• mean_chi_square - Mean Chi-Square statistic across all cells.

• median_chi_square - Median Chi-Square statistic across all cells.

• threshold - The significance threshold used.

• n_batches - Number of batches in the data.

References

Büttner, et al., Nat. Methods, 2019

Examples

# kBET rejection rate over three batches
sc <- demo_single_cells(
  syn_data_params = params_sc_synthetic_data(
    n_cells = 600L, n_genes = 50L, n_batches = 3L
  )
)
calculate_kbet_sc(sc, batch_column = "batch_index", .verbose = FALSE)

unlink(sc@dir_data, recursive = TRUE, force = TRUE)

Description

Computes the Local Inverse Simpson's Index (LISI) on the kNN graph: the effective number of labels in each cell's neighbourhood. On batch labels this is iLISI, where higher means better mixing. On cell type labels it is cLISI, where lower means cell types stay apart. Unlike kBET, LISI does not compare against global proportions, so it also works on graph-based corrections like BBKNN.

The normalised score follows scIB and lands in ⁠[0, 1]⁠, higher is better for both: iLISI as (median - 1) / (n - 1), cLISI as (n - median) / (n - 1).

Usage

calculate_lisi_sc(
  object,
  label_column,
  type = c("batch", "cell_type"),
  weighted = FALSE,
  perplexity = 30,
  .verbose = TRUE
)

Arguments

Value

A LisiScores object with the following elements

• per_cell - Per-cell LISI scores in ⁠[1, n_labels]⁠.

• mean_lisi - Mean LISI across all cells.

• median_lisi - Median LISI across all cells.

• lisi_norm - The normalised score in ⁠[0, 1]⁠, higher is better.

• n_labels - Number of distinct labels.

• type - "batch" (iLISI) or "cell_type" (cLISI).