multiomeR 1.2.0

multiomeR 1.2.0 revises the statistical methods after an audit: GEX normalization, WNN modality weights, clustering, consensus peaks, doublet detection, and the differential, gene-set, composition, genetic-enrichment and peak–gene tests. Every aggregation's results change, from QC and cell-type labels to all module results, so plan a full rebuild after updating. The release also adds an optional WNN-specific marker list, per-cluster composition tables, lower memory use for CellBender input, targets stores that can move between disks, and the public configuration, figures and results of the multiomeR manuscript.

What you need to do

1. Update to this release. Releases tag stable, validated commits, so update to the tag rather than pulling main, which can carry development work that is not yet validated. Commit or stash local changes first; resolve any merge conflicts as usual. While no other R process uses the environment, run:

   cd multiomeR  # your clone
   git fetch origin tag v1.2.0
   git merge --ff --no-edit v1.2.0
   pixi install --locked --run-post-link-scripts
   pixi run --use-environment-activation-cache install-r-github-packages

   This installs Seurat 5.6.0 from GitHub, with its build dependencies r-pbapply and r-rcppthread from the lock file.

2. Give your configuration directory its own crew_controllers.R. The root crew_controllers.R is gone: each configuration directory now carries controllers that suit its aggregations, and the default moved to configuration/crew_controllers.R. If you selected your own directory in configuration.local, copy configuration/crew_controllers.R into it and adapt it, moving in any changes you made to the former root file.

3. Update moved paths in your configuration and scripts:

   • cfg_pipeline_parameters.tsv, QC_checkpoint_manifest.tsv, QC_metric_manifest.tsv and the demo download lists are in manifests/.
   • reference_metadata/ is resources/reference_metadata/; add new Cell Ranger ARC references there.
   • module_<name>/ is extra_targets/module_<name>/.
   • The example donor metadata tables formerly in example_data/ are in configuration/ and configuration_dev/, and the download script is in scripts/. example_data/ now only receives downloads; three unused tables were removed.
   • mixed_human_31x and the nested comparison_1x to comparison_20x aggregations moved to the new configuration_dev/, with their donor metadata and GWAS selection; select it to rebuild them. manuscript_figures/ is replaced by pipeline_manuscript_figures/.

4. Check your configuration against the changed defaults and checks:

   • The cell-cycle score is no longer regressed out by default. Set aggregation_SCT_regress_vars: [CC.Difference] to keep it.
   • A differential-analysis formula that uses cluster without cell_type_formula or random_effect now stops with an error: it treated each donor's cell types as independent samples. Rewrite such formulas.
   • CellBender counts are used whenever a GEM well gives a CellBender H5 file; GEM_well_add_cellbender is no longer needed and a leftover column is ignored. Set GEM_well_cellbender_h5_file to NA to keep Cell Ranger counts.
   • Marker symbols that several genes share now stop marker validation, as missing symbols do, because they were scored on only one of the genes. Replace them.
   • The example configurations now exclude nuclei with more than 5% of their peak counts in blacklisted regions (was 1%, which removed cell types unevenly) and AMULET multiplets at the authors' q < 0.01 (was 0.05). Configurations copied from the old examples keep the old limits; change them to follow the examples.

5. Update code that reads pipeline outputs:

   • Each unassigned cluster is labeled Unassigned_cluster_<id> instead of one merged Unassigned.
   • AMULET's fragment-count column is amulet_nFrags instead of nFrags, like the other AMULET columns.

6. Demo users: run pixi run --use-environment-activation-cache setup-demo again to download controller_PBMC_human.

What happens to your aggregations

