Gene Program Discovery with scConvexNMF
2026-10-09
Overview
scConvexNMF is an R package for identifying gene programs (GPs) in single-cell RNA-seq data using a convex non-negative matrix factorization framework. The model decomposes a raw count matrix into:
- W — GP activity scores per cell
- G — gene loadings per GP
- H — contribution of each GP to reconstructing each gene’s expression
A key feature is the ability to control for group-level effects (e.g., clonal subgroups, samples) during factorization, isolating shared transcriptional programs from confounders. A separate rescoring step then lets you quantify GP activity without removing that variation, which is useful for cross-sample comparisons.
This vignette demonstrates the full workflow on a subset of malignant cells from a renal cell carcinoma (RCC) dataset.
Installation
# CRAN dependencies
install.packages(c("Rcpp", "Matrix", "quadprog"))
# Bioconductor dependencies
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(c("SingleCellExperiment", "ComplexHeatmap"))
# scConvexNMF (from GitHub)
# install.packages("remotes")
remotes::install_github("csglab/scConvexNMF")
# Optional: for pathway annotation
install.packages("msigdbr")1. Load Libraries
library(Rcpp)
library(Matrix)
library(quadprog)
library(ComplexHeatmap)
library(circlize)
library(ggplot2)
library(SingleCellExperiment)
library(scConvexNMF)
set.seed(43)2. Load Data
This tutorial will use RCC data from our study. The data is accessible in GEO accession number GSE328575.
sce <- readRDS("GSE328575_sce.rds")
sce## class: SingleCellExperiment
## dim: 30462 85563
## metadata(0):
## assays(1): counts
## rownames(30462): A1BG A1BG-AS1 ... TRAV11 ZNF729
## rowData names(8): ID Symbol ... percentage_gene_gc_content
## unique_gene_Symbol
## colnames(85563): PID003_TP1_T1_1_1...AAACCCACACTTGTCC-1
## PID003_TP1_T1_1_1...AAACCCAGTCTTAGTG-1 ...
## PID022_TM1_P1_1_1...TTTGTTGCAGACGGAT-1
## PID022_TM1_P1_1_1...TTTGTTGGTATGTCCA-1
## colData names(33): Barcode_Sample Sample ... GP_malignant GP_mac
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):3. Prepare Input Matrix
We subset to malignant cells from three samples and filter lowly expressed genes to keep the example tractable.
meta <- colData(sce)
# Keep only malignant cells from three selected samples
meta <- meta[meta$cell_type == "Malignant", ]
meta <- meta[meta$Sample %in% c("PID019_TP1_T1_1_1",
"PID024_TP1_T1_1_1",
"PID023_TP1_T1_1_1"), ]
# Remove cells without an assigned clone (required for group-level modeling)
meta <- meta[!is.na(meta$numbat_clone), ]
cat("Cells per sample:\n")## Cells per sample:print(table(meta$Sample))##
## PID019_TP1_T1_1_1 PID023_TP1_T1_1_1 PID024_TP1_T1_1_1
## 349 224 294# Extract raw count matrix and filter genes expressed in > 5 cells (count > 3)
rc <- assay(sce, "counts")[, rownames(meta)]
genes_use <- rowSums(rc > 3) > 5
cat("\nGenes retained:", sum(genes_use), "of", nrow(rc), "\n")##
## Genes retained: 6080 of 30462rc <- rc[genes_use, ]4. Select the Number of Gene Programs (K)
Rationale
The model fitted values combine:
- A WH decomposition — the gene program structure
- Cell- and gene-specific intercepts — a baseline that captures per-cell library size and per-gene mean expression
We call the prediction based on intercepts alone the offset. To choose K, we evaluate the added explanatory power of the GP structure by comparing the Poisson log-likelihood of the full model to the offset-only model on held-out (masked) entries.
vec_K <- seq(4, 16, 2)
cat("Clone distribution per sample:\n")## Clone distribution per sample:print(table(meta$Sample, meta$numbat_clone))##
## PID019_TP1_T1_1_1...3 PID019_TP1_T1_1_1...4
## PID019_TP1_T1_1_1 257 92
## PID023_TP1_T1_1_1 0 0
## PID024_TP1_T1_1_1 0 0
##
## PID023_TP1_T1_1_1...2 PID024_TP1_T1_1_1...2
## PID019_TP1_T1_1_1 0 0
## PID023_TP1_T1_1_1 224 0
## PID024_TP1_T1_1_1 0 294list_models <- list()
list_fit_stats <- list()
for (K in vec_K) {
cat("Fitting model for K =", K, "\n")
model <- new("scConvexNMF")
model$setup(
M = rc, # Raw count matrix (genes × cells)
K = K, # Number of gene programs
Groups = meta$numbat_clone, # Group-level effects to regress out
maskingRatio = 0.01 # Fraction of entries held out for CV
)
model$initialize_LVs()
model$optimize(20) # Number of optimization iterations
masked_stats <- getMaskedStats(model)
list_models[[as.character(K)]] <- model
list_fit_stats[[as.character(K)]] <- data.frame(
K = K,
corr_full_model = masked_stats$corr_full_model,
corr_only_offset = masked_stats$corr_only_offset,
meanLL_full_model = masked_stats$meanLL_full_model,
meanLL_only_offset = masked_stats$meanLL_only_offset
)
}## Fitting model for K = 4## 0.8180118## Fitting model for K = 6## 0.8180118## Fitting model for K = 8## 0.8180118## Fitting model for K = 10## 0.8180118## Fitting model for K = 12## 0.8180118## Fitting model for K = 14## 0.8180118## Fitting model for K = 16## 0.8180118Cross-Validation Results
Plot the gain in mean Poisson log-likelihood (full model vs. offset-only) across K values. Look for the point at which added K yields diminishing returns.
df_res <- do.call(rbind, list_fit_stats)
ggplot(df_res, aes(x = K, y = meanLL_full_model - meanLL_only_offset)) +
geom_point(size = 2) +
geom_line() +
labs(
x = "Number of gene programs (K)",
y = "ΔMean Poisson log-likelihood\n(full model − offset only)",
title = "Cross-validation: gain from GP structure"
) +
theme_classic()
5. Inspect the Selected Model
Based on the cross-validation plot, we select K = 14.
K <- 14
model <- list_models[[as.character(K)]]Extract W, G, and H
WGH <- getWGH_scCNMF(model)
W <- WGH$W[rownames(meta), ] # Cell × GP activity matrix
G <- WGH$G # Gene × GP loading matrix- W (
cells × K): how active each GP is in each cell - G (
genes × K): how much each gene contributes to each GP - H (
K × genes): how much each GP contributes to reconstructing each gene’s expression
GP–GP Correlation
Correlations between GPs reveal co-activation patterns.
draw(
Heatmap(cor(W),
name = "Pearson r",
column_title = "Gene program correlation",
)
)
6. Pathway Enrichment Analysis
Load Hallmark Gene Sets
library(msigdbr)
geneset <- data.frame(
msigdbr(species = "human", collection = "H")
)[, c("gene_symbol", "gs_name")]
# Binary gene × pathway membership matrix
path <- (table(geneset$gene_symbol, geneset$gs_name) > 0) * 1LWeighted Jaccard Overlap
We compute a weighted Jaccard similarity between each GP’s gene weights (G) and each Hallmark pathway, then visualize the result as a heatmap.
wJ_G <- weightedJaccard_G(model, path)
wJ_G_plot <- visualize_jaccard(wJ_G)
draw(wJ_G_plot$heatmap)
GPs 6, 12, and 13 show the strongest association with the Hallmark Hypoxia pathway.
7. Visualize Hypoxia-Associated Gene Programs
We focus on GPs 6, 12, and 13 and visualize their top contributing genes alongside cell metadata (sample and tumor site).
vec_lvs <- c("LV6", "LV12", "LV13")
meta_heat <- meta[, c("tumor_site", "Sample")]
ht <- visualize_Yp(
model,
vec_lvs = vec_lvs,
n_genes = 10,
meta_heat = meta_heat,
wj = wJ_G,
path = path
)## 'magick' package is suggested to install to give better rasterization.
##
## Set `ht_opt$message = FALSE` to turn off this message.draw(ht)
8. Rescore GP Activities Without Group-Level Correction
The original model removes clone-level variation when estimating GP activities. Rescoring reuses the fitted gene loadings (G) but re-estimates W without regressing out group effects — making activities directly comparable across samples.
rescoringModel <- new("scConvexNMF")
rescoringModel$setup(M = rc, reuseG = model)## 0.8180118rescoringModel$initialize_LVs()
rescoringModel$optimize(20)
WGH_rescored <- getWGH_scCNMF(rescoringModel)
W_rescored <- WGH_rescored$W[rownames(meta), ]Compare Original vs. Rescored Activities for GP12
# Align row order
W <- WGH$W[rownames(meta), ]
W_rescored <- W_rescored[rownames(meta), ]
ggdat <- data.frame(
GP12_original = W[, "LV12"],
GP12_rescored = W_rescored[, "LV12"],
meta
)
p1 <- ggplot(ggdat, aes(x = Sample, y = GP12_original, fill = Sample)) +
geom_violin(trim = FALSE) +
geom_boxplot(width = 0.2, color = "grey40", outlier.shape = NA) +
labs(title = "GP12 — Clone-corrected (original)", y = "GP12 activity") +
theme_classic() +
theme(legend.position = "none", axis.text.x = element_text(angle = 30, hjust = 1))
p2 <- ggplot(ggdat, aes(x = Sample, y = GP12_rescored, fill = Sample)) +
geom_violin(trim = FALSE) +
geom_boxplot(width = 0.2, color = "grey40", outlier.shape = NA) +
labs(title = "GP12 — Rescored (no correction)", y = "GP12 activity") +
theme_classic() +
theme(legend.position = "none", axis.text.x = element_text(angle = 30, hjust = 1))
gridExtra::grid.arrange(p1, p2, ncol = 2)
Sample PID019 shows higher GP12 activit.
9. Save Results
saveRDS(model, file = "model.rds")
saveRDS(rescoringModel, file = "rescoringModel.rds")Session Info
sessionInfo()## R version 4.4.3 (2025-02-28)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 20.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_CA.UTF-8 LC_COLLATE=en_CA.UTF-8
## [5] LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8
## [7] LC_PAPER=en_CA.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/Toronto
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] msigdbr_26.1.0 scConvexNMF_0.0.0.9000
## [3] SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
## [5] Biobase_2.66.0 GenomicRanges_1.58.0
## [7] GenomeInfoDb_1.42.1 IRanges_2.40.0
## [9] S4Vectors_0.44.0 BiocGenerics_0.52.0
## [11] MatrixGenerics_1.18.1 matrixStats_1.4.1
## [13] ggplot2_3.5.1 circlize_0.4.17
## [15] ComplexHeatmap_2.22.0 quadprog_1.5-8
## [17] Matrix_1.7-2 Rcpp_1.0.13
##
## loaded via a namespace (and not attached):
## [1] babelgene_22.9 gtable_0.3.6 shape_1.4.6.1
## [4] rjson_0.2.23 xfun_0.48 bslib_0.8.0
## [7] GlobalOptions_0.1.3 lattice_0.22-5 Cairo_1.6-5
## [10] vctrs_0.6.5 tools_4.4.3 generics_0.1.3
## [13] curl_5.2.3 parallel_4.4.3 tibble_3.2.1
## [16] fansi_1.0.6 highr_0.11 cluster_2.1.8
## [19] pkgconfig_2.0.3 RColorBrewer_1.1-3 assertthat_0.2.1
## [22] lifecycle_1.0.4 GenomeInfoDbData_1.2.13 farver_2.1.2
## [25] compiler_4.4.3 munsell_0.5.1 codetools_0.2-19
## [28] clue_0.3-68 htmltools_0.5.8.1 sass_0.4.9
## [31] yaml_2.3.10 pillar_1.9.0 crayon_1.5.3
## [34] jquerylib_0.1.4 cachem_1.1.0 DelayedArray_0.32.0
## [37] iterators_1.0.14 abind_1.4-8 foreach_1.5.2
## [40] tidyselect_1.2.1 digest_0.6.37 dplyr_1.1.4
## [43] labeling_0.4.3 fastmap_1.2.0 SparseArray_1.6.2
## [46] colorspace_2.1-1 cli_3.6.3 magrittr_2.0.3
## [49] S4Arrays_1.6.0 utf8_1.2.4 withr_3.0.1
## [52] scales_1.3.0 UCSC.utils_1.2.0 rmarkdown_2.28
## [55] XVector_0.46.0 httr_1.4.7 gridExtra_2.3
## [58] png_0.1-8 GetoptLong_1.1.0 evaluate_1.0.1
## [61] knitr_1.48 doParallel_1.0.17 rlang_1.1.4
## [64] glue_1.8.0 jsonlite_1.8.9 R6_2.5.1
## [67] zlibbioc_1.52.0