4  Algorithmic implementations, deviations and validation

multiomeR reimplements a small number of reference algorithms so they can operate on the workflow’s native matrices and graph state. This page is the maintained record of what those implementations preserve, where they deliberately differ, and what the executable validation establishes.

The evidence labels used below are intentionally narrow:

Passing these fixtures does not validate every dataset, parameter regime, approximate-neighbor realization, biological interpretation, or downstream target. The test suite distinguishes fast unit tests for isolated data contracts from slower parity and integration tests that load the project runtime or compare external reference implementations. The CI workflow runs the complete suite when tests, relevant helpers, or the Pixi environment change.

Run the complete suite with pixi run test. The narrower pixi run test-algorithm-validation task runs only the slow UCell, AMULET, WNN, and SCAVENGE parity and acceptance tests.

Implementation Evidence status Maintained reference Current fixed-fixture result
BPCells-native UCell Reference-parity tested UCell 2.14.0 Exact values, dimensions, and dimnames
BPCells-native AMULET Reference-parity tested scDblFinder 1.24.0 Exact metrics and multi-chromosome loci, including order
Native WNN Reference-similarity tested Seurat 5.5.0 Small-SNN pilot: weight Spearman 0.989; mean neighbor overlap 0.996
Sparse SCAVENGE propagation Algorithmically derived and reference-parity tested SCAVENGE 1.0.2 at 8ee8b173d965 Closed-form delta 8.61e-13; pinned-reference propagation delta 1.11e-16; exact streamed exceedance counts and significant-cell calls

4.1 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. The maintained comparison also covers UCell::AddModuleScore_UCell().

Reason for reimplementation. The workflow keeps gene-by-cell counts in BPCells-backed matrices. Materializing the complete matrix in memory 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.

Behavior preserved. The implementation preserves UCell signature syntax (+ and - suffixes), descending per-cell ranks, configurable tie handling, maxRank capping, impute/skip behavior for missing genes, negative-signature weighting, lower-bound clipping, signature names, and cell order.

Deliberate deviations and consequences. Matrix materialization is limited to one cell chunk at a time, and optional fork workers operate across chunks. This changes memory and execution behavior but not the tested score values. The helper returns a data frame instead of mutating a Seurat object. The target-level marker validator rejects configured genes missing from the Cell Ranger reference before normal pipeline scoring, whereas the lower-level helper still exposes UCell’s impute/skip modes for explicit use.

Implementation and wiring. The scorer is calculate_BPCells_UCell_scores_from_matrix() in R/processing_GEX_helpers.R. extra_targets/general_aggregation_targets.R validates UCell_GEX_marker_genes_list, and extra_targets/GEX_graph_and_cluster_targets.R computes cell-level scores before GEX cell-type assignment. ATAC and WNN cell-type targets reuse those metadata scores.

Validation. tests/testthat/test-scoring-parity.R creates a deterministic 500-gene by 37-cell matrix, writes the project input as BPCells, and compares signed signatures with imputed and skipped missing genes. It requires identical() values, dimensions, and dimnames against UCell 2.14.0. It also retains exact Seurat AddModuleScore and cell-cycle helper checks plus a metadata-join contract.

The UCell 2.14.0 reference call with imputed missing genes can emit non-fatal R stack-imbalance warnings under the repository’s R 4.5 environment. A narrow diagnostic reproduced them in the UCell reference call but not in the repository scorer. The validation therefore runs the reference call in a disposable callr process and compares its returned matrix in the clean parent session; this isolates the package warning without weakening the equality assertion or changing the runtime.

Status and rerun. Reference-parity tested against UCell 2.14.0 with exact equality; passing in the current locked environment.

pixi run test-scoring-parity

4.2 BPCells-native AMULET

Reference algorithm. scDblFinder::amulet() detects likely scATAC-seq doublets from the number of genomic loci covered by more than two fragments. Its underlying getFragmentOverlaps() implementation filters fragment sizes and excluded regions, calculates per-barcode fragment and overlap counts, removes 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 original fragment TSV to scDblFinder materializes chromosome-scale GRanges objects and previously requested six cores and 60 GB. The local implementation streams the existing BPCells fragment target and retains only one chromosome’s selected fragments while calculating coverage runs.

