multiomeR 1.3.0

Update to multiomeR 1.3.1 instead. It adds fixes for Harmony with several covariates and for the differential analyses. Follow its notes, which point back here for the steps below.

multiomeR 1.3.0 counts ATAC fragments instead of insertions, makes clustering and UMAPs reproducible for any thread count, and roughly halves the time of a full build. A 91-sample analysis with all modules now takes an estimated 5.0 hours on the benchmark's workers, against 12.9 hours before. Genetic enrichment adds S-LDSC heritability enrichment per cell type, and its cell-type chromVAR test now takes P-values from moved copies of each GWAS's loci. Modules can now analyse chosen datasets only, donors can give several samples, and an aggregation can read several donor tables. Every aggregation's results change, so plan a full rebuild after updating. Four parameters and two ATAC QC metrics are renamed or removed, so update your configuration first.

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. This release needs Pixi 0.82.1 or later, so run pixi self-update first if pixi --version shows an older one. Commit or stash local changes first, and resolve any merge conflicts as usual. While no other R process uses the environment, run:

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

   The environment changes only its CA certificate bundle, which no longer trusts eight compromised eMudhra certificate authorities. Native kernels now compile once into your user cache folder instead of in every worker.

2. Split the cluster-annotation margin. aggregation_cluster_annotation_min_advantage becomes aggregation_GEX_cluster_annotation_min_advantage, aggregation_ATAC_cluster_annotation_min_advantage and aggregation_WNN_cluster_annotation_min_advantage, each 0.05 by default. Set all three to your old value, for example:

   aggregation_GEX_cluster_annotation_min_advantage: 0.02
   aggregation_ATAC_cluster_annotation_min_advantage: 0.02
   aggregation_WNN_cluster_annotation_min_advantage: 0.02

   With equal values, no target command changes. The WNN margin sets only the final labels, so you can later tune it without rebuilding the GEX and ATAC processing.

3. Remove three retired parameters: aggregation_ATAC_peak_matrix_mode, aggregation_UMAP_GEX_PCs and aggregation_UMAP_ATAC_PCs. Peaks now always count fragments. UMAPs reuse the clustering neighbours, so they use aggregation_GEX_data_PCs and aggregation_ATAC_data_PCs. Their defaults, 1:30 and 2:30, match the old UMAP defaults. A configuration that still sets a parameter of step 2 or 3 stops when the graph is built, with "defines parameter(s) absent from the manifest".

4. Move your ATAC QC filters to the fragment metrics. atac_peak_counts_frac and atac_peak_counts_blacklist_frac are gone. A filter that names them stops QC_excluded_BCs_list.ATAC with "object not found" when that target runs. The graph still builds, so check aggregation_QC_exclude_list_combined_object before you run. nCount_ATAC remains, but it now counts fragments, so an unchanged threshold on it removes more nuclei. These replacements of the former example filters exclude nearly the same nuclei (Jaccard 0.93–0.95 in three aggregations):

   Before, with the default insertion counting	After
   nCount_ATAC < 1000	atac_fragments_in_peaks < 590
   atac_peak_counts_frac < 0.1	atac_fragments_in_peaks_frac < 0.12
   atac_peak_counts_blacklist_frac > 0.05	atac_fragments_in_blacklist_frac > 0.012

   In another calibration, nCount_ATAC < 200 and atac_peak_counts_frac < 0.02 became atac_fragments_in_peaks < 125 and atac_fragments_in_peaks_frac < 0.025. Count limits thus fall to about 0.6 times their old value, and fraction limits rise about 1.2-fold. Confirm other values in the checkpoint 4 QC plots. If you counted with aggregation_ATAC_peak_matrix_mode: overlaps, which counted a fragment about once, keep about the same count and in-peak fraction limits. The blacklist fraction now divides by all fragments rather than by the peak counts, so for human data 0.05 becomes about 0.012 in either mode. It also depends on the genome: the mm10 blacklist covers 239 Mb against 72 Mb for GRCh38, so a typical mouse nucleus has 1.7% instead of 0.03% of its fragments there. The mouse demo therefore uses 0.025. Set a mouse limit from its distribution, never the human value. Also replace the old names in aggregation_continuous_vars and in your own scripts.

