Learning goals

This tutorial explains two distinct applications of NMF in Yan et al. (2026): one to patient-level cancer-cell pseudo-bulk profiles and a second to individual cancer cells within each patient. You will:

Build a genes × patients expression matrix.
Interpret its two NMF outputs: gene-to-program matrix W and program-to-patient matrix H.
Reconstruct observed expression using matrix multiplication.
Implement NMF from scratch in base R.
Aggregate single-cell counts into patient-level pseudo-bulk profiles.
Assign each patient to its dominant archetype and inspect biological signatures.
Perform per-patient NMF on individual cancer cells and compare discovered programs across patients.
Interpret how recurrent programs are combined into the 13 metaprograms in the paper.
Understand which steps reflect the actual Yan et al. (2026) methodology.

Important distinction: The paper performs NMF twice at different observational levels: (1) pseudo-bulk cancer cells across patients gives 4 archetypes; (2) single cancer cells within each patient gives numerous transcriptional programs that are consolidated into 13 metaprograms. 8 ecotypes come later from a different co-occurrence network analysis, not either NMF rank.

Paper: Yan, Y. et al. Ecotypes of triple-negative breast cancer in response to chemotherapy. Nature (2026). https://doi.org/10.1038/s41586-026-10469-9 (particularly Fig. 2 and Methods sections “Creating pseudo-bulk RNA-seq data,” “Performing NMF,” and “Identifying archetypes”).

Part I — NMF fundamentals

1. What is NMF?

NMF assumes that the observed expression matrix can be approximately written as the product of two non-negative matrices:

\[ \underbrace{V}_{\text{genes} \times \text{patients}} \approx \underbrace{W}_{\text{genes} \times \text{programs}} \underbrace{H}_{\text{programs} \times \text{patients}}. \]

V: Input expression values for each gene in each patient.
W: How much each gene contributes to each latent expression program.
H: How strongly each patient uses each latent program.
K: Number of latent programs. Here we start with \(K=2\) for learning; the study selected \(K=4\).

All values in \(V\), \(W\), and \(H\) are non-negative. Unlike a hard cluster assignment, each patient may have non-zero loadings on several programs. A dominant-program label is obtained by selecting the program with the largest loading for a patient.

2. An expression matrix you can understand by inspection

Imagine measuring four genes in four TNBC patients. Two genes represent proliferation and two represent interferon response. The following values are synthetic, not data from the paper.

V <- matrix(
    c(6, 3, 0, 0, 4, 2, 2, 1, 0, 0, 6, 3, 2, 1, 4, 2),
    nrow = 4,
    dimnames = list(
        c("MKI67", "TOP2A", "ISG15", "IFIT1"),
        c("P1", "P2", "P3", "P4")
    )
)
V
##       P1 P2 P3 P4
## MKI67  6  4  0  2
## TOP2A  3  2  0  1
## ISG15  0  2  6  4
## IFIT1  0  1  3  2

Patient 1 has high proliferation-gene expression, Patient 3 has high interferon-gene expression, and Patients 2 and 4 express both programs.

heatmap(
    V,
    Rowv = NA,
    Colv = NA,
    scale = "none",
    col = colorRampPalette(c("white", "skyblue", "steelblue4"))(40),
    margins = c(6, 7),
    main = "Observed expression matrix V"
)

3. Understand W and H before fitting anything

In this specially constructed example, an exact decomposition exists:

\[ W = \begin{bmatrix} 2 & 0 \\ 1 & 0 \\ 0 & 2 \\ 0 & 1 \end{bmatrix},\qquad H = \begin{bmatrix} 3 & 2 & 0 & 1 \\ 0 & 1 & 3 & 2 \end{bmatrix}. \]

Here columns of W are the shared expression patterns:

Program 1: MKI67 and TOP2A (proliferation).
Program 2: ISG15 and IFIT1 (interferon).

Rows of H give program contributions across patients.

W_known <- matrix(
    c(2, 1, 0, 0, 0, 0, 2, 1),
    nrow = 4,
    dimnames = list(rownames(V), c("Proliferation", "Interferon"))
)
H_known <- matrix(
    c(3, 0, 2, 1, 0, 3, 1, 2),
    nrow = 2,
    dimnames = list(colnames(W_known), colnames(V))
)
W_known
##       Proliferation Interferon
## MKI67             2          0
## TOP2A             1          0
## ISG15             0          2
## IFIT1             0          1
H_known
##               P1 P2 P3 P4
## Proliferation  3  2  0  1
## Interferon     0  1  3  2