• Everything after the per-GEM-well inputs rebuilds. Every GEM well reruns amulet_metrics_tibble and full_metadata_tibble, and every aggregation target rebuilds. Cell Ranger, CellBender, cellsnp-lite, Vireo and the other per-GEM-well inputs stay current. The first run downloads the AMULET excluded-region list.
• Cell-type labels may change. GEX normalization, WNN weights, clustering and doublet calls all change, and each unassigned cluster keeps its own label. Review the labels at checkpoints 3, 7 and 8 before the modules run; the new aggregation_WNN_marker_genes lets you tune the final labels without rebuilding the GEX and ATAC processing.
• Peaks change. Each group's peaks must be called in both random halves of its nuclei, so small groups lose their weakly supported peaks.
• Module results change: differential expression and accessibility, gene-set and composition tests, genetic enrichment, and peak–gene links, which are now calibrated and fewer.
• Targets that rerun and keep their results. Targets that open BPCells matrices or fragments now record project-relative paths, so they rerun with unchanged data and their consumers rerun with them: GEX_counts_BPCells_matrix, fragments_w_prefix_bpcells, aggregated_counts_BPCells_matrix.GEX, consensus_peak_BPCells_matrix.ATAC, motif_family_accessibility_BPCells_matrix.ATAC, gene_score_archr_BPCells_matrix.ATAC and the pseudobulk count matrices. GEM_well_config_tsv, GEM_well_metadata_tibble, GEM_well_annotation_metadata_tibble, GEX_counts_BPcells_matrix_dir, chromHMMs_list_proj.ATAC, the five downloaded reference-file targets, the file targets of the moved manifests, reference metadata and donor metadata, and each GEM well's Cell Ranger reference lookup (cellranger_reference_json_file) rerun and reproduce their results. The exclusion lists (excluded_barcodes_by_type_list, excluded_cellranger_only_barcodes_by_type_list) rerun for GEM wells whose QC list lacks not_found_in_GEX_matrix, which includes every GEM well of the public configurations; their changes stop at the checkpoint 1 plots and retention tables. The new composition tables build without rerunning anything else.
• Demo: its aggregation targets rebuild with the new controller_PBMC_human well.
• In configuration_dev/, the Figure 1 targets rebuild. Supplementary Figure S2 reads the recorded runtimes of the compared workflows from the store, so rebuild its data before rerunning targets of the comparison_ aggregations.

What changed

Every change under Processing, Quality control, Differential analyses, Genetic enrichment and Peak–gene links alters results.

Processing

• The BPCells backend's Pearson residuals use one overdispersion, theta = 100, for every gene and are clipped at ±√N, the analytic Pearson residuals of Lause et al. (2021). The per-gene method-of-moments theta counted differences between cell types as overdispersion and, with the ±10 clip, flattened marker variation: on a multi-tissue test set, 24 of 27 canonical markers now enter the 3,000 variable genes, against 10.
• aggregation_SCT_regress_vars is unset by default, so the cell-cycle score, mostly noise in post-mitotic nuclei, is no longer regressed out; CC.Difference remains an example value.
• WNN modality weights are estimated within seeded random subsets of at most 10,000 nuclei, so they no longer drift towards 0 or 1 as an aggregation grows. Neighbours are still chosen among all nuclei, and smaller aggregations are unchanged. Each nucleus is now first in its own neighbour list.
• Leiden clustering runs for ten iterations, as Seurat does, instead of igraph's default two.
• Consensus peaks are built as ArchR does. Each group's discovery nuclei are split into two random halves, MACS3 calls peaks in each half, and only peaks called in both are kept; each group then keeps its strongest peaks up to 500 per nucleus and 150,000 in all, and overlapping peaks across groups are ranked by the quantile of their significance within their group. Deep groups no longer dominate, and small groups lose their weakly supported tail (a 1,200-nucleus group keeps 54% of its peaks). A half takes about half the time and memory of a full call, so MACS3 workers request 16 GB. Consensus peaks keep their group in cluster.
• The JASPAR familial root motifs lose flank columns with less than half of the largest column count; these low-support columns restricted the matches of 56 of 233 families to a sub-family.
• Each unassigned cluster is labeled Unassigned_cluster_<id> instead of all being merged into one Unassigned group, which peak calling, pseudobulks, composition and genetic enrichment analysed as a real cell type.

Quality control