5. Turn on ATAC gene activity scores if you use them. They are now off by default. This drops their only plot, marker_gene_activity_dot_plot.7_ATAC_QC, and saves about 2.5 hours of worker time per 740,000-nucleus build. Set aggregation_ATAC_gene_activity_scores: true to keep the plot. The Seurat/Signac export does not use the scores.

6. Update code that reads pipeline outputs:

   • The cell metadata tables, metadata_w_clusters_tibble.* and metadata_w_cell_types_tibble.*, no longer carry UMAP coordinates. Join UMAP_embeddings_tibble.GEX, .GEX_non_harmony, .ATAC or .WNN, or read metadata_w_cell_types_analysis_tibble.*, which the plots use. The Seurat exports keep their UMAPs.
   • The ATAC QC columns are renamed as in step 4, and atac_peak_counts_blacklist is gone. atac_fragments_in_peaks is new. atac_fragments_in_blacklist_frac is part of each GEM well's metadata from checkpoint 1.
   • nCount_ATAC and the ATAC peak matrices count fragments, so their values are lower than the former insertion counts.
   • GEX count folders store unsigned integers instead of doubles. Targets that open them still return doubles. If you open the folders yourself, open_BPCells_counts_as_double() does the same.
   • The cell-type chromVAR deviation tables take z, z_p and z_q from the locus-shift test and gain background_z, null_mean and null_sd. Their support_label is gone. The peak, variant and locus contribution tables now split the deviation only, so they lose z_contribution, their z and support_label.

7. Give genetic-enrichment workers network access to Zenodo and the GWAS Catalog. S-LDSC runs in every aggregation with the genetic_enrichment module. It downloads a 1.5 GB LD reference once, and each GWAS's summary statistics once for all aggregations. To use other summary statistics, name a path, URL or GWAS Catalog accession in the study's optional summaryStatistics field. Studies without summary statistics stay as grey rows.

8. Coming from 1.1.0? Also follow steps 2 to 6 of the 1.2.0 notes.

What happens to your aggregations

• Every aggregation rebuilds, and its results change. The rebuild starts at each GEM well's full_metadata_tibble, which gains the blacklist fraction, and at the per-well peak matrices, which now count fragments. Review the labels at checkpoints 3, 7 and 8 before the modules run.
• What changes:
  • ATAC QC exclusions, as step 4 describes.
  • Peaks: each peak regains its first base, MACS3's maximum gap is fixed, and each peak-calling group draws its random halves anew once. 98.1% of consensus peaks are unchanged by the fixed gap. The new halves change the consensus peaks as much as a new seed would (Jaccard 0.93).
  • Clusters, UMAPs and cell-type labels: neighbour indexes are built on one thread, WNN draws one set of modality weights instead of two, and the GEX and ATAC UMAPs start from a PCA layout. Before, two identical rebuilds already gave different clusters (adjusted Rand index 0.85 for GEX and 0.77 for ATAC). Rebuilds now reproduce for any thread count.
  • All module results, through these inputs. The differential analyses and peak–gene tests themselves compute the same values faster. Log fold changes agree within 2.1e-9 and peak–gene P-values within 9.5e-14.
  • Cell-type chromVAR z-scores, P-values, q-values and heatmaps, as Genetic enrichment describes. Their deviations change only through their inputs.
• Targets added: S-LDSC and the locus-shift test build for every aggregation with the genetic_enrichment module. The first run downloads the LD reference and the summary statistics: on the benchmark, the LD reference took 2 minutes and the 26 GWAS sources 38 minutes of work in parallel branches. S-LDSC's own steps then took 33 minutes of work at 31 samples and 45 minutes at 91.
• Per-GEM-well targets that rerun and reproduce their results:
  • cellsnp-lite reruns once for every genotyped well, about 1 to 1.5 hours per large well on 6 cores. Vireo then reruns once with identical calls.
  • AMULET reruns once, since its native kernel now loads from the user cache.
  • Each well's GEX counts are rewritten as unsigned integers.
  • Cell Ranger, CellBender and the downloads stay current.