Example calculation: reconstruct Patient 2

Patient 2 has loadings (2, 1), so the expression vector is

\[ 2\begin{bmatrix}2\\1\\0\\0\end{bmatrix} +1\begin{bmatrix}0\\0\\2\\1\end{bmatrix} =\begin{bmatrix}4\\2\\2\\1\end{bmatrix}. \]

This exactly matches the second column of V.

p2_reconstructed <- W_known %*% H_known[, "P2", drop = FALSE]
p2_reconstructed
##       P2
## MKI67  4
## TOP2A  2
## ISG15  2
## IFIT1  1

V_reconstructed <- W_known %*% H_known
V_reconstructed
##       P1 P2 P3 P4
## MKI67  6  4  0  2
## TOP2A  3  2  0  1
## ISG15  0  2  6  4
## IFIT1  0  1  3  2
all.equal(unname(V), unname(V_reconstructed))
## [1] TRUE

Takeaway: NMF represents each patient as an additive mixture of shared programs. Non-negativity makes negative cancellations impossible.

Experiment with patient loadings

Change the two values below and rerun the code to see how the resulting gene expression changes. If you use RStudio/Quarto, you can edit these parameters directly.

proliferation_loading <- 2
interferon_loading <- 1

new_H <- matrix(c(proliferation_loading, interferon_loading), ncol = 1)
rownames(new_H) <- colnames(W_known)
colnames(new_H) <- "NewPatient"
new_expression <- W_known %*% new_H

barplot(
    new_expression[, 1],
    names.arg = rownames(V),
    col = c("#447bb8", "#447bb8", "#ea9757", "#ea9757"),
    ylab = "Reconstructed expression",
    main = "Change program contributions, then rerun"
)

new_expression
##       NewPatient
## MKI67          4
## TOP2A          2
## ISG15          2
## IFIT1          1

4. How does the algorithm learn W and H?

In reality, we only observe V. The algorithm begins with non-negative guesses for W and H, then adjusts them to minimize the discrepancy between V and WH. One common objective is

\[ \min_{W,H\geq 0}\,\lVert V-WH\rVert_F^2 =\min_{W,H\geq 0}\sum_{g,p}(V_{gp}-(WH)_{gp})^2. \]

The following multiplicative updates are a simple educational way to solve the problem (they are not the exact solver used in the paper):

\[ H \leftarrow H\odot\frac{W^\mathsf{T}V}{W^\mathsf{T}WH+\epsilon}, \qquad W \leftarrow W\odot\frac{VH^\mathsf{T}}{WHH^\mathsf{T}+\epsilon}. \]

The products, fractions, and updates indicated by \(\odot\) and the division bars are applied element-wise where appropriate. \(\epsilon\) prevents division by zero.

Fit NMF from scratch with base R

set.seed(123)
K <- 2
n_iter <- 3000
eps <- 1e-10

W <- matrix(runif(nrow(V) * K, 0.1, 1), nrow = nrow(V))
H <- matrix(runif(K * ncol(V), 0.1, 1), nrow = K)
loss <- numeric(n_iter)

for (iter in seq_len(n_iter)) {
    # Update the patient loading matrix H
    H <- H * ((t(W) %*% V) / (t(W) %*% W %*% H + eps))

    # Update gene-program matrix W
    W <- W * ((V %*% t(H)) / (W %*% H %*% t(H) + eps))

    # Measure reconstruction error
    V_hat <- W %*% H
    loss[iter] <- sum((V - V_hat)^2)
}

# Optional rescaling removes an arbitrary factor scale.
# The product WH is unchanged by this operation.
factor_scale <- colSums(W)
W <- sweep(W, 2, factor_scale, "/")
H <- sweep(H, 1, factor_scale, "*")

rownames(W) <- rownames(V)
colnames(W) <- paste0("Program", seq_len(K))
rownames(H) <- colnames(W)
colnames(H) <- colnames(V)

round(W, 3)
##       Program1 Program2
## MKI67    0.667    0.000
## TOP2A    0.333    0.000
## ISG15    0.000    0.667
## IFIT1    0.000    0.333
round(H, 3)
##             P1 P2    P3 P4
## Program1 9.001  6 0.000  3
## Program2 0.000  3 9.001  6
round(W %*% H, 3)
##          P1 P2    P3 P4
## MKI67 6.000  4 0.001  2
## TOP2A 3.000  2 0.000  1
## ISG15 0.001  2 6.000  4
## IFIT1 0.000  1 3.000  2
tail(loss, 1)
## [1] 1.972334e-06

