NMF Meta-Program (MP) Analysis

Overview

NMF meta-program analysis asks a different question from clustering. Clustering groups cells
into discrete states; NMF factorises a cell-state expression matrix V ≈ W·H into

• W — a gene loading matrix (genes x rank): each column is one non-negative "program"
  (a set of genes that go up together), scored by contribution to the factor;
• H — a cell loading matrix (rank x cells): how strongly each cell expresses each program.

Programs are reusable, partly overlapping axes of variation (cell cycle, stress, hypoxia,
interferon, lineage programs…), not exclusive clusters — a cell can score high on several.
Because a single sample's programs are noisy, the standard workflow (Tirosh-lab lineage,
Gavish et al., Nature 2023, doi:10.1038/s41586-023-06130-4 — 41 consensus MPs across
1,163 tumours, K = 4–9 per sample) is:

1. run NMF per sample over a range of ranks, giving many candidate programs;
2. keep only robust programs (recovered at several ranks / overlapping programs);
3. cluster the surviving programs by their shared top genes (Jaccard/overlap) and merge
   each cluster into one meta-program (MP) = the genes most frequently shared inside it.

Upstream package: NMF (Gaujoux & Seoighe 2010) — <https://cran.r-project.org/package=NMF>.

Inputs

The reference run in this container: neutrophil.rds, 34,058 genes x 1,449 cells,
samples NP10-AMZ-2 (72 cells), NP11-DP-1 (68), NP12-AMZ-1 (1,309); meta.data columns
include sample, percent.mt, sub, group, response, and reductions pca, harmony, umap.

Steps

Reusable, fully command-line driven script: scripts/nmf_metaprogram.R.txt
(no hard-coded paths, no source() of external helper files — the palette and the
robust-program filter are inlined).