• Targets built once: gene annotations, blacklists and chromatin states per reference genome, in about 9 minutes.
• Targets removed: the four ATAC gene-activity targets of each aggregation unless you enable them (step 5), reference_Ensembl_annotations_GRanges_list and cross_UMAP_barcodes.ATAC. prune_orphaned_plots() lists the images that removed or renamed targets leave behind.
• Mouse demo: brain_mouse excludes nuclei whose blacklist fraction exceeds 0.025, about 1.5 times its median. It keeps 3,532 instead of 2,751 nuclei, because the former limit sat close to the mouse median and removed more.
• Demo: immune_human_2x keeps clusters of 50 nuclei or more, so its MAIT cluster, now 99 nuclei, stays. A WNN-only Cytotoxic_T set labels a 918-nucleus cluster of effector CD8 and other cytotoxic T cells. The demo keeps 7,887 instead of 7,838 nuclei, and 426 instead of 725 stay unassigned.

What changed

The changes under ATAC counting and QC, Clustering, UMAP and WNN, Peak calling, and Genetic enrichment alter results. The other sections keep results unless they say otherwise.

ATAC counting and QC

• The consensus-peak matrix counts each fragment once per peak that holds one of its cut sites. This is Paired Insertion Counting (Miao and Kim, 2024), BPCells' "fragments" mode and Signac's default. Insertion counting counted a fragment with both cut sites in one peak twice. aggregation_ATAC_peak_matrix_mode is removed.
• Amounts now count fragments everywhere. Positions, such as peak-calling pileups, TSS enrichment and coverage tracks, keep using cut sites.
• Peaks reach the matrix as genomic ranges. The data frame they came in before made BPCells read their 1-based starts as 0-based, so each peak lost its first base.
• The QC metrics count each fragment once from the stored fragments, however many regions it touches. atac_fragments_in_peaks and atac_fragments_in_peaks_frac are counted per aggregation. atac_fragments_in_blacklist_frac is counted per GEM well, so it no longer depends on the aggregation.
• nCount_ATAC stays the matrix column sum, the Seurat-style depth of exports and summaries.
• The example configurations filter on atac_fragments_in_peaks < 590, atac_fragments_in_peaks_frac < 0.12 and atac_fragments_in_blacklist_frac > 0.012, or > 0.025 for the mouse demo.

Clustering, UMAP and WNN

• One neighbour search per modality, in the new knn.GEX and knn.ATAC targets, serves clustering, the main UMAP and the full-dimension sweep (recall 0.9989 against exact neighbours). A resolution change now reruns only Leiden.
• The HNSW index is built on one thread and searched on all threads. Threaded builds made identical rebuilds disagree. Single-threaded builds reproduce for any thread count, at about 3 more minutes per clustering step.
• UMAPs run uwot's batch optimisation on all threads. It also reproduces for any thread count, and takes 1.3 instead of 7 minutes per UMAP.
• The GEX and ATAC UMAPs start from a PCA layout, which scored better on every layout measure. For ATAC, global structure rose from 0.13 to 0.72.
• The WNN UMAP keeps its spectral start, now solved for the largest eigenvalues of the normalised adjacency. Its layouts match uwot's own start (correlation 0.999), in less time.
• GEX PCA and LSI read their input in 24 fixed chunks, so their results no longer depend on the thread count.
• GEX PCA and LSI compute at least 50 components. The elbow and loading plots then show whether further components carry structure, and raising the dimensions up to 50 does not recompute them.
• WNN keeps twice the neighbours in WNN_results_raw and restricts them to the retained nuclei instead of recomputing (overlap 0.993). Nuclei that the restriction leaves short of neighbours are reselected by exact search. WNN clusters and the final graph share one draw of modality weights instead of two.
• The UMAP sweeps use the same start and input as the main UMAPs. The ATAC sweep now runs on all nuclei, like the main ATAC UMAP. For both modalities, the sweep panel at the main settings equals the main UMAP.
• Without Harmony covariates, the non-Harmony GEX UMAP reuses the Harmony layout instead of computing it again.

Peak calling