Behavior preserved. calculate_amulet_metrics_BPCells() preserves barcode selection, minimum-fragment thresholds, maximum fragment size, excluded regions, nFrags, uniqFrags, nAbove2, total.nAbove2, p-values, q-values, and high-overlap-site removal. The lower-level loci return also preserves scDblFinder’s cell-major, chromosome, and coordinate ordering. Cell Ranger’s inclusive end-insertion convention is shifted back by one base before calculation so the BPCells representation matches scDblFinder’s BED import.

Deliberate deviations and consequences. Only unique-fragment operation is supported. BPCells fragment objects do not retain Cell Ranger’s PCR-duplicate count column, so requesting non-unique expansion fails explicitly instead of silently changing nFrags. BPCells does not export its fragment iterator header; the native helper therefore mirrors that private C++ interface, verifies the exact project-pinned BPCells commit 28759cdd5125 before use, and compiles a small shared library in each worker’s temporary directory. A BPCells upgrade must revalidate this interface and the exact parity fixture before updating the pin. The target is single-threaded and requests the standard 16-GB worker tier. A native-only probe on the stored healthy_PBMC_human fragments processed 2,711 selected cells in 15.4 seconds with 0.64 GB peak RSS, including R startup and native compilation; this supports the reduced allocation for the public fixture but is not a memory guarantee for larger datasets.

Implementation and wiring. Native iteration and coverage-run calculation are implemented in src/amulet_bpcells.cpp. The R wrapper, ABI check, high-overlap filtering, and AMULET statistics are in R/amulet_BPCells_helpers.R. amulet_metrics_tibble in extra_targets/per_GEM_well_targets.R consumes the existing prefixed BPCells fragments, restores unprefixed barcode keys, and preserves the existing downstream metrics shape.

Validation. tests/testthat/test-amulet-parity.R compares against scDblFinder 1.24.0 in disposable callr processes. It requires identical() results for the bundled fragment-file metrics, prefixed pipeline barcodes, a deterministic 12,000-fragment multi-chromosome loci fixture, and the corresponding full AMULET metrics. It also requires an explicit error for unsupported PCR-duplicate expansion. The disposable reference processes isolate stack-imbalance warnings emitted by the current scDblFinder reference under R 4.5 without weakening the returned-object comparison.

Status and rerun. Reference-parity tested against scDblFinder 1.24.0 with exact equality; passing in the current locked environment.

pixi run test-amulet-parity

4.3 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 has aligned RNA PCA/Harmony and ATAC LSI/Harmony matrices 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 optional Seurat/Signac export.

Behavior preserved. weighted_nearest_neighbors_BPCells() preserves row-wise L2 normalization, per-modality nearest neighbors, nearest-distance correction, Seurat’s small-SNN far-neighbour kernel bandwidth, within/cross prediction kernels, capped modality affinity ratios, normalized cell-specific modality weights, candidate-set union, weighted neighbor ranking, and Seurat’s transformation from weighted affinity to neighbor distance.

Deliberate deviations and consequences. BPCells HNSW replaces Seurat’s Annoy search, so approximate candidate sets need not be identical. The helper does not expose Seurat’s optional smoothing or cross-constant list, and BPCells builds downstream SNN state rather than storing Seurat Neighbor and Graph objects. Its current seed argument is not consulted by the HNSW calls, so it must not be interpreted as controlling neighbor-search randomness. These choices can change weights, selected neighbors, SNN edges, clusters, and UMAP coordinates; correlation and overlap are therefore the validation contract, not exact equality.

Implementation and wiring. The implementation and graph consumers are in R/processing_multimodal_helpers.R; the tracked project-owned small-SNN kernel is in src/wnn_snn_bandwidth.cpp. extra_targets/WNN_targets.R aligns modality embeddings, creates WNN_results_raw, filters small clusters, optionally recomputes WNN_results, and wires that state into WNN UMAP, clustering, metadata, and cell-type targets.

