Primary-module methods

This page describes the analyses of the primary module in the order of the checkpoints in Run your own analysis. The same descriptions form the Supplementary Methods of the multiomeR manuscript. Configurable settings and their defaults are listed in the parameter browser, and the reimplemented algorithms are compared with their references in Algorithm validation.

Each section ends with its source files and a simplified view of its target graph. The views keep real target names but omit lower-level targets; Change the pipeline explains how they are made.

Preprocessing

Checkpoint: 1: Pre-aggregation QC

Each GEM well supplies a Cell Ranger ARC count directory. Gene-expression counts are imported into BPCells from the filtered feature-barcode matrix, or from a CellBender output when configured (Parks and Greenleaf 2025; Fleming et al. 2023). ATAC fragments are imported from the Cell Ranger fragment file. Both carry a GEM-well-prefixed barcode, so nuclei remain distinct across wells. Per-well quality-control metrics are calculated for the barcodes that Cell Ranger called as cells. BPCells calculates ATAC metrics such as TSS enrichment and nucleosome signal. The blacklist fraction is the share of a nucleus’s fragments with at least one cut site in a blacklisted region of the genome (Amemiya et al. 2019; Ogata et al. 2023). Gene-expression metrics, such as the percentage of counts from mitochondrial genes (names beginning with MT-, ignoring case), are calculated from the imported matrix. Optional genotype demultiplexing runs cellsnp-lite on the ATAC alignments, by default with a minimum minor-allele frequency of 0.1 and 20 reads per SNP (Huang and Huang 2021). Vireo then demultiplexes the barcodes with the supplied donor genotypes (Huang et al. 2019). The most likely singlet donor replaces the GEM-well donor identifier of every barcode. Vireo’s singlet, doublet and unassigned calls and probabilities are recorded as metadata for the exclusions. A BPCells-based implementation of AMULET reads the stored fragments of up to 1,000 bp one by one, without loading them all into memory (Thibodeau et al. 2021). This avoids the chromosome-scale GRanges objects that scDblFinder builds from the text fragment file (Germain et al. 2022). It skips fragments on chrM, chrX and chrY and, for GRCh38, fragments in the AMULET authors’ list of repeats, segmental duplications and blacklisted regions, where multi-mapping reads cause spurious overlaps. It reproduces the overlap metrics and q-values of scDblFinder::amulet() but applies no automatic filter. Because AMULET performs well only on deeply sequenced nuclei, filtering on its q-values is left to the user, and its authors recommend q < 0.01 (Thibodeau et al. 2021).

Configured GEM-well exclusions are filter expressions applied to the called cells of each well. Demultiplexing and AMULET results can therefore be used as exclusion criteria before aggregation-level analysis. Called cells missing from the imported gene-expression matrix are always excluded. The remaining nuclei of the selected wells are combined without further filtering.

GEM-well settings are described in GEM well table, and the AMULET implementation is compared with scDblFinder in Algorithm validation.

Source and target graph

Source: extra_targets/per_GEM_well_targets.R, extra_targets/general_aggregation_targets.R, R/parallel_GEM_well_preprocessing_helpers.R, R/amulet_BPCells_helpers.R, R/processing_GEX_helpers.R, R/processing_ATAC_helpers.R.

This view covers per-GEM-well processing and QC, including optional ambient RNA correction, donor demultiplexing, doublet detection, barcode filtering, and handoffs into aggregation-level GEX and ATAC objects.