• Cut sites are written once per GEM well and half, one chromosome at a time, and each MACS3 branch reads its half's files from all wells. On a 740,000-nucleus aggregation, this took 1.2 instead of 27.6 hours of work and 40 instead of 73 GB of store. A 24,000-nucleus well needs 2.9 instead of 27 GB of memory.
• MACS3's maximum gap is fixed at 150 bp, the extension size. MACS3 had set it from the first ten input lines, so peaks depended on line order (52–267 bp). 98.1% of consensus peaks are unchanged (base-pair Jaccard 0.985).
• Each group's random halves are seeded by the group's name. Dropping or reordering a group now keeps the other groups' halves, where before only 0 to 34 of 72 groups kept theirs.
• MACS3 writes its temporary pileups, at least 1.4 GB per half, to the worker's temporary folder (TMPDIR) instead of /tmp.

Genetic enrichment

• S-LDSC heritability enrichment: stratified LD score regression tests whether the SNPs in each WNN cell type's accessible peaks explain more heritability per SNP than SNPs in peaks in general. It uses the genome-wide polygenic signal, so it complements chromVAR's fine-mapped loci.
  • Each SNP in a consensus peak carries the peak's share of the cell type's pseudobulk accessibility. A gene annotation marks gene body ±100 kb of the cell type's 10% most specifically expressed protein-coding genes.
  • Each cell type is fitted with an all-peaks annotation and the 97 baseline-LD v2.2 annotations, the LDSC-SEG design, and also with genes and with both jointly. LD comes from 1000 Genomes Europeans on GRCh38. The heatmap marks GWAS of other ancestries.
  • Results give tau*, z, P, q and each GWAS's top set of cell types. Further plots compare S-LDSC with chromVAR, show how much of each coefficient its most influential blocks carry, and compare annotations of matched size.
  • The R implementation reproduces ldsc 1.0.1: LD scores to its three-decimal output, coefficients and standard errors to six digits.
  • Summary statistics come from the GWAS Catalog for a GCST sourceId, preferring harmonised GWAS-SSF files, or from the study's summaryStatistics field.
  • A cell type that expresses too few genes for a genome-wide annotation gets no gene or joint model, and its peak model still runs.
• Cell types with at least 100 nuclei: the cell-type chromVAR analyses and S-LDSC profile only cell types with at least 100 nuclei. Single unassigned nuclei, which an aggregation with a smaller minimum cluster size keeps, no longer distort the equal-weight background and the S-LDSC annotations, or make the cell-type figures unreadable. Aggregations with the default minimum cluster size of 100 are unaffected.
• Locus-shift test: the cell-type chromVAR z-score took its P-value from the normal tail. GWAS place several correlated peaks at each locus, and sums over a few cell-type-specific peaks are skewed, so that tail was not safe. Each GWAS is now scored against 10,000 copies whose credible sets move 0.5–2 Mb along their chromosome, keeping their spacing and fine-mapping probabilities. z_p is the empirical one-sided P-value. Held-out copies gave the nominal shares of P-values below 0.05, 0.01 and 0.001.
• Background variance: betterChromVAR 1.1.10 divides the background mean by the bin size twice, so its variance approaches the mean squared count. This shrank the z-scores of deep pseudobulks. The cell-type z-scores now use the variance of the sampled background peaks, which these sets match. Deviations are unchanged.
• Together, on the 31-sample example aggregation, 38 instead of 16 cell type–GWAS pairs reach q < 0.05, for 21 instead of 12 studies.
• Top sets: the chromVAR and S-LDSC heatmaps outline each GWAS's top set of cell types. Its anchor, outlined in black, is the significant cell type with the largest estimate. Cell types not significantly below it are outlined in grey. By the Bonferroni inequality, the set holds the cell type with the largest true value with at least 95% probability. One star marks P < 0.05 and two mark q < 0.05.

Speed and memory

Measured on the release candidate with the benchmark aggregations of configuration_dev/. Run times are estimated on the benchmark's fixed workers, as the earlier values were. Those came from a benchmark without S-LDSC, so the gains are understated:

• Like for like on 31 samples (376,000 nuclei), multiomeR with its Seurat exports has a critical path of 1.5 hours, against 2.2 hours before.
• A 91-sample analysis (1.16 million nuclei) with all modules takes an estimated 5.0 hours on the benchmark's workers, against 12.9 hours. It uses 51 instead of 115 CPU-hours, at most 88 instead of 104 GB in one worker, and 432 instead of 567 GB of disk.
• A 31-sample analysis builds in an estimated 1.1 instead of 1.8 hours. Adding a sample reruns 92% of a full build, against 102% before.