Validation. tests/testthat/test-wnn-parity.R compares two deterministic RNA/ATAC fixtures with Seurat 5.5.0. The headline production-like fixture uses 400 cells, 12 dimensions, k = 30, and candidate range 200, matching the pipeline’s configured neighbor count and native WNN candidate range. The pre-migration small-SNN pilot gave modality-weight Spearman 0.989 and mean neighbour-set overlap 0.996 on this fixture. On the 26,667-cell production object, it gave RNA/ATAC weight Spearman 0.988/0.988 and mean neighbour overlap 0.979. These figures supported the migration, but the maintained validation was deliberately not rerun as part of the production commit. The existing thresholds remain the acceptance contract. The current locked BPCells 0.3.1 build is pinned at 28759cdd5125.

Status and rerun. The small-SNN default is based on the reference-similarity pilot above; post-migration validation remains pending. This does not assert exact equality of selected neighbours, SNN weights, clustering, or UMAP.

pixi run test-algorithm-validation

4.4 Sparse SCAVENGE propagation and significance

Reference algorithm. SCAVENGE 1.0.2 at commit 8ee8b173d965 selects high chromVAR Z-score seed cells, constructs a binary mutual-nearest-neighbor adjacency 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 last commit and dependency stack predate the pipeline’s current R/Bioconductor environment. multiomeR needs sparse propagation over native RNA PCA, ATAC LSI, and multimodal WNN SNN matrices and must avoid materializing a cell-by-permutation score matrix for large cell sets.

Behavior preserved. The implementation converts nonzero graph support to binary adjacency before analysis and preserves the one-sided Z-score seed threshold and top-percent cap, column-normalized transition matrix, equal seed restart mass, iterative random walk, 0.95 propagation-score cap, min-max scaling, Z-score scale factor, sequential base-R degree-matched seed sampling, and strict per-cell comparison with permuted propagation scores. The sampled seed-index lists are retained, but the cell-by-permutation score matrix is not: a native worker streams random walks and accumulates only per-cell exceedance counts and the cluster medians needed downstream. Random walks, rather than random-number generation, are parallelized, so the sampled null is invariant to the requested core count.

Deliberate deviations and consequences. The reference workflow constructs a mutual-kNN graph, whereas multiomeR uses the binary support of its BPCells-derived PCA, LSI, or WNN SNN graph; edge weights are discarded, but graph topology can still differ. Seed and scale-factor helpers guarantee at least one selected cell for small inputs. The degree sampler also handles a one-cell candidate stratum explicitly, avoiding base R’s special interpretation of sample(x, 1) when x is one positive integer. The random walk validates graph inputs and has a maximum-iteration guard. Cell-level empirical P-values and significance calls follow the reference exceedance fraction and threshold. Cluster-level permutation medians, add-one P-values, and Benjamini–Hochberg adjustment within each grouping column are pipeline extensions.

Implementation and wiring. Seed selection, sparse random walk, streaming degree-matched permutations, TRS construction, and cluster-level null statistics are in R/SCAVENGE_helpers.R. module_genetic_enrichment/SCAVENGE_graph_targets.R constructs each graph and maps chromVAR Z-score records into SCAVENGE_result_records, from which cell-level TRS and cluster-level summaries are extracted; SCAVENGE_group_targets.R combines summaries and plots.

Validation. tests/testthat/test-scavenge-parity.R uses a deterministic 60-cell fixture with repeated heterogeneous-degree graph blocks, nonuniform input edge weights, and three enriched seeds. The production helper receives the weighted graph, so the fixture also tests conversion to binary adjacency. First, the iterative sparse random walk is compared with the closed-form solution

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

with a maximum absolute tolerance of 1e-10; the current delta is 8.61e-13. Second, compact local reference functions reproduce the relevant SCAVENGE 1.0.2 code at the pinned commit without installing its historical dependency stack. The random-walk delta against that reference is 1.11e-16, the transformed-score delta is 3.33e-16, and all 199 fixed-RNG degree-matched seed samples, streamed per-cell exceedance counts, empirical P-values, and significant-cell calls are identical. One- and two-core native results are also identical.

Status and rerun. Random-walk propagation is algorithmically derived against the closed form. Seed selection, binary propagation, transformed scores, sequential permutation sampling, streamed exceedance counts, empirical P-values, and significant-cell calls are reference-parity tested against the pinned source calculation; cluster summaries are documented pipeline extensions. This does not establish parity of mutual-kNN versus pipeline graph construction, chromVAR inputs, or biological interpretation.

pixi run test-algorithm-validation

Run the complete maintained test suite, including the fast helper contracts, with:

pixi run test