Algorithm validation

multiomeR reimplements a small number of reference algorithms so they can operate on the workflow’s native matrices and graph state. The algorithms themselves are described on the methods pages. This page records why each was reimplemented, where it deliberately differs from its reference, and what the executable validation establishes.

The evidence labels are intentionally narrow:

Passing these fixtures does not validate every dataset, parameter regime, approximate-neighbor realization, biological interpretation, or downstream target. Apart from a few checks of configuration and input handling, the test suite contains only such reference comparisons. Tests that compare with an installed package assert the version they were written against, and pixi.lock locks the versions. Run the complete suite locally with pixi run --use-environment-activation-cache test after changing tests, tested helpers, native sources, or the Pixi environment; cloud CI runs only checks that need no environment. The narrower pixi run --use-environment-activation-cache test-algorithm-validation task runs only the slow UCell, AMULET, WNN, and SCAVENGE tests.

Implementation Evidence status Reference Validation contract
BPCells-native UCell Reference-parity tested UCell Identical values, dimensions, and dimnames
BPCells-native AMULET Reference-parity tested scDblFinder Identical metrics and loci, including order
BPCells-backed ATAC scDblFinder feature aggregation Not reference-tested scDblFinder None; scores and calls can differ from its internal aggregation
Native WNN Reference-similarity tested Seurat Modality-weight Spearman and neighbor-overlap thresholds per fixture
Sparse SCAVENGE propagation Algorithmically derived and reference-parity tested SCAVENGE source at 8ee8b173d965 Closed-form propagation; identical seeds; propagation and trait relevance scores within 1e-12
Peak–gene donor-slope REML and Kenward–Roger kernels Reference-parity tested lme4 and pbkrtest Coefficients, df and P-values within 1e-6; identical fit statuses
Peak–gene compact BH breakpoints Reference-parity tested stats::p.adjust() Identical FDR per chromosome slice
S-LDSC LD scores and regression Not reference-tested ldsc 1.0.1 None in the suite; a development comparison is described below

The peak–gene rows belong to the analyses in Peak–gene correlation methods; their tests are test-peak-gene-hierarchical-parity.R and test-peak-gene-correlation-parity.R.

BPCells-native UCell scoring

Reference algorithm. UCell::ScoreSignatures_UCell() calculates per-cell signature scores from descending feature ranks, caps ranks at maxRank, combines positive and negative signatures, and clips negative combined scores to zero.

Reason for reimplementation. The workflow keeps gene-by-cell counts in BPCells-backed matrices. Materializing the complete matrix or building a Seurat object solely for marker scoring would discard that storage contract, so multiomeR ranks bounded cell chunks and returns metadata-ready scores directly.

Deliberate deviations and consequences. Only one cell chunk is materialized at a time, and optional fork workers operate across chunks; this changes memory and execution behavior but not the tested values. The helper returns a data frame instead of mutating a Seurat object. The target-level marker validator rejects configured genes missing from the reference, whereas the lower-level helper still exposes UCell’s impute and skip modes. The production annotation ranks each cell from its nonzero counts in a compiled kernel, because zero counts tie below the rank cap, and scores signed signatures per cell before aggregation, which costs more computation than a positive-only rank-summary shortcut. It adds the chunk sums in the order of the R implementation, so its results are identical to it.

Implementation. calculate_BPCells_UCell_scores_from_matrix() and rank_UCell_count_chunk() in pipeline_manuscript_figures/algorithm_parity_helpers.R score single cells for the parity comparisons; the cluster annotation is in R/cluster_annotation_helpers.R with its kernel src/cluster_UCell_chunk.cpp, and is described in Cell-type annotation.

Validation. tests/testthat/test-scoring-parity.R compares signed signatures, with imputed and skipped missing genes, on a deterministic BPCells fixture and requires identical() values, dimensions, and dimnames. It also runs the production annotation path on unsigned, signed and negative-only signatures and compares per-cell scores, cluster means and matched-control summaries with reference scores averaged within clusters, within 1e-12, and checks that chunking and fork workers leave them unchanged. The production cell-cycle scorer is compared with Seurat::CellCycleScoring(): phases are identical and scores agree within 1e-6, because BPCells normalizes counts at lower floating-point precision.

pixi run --use-environment-activation-cache test-scoring-parity

BPCells-native AMULET

Reference algorithm. scDblFinder::amulet() detects likely scATAC-seq doublets from the number of genomic loci covered by more than two fragments, after filtering fragment sizes and excluded regions and removing loci recurrently covered across many cells, and derives Poisson P-values with Benjamini–Hochberg correction.