Inspect convergence

plot(
    loss,
    type = "l",
    lwd = 2,
    col = "#3b82bd",
    xlab = "Iteration",
    ylab = "Sum of squared reconstruction errors",
    main = "NMF convergence"
)

Interpretation: As iterations proceed, the reconstruction error should generally fall toward a low value. Factor indices can swap between runs (e.g., Program 1 and Program 2 exchanging order), and the W/H scaling is not unique. The learned matrices need not exactly match the numbers in the hand-worked example.

# Label each patient by its dominant learned factor
learned_labels <- apply(H, 2, which.max)
names(learned_labels) <- colnames(H)
learned_labels
## P1 P2 P3 P4 
##  1  1  2  2

Caution: NMF can reach different local solutions from different starting points. In practice, use multiple initializations and check stability, not just reconstruction error.

Part II — Patient-level NMF: discover differences between patients

5. Go from single cancer cells to patient pseudo-bulk

The study used cancer cells only, summed UMI counts gene by gene across the cancer cells of each patient, and derived a patient-level expression matrix. This avoids confounding from immune and stromal cells in conventional whole-tissue bulk expression.

For gene \(g\) and patient \(p\):

\[ X_{g,p} = \sum_{c\in\text{cancer cells from patient }p} X_{g,c}. \]

Here is an intentionally tiny two-patient example:

counts <- matrix(
    c(3, 1, 0, 0, 3, 2, 0, 0, 2, 1, 1, 0, 2, 1, 1, 1),
    nrow = 4,
    dimnames = list(
        c("MKI67", "TOP2A", "ISG15", "IFIT1"),
        c("P1_C1", "P1_C2", "P2_C1", "P2_C2")
    )
)
patient <- c("P1", "P1", "P2", "P2")

pseudobulk <- sapply(
    unique(patient),
    function(p) rowSums(counts[, patient == p, drop = FALSE])
)
counts
##       P1_C1 P1_C2 P2_C1 P2_C2
## MKI67     3     3     2     2
## TOP2A     1     2     1     1
## ISG15     0     0     1     1
## IFIT1     0     0     0     1
pseudobulk
##       P1 P2
## MKI67  6  4
## TOP2A  3  2
## ISG15  0  2
## IFIT1  0  1

What happens in a real project? Use raw UMI counts when aggregating. Normalize the resulting patient-level pseudo-bulk expression before NMF. The authors describe using a variance-stabilizing transformation (DESeq2 vst) when constructing pseudo-bulk RNA-seq data.

Paper-specific expression processing

The Methods describe the following steps for patient archetypes:

Construct cancer-cell-only patient pseudo-bulk profiles.
Retain the top 75% of genes based on mean- and variance-based expression criteria.
Prepare the normalized, log2-scaled expression matrix.
Center each gene across patients; replace values below zero with zero to make the input non-negative.
Fit NMF using RcppML v0.3.7 with tolerance 1e-5, L1 sparsity parameter 0.05, and candidate ranks K = 2 through 10.
Evaluate ranks using silhouette scores calculated from the patient-to-patient Spearman expression-correlation matrix and archetype assignments; select K = 4.
Use each patient’s highest H loading to assign ARC1–ARC4.

The exact gene filter and normalized values require the underlying data and complete computational setup. The next snippet demonstrates only the centering/truncation principle using our toy data; it is not a full reproduction of the paper.

# Educational preprocessing demonstration only:
# assume V has already been appropriately normalized/log-scaled
logV <- log2(V + 1)
centered <- sweep(logV, 1, rowMeans(logV), FUN = "-")
nonnegative_V <- pmax(centered, 0)

round(logV, 2)
##         P1   P2   P3   P4
## MKI67 2.81 2.32 0.00 1.58
## TOP2A 2.00 1.58 0.00 1.00
## ISG15 0.00 1.58 2.81 2.32
## IFIT1 0.00 1.00 2.00 1.58
round(centered, 2)
##          P1    P2    P3    P4
## MKI67  1.13  0.64 -1.68 -0.09
## TOP2A  0.85  0.44 -1.15 -0.15
## ISG15 -1.68 -0.09  1.13  0.64
## IFIT1 -1.15 -0.15  0.85  0.44
round(nonnegative_V, 2)
##         P1   P2   P3   P4
## MKI67 1.13 0.64 0.00 0.00
## TOP2A 0.85 0.44 0.00 0.00
## ISG15 0.00 0.00 1.13 0.64
## IFIT1 0.00 0.00 0.85 0.44
stopifnot(all(nonnegative_V >= 0))