• AMULET runs with chrM, chrX and chrY excluded in full and, for GRCh38, the authors' list of repeats, segmental duplications and blacklisted regions, downloaded once into the store and checked against its SHA-256 checksum. The example configurations exclude nuclei at the authors' q < 0.01 instead of 0.05.
• In insertions mode, the peak fraction is divided by twice the fragment count, so atac_peak_counts_frac is a fraction and filters on it take effect.
• ATAC scDblFinder groups peaks on singular-value-scaled LSI loadings, as scDblFinder itself does.
• The example configurations exclude nuclei with more than 5% of their peak counts in blacklisted regions instead of 1%, which removed cell types unevenly.
• Marker symbols that several genes share are rejected; they were scored on one of the genes.

Differential analyses

• GEX markers and ATAC differential peaks are tested on counts per 10,000 instead of raw counts, with Seurat-style fold changes, so they no longer follow each group's sequencing depth.
• Motif-family pseudobulks are tested on bias-corrected chromVAR deviations without centring or quantile normalization, which collapsed cell-type-defining families and produced false positives.
• Contrasts between annotation classes in the paired route are fitted on each donor's combined response, weighted by its precision weights, instead of using one residual correlation between classes for all features, which made low-expression features anti-conservative.
• The Hallmark and Reactome tests run with camera, estimating each set's inter-gene correlation, instead of cameraPR's default of 0.01; in permutations of a donor covariate, the old test called 17% of sets at P < 0.05. Sets need at least ten tested genes.
• The random-effect route uses edgeR's voomLmFit(), so precision weights account for the correlation block.
• A joint formula that uses cluster without cell_type_formula or random_effect stops with an error. The example formula now tests one annotation class.
• Composition models test by default only assigned cell types present in at least half of the donors.

Genetic enrichment

• The cell-type pseudobulks are summed by the final WNN labels, as SCAVENGE, the differential analyses and the peak–gene module do, instead of by the GEX labels carried into the WNN metadata, so labels set only at the WNN level reach the pseudobulk enrichment too.
• gchromVAR's background variance is scaled with squared peak weights, the variance of a weighted sum, in both the cell-type contributions and the single-nucleus z-scores. Summed weights deflated the z-scores by a trait-dependent factor whenever fine-mapping probabilities were below one.

Peak–gene links

• The model includes up to ten principal components of each cell type's residual peak accessibility, one per ten residual degrees of freedom when fewer, from the new target peak_gene_correlation_latent_factors.WNN, so peaks and genes that vary along the same chromatin states no longer appear linked.
• Each candidate pair is calibrated against its own gene, following the idea of Signac's LinkPeaks: up to 200 of the branch's candidate peaks more than 10 Mb from the gene's TSS are scored with the same model, the pair's statistic is standardized by this background (hierarchical_gene_z) and rescaled within peak-detection quintiles, and its empirical P-value against the cell type's background is Benjamini–Hochberg-corrected. Links also need a positive gene z. Before, the model P-values were corrected directly, and random pairs on different chromosomes passed the FDR cutoff almost as often as cis pairs. On held-out background pairs, the new P-values fall below 0.05 and 0.01 within 5% of the nominal rates.
• The random-slope test requires within-donor peak variation in at least eight donors and leaves unassigned clusters out.
• SuSiE's purity filter (minimum absolute correlation 0.5) is applied, so credible sets no longer span most candidate peaks.

Options and outputs

• The optional aggregation_WNN_marker_genes replaces, adds or, given as ~, drops marker sets for the WNN labels only, with their own matched controls. The final cell types can then be tuned without rebuilding the GEX and ATAC processing, whose doublet detection and peak calling use the GEX labels.
• cluster_composition_table.tsv is written at checkpoints 3, 7 and 8: for each cluster and each categorical variable, GEM_well_ID and donor_id, the cluster size, its most frequent category and that category's share, and every category above 5%.
• The peak–gene results gain the gene-standardized statistic hierarchical_gene_z, the calibrated hierarchical_background_pvalue and each cell type's n_background_pairs.
• crew_controllers.R is read from the selected configuration directory, and CellBender counts follow the H5 file alone (see What you need to do).
• The new configuration_dev/ builds the manuscript's results in its own targets store (TAR_PROJECT=dev, store outputs_dev), with Slurm workers on one node. configuration/ keeps the demo and the small examples, and validation covers whichever validation aggregations the selected directory defines.

Performance and fixes