1. Inspect the object first (column names are the #1 failure source):

r srt <- readRDS("neutrophil.rds"); print(colnames(srt@meta.data)) print(names(srt@reductions)); print(table(srt$sample))

2. Run the pipeline:

bash Rscript scripts/nmf_metaprogram.R.txt \ --input-rds /path/to/neutrophil.rds \ --output-dir /path/to/out \ --subset-col sub --subset-value Neutrophils \ --sample-col sample --min-cells 30 \ --ranks 4:9 --nrun 10 --seed 123 --cores 8 \ --top-genes 50 --reduction umap --group-col group

The per-sample core (per sample: subset → NormalizeData → ScaleData(hvg) → clip → NMF):

r hvg <- VariableFeatures(srt)[1:2000] # after FindVariableFeatures(method="vst") sub <- NormalizeData(sub, normalization.method = "LogNormalize", scale.factor = 10000) sub <- ScaleData(sub, features = hvg, vars.to.regress = "percent.mt") d <- as.matrix(SeuratObject::LayerData(sub, "scale.data")) d[!is.finite(d)] <- 0; d[d < 0] <- 0 # NMF refuses negatives d <- d[rowSums(d) > 0, ] # NMF refuses all-zero rows (empty program) set.seed(123) fit <- NMF::nmf(d, rank = 4:9, method = "snmf/r", nrun = 10) # -> NMF.rank object W <- lapply(fit$fit, NMF::basis) # fit$fit[["4"]] is NMFfitX1; W = genes x 4

3. Programs → MPs (what the script does after the loop): top-50 genes per program
   (uppercased), robust filter, overlap matrix, hclust on 2*top - overlap, reorder columns,
   then a greedy expansion — start a cluster from the most-connected program, add programs with
   >= --min-intersect-cluster shared genes, redefine the cluster gene set as the 50 most
   frequently shared genes (border ties broken by raw W scores), repeat while the best
   remaining program still has > --min-group-size partners.

r nmf_programs <- lapply(out1, function(x) apply(x, 2, function(y) toupper(names(sort(y, decreasing = TRUE))[1:50]))) sel <- robust_nmf_programs(nmf_programs, intra_min = 35, intra_max = 10, inter_filter = TRUE, inter_min = 10) inter <- crossprod(program_indicator) # == pairwise length(intersect(a, b)) hc <- hclust(as.dist(2 * 50 - inter), method = "average")

4. Read the results (see Outputs) and, if too few/many MPs, relax or tighten the thresholds
   in proportion to --top-genes (see Pitfalls).

5. Map MPs onto cells (optional, on by default): AddModuleScore with each MP gene set →
   Stat/mp_cell_scores.csv + per-MP UMAP pages + violin per --group-col.

Parameters

Outputs

<out>/Rdata/<sample>_SP_NMF.RData          # NMF.rank object per sample (object name: NMFs_per_sample)
<out>/Rdata/<prefix>.Genes.nmf.All.RData   # out1: named list of W matrices, genes x programs
<out>/Rdata/<prefix>.MP_list.RData         # MP_list (MP -> genes), Cluster_list (MP -> programs), params
<out>/Stat/<prefix>.MPs.csv                # MP x top genes (wide, NA-padded)   <- the deliverable table
<out>/Stat/mp_gene_list.csv                # long: MP, rank_in_MP, gene
<out>/Stat/program_cluster_membership.csv  # MP, program (<sample>_rank<k>_<i>), sample
<out>/Stat/jaccard_intersect_matrix.csv    # programme x program shared-gene counts (ordered by cluster)
<out>/Stat/jaccard_long.csv                # long + jaccard_pct = 100*x/(2*top - x)
<out>/Stat/rank_measures.csv               # per sample x rank: evar, cophenetic, residuals, silhouette.consensus…
<out>/Stat/mp_sample_composition.csv       # how many programs of each sample build each MP (confounding check)
<out>/Stat/mp_cell_scores.csv              # cell, sample, [subset/group cols], MP scores (AddModuleScore)
<out>/Stat/nmf_metaprogram.log, sessionInfo.txt
<out>/Img/<prefix>_cells-2MetaProgram.pdf  # program similarity heat map (magma palette, Jaccard %)
<out>/Img/<prefix>_rank_survey.pdf         # cophenetic / evar / silhouette.consensus vs rank
<out>/Img/<prefix>_MP_sample_heatmap.pdf   # MP x sample program counts
<out>/Img/<prefix>_MP_umap.pdf             # one page per MP: cells coloured by MP score
<out>/Img/<prefix>_MP_score_by_<group>.pdf # violin of MP score per group (when --group-col given)

Structure of the per-sample object (this is the layout the original NP*_SP_NMF.RData files have):

load("<out>/Rdata/NP10-AMZ-2_SP_NMF.RData")        # NMFs_per_sample
class(NMFs_per_sample)          # "NMF.rank"  (from nmf(rank = c(4,...))) ; names: measures, consensus, fit
NMFs_per_sample$measures        # matrix, one row per rank: evar, cophenetic, residuals, silhouette.consensus, nrun…
NMFs_per_sample$fit$`4`         # NMFfitX1 (consensus of the 10 runs); slots: fit(NMF), method, seed, rng, nrun
W <- NMFs_per_sample$fit$`4`@fit@W   # == NMF::basis(fit), genes x 4, rownames = genes, non-negative
H <- NMFs_per_sample$fit$`4`@fit@H   # == NMF::coef(fit),  4 x cells

out1 (from <prefix>.Genes.nmf.All.RData) is the concatenation of all bases per sample with
column names <sample>_rank<k>_<i> (e.g. NP10-AMZ-2_rank4_1), 39 columns for ranks 4:9.
Those column names are what the merge step parses (sub("_rank.*$", "") → sample).

Pitfalls

• Randomness and seed. nrun > 1 runs multiple random starts, so programs shift with the RNG.
  Always fix --seed and record it (the run config is saved into <prefix>.MP_list.RData$params).
  Two runs that differ only by seed produce program sets that are highly similar but not identical;
  if you compare groups/conditions, use the same seed and the same --ranks.
• Rank: the scan is for candidate generation, not model selection. The workflow concatenates the
  programs of every rank and merges them, so a bigger rank gives more (finer) candidates rather
  than "the right K". Use Stat/rank_measures.csv / rank_survey.pdf as diagnostics: evar
  increases monotonically with rank (here 0.139 → 0.240 for ranks 4 → 9 — useless for selection),
  the NMF vignette's rule of thumb is cophenetic >= 0.9 for a reproducible rank (this dataset scores
  0.987–0.995 at every rank, so it does not discriminate here), and silhouette.consensus peaks at
  the most reproducible rank (rank 7, 0.867, in the reference run). Also remember the overlap
  thresholds are counts of shared genes, so larger ranks need proportionally larger thresholds.
• --top-genes and the overlap thresholds are coupled. 35 / 15 / 10 were tuned for
  --top-genes 50 on a multi-sample dataset. Lowering --top-genes to 20 while keeping them makes
  an MP impossible to form (max possible overlap is 20) — the script then falls back to
  same-sample-only filtering and warns, and may return 0 MPs. Scale the thresholds with --top-genes.