Reason for reimplementation. The per-GEM-well workflow already stores Cell Ranger ATAC fragments as compressed BPCells directories. Passing the fragment file to scDblFinder materializes chromosome-scale GRanges objects; the local implementation streams the BPCells fragments and retains only one chromosome’s selected fragments while calculating coverage runs.

Deliberate deviations and consequences. Only unique-fragment operation is supported: BPCells fragments do not retain Cell Ranger’s PCR-duplicate counts, so requesting non-unique expansion fails explicitly. Cell Ranger’s inclusive end coordinate is shifted back by one base to match scDblFinder’s BED import. BPCells does not export its fragment iterator header, so the native helper mirrors that private interface, verifies the pinned BPCells commit before use, and compiles a small shared library in each worker. A BPCells upgrade must revalidate this interface and the parity fixture before updating the pin. The parity fixtures use scDblFinder’s default excluded regions; production excludes the mitochondrial and sex chromosomes in full and, for GRCh38, the AMULET authors’ repeat list.

Implementation. src/amulet_bpcells.cpp iterates fragments and calculates coverage runs; calculate_amulet_metrics_BPCells() in R/amulet_BPCells_helpers.R holds the interface check, high-overlap filtering and statistics; the amulet_metrics_tibble target consumes the prefixed BPCells fragments.

Validation. tests/testthat/test-amulet-parity.R requires identical() results against scDblFinder for the bundled fragment-file metrics, the production call with prefixed Cell Ranger barcodes, a deterministic multi-chromosome loci fixture, and the corresponding full AMULET metrics, and requires an explicit error for PCR-duplicate expansion.

pixi run --use-environment-activation-cache test-amulet-parity

BPCells-backed ATAC scDblFinder feature aggregation

Reference algorithm. With aggregateFeatures = TRUE, scDblFinder::scDblFinder() performs its own TF-IDF-based feature clustering and sums peaks into feature groups before artificial-doublet classification.

Reason for adaptation. Materializing and transforming every peak within each GEM well can exhaust worker memory before classification. The pipeline instead derives feature groups once from the aggregation’s LSI loadings, sums the disk-backed peak matrix with BPCells, and passes the compact matrix to scDblFinder with aggregateFeatures = FALSE.

Deliberate deviations and consequences. This is not a reimplementation of the classifier, which runs unchanged. Global LSI-derived groups replace scDblFinder’s per-GEM-well feature groups, so the aggregated matrix, and therefore scores and calls, can differ; exact parity is not expected. The GEX path calls scDblFinder::scDblFinder() directly through the same memory-bounding per-GEM-well wrapper.

Implementation. get_feature_groups_from_LSI() and aggregate_BPCells_rows_by_group() in R/processing_GEX_helpers.R, called from extra_targets/ATAC_targets.R.

Validation scope. The reference-parity fixtures do not establish equivalence for this path. The method is described in Filtering and doublet removal.

Native weighted nearest neighbors

Reference algorithm. Seurat::FindMultiModalNeighbors() constructs cell-specific modality weights from within- and cross-modality neighborhood prediction, collects candidate neighbors across modalities, and selects a weighted multimodal neighbor set.

Reason for reimplementation. The pipeline already holds aligned RNA and ATAC embeddings and needs reusable neighbor indices, distances, and modality weights without creating a Seurat object. Native graph state also feeds UMAP, Leiden clustering, SCAVENGE, and the Seurat/Signac export.

Deliberate deviations and consequences. HNSW search with RcppHNSW replaces Seurat’s Annoy search, so candidate sets need not be identical. The helper does not expose Seurat’s optional smoothing or cross-constant list, and BPCells builds the downstream SNN graph rather than Seurat Neighbor and Graph objects. Modality weights are estimated within random subsets of at most 10,000 cells, whose assignment the seed argument fixes, because Seurat’s weights become nearly binary in large datasets; the validation fixtures are smaller, so this deviation is not tested. The seed does not control neighbor-search randomness. These choices can change weights, selected neighbors, SNN edges, clusters, and UMAP coordinates, so correlation and overlap, not exact equality, are the validation contract.

Implementation. weighted_nearest_neighbors_BPCells() in R/processing_multimodal_helpers.R with the small-SNN bandwidth kernel in src/wnn_snn_bandwidth.cpp; extra_targets/WNN_targets.R aligns the embeddings and wires the result into clustering, UMAP and metadata targets.

Validation. tests/testthat/test-wnn-parity.R compares deterministic RNA/ATAC fixtures with Seurat at production search settings: a 400-cell fixture with candidate range 200 and k of 20 and 50, and a 160-cell stress fixture with k = 15 and candidate range 50. Each case must exceed its thresholds for modality-weight Spearman correlation and mean and lower-quartile neighbor-set overlap, and one- and two-thread results must be identical. This does not assert equality of selected neighbors, SNN weights, clustering, or UMAP.