• CellBender matrices are read in blocks of 2,000 barcodes, keeping only the gene features, so a 34,498-nucleus GEM well peaks at 5.2 GB instead of 17.1 GB. The written matrices are byte-identical.
• Seurat 5.6.0, installed from GitHub, has threaded C++ kernels for normalization, PCA, neighbour search and UMAP; they get the cores of each Slurm worker through Seurat.nthreads.
• Opened BPCells matrices and fragments and downloaded reference files are recorded by unresolved, project-relative paths, so a targets store keeps working after it moves, for example from node-local disk to shared storage.
• Roadmap chromHMM annotations load per aggregation, so activating another aggregation or selecting another configuration directory no longer rebuilds them.
• Nuclei missing from the GEX count matrix are always excluded at checkpoint 1, so the checkpoint 1 plots and retention tables report them instead of the checkpoint 3 retention table counting them as small GEX clusters.
• All unassigned clusters are drawn as one track in the ATAC coverage plots.

Demo

• controller_PBMC_human, 10x Genomics' 10k PBMC set from the Chromium Controller, replaces unsorted_PBMC_human. The unsorted 3k well came from the same donor as healthy_PBMC_human, so the demo's donors pbmc1 and pbmc2 are now two people, and its categorical UMAPs and composition tables show the donor.
• Harmony corrects for the GEM well, which now separates donors and runs, and cell types need a 0.02 margin over the marker background instead of 0.05, which the T-cell clusters, inflated by the larger well's T cells, did not reach. The NK markers drop KLRD1, which effector CD8 T cells share, for SH2D1B, and a MAIT set (SLC4A10, ZBTB16, KLRB1, NCR3) labels the MAIT cells. With the revised methods, three small WNN clusters (725 of 7,838 nuclei) match no marker set above background and stay unassigned.
• The download grows from about 1.3 to 4.6 GB, the stored results from about 3.3 to 6.6 GB and the run from about 20 to 30 minutes with 16 threads.

Manuscript

• The six 10x Genomics PBMC GEM wells of mixed_human_31x are labeled with their two donors: pbmc1 for the four from one woman and pbmc2 for the two from one man, as 10x Genomics' dataset descriptions, sex markers and mitochondrial genotypes agree. They were six donors before. The comparison_ copies keep their labels, since those aggregations run no donor-level analyses.
• pipeline_manuscript_figures/ builds Figure 1, Supplementary Figure S2 and Supplementary Table S4 of the manuscript from public data in configuration_dev/, with the results and their data under results/.
• Figure 1, Supplementary Figure S2 and Supplementary Table S4 are rebuilt with the revised methods. Figure 1A–B now shows a peak 72 kb from SGPL1 in NK cells, its trait panels show assigned cell types only, and figure_1_data writes the values the manuscript quotes. Panels C to E spell cell types alike, and panel D shows protein-coding genes with labels that no longer overlap.
• Supplementary Figure S2 compares multiomeR with a current Seurat v5/Signac v2 workflow on BPCells-backed assays, counting only the analyses both workflows perform. Its peak quantification uses fragtk, which pixi run --use-environment-activation-cache install-fragtk installs.

Documentation

• The In depth pages and the shared methods text describe the revised methods, and the checkpoint 8 guidance mentions the WNN marker list.
• New guidance covers aggregations that combine tissues, previewing the per-GEM-well targets of checkpoint 1, and leaving the costly Seurat exports out until the cell-type labels are accepted.
• The manual describes the per-cluster composition tables.
• The workflow overview figure no longer shows the dropped HC3 scan.
• Bibliography entries are complete and case-protected.
• CITATION.cff is added; GitHub shows it as "Cite this repository", and Zenodo uses it for the authors of archived releases. The README and the manual's home page ask users to cite multiomeR, targets and BPCells, and the methods behind the results they report; CITATION.cff lists targets and BPCells as references.

Development

• The WNN and cell-cycle reference tests require the pinned Seurat 5.6.0, and the peak–gene parity test scores background pairs.
• validation/ and RELEASES.md moved into dev/, release validation reports are kept in dev/validation/reports/, the scripts that make Cell Ranger ARC test subsets moved to dev/, and the unused air.toml was removed.