Centering reveals expression above each gene’s across-patient average. Truncating negatives prevents below-average expression from becoming negative NMF inputs. Note that this transformation changes the interpretation relative to raw counts.

6. Assign four archetypes from H

In the study, each NMF factor was interpreted as an archetype. For this hypothetical example, consider the following illustrative normalized H matrix with four factors:

H4 <- rbind(
    ARC1 = c(0.80, 0.10, 0.05, 0.10),
    ARC2 = c(0.10, 0.70, 0.10, 0.10),
    ARC3 = c(0.05, 0.15, 0.75, 0.20),
    ARC4 = c(0.05, 0.05, 0.10, 0.60)
)
colnames(H4) <- paste0("P", 1:4)
H4
##        P1   P2   P3  P4
## ARC1 0.80 0.10 0.05 0.1
## ARC2 0.10 0.70 0.10 0.1
## ARC3 0.05 0.15 0.75 0.2
## ARC4 0.05 0.05 0.10 0.6

patient_labels <- apply(H4, 2, function(h) rownames(H4)[which.max(h)])
patient_labels
##     P1     P2     P3     P4 
## "ARC1" "ARC2" "ARC3" "ARC4"

Patient P3 has the largest loading on ARC3 (0.75), so P3 receives the ARC3 label. This is not a claim that the patient contains only that program.

barplot(
    H4,
    col = c("#dc6974", "#5689cb", "#55a786", "#e5a14e"),
    legend.text = rownames(H4),
    args.legend = list(x = "topright", bty = "n"),
    ylab = "Relative archetype loading",
    main = "Illustrative patient loading matrix H"
)

Why are these called archetypes? NMF only discovers unnamed factors. Researchers use the genes in W, differential-expression results, and pathway enrichment to interpret the biology and name the factors. The illustrative four-patient matrix above is not fitted from real data or the earlier two-factor toy V.

7. From NMF to differential expression and pathways

The complete conceptual chain is:

Cancer cells → sum raw counts within each patient → normalize/filter → NMF → dominant archetype from H → determine factor-marker genes (W) and DEGs → interpret gene signatures using pathway enrichment.

NMF itself does not perform DEG analysis or gene set enrichment. Those are downstream analyses.

Identify factor-associated genes

For factor \(k\), rank genes by \(W_{gk}\). A larger value means that gene contributes more strongly to the factor, though gene scaling and preprocessing affect interpretation.

# Inspect strongest genes in each program in the teaching example
for (k in colnames(W)) {
    cat("\n", k, "\n", sep = "")
    print(sort(W[, k], decreasing = TRUE))
}
## 
## Program1
##        MKI67        TOP2A        ISG15        IFIT1 
## 6.665681e-01 3.332841e-01 9.852865e-05 4.926513e-05 
## 
## Program2
##        ISG15        IFIT1        MKI67        TOP2A 
## 6.665678e-01 3.332839e-01 9.882438e-05 4.941043e-05

Differential-expression analysis: the idea

Compare patients assigned to one archetype against patients assigned to other archetypes to identify genes whose expression differs across the groups. With real pseudo-bulk counts and sufficient replicates, a count-aware method such as DESeq2 is often appropriate; modeling choices should reflect cohort design and sample sizes.

The study’s archetype signatures were based on shared genes between NMF factor marker genes and archetype DEGs. The published Methods describe their marker-gene selection and statistical criteria. A toy dataset with four patients is too small to provide a meaningful statistical DEG analysis, so no fabricated significance values are reported here.

Gene set enrichment analysis (GSEA): interpret the signatures

Once you have ranked genes or an archetype signature, compare them with curated pathway gene sets to see which biological processes are over-represented. Yan et al. examined MSigDB hallmark gene sets (Fig. 2e).