These changes keep their results:

• LSI holds its input in memory and takes 12–13 instead of 61 minutes on 786,000 nuclei. It picks the fastest layout whose projected memory fits its budget. GEX normalization takes 306 instead of 690 seconds.
• Critical path: cell typing, doublet calls, peak calling and the modules no longer wait for the UMAPs, about 31 minutes on 740,000 nuclei. Cell-cycle scores for the Seurat exports have their own target, which takes about 100 seconds off the GEX normalization.
• WNN: a compiled neighbour selection takes 9 instead of 30 minutes and 9 instead of 25 GB. Restricting to the retained nuclei takes 6 seconds instead of a 26-minute recomputation. Leiden graphs build in 3 instead of 35 seconds.
• Plots: continuous UMAPs read one slim input per family, which cuts their work from 6.6 to 0.85 hours on 740,000 nuclei. Each plot is built once, so large UMAPs save in about 4 instead of 7 seconds. UMAP sweeps compute their layouts in their own targets, so a label change no longer reruns 3.2 hours of sweep layouts.
• Metadata: stored metadata columns drop their repeated barcode names. The tables shrink from 363 to 193 MB, and they load in 2.3 instead of 16 seconds.
• Differential analyses: a closed-form, vectorised solution replaces limma's per-feature blocked fits, and contrasts are fitted together. results_tibble takes 223 instead of 1,879 seconds, and feature_matrix_fit 282 instead of 548 seconds.
• Genetic enrichment: locus summaries take 13 instead of 369 seconds. Single-nucleus z-scores run in 39 instead of 1,833 branches, and the variant detail plots take 1.15 instead of 2.79 hours. SCAVENGE and the variant details read slim cell tables.
• Peak–gene links: pseudobulks are summed in one pass, 0.3 instead of 4.5 hours of work. The latent factors come from per-chromosome cross products on a 16 GB instead of a 32 GB worker. Focal coverage for all top links takes 4.5 instead of 93 minutes of work. The REML fit and the Kenward–Roger test share one kernel, about twice as fast.
• Seurat/Signac export: each GEM well builds its fragment object once for all aggregations, which saved 2.8 hours of work on an 83-well aggregation.
• Genotyping: cellsnp-lite streams each chromosome (-T) and uses -R only at positions that the VCF lists more than once. This is 1.2–1.3 times faster with 22% less CPU, and the counts are byte-identical. Vireo runs one process on the 32 GB tier, since it ran single-threaded anyway.
• Storage: GEX counts are stored as unsigned integers, 4.8 times smaller on 50,000 nuclei.
• Workers: native kernels compile once per user. Before, they compiled in every worker, for 40–139 seconds each. Workers skip the driver-only start-up, and helpers travel without source references, 3 instead of 13 MB per task. No full garbage collection runs before each target, which took 0.2 of the 0.36 seconds of overhead per target.

New options

• Modules on chosen datasets: differential_analyses_GEM_well_datasets and peak_gene_correlation_GEM_well_datasets give a module only the nuclei and donor pseudobulks of the listed GEM_well_dataset values. Background tissues can then shape clustering and genetic enrichment without entering these analyses. Each donor's nuclei must come from the listed datasets or from none of them.
• Several samples per donor: an optional cfg_GEM_well_samples.tsv in the configuration directory maps each donor of a GEM well to a sample (GEM_well_ID, donor_id, sample_id). Samples become the replicates of the donor table, pseudobulks and differential analyses, while demultiplexing keeps the genetic donors. Each sample inherits its donor's row of the donor table. Add sample rows only for variables that differ, such as a condition. Sample IDs may not contain underscores.
• Several donor tables: aggregation_donor_id_metadata_tsv takes a list of files, for example one per cohort, and joins them into one table. Conflicting values for a donor stop with an error that names the files, donors and columns. A single file is read as before.
• cellsnp-lite thresholds per GEM well: the optional cfg_GEM_wells.tsv columns GEM_well_cellsnp_min_count and GEM_well_cellsnp_min_MAF default to the previous 20 reads and MAF 0.1. With donor genotypes supplied, 1 read and MAF 0 cut unassigned nuclei from over 90% to 10% in a well with weak ATAC enrichment.
• Annotation margins per modality (step 2) and optional ATAC gene activity scores (step 5).
• GWAS sources: a sourceId ending in .parquet names a local file, and any other ID names an Open Targets study. FinnGen studies such as FINNGEN_R12_T2D can therefore be configured.