pixi run --use-environment-activation-cache test-algorithm-validation

Sparse SCAVENGE propagation

Reference algorithm. SCAVENGE at commit 8ee8b173d965 selects high chromVAR z-score seed cells, constructs a binary mutual-nearest-neighbor graph, performs a column-normalized random walk with restart, caps and rescales the propagation score into a trait relevance score (TRS), and uses degree-matched seed permutations to identify significant cells.

Reason for reimplementation. The reference package’s dependency stack predates the pipeline’s R/Bioconductor environment. multiomeR needs sparse propagation over its native SNN graphs.

Deliberate deviations and consequences. The reference builds a mutual-kNN graph, whereas multiomeR uses the binary support of its BPCells-derived SNN graph; edge weights are discarded, but topology can differ. Seed and scale-factor helpers guarantee at least one selected cell for small inputs, and the random walk has a maximum-iteration guard. The degree-matched seed permutations and significant-cell calls are omitted: multiomeR reports trait relevance scores and their group summaries without P-values.

Implementation. R/SCAVENGE_helpers.R; extra_targets/module_genetic_enrichment/SCAVENGE_graph_targets.R builds the graph and cell-level trait relevance scores, and SCAVENGE_group_targets.R summarizes and plots them by cluster.

Validation. tests/testthat/test-scavenge-parity.R uses a deterministic fixture with heterogeneous-degree graph blocks and nonuniform edge weights, so it also tests conversion to binary adjacency. The iterative random walk must match the closed-form solution

[ s = r(I - (1-r)P)^{-1}p_0 ]

within 1e-10. Compact local reference functions reproduce the relevant SCAVENGE code at the pinned commit without installing its dependency stack; seed selection must be identical, and propagation and transformed scores must agree within 1e-12. This does not establish parity of graph construction, chromVAR inputs, or biological interpretation.

pixi run --use-environment-activation-cache test-algorithm-validation

Stratified LD score regression

Reference algorithm. ldsc 1.0.1 cleans summary statistics with munge_sumstats.py, computes annotation LD scores from PLINK reference genotypes with ldsc.py --l2, and fits stratified LD score regression with a block jackknife with ldsc.py --h2 (Bulik-Sullivan et al. 2015; Finucane et al. 2015).

Reason for reimplementation. ldsc needs Python 2.7, which the workflow’s environment does not provide, and computes LD scores for every reference variant. The R implementation reads the PLINK genotypes directly, computes LD scores only at the HapMap3 regression variants, and fits all cell types of a GWAS from shared block cross-products.

Deliberate deviations and consequences. Regression weights use the total LD score of the annotations shared by all models, the baseline-LD and all-peaks annotations, rather than that of each cell type’s model. All cell types of a GWAS therefore share one set of block cross-products. In the development comparison below, this changed coefficient z-scores by at most 0.003.

Implementation. munge_GWAS_sumstats(), get_SLDSC_LD_score_record(), fit_SLDSC_candidate_regressions() and jackknife_LDSC_blocks() in R/SLDSC_helpers.R; the targets are in extra_targets/module_genetic_enrichment/SLDSC_targets.R, and the method is described in Heritability enrichment with S-LDSC.

Validation scope. No fixture in the test suite runs ldsc, because the environment does not provide it. During development, the implementation was compared once with ldsc 1.0.1 on the peak annotations and GWAS of mixed_human_31x. Chromosome 22 LD scores agreed within 0.0005, the rounding of ldsc’s output, and munging kept the same variants. The coefficient z-score, total heritability and intercept agreed to six digits. On the annotations and eight UK Biobank GWAS of Varshney et al., fitted with their 37-annotation baseline, the implementation reproduced their published z-scores with a correlation of 0.99985 (Varshney et al. 2024). These comparisons are not rerun by the suite.

Bulik-Sullivan, Brendan K., Po-Ru Loh, Hilary K. Finucane, et al. 2015. “LD Score Regression Distinguishes Confounding from Polygenicity in Genome-Wide Association Studies.” Nature Genetics 47: 291–95. https://doi.org/10.1038/ng.3211.
Finucane, Hilary K., Brendan Bulik-Sullivan, Alexander Gusev, et al. 2015. “Partitioning Heritability by Functional Annotation Using Genome-Wide Association Summary Statistics.” Nature Genetics 47: 1228–35. https://doi.org/10.1038/ng.3404.
Varshney, Arushi, Nandini Manickam, Peter Orchard, et al. 2024. “Population-Scale Skeletal Muscle Single-Nucleus Multi-Omic Profiling Reveals Extensive Context Specific Genetic Regulation.” bioRxiv, ahead of print. https://doi.org/10.1101/2023.12.15.571696.