Archetype	Interpretation	Reported enriched processes
ARC1	Luminal secretory-like	Apoptosis, TNFα/NFκB signaling, TP53 signaling, EMT
ARC2	Basal-like	EMT, myogenesis
ARC3	Interferon-responsive	Interferon-α and interferon-γ responses, inflammatory response
ARC4	Androgen receptor-enriched	Androgen response, early oestrogen response

Important: Differential expression and enrichment do not establish that a pathway is causally driving a phenotype; they characterize expression associations.

Part III — Single-cell NMF: discover heterogeneity within each patient

8. Why run NMF a second time?

Patient-level pseudo-bulk NMF collapses all cancer cells from one patient into one observation. That is useful for comparing patients, but it hides differences between cancer cells inside a tumour.

The second application retains individual cancer cells. Crucially, the authors performed NMF separately for each patient/sample, not one global cell-level NMF across all patients. Each per-sample NMF creates gene programs describing heterogeneity within that sample. They then compared program marker-gene sets across samples to identify recurrent programs and merged related programs into 13 metaprograms (Fig. 3 and Methods, “Identifying metaprograms of cancer cells”).

Question	Patient-level NMF (Part II)	Single-cell NMF (Part III)
What is one observation (one column of V)?	One patient’s aggregated cancer-cell profile	One individual cancer cell within one sample
Matrix dimensions	Genes × patients	Genes × cancer cells (separate matrix for each sample)
What do W columns mean?	Patient-level archetype-associated gene programs	Within-sample cancer-cell transcriptional programs
What do H columns mean?	Program loadings for each patient	Program loadings for each individual cancer cell
How are factors combined?	Choose a dominant archetype per patient	Compare gene marker sets across per-sample runs; merge recurrent programs
Paper’s result	4 archetypes, ARC1–ARC4	13 recurrent metaprograms, M1–M13
Main biological question	Why do patients differ?	How do cancer cells within a tumour differ, and which programs recur across tumours?

Visual intuition: patient pseudo-bulk converts 100 cancer cells from Patient A into one column; cell-level NMF keeps 100 columns for Patient A and decomposes their heterogeneity.

9. A single-cell example with three tumours

We simulate 3 patients, each having 24 cancer cells. Cells can show mixtures of three programs: cell cycle (MKI67/TOP2A), interferon (ISG15/IFIT1), and hypoxia (VEGFA/NDRG1). These are teaching data, not sequencing data from the paper. For simplicity, we create a small positive expression matrix directly; real scRNA-seq preprocessing is more complex.

set.seed(2026)
genes_sc <- c("MKI67", "TOP2A", "ISG15", "IFIT1", "VEGFA", "NDRG1")
# W_true has one highly specific 2-gene signature per program
W_true <- matrix(
    c(
        2,
        1,
        0,
        0,
        0,
        0,
        0,
        0,
        2,
        1,
        0,
        0,
        0,
        0,
        0,
        0,
        2,
        1
    ),
    nrow = 6,
    dimnames = list(genes_sc, c("Cycle", "IFN", "Hypoxia"))
)

make_sample <- function(patient, dominant) {
    n <- 24
    H_true <- matrix(rexp(3 * n, rate = 2), nrow = 3)
    H_true[dominant, ] <- H_true[dominant, ] +
        rep(c(2, 0.2, 0.2), length.out = n)
    X <- W_true %*% H_true + matrix(runif(6 * n, 0, 0.12), nrow = 6)
    rownames(X) <- genes_sc
    colnames(X) <- paste0(patient, "_cell", seq_len(n))
    X
}

sc_samples <- list(
    P1 = make_sample("P1", 1),
    P2 = make_sample("P2", 2),
    P3 = make_sample("P3", 3)
)

# Observe the difference in granularity
sapply(sc_samples, dim) # each: 6 genes x 24 cells
##      P1 P2 P3
## [1,]  6  6  6
## [2,] 24 24 24

par(mfrow = c(1, 3), mar = c(5, 5, 3, 1))
for (p in names(sc_samples)) {
    image(
        t(sc_samples[[p]][nrow(sc_samples[[p]]):1, ]),
        axes = FALSE,
        main = p,
        xlab = "Cells",
        ylab = "Genes",
        col = hcl.colors(30, "YlOrRd", rev = TRUE)
    )
    axis(
        2,
        at = seq(0, 1, length.out = length(genes_sc)),
        labels = rev(genes_sc),
        las = 2,
        cex.axis = 0.7
    )
}

par(mfrow = c(1, 1))