Outputs and diagnostics

• Resolution sweep: checkpoints 3 and 7 add resolution_sweep_table.tsv and resolution_sweep_plot.png. They show cluster counts and sizes at 0.5, 1 and 2 times the configured resolution, from the cached neighbours.
• Components: the elbow, loading, spread and association plots show all 50 or more computed components.
• Genotyping power: the demultiplexing summary gives each genotyped well's SNPs that passed cellsnp-lite's filters, and the median number covered per nucleus. Wells with too little genotype data then stand out before pooling.
• cellsnp-lite stops when no SNP passes its filters and names the likely cause, contig names or genome build differing between VCF and BAM. Before, it passed and Vireo failed later.
• Plot files are tracked by output location, so a renamed branch replaces its predecessor's files. The new prune_orphaned_plots() lists the images under the store's plots/ that no current target owns. It deletes them only with delete = TRUE. In one shared store it found 20,003 of 25,233 images (7.35 GB).
• Run checks: list_unfinished_targets() tells whether a run built its selection, since tar_make() returns normally when targets error or workers crash. list_errored_targets() keeps the stored errors of the current graph, not only those of the latest run.
• Large plots: a PNG whose side would pass ragg's 50,000-pixel limit is saved at a lower resolution instead of failing. The faceted UCell advantage pages are saved at 150 dpi, and every page shows all three role keys.
• Start-up check: the project runtime stops with a repair hint when R's build configuration (Makeconf, ldpaths) is empty. Otherwise, every native kernel failed with an unrelated-looking error.
• Reference annotations are built once per reference instead of once per aggregation. An EnsDb loads only after its MD5 checksum matches, so a truncated download fails at once.

Fixes

• Marker tests no longer fail when a cell-type or cluster level loses all its nuclei to later filters.
• scDblFinder leaves nuclei without counts in the doublet features unscored, instead of stopping the GEM well's doublet detection.
• The motif heatmaps by GEX cluster and cell type leave out nuclei without finite chromVAR scores instead of stopping.
• The AnnotationHub cache takes a directory lock only. On NFS, two targets could hold the file lock and the fallback lock at once, so one read an annotation that the other was still downloading.

Documentation

• The methods describe fragment counting, with Paired Insertion Counting cited, and the shared neighbour searches, with RcppHNSW cited. They also cover the UMAP starts, the chunked LSI and GEX PCA, the closed-form blocked fits, and the WNN restriction with reselection.
• The minimum cluster size is described as the outlier filter it is. In a 740,000-nucleus aggregation, 7,693 of 7,810 raw clusters held fewer than 100 nuclei. 7,282 of them held one nucleus, and together they held 1.1% of the nuclei.
• The genetic-enrichment page tabulates the columns multiomeR reads from a local credible-set file, and the parameter description says which sourceId values name local files.
• The donor metadata page documents cfg_GEM_well_samples.tsv. The troubleshooting page covers prune_orphaned_plots() and list_unfinished_targets().
• The genetic-enrichment methods describe S-LDSC, the locus-shift test, the corrected background variance and the top sets. Algorithm validation records how S-LDSC was compared with ldsc 1.0.1 and with published z-scores.
• The graph views, demo outputs, output gallery, defaults snapshot and parameter browser are refreshed.

Development

• The agent skills share their wording with the private repository's. A new rebuild-impact skill holds the rerun rules, and other skills drop generic guidance. The targets skills record what agents kept getting wrong. Claude Code finds the skills through .claude/skills.
• The validation preflight checks that the pinned S-LDSC reference on Zenodo answers.
• The peak–gene parity test checks the fused kernel against lme4 and pbkrtest.
• The release validation passed every check on the release commit. It finished after publication, so its report, VALIDATION-1.3.0.md, is on main rather than in the tag.