• Sample/patient confounding (most common scientific error). Every program carries its sample of
  origin, and a sample with 10x more cells contributes 10x more programs. Always check
  mp_sample_composition.csv: an MP built exclusively from one sample is sample-specific, not a
  consensus program. Long-tail sample sizes (e.g. 72 vs 1,309 cells) can be equalised with
  --max-cells-per-sample.
• Batch regression eats the signal. --regress percent.mt (the source default) removes
  mitochondrial-content-driven covariance; regression is per sample, and regressing out too many
  covariates flattens the programs. Use --regress none to keep raw scaled expression.
• Scaling + clipping = "up-regulation only". After z-scoring per gene, negative values are set to
  0, so a program describes genes above the sample mean; a gene below the mean contributes nothing.
  Genes with constant expression in a sample become all-zero rows and are dropped — the reference run
  keeps 1,221 of the 2,000 HVGs for the 72-cell sample. That is expected, not a bug.
• nmf() compatibility traps (all verified in this container, NMF 0.28):
  sparse dgCMatrix is rejected (no inherited method for signature x = "dgCMatrix") → always
  as.matrix(); a matrix with negatives aborts the fit (input matrix x contains some negative entries); all-zero rows/columns yield degenerate programs.
• Seurat v5 API traps. VariableFeatures(srt, selection.method = "vst") is defunct in
  SeuratObject >= 5.0.0 (the original script was written against an older version) → call
  FindVariableFeatures(...) then VariableFeatures(srt)[1:n]. Read layers with
  SeuratObject::LayerData(obj, "scale.data"), not @scale.data. A stale scale.data layer inherited
  from the input object triggers Different features in new layer data than already exists for scale.data; the script drops it before re-scaling.
• AddModuleScore can fail on small matrices: control genes are sampled without replacement inside
  expression bins, so ctrl = 100, nbin = 24 may be impossible when few genes are available; the script
  retries with decreasing ctrl/nbin. MP gene sets smaller than 5 matched genes are skipped (usually a
  species / gene-symbol naming mismatch, e.g. mouse Xkr4 vs human XKR4).
• The MP table is uppercased. Genes are upper-cased before comparing programs, so MPs.csv contains
  uppercase symbols (HES4, ISG15); map back with toupper(rownames(obj)) when scoring or plotting.
• Cost. length(ranks) * nrun NMF fits per sample, each on genes x cells. In the reference run
  the recorded per-rank cpu.all for the 72-cell sample (rank 4:9, nrun 10) sums to ~198 s
  (~3.3 min), and the 1,309-cell sample is far heavier — which is why the source's
  NP12-AMZ-1_SP_NMF.RData is 167 MB, ~20x the 72-cell object. Use --cores for the restarts and
  --max-cells-per-sample when a sample is much larger than the others.

Environment

• Interpreter used here: /opt/conda/bin/Rscript (R 4.5.3); libraries
  /opt/user_R_libs, /opt/public_R_libs, /opt/conda/lib/R/library.
• Installed and verified: NMF 0.28, Seurat 5.5.1, SeuratObject 5.4.0,
  SingleCellExperiment 1.32.0, ComplexHeatmap 2.26.1, tidyverse 2.0.0, ggplot2 4.0.3,
  reshape2 1.4.5, viridis 0.6.5, RColorBrewer 1.1.3, matrixStats 1.5.0, AnnotationDbi 1.72.0,
  future 1.75.0, doParallel 1.0.17, ggsci 5.2.0, tidydr 0.0.6.
• Not installed: doFuture — future::plan() does not parallelise NMF runs unless the foreach
  backend is registered and nmf(.pbackend = NULL) is used; the script instead uses NMF's own
  doParallel backend (--cores, .options = "P<cores>").
• The original analysis ran in a Singularity/Apptainer base image with conda activate R.4.4.3
  (per that project's README); the two helper files it source()d
  (custom_magma.R, robust_nmf_programs.R) are inlined in this script.

> Packaging note. Scripts and tables in this skill carry a .txt suffix on their
> original extension (e.g. .R.txt, .tsv.txt) because the SkillHub upload allowlist
> does not accept .R or .tsv. This is cosmetic: Rscript is extension-agnostic, so
> Rscript scripts/<name>.R.txt ... behaves exactly like .R, and tab-separated data
> parses identically from .tsv.txt. Keep the suffix when repackaging.