Each column in the three heatmaps is a cell. A patient can contain cells strongly expressing different programs; patient-level averaging loses this distributional information.

10. Run NMF separately on each patient’s cells

We will reuse the multiplicative-update principle from Part I. This simplified implementation fixes rank \(K=3\) in every sample to make results easy to inspect. The paper actually tested K = 2, 4, …, 30 within each sample, so the toy analysis is not an exact reproduction.

nmf_toy <- function(V, K = 3, n_iter = 1500, seed = 42) {
    set.seed(seed)
    stopifnot(all(V >= 0), K <= min(dim(V)))
    eps <- 1e-10
    W <- matrix(runif(nrow(V) * K, 0.2, 1), nrow(V), K)
    H <- matrix(runif(K * ncol(V), 0.2, 1), K, ncol(V))
    for (it in seq_len(n_iter)) {
        H <- H * ((t(W) %*% V) / (t(W) %*% W %*% H + eps))
        W <- W * ((V %*% t(H)) / (W %*% H %*% t(H) + eps))
    }
    # Resolve arbitrary factor scaling (without changing W %*% H)
    scales <- colSums(W)
    W <- sweep(W, 2, scales, "/")
    H <- sweep(H, 1, scales, "*")
    rownames(W) <- rownames(V)
    colnames(W) <- paste0("Factor", seq_len(K))
    rownames(H) <- colnames(W)
    colnames(H) <- colnames(V)
    list(W = W, H = H, error = sum((V - W %*% H)^2))
}

fits_sc <- lapply(sc_samples, nmf_toy, K = 3)
sapply(fits_sc, function(fit) fit$error)
##         P1         P2         P3 
## 0.07738055 0.06744504 0.10160884
round(fits_sc$P1$W, 2) # genes x within-P1 programs
##       Factor1 Factor2 Factor3
## MKI67    0.65    0.03    0.13
## TOP2A    0.33    0.02    0.07
## ISG15    0.01    0.61    0.01
## IFIT1    0.01    0.31    0.01
## VEGFA    0.00    0.02    0.52
## NDRG1    0.00    0.02    0.26
round(fits_sc$P1$H[, 1:6], 2) # programs x first six P1 cells
##         P1_cell1 P1_cell2 P1_cell3 P1_cell4 P1_cell5 P1_cell6
## Factor1     6.24     0.44     0.00     6.30     0.00     5.68
## Factor2     0.00     0.12     0.09     2.27    10.68     1.72
## Factor3     3.05     7.83    10.39     2.18     1.59     1.98

Interpretation: In cell-level analysis, H[, "P1_cell1"] gives the three program loadings for one cell, not for the entire patient. The factor columns in W are unlabeled, and their order can change between samples; Factor1 in P1 is not necessarily Factor1 in P2.

Label the cells by their highest loading (for visualization only)

P1_fit <- fits_sc$P1
cell_factor <- apply(P1_fit$H, 2, which.max)
table(cell_factor)
## cell_factor
##  1  2  3 
## 12  6  6
barplot(
    P1_fit$H[, 1:12, drop = FALSE],
    col = c("#4D88CC", "#E59655", "#53AA87"),
    border = NA,
    las = 2,
    cex.names = 0.6,
    ylab = "NMF loading",
    main = "P1: mixed programs in individual cells"
)

Important: Choosing a dominant factor here is a didactic visualization. The paper derived recurrent metaprograms from factor marker-gene overlap across samples; it did not obtain 13 metaprograms simply by pooling these per-cell dominant labels.

11. How do per-sample NMF factors become recurrent metaprograms?

Imagine P1’s Factor2 and P2’s Factor1 are both driven by ISG15 and IFIT1. Even though factor numbers differ, their marker genes overlap. The authors use Jaccard similarity to compare marker-gene sets:

\[ J(A,B)=\frac{|A\cap B|}{|A\cup B|}. \]

A Jaccard value of 1 means identical gene sets; 0 means no overlap.

Small hands-on example: detect recurrence across tumours

For this teaching demo, define marker genes as the two highest-weight genes per factor. The paper instead used a special W-based marker selection rule with an allowance for limited overlap; top-two is only a simplified substitute.