%%{init: {
  'theme': 'base',
  'themeVariables': {
    'primaryTextColor': '#1f2937',
    'lineColor': '#4b5563',
    'clusterBkg': '#f8fafc',
    'clusterBorder': '#cbd5e1',
    'edgeLabelBackground': '#ffffff00',
    'fontFamily': 'Helvetica, Arial, sans-serif'
  },
  'flowchart': {
    'defaultRenderer': 'elk'
  },
  'themeCSS': '.cluster rect { fill: #f8fafc !important; stroke: #cbd5e1 !important; } .cluster-label text, .cluster-label span { font-weight: 700 !important; font-size: 16px !important; fill: #475569 !important; color: #475569 !important; }'
}}%%
%%{init: {'flowchart': {'defaultRenderer': 'elk'}}}%%
flowchart TB
  cellranger_barcodes_tsv.controller_PBMC_human(["cellranger_barcodes_tsv.&lt;GEM_well_ID&gt;"]):::none --> amulet_metrics_tibble.controller_PBMC_human(["amulet_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none
  cellranger_barcodes_tsv.controller_PBMC_human(["cellranger_barcodes_tsv.&lt;GEM_well_ID&gt;"]):::none --> ATAC_blacklist_metrics_tibble.controller_PBMC_human(["ATAC_blacklist_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none
  cellranger_barcodes_tsv.controller_PBMC_human(["cellranger_barcodes_tsv.&lt;GEM_well_ID&gt;"]):::none --> cellsnp_dir.controller_PBMC_human(["cellsnp_dir.&lt;GEM_well_ID&gt;"]):::none
  full_metadata_tibble.controller_PBMC_human(["full_metadata_tibble.&lt;GEM_well_ID&gt;"]):::none --> excluded_cellranger_only_barcodes_by_type_list.controller_PBMC_human(["excluded_cellranger_only_barcodes_by_type_list.&lt;GEM_well_ID&gt;"]):::none
  ATAC_blacklist_metrics_tibble.controller_PBMC_human(["ATAC_blacklist_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none --> full_metadata_tibble.controller_PBMC_human(["full_metadata_tibble.&lt;GEM_well_ID&gt;"]):::none
  ATAC_qc_metrics_tibble.controller_PBMC_human(["ATAC_qc_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none --> full_metadata_tibble.controller_PBMC_human(["full_metadata_tibble.&lt;GEM_well_ID&gt;"]):::none
  vireo_donor_ids_tibble.controller_PBMC_human(["vireo_donor_ids_tibble.&lt;GEM_well_ID&gt;"]):::none --> full_metadata_tibble.controller_PBMC_human(["full_metadata_tibble.&lt;GEM_well_ID&gt;"]):::none
  amulet_metrics_tibble.controller_PBMC_human(["amulet_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none --> full_metadata_tibble.controller_PBMC_human(["full_metadata_tibble.&lt;GEM_well_ID&gt;"]):::none
  cellbender_h5_file.controller_PBMC_human(["cellbender_h5_file.&lt;GEM_well_ID&gt;"]):::none --> GEX_counts_BPcells_matrix_dir.controller_PBMC_human(["GEX_counts_BPcells_matrix_dir.&lt;GEM_well_ID&gt;"]):::none
  cellsnp_dir.controller_PBMC_human(["cellsnp_dir.&lt;GEM_well_ID&gt;"]):::none --> vireo_donor_ids_tibble.controller_PBMC_human(["vireo_donor_ids_tibble.&lt;GEM_well_ID&gt;"]):::none
  cellranger_barcodes_tsv.controller_PBMC_human(["cellranger_barcodes_tsv.&lt;GEM_well_ID&gt;"]):::none --> vireo_donor_ids_tibble.controller_PBMC_human(["vireo_donor_ids_tibble.&lt;GEM_well_ID&gt;"]):::none
  fragments_w_prefix_bpcells_dir.controller_PBMC_human(["fragments_w_prefix_bpcells_dir.&lt;GEM_well_ID&gt;"]):::none --> amulet_metrics_tibble.controller_PBMC_human(["amulet_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none
  fragments_w_prefix_bpcells_dir.controller_PBMC_human(["fragments_w_prefix_bpcells_dir.&lt;GEM_well_ID&gt;"]):::none --> ATAC_blacklist_metrics_tibble.controller_PBMC_human(["ATAC_blacklist_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none
  fragments_w_prefix_bpcells_dir.controller_PBMC_human(["fragments_w_prefix_bpcells_dir.&lt;GEM_well_ID&gt;"]):::none --> ATAC_qc_metrics_tibble.controller_PBMC_human(["ATAC_qc_metrics_tibble.&lt;GEM_well_ID&gt;"]):::none
  fragments_w_prefix_bpcells_dir.controller_PBMC_human(["fragments_w_prefix_bpcells_dir.&lt;GEM_well_ID&gt;"]):::none --> combined_BPCells_fragment_obj.ATAC.ENCODE_heart_LV_6x(["combined_BPCells_fragment_obj.ATAC.&lt;aggregation_name&gt;"]):::none
  GEX_counts_BPcells_matrix_dir.controller_PBMC_human(["GEX_counts_BPcells_matrix_dir.&lt;GEM_well_ID&gt;"]):::none --> full_metadata_tibble.controller_PBMC_human(["full_metadata_tibble.&lt;GEM_well_ID&gt;"]):::none
  GEX_counts_BPcells_matrix_dir.controller_PBMC_human(["GEX_counts_BPcells_matrix_dir.&lt;GEM_well_ID&gt;"]):::none --> aggregated_GEX_BPCells_matrix_dir.GEX.ENCODE_heart_LV_6x(["aggregated_GEX_BPCells_matrix_dir.GEX.&lt;aggregation_name&gt;"]):::none
  donor_id_metadata_tibble.ENCODE_heart_LV_6x(["donor_id_metadata_tibble.&lt;aggregation_name&gt;"]):::none

GEX

Checkpoints: 2: GEX dimension reduction, 3: GEX clusters and cell types

Normalization and PCA

For each aggregation, the quality-controlled gene-expression matrices are combined with their GEM-well-prefixed barcodes. Genes with 50 or fewer total counts are excluded before dimensional reduction. Normalization and dimensional reduction use either the disk-backed BPCells Pearson-residual workflow or Seurat SCTransform (Parks and Greenleaf 2025; Hafemeister and Satija 2019). With the BPCells backend, analytic Pearson residuals are computed with one overdispersion parameter (θ = 100) for all genes and clipped to ±√N for N nuclei (Lause et al. 2021). Per-gene estimates of θ would absorb differences between cell types and flatten marker-gene residuals, and the √N clip keeps markers of rare populations (Lause et al. 2021). Configured cell-level covariates can be regressed out. The 3,000 genes with the highest residual variance are kept, and a truncated singular-value decomposition gives at least 50 principal components, of which the configured ones are used downstream. It reads these genes’ counts from memory in a fixed number of blocks, so its result does not depend on the number of threads. The Seurat backend applies SCTransform with its v2 regularization to the same cells and genes (Choudhary and Satija 2022). It keeps the same number of variable features and computes the principal components from the residual matrix. A cell-cycle difference score can be regressed out but is not by default. Its S and G2M scores are computed from log-normalized counts with the 2019 Seurat cell-cycle gene sets (Tirosh et al. 2016).

Batch correction and clustering

Harmony correction can be applied separately to the selected GEX principal components and ATAC LSI dimensions, with at most 25 iterations and lambda = 1 (Korsunsky et al. 2019; Patikas et al. 2026). It is off by default, since correcting a variable can also remove biology confounded with it (Luecken et al. 2022). Harmony corrects each configured metadata column as a separate batch variable, and the ATAC correction can add its own columns. Nuclei with missing covariate values are removed with a warning.

For each modality, the nearest neighbors of every nucleus are found once on the selected dimensions, by hierarchical navigable small world (HNSW) search with cosine distance (Malkov and Yashunin 2020; Melville 2024). The search index is built on one thread, so a rebuild finds the same neighbors, while the search itself runs in parallel. Clustering and UMAP share these neighbors. They are converted to a shared-nearest-neighbor (SNN) graph with Jaccard weights, removing edges with weights below 1/15. The graph is clustered with the modularity objective of the Leiden algorithm in ten iterations (Parks and Greenleaf 2025; Traag et al. 2019). Clusters are renumbered by size, and clusters of fewer than 100 nuclei by default are removed from the GEX, ATAC and WNN clusterings. UMAP embeddings are computed with uwot from the same neighbors (McInnes et al. 2018; Melville 2026). GEX and ATAC layouts start from the first two principal components of the clustering dimensions, and WNN layouts from a spectral embedding of the neighbor graph. uwot’s batch optimization gives the same layout for any number of threads.

Most clusters below the size threshold hold a single nucleus or a few nuclei that share too few neighbors with any larger cluster, so the filter mainly removes such outlier nuclei. Rare populations smaller than the threshold are removed as well.

Cell-type annotation

Cell-type annotation uses GEX marker signatures supplied by the user. Signatures list positive markers and, optionally, genes expected to be absent. Genes missing from the count matrix stop the run with an error before scoring. Scores are computed with a BPCells-based implementation of UCell scoring (Andreatta and Carmona 2021; Parks and Greenleaf 2025). It processes the raw GEX counts in blocks of limited size and reproduces UCell’s per-cell descending ranks, rank truncation at 1,500 genes, positive and negative signatures and lower-bound clipping. Working in blocks avoids loading the whole matrix or building a Seurat object only for scoring. The pipeline does not use the UCell package itself. Validation instead requires agreement with the UCell reference to within 10−12. Signed signatures are clipped at zero per cell before averaging. The same scoring annotates the GEX, ATAC and WNN clusters. WNN labels can use their own marker sets, whose matched controls come from the same GEX reference nuclei.

Raw UCell scores are not comparable between signatures, because some gene sets score highly by chance. Each label’s cluster-level mean score is therefore compared with 999 random control signatures matched on gene abundance and detection. The control reference samples up to 50 nuclei per GEM well from the GEX metadata before doublet filtering. The GEX and ATAC annotations thus share the same controls, and the WNN annotation draws its own from the same reference, whatever nuclei each stage keeps. For each control replicate, all marker genes of the panel are taken in random order. Each is replaced by a gene drawn from its 50 nearest eligible candidates not yet used in that replicate. Candidates exclude all marker genes and undetected genes. A label’s adjusted score in a cluster is its observed mean minus the 95th percentile of its matched controls. The label with the highest adjusted score is the candidate. Its advantage is the smaller of its lead over zero and its lead over the runner-up. The candidate is assigned when its advantage is positive and at least a minimum configured separately for the GEX, ATAC and WNN clusters (0.05 by default). Otherwise, the cluster remains unassigned, with the candidate and reason kept. Requiring a lead over both the background and the runner-up leaves ambiguous clusters unassigned rather than forcing a label. Raising the minimum can therefore only make assigned clusters unassigned, never relabel them. Each unassigned cluster keeps its own label, Unassigned_cluster_<id>, so doublet detection, pseudobulks and per-cell-type tests never pool unrelated unassigned clusters. Because annotation depends on the supplied signatures and their level of detail, labels can represent either cell types or broader groups of cells, such as a lineage or a tissue of origin.

Annotation outputs also report, for each cluster, how many positive markers are detected in at least 10% of its nuclei.

The BPCells-native UCell scorer is compared with UCell in Algorithm validation.

Doublet removal

After GEX clustering, scDblFinder is applied per GEM well, with dbr.sd = 1, to the gene-expression counts of the kept nuclei (Germain et al. 2022). These counts are CellBender-corrected where configured. Doublets form only within a GEM well, and dbr.sd = 1 sets the doublet threshold from misclassification alone rather than from an expected doublet rate. scDblFinder receives cluster labels from the annotation: the cell-type label of an assigned cluster, or otherwise its annotation status and identifier. A configured map can optionally collapse these labels. By default, nuclei called as doublets are removed, as are all nuclei of Leiden clusters whose doublet fraction exceeds 0.5.

Source and target graph

Source: extra_targets/GEX_merge_and_dim_reduc_targets.R, extra_targets/GEX_graph_and_cluster_targets.R, R/processing_GEX_helpers.R, R/processing_ATAC_helpers.R, R/cluster_annotation_helpers.R, src/cluster_UCell_chunk.cpp.

This view runs from the aggregated counts through PCA, Harmony correction, clustering and UMAP to cell-type annotation and GEX doublet calls.

%%{init: {
  'theme': 'base',
  'themeVariables': {
    'primaryTextColor': '#1f2937',
    'lineColor': '#4b5563',
    'clusterBkg': '#f8fafc',
    'clusterBorder': '#cbd5e1',
    'edgeLabelBackground': '#ffffff00',
    'fontFamily': 'Helvetica, Arial, sans-serif'
  },
  'flowchart': {
    'defaultRenderer': 'elk'
  },
  'themeCSS': '.cluster rect { fill: #f8fafc !important; stroke: #cbd5e1 !important; } .cluster-label text, .cluster-label span { font-weight: 700 !important; font-size: 16px !important; fill: #475569 !important; color: #475569 !important; }'
}}%%
%%{init: {'flowchart': {'defaultRenderer': 'elk'}}}%%
flowchart TB
  PCA_BPCells.GEX.ENCODE_heart_LV_6x(["PCA_BPCells.GEX.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none --> knn.GEX.ENCODE_heart_LV_6x(["knn.GEX.&lt;aggregation_name&gt;"]):::none
  PCA_BPCells.GEX.ENCODE_heart_LV_6x(["PCA_BPCells.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_tibble.GEX.ENCODE_heart_LV_6x(["metadata_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  scDblFinder_results_df.GEX.ENCODE_heart_LV_6x(["scDblFinder_results_df.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  GEX_normalization.GEX.ENCODE_heart_LV_6x(["GEX_normalization.GEX.&lt;aggregation_name&gt;"]):::none --> PCA_BPCells.GEX.ENCODE_heart_LV_6x(["PCA_BPCells.GEX.&lt;aggregation_name&gt;"]):::none
  knn.GEX.ENCODE_heart_LV_6x(["knn.GEX.&lt;aggregation_name&gt;"]):::none --> PCA_clusters.GEX.ENCODE_heart_LV_6x(["PCA_clusters.GEX.&lt;aggregation_name&gt;"]):::none
  PCA_BPCells.GEX.ENCODE_heart_LV_6x(["PCA_BPCells.GEX.&lt;aggregation_name&gt;"]):::none --> UMAP_embeddings_tibble.GEX_non_harmony.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.GEX_non_harmony.&lt;aggregation_name&gt;"]):::none
  UMAP_embeddings_tibble.GEX.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> UMAP_embeddings_tibble.GEX_non_harmony.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.GEX_non_harmony.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none --> UMAP_embeddings_tibble.GEX.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  knn.GEX.ENCODE_heart_LV_6x(["knn.GEX.&lt;aggregation_name&gt;"]):::none --> UMAP_embeddings_tibble.GEX.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  aggregated_GEX_BPCells_matrix_dir.GEX.ENCODE_heart_LV_6x(["aggregated_GEX_BPCells_matrix_dir.GEX.&lt;aggregation_name&gt;"]):::none --> GEX_normalization.GEX.ENCODE_heart_LV_6x(["GEX_normalization.GEX.&lt;aggregation_name&gt;"]):::none
  aggregated_GEX_BPCells_matrix_dir.GEX.ENCODE_heart_LV_6x(["aggregated_GEX_BPCells_matrix_dir.GEX.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.GEX.ENCODE_heart_LV_6x(["scDblFinder_results_df.GEX.&lt;aggregation_name&gt;"]):::none
  aggregated_GEX_BPCells_matrix_dir.GEX.ENCODE_heart_LV_6x(["aggregated_GEX_BPCells_matrix_dir.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  PCA_clusters.GEX.ENCODE_heart_LV_6x(["PCA_clusters.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  metadata_tibble.GEX.ENCODE_heart_LV_6x(["metadata_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  PCA_clusters.GEX.ENCODE_heart_LV_6x(["PCA_clusters.GEX.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.GEX.ENCODE_heart_LV_6x(["scDblFinder_results_df.GEX.&lt;aggregation_name&gt;"]):::none
  metadata_tibble.GEX.ENCODE_heart_LV_6x(["metadata_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.GEX.ENCODE_heart_LV_6x(["scDblFinder_results_df.GEX.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.GEX.ENCODE_heart_LV_6x(["scDblFinder_results_df.GEX.&lt;aggregation_name&gt;"]):::none
  donor_id_metadata_tibble.ENCODE_heart_LV_6x(["donor_id_metadata_tibble.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none
  donor_id_metadata_tibble.ENCODE_heart_LV_6x(["donor_id_metadata_tibble.&lt;aggregation_name&gt;"]):::none --> GEX_normalization.GEX.ENCODE_heart_LV_6x(["GEX_normalization.GEX.&lt;aggregation_name&gt;"]):::none
  metadata_tibble.GEX.ENCODE_heart_LV_6x(["metadata_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none

ATAC

Checkpoints: 4: Peak QC, 5: ATAC filtering, 6: ATAC dimension reduction, 7: ATAC clusters and motifs

Peaks and LSI

Fragments on the standard chromosomes are used for peak calling. Peak-calling groups are defined by a configurable metadata column of the nuclei kept after GEX processing, by default their GEX clusters. GEX clusters serve because ATAC clustering needs the peak matrix. Calling peaks per group keeps those of rare populations, which a pooled call would dilute (Stuart et al. 2021; Granja et al. 2021). Groups larger than 50,000 nuclei are randomly downsampled for peak discovery, which limits the runtime and memory of peak calling. Each group’s discovery nuclei are split into two random halves. Both draws are seeded by the group’s name, so a group’s halves do not change when other groups are added or removed. Peaks are called in each half with either MACS3 or the BPCells tile-based caller (Zhang et al. 2008; Parks and Greenleaf 2025). MACS3 runs on each half’s fragment left ends, written in one pass per GEM well, with a −75-bp shift, a 150-bp extension, a 150-bp maximum gap and summit calling. Only one of each fragment’s two insertions thus enters the pileup, as when MACS3 reads a fragment file. The gap is fixed because MACS3’s default, the mean length of the first ten input lines, would be 1 bp for these single-base cut sites. The tile-based caller runs on the stored fragments, using 500-bp tiles and an FDR of 0.01. Consensus peaks follow the ArchR strategy of iterative fixed-width peak construction (Granja et al. 2021). Summits are extended by 250 bp to 500-bp peaks, and peaks overlapping any blacklisted region of the reference genome are removed (Amemiya et al. 2019; Ogata et al. 2023). Overlapping peaks are collapsed by summit significance, first within each half and then across the two halves. Only peaks called in both halves are kept. Like ArchR’s pseudo-replicate rule (Granja et al. 2021), this removes weak peaks called in only one half. Each group keeps its strongest peaks, up to 500 per discovery nucleus and 150,000 in all, as ArchR does by default (Granja et al. 2021). The group-level sets are then combined into one non-overlapping consensus set. Here, overlapping peaks are ranked by the quantile of their summit significance within their group. The caps and within-group quantiles keep large or deeply sequenced groups from dominating the consensus (Corces et al. 2018). For all nuclei kept after GEX processing, BPCells counts each fragment once in every peak that holds at least one of its two Tn5 cut sites. This paired-insertion counting, also Signac’s default, counts a fragment with both cut sites in a peak once rather than twice, and ignores fragments that only span a peak (Miao and Kim 2024; Stuart et al. 2021). Signals that describe positions, such as peak-calling pileups, TSS enrichment and coverage tracks, still use cut sites.

ATAC dimensional reduction starts from the consensus peak matrix of the nuclei that pass the aggregation-level ATAC exclusions. Signac’s TF-IDF method 1, with a scale factor of 10,000, normalizes this matrix (Stuart et al. 2021). A BPCells singular-value decomposition of the result gives the latent semantic indexing (LSI) embeddings (Parks and Greenleaf 2025). It reads the matrix in a fixed number of blocks, so its result does not depend on the number of threads. When memory allows, the transformed matrix or its counts are held in memory during the decomposition. At least 50 components are computed, of which the configured ones are used downstream. The pipeline does not test them for correlation with sequencing depth. Instead, the default dimension settings leave out the first component, which often reflects depth (Stuart et al. 2021).

Filtering and doublet removal

ATAC clusters are built from the LSI dimensions as described in Batch correction and clustering.

ATAC processing starts from the nuclei kept after GEX processing. The peak-calling groups and the consensus peak matrix are built on this cell set. Per-nucleus quality-control metrics are checked against configured aggregation-level exclusions. The fragments in peaks and their fraction of all fragments are counted from the stored fragments, each fragment once, even when its two cut sites fall in different peaks. After ATAC clustering, a separate ATAC-based scDblFinder analysis follows per GEM well, with dbr.sd = 1 (Germain et al. 2022). Instead of scDblFinder’s internal per-well feature aggregation, peaks are grouped once into 50 feature groups. The groups come from k-means clustering of the peaks’ aggregation-wide LSI coordinates on the configured dimensions within 2–20 (loadings scaled by the singular values, as in scDblFinder). The peaks of each group are summed with BPCells before scDblFinder’s unchanged classifier runs. This change limits memory use, but its scores are not expected to equal those from the internal aggregation. That aggregation transforms every peak within each well, which can exhaust memory before classification. ATAC doublets are removed at the nucleus and cluster level as for GEX.

Algorithm validation describes how this feature aggregation differs from scDblFinder’s own.

Motif accessibility

Transcription-factor motifs are represented by the 233 sequence-similarity families of the JASPAR 2026 CORE vertebrate collection (Ovek Baydar et al. 2026). Families are used because related transcription factors share near-identical motifs, whose matches and accessibility cannot be told apart. The pipeline includes a copy of the familial root motifs and the family membership table. Each root motif is trimmed of flank columns whose total count is less than half of its largest column total. These columns come from a minority of the family’s aligned member motifs, so their counts are sparse and can look specific by chance. Scanning with them would restrict matches to that minority’s sub-family. The trimmed root motifs are matched to the consensus peaks with motifmatchr at its default threshold of P < 5 × 10−5 (Schep 2025). Configured transcription factors are matched to families by name, or by name and motif identifier when a symbol belongs to more than one family. betterChromVAR applies the chromVAR workflow with GC-bias correction to compute analytic motif-family deviations and z-scores for individual nuclei (Schep et al. 2017; Germain 2026). All families are tested per ATAC annotation class with one-versus-rest Wilcoxon tests, with Benjamini–Hochberg correction (Benjamini and Hochberg 1995) across all family–class tests. They are also summarized as cell-weighted means per ATAC cluster, GEX cluster and GEX cell type. The configured transcription factors select only the families shown on embeddings. These values measure accessibility associated with a motif family and are not direct measurements of transcription-factor activity.

Source and target graph

Source: extra_targets/ATAC_targets.R, R/processing_ATAC_helpers.R, R/processing_GEX_helpers.R, R/celltype_labeling_helpers.R, resources/.

This view runs from the combined fragments and the GEX-annotated nuclei through consensus peaks and LSI to ATAC clusters, doublet calls, motif-family accessibility and coverage tracks.

%%{init: {
  'theme': 'base',
  'themeVariables': {
    'primaryTextColor': '#1f2937',
    'lineColor': '#4b5563',
    'clusterBkg': '#f8fafc',
    'clusterBorder': '#cbd5e1',
    'edgeLabelBackground': '#ffffff00',
    'fontFamily': 'Helvetica, Arial, sans-serif'
  },
  'flowchart': {
    'defaultRenderer': 'elk'
  },
  'themeCSS': '.cluster rect { fill: #f8fafc !important; stroke: #cbd5e1 !important; } .cluster-label text, .cluster-label span { font-weight: 700 !important; font-size: 16px !important; fill: #475569 !important; color: #475569 !important; }'
}}%%
%%{init: {'flowchart': {'defaultRenderer': 'elk'}}}%%
flowchart TB
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> ATAC_coverage_marker_genes.ENCODE_heart_LV_6x(["ATAC_coverage_marker_genes.&lt;aggregation_name&gt;"]):::none
  JASPAR_motif_family_members_tibble(["JASPAR_motif_family_members_tibble"]):::none --> ATAC_marker_motif_families_list.ENCODE_heart_LV_6x(["ATAC_marker_motif_families_list.&lt;aggregation_name&gt;"]):::none
  metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none --> chromVAR_obj.ATAC.ENCODE_heart_LV_6x(["chromVAR_obj.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none --> chromVAR_obj.ATAC.ENCODE_heart_LV_6x(["chromVAR_obj.ATAC.&lt;aggregation_name&gt;"]):::none
  chromVAR_obj.ATAC.ENCODE_heart_LV_6x(["chromVAR_obj.ATAC.&lt;aggregation_name&gt;"]):::none --> chromVAR_TF_motif_family_matrix.ATAC.ENCODE_heart_LV_6x(["chromVAR_TF_motif_family_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  peak_TF_motif_family_matrix.ATAC.ENCODE_heart_LV_6x(["peak_TF_motif_family_matrix.ATAC.&lt;aggregation_name&gt;"]):::none --> chromVAR_TF_motif_family_matrix.ATAC.ENCODE_heart_LV_6x(["chromVAR_TF_motif_family_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none
  within_clusters_collapsed_peaks_per_cluster_GRanges.ATAC.ENCODE_heart_LV_6x["within_clusters_collapsed_peaks_per_cluster_GRanges.ATAC.&lt;aggregation_name&gt;"]:::none --> consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none
  LSI_BPCells.ATAC.ENCODE_heart_LV_6x(["LSI_BPCells.ATAC.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  JASPAR2026_vertebrate_familial_root_motifs_tf(["JASPAR2026_vertebrate_familial_root_motifs_tf"]):::none --> JASPAR_familial_root_motif_matrix_list(["JASPAR_familial_root_motif_matrix_list"]):::none
  JASPAR_familial_root_motif_matrix_list(["JASPAR_familial_root_motif_matrix_list"]):::none --> JASPAR_motif_family_members_tibble(["JASPAR_motif_family_members_tibble"]):::none
  JASPAR2026_vertebrate_motif_families_tsv(["JASPAR2026_vertebrate_motif_families_tsv"]):::none --> JASPAR_motif_family_members_tibble(["JASPAR_motif_family_members_tibble"]):::none
  harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none --> knn.ATAC.ENCODE_heart_LV_6x(["knn.ATAC.&lt;aggregation_name&gt;"]):::none
  knn.ATAC.ENCODE_heart_LV_6x(["knn.ATAC.&lt;aggregation_name&gt;"]):::none --> LSI_clusters.ATAC.ENCODE_heart_LV_6x(["LSI_clusters.ATAC.&lt;aggregation_name&gt;"]):::none
  scDblFinder_results_df.ATAC.ENCODE_heart_LV_6x(["scDblFinder_results_df.ATAC.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  JASPAR_familial_root_motif_matrix_list(["JASPAR_familial_root_motif_matrix_list"]):::none --> peak_TF_motif_family_matrix.ATAC.ENCODE_heart_LV_6x(["peak_TF_motif_family_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none --> peak_TF_motif_family_matrix.ATAC.ENCODE_heart_LV_6x(["peak_TF_motif_family_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  peak_calling_cut_sites_per_GEM_well.fragments.ATAC.ENCODE_heart_LV_6x["peak_calling_cut_sites_per_GEM_well.fragments.ATAC.&lt;aggregation_name&gt;"]:::none --> peaks_per_cluster_narrowPeaks.peaks.ATAC.ENCODE_heart_LV_6x["peaks_per_cluster_narrowPeaks.peaks.ATAC.&lt;aggregation_name&gt;"]:::none
  knn.ATAC.ENCODE_heart_LV_6x(["knn.ATAC.&lt;aggregation_name&gt;"]):::none --> UMAP_embeddings_tibble.ATAC.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none --> UMAP_embeddings_tibble.ATAC.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  peaks_per_cluster_narrowPeaks.peaks.ATAC.ENCODE_heart_LV_6x["peaks_per_cluster_narrowPeaks.peaks.ATAC.&lt;aggregation_name&gt;"]:::none --> within_clusters_collapsed_peaks_per_cluster_GRanges.ATAC.ENCODE_heart_LV_6x["within_clusters_collapsed_peaks_per_cluster_GRanges.ATAC.&lt;aggregation_name&gt;"]:::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> peak_calling_cut_sites_per_GEM_well.fragments.ATAC.ENCODE_heart_LV_6x["peak_calling_cut_sites_per_GEM_well.fragments.ATAC.&lt;aggregation_name&gt;"]:::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> peaks_per_cluster_narrowPeaks.peaks.ATAC.ENCODE_heart_LV_6x["peaks_per_cluster_narrowPeaks.peaks.ATAC.&lt;aggregation_name&gt;"]:::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> within_clusters_collapsed_peaks_per_cluster_GRanges.ATAC.ENCODE_heart_LV_6x["within_clusters_collapsed_peaks_per_cluster_GRanges.ATAC.&lt;aggregation_name&gt;"]:::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> LSI_BPCells.ATAC.ENCODE_heart_LV_6x(["LSI_BPCells.ATAC.&lt;aggregation_name&gt;"]):::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.ATAC.ENCODE_heart_LV_6x(["scDblFinder_results_df.ATAC.&lt;aggregation_name&gt;"]):::none
  chromVAR_obj.ATAC.ENCODE_heart_LV_6x(["chromVAR_obj.ATAC.&lt;aggregation_name&gt;"]):::none --> motif_family_accessibility_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["motif_family_accessibility_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.GEX.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.GEX.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  LSI_clusters.ATAC.ENCODE_heart_LV_6x(["LSI_clusters.ATAC.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.ATAC.ENCODE_heart_LV_6x(["scDblFinder_results_df.ATAC.&lt;aggregation_name&gt;"]):::none
  LSI_clusters.ATAC.ENCODE_heart_LV_6x(["LSI_clusters.ATAC.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.ATAC.ENCODE_heart_LV_6x(["scDblFinder_results_df.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.ATAC.ENCODE_heart_LV_6x(["scDblFinder_results_df.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.ATAC.ENCODE_heart_LV_6x(["scDblFinder_results_df.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none --> consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none --> chromVAR_obj.ATAC.ENCODE_heart_LV_6x(["chromVAR_obj.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none --> peak_TF_motif_family_matrix.ATAC.ENCODE_heart_LV_6x(["peak_TF_motif_family_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none --> LSI_BPCells.ATAC.ENCODE_heart_LV_6x(["LSI_BPCells.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["consensus_peak_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  consensus_peak_GRanges.ATAC.ENCODE_heart_LV_6x(["consensus_peak_GRanges.ATAC.&lt;aggregation_name&gt;"]):::none --> LSI_BPCells.ATAC.ENCODE_heart_LV_6x(["LSI_BPCells.ATAC.&lt;aggregation_name&gt;"]):::none
  chromVAR_TF_motif_family_matrix.ATAC.ENCODE_heart_LV_6x(["chromVAR_TF_motif_family_matrix.ATAC.&lt;aggregation_name&gt;"]):::none --> motif_family_accessibility_BPCells_matrix_dir.ATAC.ENCODE_heart_LV_6x(["motif_family_accessibility_BPCells_matrix_dir.ATAC.&lt;aggregation_name&gt;"]):::none
  LSI_BPCells.ATAC.ENCODE_heart_LV_6x(["LSI_BPCells.ATAC.&lt;aggregation_name&gt;"]):::none --> scDblFinder_results_df.ATAC.ENCODE_heart_LV_6x(["scDblFinder_results_df.ATAC.&lt;aggregation_name&gt;"]):::none

WNN

Checkpoint: 8: WNN integration

Weighted nearest neighbors

GEX and ATAC representations are combined with a BPCells-based implementation of the Seurat weighted-nearest-neighbor (WNN) strategy (Hao et al. 2021; Parks and Greenleaf 2025). The native implementation returns neighbor indices, distances and modality weights that clustering, UMAP and SCAVENGE (Yu et al. 2022) reuse, without building a Seurat object. The GEX and ATAC embeddings, after optional Harmony correction, are aligned by barcode and L2-normalized. For each modality separately, 200 candidate neighbors are then found with the same reproducible hierarchical navigable small world (HNSW) search as for clustering, here with Euclidean distance on the normalized embeddings (Malkov and Yashunin 2020). multiomeR sets each cell’s kernel width in each modality from the configured neighbor count, following Seurat’s small-SNN bandwidth strategy. Per-cell modality weights are computed from these kernels by comparing within-modality with cross-modality predictions. The candidates of both modalities are then pooled, and those with the highest weighted kernel scores become the final neighbors. In our tests, Seurat’s weights approached 0 or 1 as the number of nuclei grew, because the kernel widths shrink with density. The weights are therefore estimated within random subsets of at most 10,000 nuclei, whereas the neighbors are selected among all nuclei. BPCells builds an SNN graph from the joint neighbors for Leiden clustering, and UMAP uses the joint neighbor lists directly. The joint lists first keep twice the configured number of neighbors. After nuclei in small WNN clusters are removed, each remaining nucleus keeps its highest-scoring remaining neighbors, so no second weight estimate is needed. Only nuclei left with too few neighbors are searched again. Nuclei left with too few are reselected among the remaining nuclei by exact search, with their original weights and kernel widths. The implementation does not use Seurat::FindMultiModalNeighbors(). Because HNSW replaces Seurat’s Annoy search, exact equality is not expected. Its validation instead requires the modality weights to correlate with Seurat’s, and the neighbors to overlap with Seurat’s, above preset thresholds.

The WNN implementation is compared with Seurat in Algorithm validation.

Final cell set

WNN analysis uses the nuclei kept after ATAC processing that are present in both embeddings. When small joint clusters are removed, each kept nucleus keeps its highest-scoring kept neighbors and its original modality weights, while the cluster labels are kept. A nucleus left with too few neighbors has them selected again among the kept nuclei, with the same weights. The UMAP is computed once, from the final neighbors. The GEX, ATAC and WNN metadata therefore describe successively smaller sets of nuclei, not interchangeable versions of the Cell Ranger cell calls.

Source and target graph

Source: extra_targets/WNN_targets.R, R/processing_multimodal_helpers.R, src/wnn_snn_bandwidth.cpp.

This view combines the Harmony-corrected GEX and ATAC embeddings into the WNN results, clusters, UMAP embeddings and final annotated metadata.

%%{init: {
  'theme': 'base',
  'themeVariables': {
    'primaryTextColor': '#1f2937',
    'lineColor': '#4b5563',
    'clusterBkg': '#f8fafc',
    'clusterBorder': '#cbd5e1',
    'edgeLabelBackground': '#ffffff00',
    'fontFamily': 'Helvetica, Arial, sans-serif'
  },
  'flowchart': {
    'defaultRenderer': 'elk'
  },
  'themeCSS': '.cluster rect { fill: #f8fafc !important; stroke: #cbd5e1 !important; } .cluster-label text, .cluster-label span { font-weight: 700 !important; font-size: 16px !important; fill: #475569 !important; color: #475569 !important; }'
}}%%
%%{init: {'flowchart': {'defaultRenderer': 'elk'}}}%%
flowchart TB
  harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none --> embedding_matrices.WNN.ENCODE_heart_LV_6x(["embedding_matrices.WNN.&lt;aggregation_name&gt;"]):::none
  metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none --> embedding_matrices.WNN.ENCODE_heart_LV_6x(["embedding_matrices.WNN.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none --> embedding_matrices.WNN.ENCODE_heart_LV_6x(["embedding_matrices.WNN.&lt;aggregation_name&gt;"]):::none
  WNN_results.ENCODE_heart_LV_6x(["WNN_results.&lt;aggregation_name&gt;"]):::none --> UMAP_embeddings_tibble.WNN.ENCODE_heart_LV_6x(["UMAP_embeddings_tibble.WNN.&lt;aggregation_name&gt;"]):::none
  embedding_matrices.WNN.ENCODE_heart_LV_6x(["embedding_matrices.WNN.&lt;aggregation_name&gt;"]):::none --> WNN_results.ENCODE_heart_LV_6x(["WNN_results.&lt;aggregation_name&gt;"]):::none
  clusters_tibble.WNN.ENCODE_heart_LV_6x(["clusters_tibble.WNN.&lt;aggregation_name&gt;"]):::none --> WNN_results.ENCODE_heart_LV_6x(["WNN_results.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.ATAC.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.ATAC.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none
  clusters_tibble.WNN.ENCODE_heart_LV_6x(["clusters_tibble.WNN.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.WNN.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.WNN.&lt;aggregation_name&gt;"]):::none
  WNN_results.ENCODE_heart_LV_6x(["WNN_results.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.WNN.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.WNN.&lt;aggregation_name&gt;"]):::none
  metadata_w_cell_types_tibble.ATAC.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.ATAC.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.WNN.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.WNN.&lt;aggregation_name&gt;"]):::none
  UCell_GEX_marker_genes_list.ENCODE_heart_LV_6x(["UCell_GEX_marker_genes_list.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.WNN.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.WNN.&lt;aggregation_name&gt;"]):::none
  harmony_embeddings_matrix.GEX.ENCODE_heart_LV_6x(["harmony_embeddings_matrix.GEX.&lt;aggregation_name&gt;"]):::none --> metadata_w_cell_types_tibble.WNN.ENCODE_heart_LV_6x(["metadata_w_cell_types_tibble.WNN.&lt;aggregation_name&gt;"]):::none
  embedding_matrices.WNN.ENCODE_heart_LV_6x(["embedding_matrices.WNN.&lt;aggregation_name&gt;"]):::none --> clusters_tibble.WNN.ENCODE_heart_LV_6x(["clusters_tibble.WNN.&lt;aggregation_name&gt;"]):::none

Amemiya, Haley M., Anshul Kundaje, and Alan P. Boyle. 2019. “The ENCODE Blacklist: Identification of Problematic Regions of the Genome.” Scientific Reports 9 (1): 9354. https://doi.org/10.1038/s41598-019-45839-z.
Andreatta, Massimo, and Santiago J. Carmona. 2021. “UCell: Robust and Scalable Single-Cell Gene Signature Scoring.” Computational and Structural Biotechnology Journal 19: 3796–98. https://doi.org/10.1016/j.csbj.2021.06.043.
Benjamini, Yoav, and Yosef Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society: Series B (Methodological) 57 (1): 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.
Choudhary, Saket, and Rahul Satija. 2022. “Comparison and Evaluation of Statistical Error Models for scRNA-seq.” Genome Biology 23 (1): 27. https://doi.org/10.1186/s13059-021-02584-9.
Corces, M. Ryan, Jeffrey M. Granja, Shadi Shams, et al. 2018. “The Chromatin Accessibility Landscape of Primary Human Cancers.” Science 362 (6413): eaav1898. https://doi.org/10.1126/science.aav1898.
Fleming, Stephen J., Mark D. Chaffin, Alessandro Arduini, et al. 2023. “Unsupervised Removal of Systematic Background Noise from Droplet-Based Single-Cell Experiments Using CellBender.” Nature Methods 20 (9): 1323–35. https://doi.org/10.1038/s41592-023-01943-7.
Germain, Pierre-Luc. 2026. betterChromVAR. Bioconductor. https://doi.org/10.18129/B9.bioc.betterChromVAR.
Germain, Pierre-Luc, Aaron Lun, Carlos Garcia Meixide, Will Macnair, and Mark D. Robinson. 2022. “Doublet Identification in Single-Cell Sequencing Data Using scDblFinder.” F1000Research 10: 979. https://doi.org/10.12688/f1000research.73600.2.
Granja, Jeffrey M., M. Ryan Corces, Sarah E. Pierce, et al. 2021. “ArchR Is a Scalable Software Package for Integrative Single-Cell Chromatin Accessibility Analysis.” Nature Genetics 53 (3): 403–11. https://doi.org/10.1038/s41588-021-00790-6.
Hafemeister, Christoph, and Rahul Satija. 2019. “Normalization and Variance Stabilization of Single-Cell RNA-Seq Data Using Regularized Negative Binomial Regression.” Genome Biology 20 (1): 296. https://doi.org/10.1186/s13059-019-1874-1.
Hao, Yuhan, Stephanie Hao, Erica Andersen-Nissen, et al. 2021. “Integrated Analysis of Multimodal Single-Cell Data.” Cell 184 (13): 3573–3587.e29. https://doi.org/10.1016/j.cell.2021.04.048.
Huang, Xianjie, and Yuanhua Huang. 2021. “cellsnp-lite: An Efficient Tool for Genotyping Single Cells.” Bioinformatics 37 (23): 4569–71. https://doi.org/10.1093/bioinformatics/btab358.
Huang, Yuanhua, Davis J. McCarthy, and Oliver Stegle. 2019. “Vireo: Bayesian Demultiplexing of Pooled Single-Cell RNA-Seq Data Without Genotype Reference.” Genome Biology 20 (1): 273. https://doi.org/10.1186/s13059-019-1865-2.
Korsunsky, Ilya, Nghia Millard, Jean Fan, et al. 2019. “Fast, Sensitive and Accurate Integration of Single-Cell Data with Harmony.” Nature Methods 16 (12): 1289–96. https://doi.org/10.1038/s41592-019-0619-0.
Lause, Jan, Philipp Berens, and Dmitry Kobak. 2021. “Analytic Pearson Residuals for Normalization of Single-Cell RNA-seq UMI Data.” Genome Biology 22 (1): 258. https://doi.org/10.1186/s13059-021-02451-7.
Luecken, Malte D., Maren Büttner, Kridsadakorn Chaichoompu, et al. 2022. “Benchmarking Atlas-Level Data Integration in Single-Cell Genomics.” Nature Methods 19 (1): 41–50. https://doi.org/10.1038/s41592-021-01336-8.
Malkov, Yu A., and D. A. Yashunin. 2020. “Efficient and Robust Approximate Nearest Neighbor Search Using Hierarchical Navigable Small World Graphs.” IEEE Transactions on Pattern Analysis and Machine Intelligence 42 (4): 824–36. https://doi.org/10.1109/tpami.2018.2889473.
McInnes, Leland, John Healy, and James Melville. 2018. “UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction.” arXiv, ahead of print. https://doi.org/10.48550/arXiv.1802.03426.
Melville, James. 2024. RcppHNSW: Rcpp Bindings for Hnswlib, a Library for Approximate Nearest Neighbors. https://doi.org/10.32614/CRAN.package.RcppHNSW.
Melville, James. 2026. uwot: The Uniform Manifold Approximation and Projection (UMAP) Method for Dimensionality Reduction. https://doi.org/10.32614/CRAN.package.uwot.
Miao, Zhen, and Junhyong Kim. 2024. “Uniform Quantification of Single-Nucleus ATAC-seq Data with Paired-Insertion Counting (PIC) and a Model-Based Insertion Rate Estimator.” Nature Methods 21 (1): 32–36. https://doi.org/10.1038/s41592-023-02103-7.
Ogata, Jonathan D., Wancen Mu, Eric S. Davis, et al. 2023. “Excluderanges: Exclusion Sets for T2T-CHM13, GRCm39, and Other Genome Assemblies.” Bioinformatics 39 (4): btad198. https://doi.org/10.1093/bioinformatics/btad198.
Ovek Baydar, Damla, Ieva Rauluseviciute, Dina R Aronsen, et al. 2026. “JASPAR 2026: Expansion of Transcription Factor Binding Profiles and Integration of Deep Learning Models.” Nucleic Acids Research 54 (D1): D184–93. https://doi.org/10.1093/nar/gkaf1209.
Parks, Benjamin, and William Greenleaf. 2025. “Scalable High-Performance Single Cell Data Analysis with BPCells.” bioRxiv, ahead of print. https://doi.org/10.1101/2025.03.27.645853.
Patikas, Nikolaos, Hongcheng Yao, Roopa Madhu, Soumya Raychaudhuri, Martin Hemberg, and Ilya Korsunsky. 2026. “Integration of Large, Complex Single-Cell Datasets with Harmony2.” bioRxiv, ahead of print. https://doi.org/10.64898/2026.03.16.711825.
Schep, Alicia. 2025. motifmatchr: Fast Motif Matching in R. https://doi.org/10.18129/B9.bioc.motifmatchr.
Schep, Alicia N., Beijing Wu, Jason D. Buenrostro, and William J. Greenleaf. 2017. “chromVAR: Inferring Transcription-Factor-Associated Accessibility from Single-Cell Epigenomic Data.” Nature Methods 14 (10): 975–78. https://doi.org/10.1038/nmeth.4401.
Stuart, Tim, Avi Srivastava, Shaista Madad, Caleb A. Lareau, and Rahul Satija. 2021. “Single-Cell Chromatin State Analysis with Signac.” Nature Methods 18 (11): 1333–41. https://doi.org/10.1038/s41592-021-01282-5.
Thibodeau, Asa, Alper Eroglu, Christopher S. McGinnis, et al. 2021. “AMULET: A Novel Read Count-Based Method for Effective Multiplet Detection from Single Nucleus ATAC-Seq Data.” Genome Biology 22 (1): 252. https://doi.org/10.1186/s13059-021-02469-x.
Tirosh, Itay, Benjamin Izar, Sanjay M. Prakadan, et al. 2016. “Dissecting the Multicellular Ecosystem of Metastatic Melanoma by Single-Cell RNA-seq.” Science 352 (6282): 189–96. https://doi.org/10.1126/science.aad0501.
Traag, V. A., L. Waltman, and N. J. van Eck. 2019. “From Louvain to Leiden: Guaranteeing Well-Connected Communities.” Scientific Reports 9 (1): 5233. https://doi.org/10.1038/s41598-019-41695-z.
Yu, Fulong, Liam D. Cato, Chen Weng, et al. 2022. “Variant to Function Mapping at Single-Cell Resolution Through Network Propagation.” Nature Biotechnology 40: 1644–53. https://doi.org/10.1038/s41587-022-01341-y.
Zhang, Yong, Tao Liu, Clifford A. Meyer, et al. 2008. “Model-Based Analysis of ChIP-Seq (MACS).” Genome Biology 9 (9): R137. https://doi.org/10.1186/gb-2008-9-9-r137.