markers <- list()
for (p in names(fits_sc)) {
    for (k in seq_len(ncol(fits_sc[[p]]$W))) {
        nm <- paste0(p, "_F", k)
        markers[[nm]] <- names(sort(fits_sc[[p]]$W[, k], decreasing = TRUE))[
            1:2
        ]
    }
}
markers
## $P1_F1
## [1] "MKI67" "TOP2A"
## 
## $P1_F2
## [1] "ISG15" "IFIT1"
## 
## $P1_F3
## [1] "VEGFA" "NDRG1"
## 
## $P2_F1
## [1] "VEGFA" "NDRG1"
## 
## $P2_F2
## [1] "MKI67" "ISG15"
## 
## $P2_F3
## [1] "ISG15" "IFIT1"
## 
## $P3_F1
## [1] "MKI67" "TOP2A"
## 
## $P3_F2
## [1] "ISG15" "IFIT1"
## 
## $P3_F3
## [1] "VEGFA" "NDRG1"

jaccard <- function(a, b) {
    length(intersect(a, b)) / length(union(a, b))
}

jaccard_mat <- outer(
    names(markers),
    names(markers),
    Vectorize(function(a, b) jaccard(markers[[a]], markers[[b]]))
)
dimnames(jaccard_mat) <- list(names(markers), names(markers))
round(jaccard_mat, 2)
##       P1_F1 P1_F2 P1_F3 P2_F1 P2_F2 P2_F3 P3_F1 P3_F2 P3_F3
## P1_F1  1.00  0.00     0     0  0.33  0.00  1.00  0.00     0
## P1_F2  0.00  1.00     0     0  0.33  1.00  0.00  1.00     0
## P1_F3  0.00  0.00     1     1  0.00  0.00  0.00  0.00     1
## P2_F1  0.00  0.00     1     1  0.00  0.00  0.00  0.00     1
## P2_F2  0.33  0.33     0     0  1.00  0.33  0.33  0.33     0
## P2_F3  0.00  1.00     0     0  0.33  1.00  0.00  1.00     0
## P3_F1  1.00  0.00     0     0  0.33  0.00  1.00  0.00     0
## P3_F2  0.00  1.00     0     0  0.33  1.00  0.00  1.00     0
## P3_F3  0.00  0.00     1     1  0.00  0.00  0.00  0.00     1

# Hierarchical clustering based on marker-gene overlap
hc <- hclust(as.dist(1 - jaccard_mat), method = "ward.D2")
plot(
    hc,
    main = "Recurrent expression programs across samples",
    xlab = "Per-sample NMF factors",
    sub = "Teaching example"
)


# Simple illustrative grouping: 3 clusters
program_group <- cutree(hc, k = 3)
print(program_group)
## P1_F1 P1_F2 P1_F3 P2_F1 P2_F2 P2_F3 P3_F1 P3_F2 P3_F3 
##     1     2     3     3     1     2     1     2     3
split(names(program_group), program_group)
## $`1`
## [1] "P1_F1" "P2_F2" "P3_F1"
## 
## $`2`
## [1] "P1_F2" "P2_F3" "P3_F2"
## 
## $`3`
## [1] "P1_F3" "P2_F1" "P3_F3"

The clusters of per-sample factors represent programs that recur in different tumours. This is the conceptual bridge from NMF to metaprograms.

The paper’s actual program-consolidation procedure

The published Methods specify the following process (this is more sophisticated than the small example):

Exclude samples with fewer than 20 cancer cells; exclude genes whose average normalized log2 expression is no greater than 0.05.
For each eligible sample, run NMF separately for K = 2, 4, …, 30 (and do not use K greater than the available number of cells).
For each factor, extract marker genes from its W contribution column using the paper’s factor-specific marker selection procedure.
Compute pairwise Jaccard similarity between factors based on marker sets. Exclude a factor as noisy if its maximum similarity to any other factor is less than 0.25.
Hierarchically cluster program similarities using hclust(..., method="ward.D2"); investigate cut counts C = 10, 12, …, 30.
Iteratively filter low-coherence groups (within-group Jaccard < 0.1) and merge compatible groups (mutual between-group Jaccard > 0.2 when also grouped together at a smaller cut), repeating until stable.
Average the gene-to-factor W contributions of programs assigned to each metaprogram, extract final marker genes, and remove technical-noise metaprograms.

The paper obtained 13 recurrent metaprograms (Fig. 3):

Metaprogram	Biological annotation	Example genes
M1	G2/M	TOP2A, MKI67, CENPF
M2	Mitochondrial	MT-ATP6, MT-CO3, MT-ND3
M3	Ribosomal	RPL26, RPS3A, RPS3
M4	Stress	FOS, JUN, FOSB
M5	Interferon	ISG15, IFI6, IFIT1
M6	HLA	HLA-A, HLA-B, CD74
M7	S phase	PCNA, ATAD2, MCM7
M8	Hypoxia	NDRG1, VEGFA, EGLN3
M9	Basal	KRT17, TAGLN, ACTA2
M10	EMT	KRT7, ANXA1, LGALS3
M11	Luminal secretory	SLPI, LTF, MMP7
M12	Cholesterol	SCD, DHCR7, HMGCR
M13	ER stress	HSPA5, DNAJB11, SDF2L1

These metaprograms are recurring cancer-cell expression states, not 13 mutually exclusive patient classes.

12. Why both analyses are needed

Consider two hypothetical tumours with identical mean expression of proliferative and interferon programs:

Tumour A: half of its cells are predominantly cycling and half are predominantly interferon-responsive.
Tumour B: every cell moderately expresses both programs.

Their patient pseudo-bulk profiles can be similar, so patient-level NMF may assign the same archetype. But the cell-level distribution is different. Single-cell NMF can expose these contrasting states, which may relate differently to treatment response or spatial niches.

# Same gene-wise average, different cell-to-cell heterogeneity
A <- cbind(c(4, 0), c(4, 0), c(0, 4), c(0, 4))
B <- matrix(rep(c(2, 2), 4), nrow = 2)
rownames(A) <- rownames(B) <- c("Cycle", "IFN")
colnames(A) <- colnames(B) <- paste0("cell", 1:4)
rowMeans(A) # 2, 2
## Cycle   IFN 
##     2     2
rowMeans(B) # 2, 2
## Cycle   IFN 
##     2     2

par(mfrow = c(1, 2), mar = c(5, 4, 3, 1))
barplot(
    A,
    col = c("#5689cb", "#e5a14e"),
    main = "Tumour A: distinct cells",
    ylab = "Program signal"
)
barplot(
    B,
    col = c("#5689cb", "#e5a14e"),
    main = "Tumour B: mixed cells",
    ylab = "Program signal"
)

par(mfrow = c(1, 1))

Relationship to the eight ecotypes

The final eight ecotypes were obtained in a third, distinct analysis: the authors measured the co-occurrence of cancer-cell metaprograms with immune and stromal cell states across patients, constructed a network, and partitioned it into eight multicellular communities (Fig. 5). Therefore:

4 archetypes = inter-patient cancer-cell profiles (patient-level NMF); 13 metaprograms = recurrent within-tumour programs (sample-wise cell-level NMF + cross-sample consolidation); 8 ecotypes = co-occurring cancer/TME cellular communities (network analysis).

13. Test yourself

If a patient’s H column is [0.1, 0.6, 0.2, 0.1], which archetype would it receive? ARC2.
Can the same patient express programs ARC1 and ARC3 even if assigned to ARC2? Yes. Its H vector can contain several non-zero contributions.
If MKI67 has a high weight in column 1 of W, what does that imply? It strongly contributes to Program 1, not necessarily that every patient has high MKI67 expression.
Why sum single-cell counts to patient pseudo-bulk? To obtain a patient-level profile, rather than having patients with many sequenced cells dominate a cell-level analysis.
Why center and truncate before NMF in this paper? To represent above-average normalized gene-expression signals with non-negative entries, as required by its NMF formulation.
Does the cell-level NMF use one giant genes × all cells matrix? No. The paper runs NMF per sample, then matches recurrent programs by marker-gene similarity.
Are M1–M13 the same as ARC1–ARC4? No. They represent a different biological and computational level.

14. Summary

The same matrix equation applies at both levels; the meaning of an observation changes:

\[ V_{\text{genes} \times \text{patients}} \approx W_{\text{genes} \times \text{programs}} H_{\text{programs} \times \text{patients}}. \]

At patient level, each V column is one patient; H assigns patient archetypes. At single-cell level, each V column is one cancer cell inside one sample; H describes per-cell program mixtures. Cross-patient comparisons of the W-derived marker genes are then used to consolidate recurrent factors into metaprograms. Biological names are assigned from the genes and subsequent annotation.

Reference

Yan Y, Lin Y, Kumar T, et al. (2026). Ecotypes of triple-negative breast cancer in response to chemotherapy. Nature. https://doi.org/10.1038/s41586-026-10469-9.