Pathway activity functions
2026-09-24
Single-sample pathway activity
A common question in transcriptomics is whether the genes belonging to a given pathway are concordantly up- or down-regulated in each sample. Several algorithms tackle this by collapsing a gene-by-sample expression matrix into a pathway-by-sample score matrix, the most widely used being GSVA and ssGSEA. biverse provides fast versions of both with some tricks under the hood to make them faster than the original implementations. Moreover, if you have multiple contrasts (different drugs vs. DMSO for example). bixverse also provides Rust-accelerated implementations from mitch.
Data
Following the GSVA vignette convention, we generate some synthetic expression data (Gaussian-distributed to mimic normalised microarray or RNA-seq log counts, and Poisson-distributed to mimic raw counts) together with a collection of random gene sets.
p <- 10000L # genes
n <- 100L # samples
# Gaussian expression matrix
X <- matrix(
rnorm(p * n),
nrow = p,
dimnames = list(paste0("g", seq_len(p)), paste0("s", seq_len(n)))
)
# Poisson count matrix
X_counts <- matrix(
rpois(p * n, lambda = 10),
nrow = p,
dimnames = list(paste0("g", seq_len(p)), paste0("s", seq_len(n)))
)
storage.mode(X_counts) <- "numeric"
# Random gene sets of varying size
gs <- as.list(sample(10:100, size = 250, replace = TRUE))
gs <- lapply(
gs,
function(n, p) paste0("g", sample(seq_len(p), size = n, replace = FALSE)),
p
)
names(gs) <- paste0("gs", seq_along(gs))GSVA
Gaussian kernel
The Gaussian kernel version of GSVA is appropriate for continuous expression values (log-CPM, microarray intensities, etc.). Running it in bixverse is a single call:
bixverse_res_gaussian <- calc_gsva(
exp = X,
pathways = gs,
kernel = "gaussian"
)
bixverse_res_gaussian[1:5, 1:5]
#> s1 s2 s3 s4 s5
#> gs1 -0.20784044 0.10323427 -0.006921166 -0.13209352 -0.084867420
#> gs2 -0.13148017 -0.33316957 0.335542436 -0.06022571 -0.066540827
#> gs3 0.04680407 0.05807882 -0.084391736 0.15393371 0.076860151
#> gs4 -0.09646102 0.15451531 -0.342093457 -0.19453201 0.014982852
#> gs5 -0.06156146 0.00818317 0.092330995 -0.06826265 -0.006438803If the original GSVA Bioconductor package is available we can verify that the results are essentially identical. Minor differences in numerical precision arise from optimisations on the Rust side, but per-pathway correlations between the two implementations are consistently above 0.99.
library(GSVA)
gsvaPar <- gsvaParam(X, gs)
gsva_res_gaussian <- as.matrix(gsva(gsvaPar, verbose = FALSE))
correlations <- diag(cor(bixverse_res_gaussian, gsva_res_gaussian))
print(sprintf(
"All pathway correlations >= 0.99: %s (min: %.4f)",
all(correlations >= 0.99),
min(correlations)
))
#> [1] "All pathway correlations >= 0.99: TRUE (min: 1.0000)"Let’s compare speed (differences will be more pronounced, the more computational fire power your system has.)
microbenchmark::microbenchmark(
gsva = {
gsvaPar <- gsvaParam(X, gs)
gsva(gsvaPar, verbose = FALSE)
},
bixverse = calc_gsva(exp = X, pathways = gs, kernel = "gaussian"),
times = 3L
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> gsva 1695.0888 1696.2992 1702.0504 1697.5097 1705.5313 1713.5529 3
#> bixverse 357.1036 358.6324 363.3861 360.1612 366.5273 372.8934 3Poisson kernel
For count data a Poisson kernel is more appropriate. The only change is setting kernel = "poisson (to note: prior to "0.4.1", this was controlled via a Boolean flag gaussian = FALSE - this still works, but will throw a deprecation warning now):
bixverse_res_poisson <- calc_gsva(
exp = X_counts,
pathways = gs,
kernel = "poisson"
)
bixverse_res_poisson[1:5, 1:5]
#> s1 s2 s3 s4 s5
#> gs1 0.01640780 -0.15821556 -0.27017533 -0.1042462692 0.008656727
#> gs2 -0.20216175 0.06584906 0.30586274 -0.0794790966 -0.076076044
#> gs3 0.07160004 0.02282096 0.12525568 0.0007774097 -0.020296061
#> gs4 -0.12638670 -0.17015543 -0.10048339 -0.2393481286 0.186163110
#> gs5 0.13422584 0.12526100 0.06837335 0.1597392880 0.086662487And vs. the original
gsvaParPoisson <- gsvaParam(X_counts, gs, kcdf = "Poisson")
gsva_res_poisson <- as.matrix(gsva(gsvaParPoisson, verbose = FALSE))
correlations <- diag(cor(bixverse_res_poisson, gsva_res_poisson))
print(sprintf(
"All pathway correlations >= 0.99: %s (min: %.4f)",
all(correlations >= 0.99),
min(correlations)
))
#> [1] "All pathway correlations >= 0.99: TRUE (min: 1.0000)"And speed:
microbenchmark::microbenchmark(
gsva = {
gsvaPar <- gsvaParam(X_counts, gs, kcdf = "Poisson")
gsva(gsvaPar, verbose = FALSE)
},
bixverse = calc_gsva(exp = X_counts, pathways = gs, kernel = "poisson"),
times = 3L
)
#> Unit: seconds
#> expr min lq mean median uq max neval
#> gsva 12.004886 12.017556 12.023936 12.03023 12.033461 12.036695 3
#> bixverse 1.680057 1.684173 1.691892 1.68829 1.697809 1.707329 3ssGSEA
ssGSEA, first described in Barbie et al., takes a similar approach but does not apply a kernel-based normalisation across samples. bixverse provides a Rust-optimised version here as well:
bixverse_res_ssgsea <- calc_ssgsea(
exp = X,
pathways = gs
)
bixverse_res_ssgsea[1:5, 1:5]
#> s1 s2 s3 s4 s5
#> gs1 -0.01806140 0.1338473 0.09748205 0.10830414 0.06392018
#> gs2 0.06450835 0.0140586 0.28922513 0.05175004 0.10593315
#> gs3 0.12308022 0.1338399 0.08946982 0.16831129 0.13581000
#> gs4 0.13075993 0.1626255 -0.08461512 -0.02482929 0.14748332
#> gs5 0.05195984 0.1300232 0.20406127 0.05773036 0.10783745Let’s compare again against the GSVA version:
ssgseaPar <- ssgseaParam(X, gs)
ssgsea_res <- as.matrix(gsva(ssgseaPar, verbose = FALSE))
correlations <- diag(cor(bixverse_res_ssgsea, ssgsea_res))
print(sprintf(
"All pathway correlations >= 0.99: %s (min: %.4f)",
all(correlations >= 0.99),
min(correlations)
))
#> [1] "All pathway correlations >= 0.99: TRUE (min: 1.0000)"And check the underlying speed differences:
microbenchmark::microbenchmark(
gsva = {
ssgseaPar <- ssgseaParam(X, gs)
gsva(ssgseaPar, verbose = FALSE)
},
bixverse = calc_ssgsea(exp = X, pathways = gs),
times = 3L
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> gsva 391.55190 397.18380 402.95502 402.81571 408.6566 414.49745 3
#> bixverse 62.71413 63.67759 64.30484 64.64104 65.1002 65.55935 3singscore
singscore takes a different angle than GSVA or ssGSEA: rather than estimating a null distribution over samples or genes, it scores each sample independently using only the ranks of the genes in the signature. This makes it well-suited to scoring single samples in isolation and to working with directional signatures (separate up- and down-regulated gene sets). bixverse provides a Rust-accelerated implementation.
Ranking
singscore operates on a ranked expression matrix rather than the raw values, so the first step is to compute per-sample ranks:
ranks <- calc_singscore_rank(exp = X)
ranks[1:5, 1:5]
#> s1 s2 s3 s4 s5
#> g1 7964 8080 1986 4418 8913
#> g2 4762 4102 1696 8781 6687
#> g3 25 1451 9200 3560 7948
#> g4 8783 228 6848 6359 895
#> g5 101 1652 2167 807 3712If your data comes from a small targeted panel (NanoString, RT-qPCR) rather than a transcriptome-wide assay, you can pass a set of stable genes to calc_singscore_rank to use the stable-gene ranking approach from Bhuva et al.:
ranks_stable <- calc_singscore_rank(
exp = X,
stable_genes = c("g1", "g2", "g3", "g4", "g5")
)The stable attribute on the rank matrix is read automatically by the downstream scoring functions to pick the appropriate bounds formula.
Scoring a single signature
For a single signature, calc_singscore returns one score and dispersion per sample. Below we use one of our random gene sets as the up-regulated set and another as the down-regulated set:
up_set <- gs[[1]]
down_set <- gs[[2]]
singscore_res <- calc_singscore(
ranks = ranks,
up_set = up_set,
down_set = down_set
)
head(singscore_res)
#> total_score total_dispersion up_score up_dispersion down_score
#> <num> <num> <num> <num> <num>
#> 1: -0.034921233 2951.120 -0.045986180 2788.775 0.011064946
#> 2: 0.089792025 3032.663 0.026824244 2656.823 0.062967781
#> 3: -0.120671525 3558.987 -0.008493127 3936.309 -0.112178398
#> 4: 0.003851381 3361.059 -0.022214909 3713.919 0.026066289
#> 5: -0.016716246 3430.000 -0.015617679 3172.769 -0.001098567
#> 6: 0.079067602 3029.698 -0.010414632 3407.020 0.089482234
#> down_dispersion sample_id
#> <num> <char>
#> 1: 3113.465 s1
#> 2: 3408.503 s2
#> 3: 3181.664 s3
#> 4: 3008.200 s4
#> 5: 3687.232 s5
#> 6: 2652.375 s6If the direction of the signature is unknown (a gene ontology term, for example), pass only up_set and set known_direction = FALSE. The scoring function then transforms ranks around their median so that genes at either extreme contribute symmetrically.
Permutation testing
To assess whether an observed score is larger than would be expected by chance for a random gene set of the same size, set n_permutations to a positive integer. The function draws random gene sets matching the size of the real signature, scores them, and returns empirical one-tailed p-values alongside the scores. The full null distribution is attached as an attribute for inspection or plotting.
singscore_perm <- calc_singscore(
ranks = ranks,
up_set = up_set,
down_set = down_set,
n_permutations = 1000L,
seed = 42L
)
head(singscore_perm)
#> total_score total_dispersion up_score up_dispersion down_score
#> <num> <num> <num> <num> <num>
#> 1: -0.034921233 2951.120 -0.045986180 2788.775 0.011064946
#> 2: 0.089792025 3032.663 0.026824244 2656.823 0.062967781
#> 3: -0.120671525 3558.987 -0.008493127 3936.309 -0.112178398
#> 4: 0.003851381 3361.059 -0.022214909 3713.919 0.026066289
#> 5: -0.016716246 3430.000 -0.015617679 3172.769 -0.001098567
#> 6: 0.079067602 3029.698 -0.010414632 3407.020 0.089482234
#> down_dispersion sample_id pval
#> <num> <char> <num>
#> 1: 3113.465 s1 1.000
#> 2: 3408.503 s2 0.001
#> 3: 3181.664 s3 1.000
#> 4: 3008.200 s4 1.000
#> 5: 3687.232 s5 1.000
#> 6: 2652.375 s6 0.001
dim(attr(singscore_perm, "null_distribution"))
#> [1] 1000 100The seed argument makes the permutations reproducible.
Scoring many signatures
When you have many signatures to score at once, calc_singscore_multi avoids the per-call overhead and parallelises across gene sets. It accepts a named list of up-regulated sets, optionally paired with a list of down-regulated sets keyed by the same names:
singscore_multi_res <- calc_singscore_multi(
ranks = ranks,
up_pathways = gs[1:10]
)
singscore_multi_res$scores[1:5, 1:5]
#> s1 s2 s3 s4 s5
#> gs1 -0.045986180 0.026824244 -0.008493127 -0.02221491 -0.015617679
#> gs2 -0.011064946 -0.062967781 0.112178398 -0.02606629 0.001098567
#> gs3 0.010112383 0.008994394 -0.024030925 0.02789075 0.016504689
#> gs4 0.002285483 0.044612363 -0.096924838 -0.07019893 0.016974970
#> gs5 -0.027956374 0.018452663 0.036204955 -0.02604292 0.001362800The returned list contains a scores matrix and a matching dispersions matrix, both with shape gene sets × samples.
Comparison with singscore
When the original singscore package is available we can confirm that the two implementations agree closely. One detail worth noting: singscore uses ties.method = "min" when ranking, while bixverse uses average ranks. For continuous data the difference is negligible, but it can introduce small discrepancies on heavily tied data.
library(singscore)
sing_ranks <- rankGenes(X)
sing_res <- simpleScore(
rankData = sing_ranks,
upSet = up_set,
downSet = down_set
)
cor(singscore_res$total_score, sing_res$TotalScore)
#> [1] 1And the speed difference:
microbenchmark::microbenchmark(
singscore = {
sing_ranks <- rankGenes(X)
simpleScore(rankData = sing_ranks, upSet = up_set, downSet = down_set)
},
bixverse = {
ranks <- calc_singscore_rank(exp = X)
calc_singscore(ranks = ranks, up_set = up_set, down_set = down_set)
},
times = 5L
)
#> Unit: milliseconds
#> expr min lq mean median uq max neval
#> singscore 56.86308 57.30997 57.44354 57.47051 57.63791 57.93622 5
#> bixverse 11.93253 12.88710 14.04361 13.58652 13.89422 17.91770 5Multi-contrast enrichment (mitch)
When you have multiple contrasts - say, differential expression results from several comparisons or different omics layers - you often want to know which pathways show coordinated changes across contrasts. The mitch method does exactly this: it ranks genes within each contrast, computes mean ranks per pathway, and then uses a MANOVA test to identify pathways that are enriched in one or more contrasts simultaneously. bixverse re-implements the core algorithm so that it slots into the same workflow as the other pathway activity functions… Also, heavily multi-threaded and designed to go brrrrrr in terms of speed.
calc_mitch expects a numeric matrix of contrast statistics (genes in rows, contrasts in columns) and a named list of gene sets. It returns a data.table containing per-pathway MANOVA p-values, FDR-adjusted p-values, and the individual contrast-level enrichment scores and p-values.
set.seed(42L)
contrast_data <- matrix(rnorm(3 * 26), nrow = 26)
colnames(contrast_data) <- sprintf("contrast_%i", 1:3)
rownames(contrast_data) <- letters
gene_sets <- list(
pathway_A = sample(letters, 4),
pathway_B = sample(letters, 5),
pathway_C = sample(letters, 6),
pathway_D = sample(letters, 7)
)
res <- calc_mitch(
contrast_mat = contrast_data,
gene_set_list = gene_sets
)
resNote that pathways smaller than the internal minimum set size are dropped from the output, and the results are sorted by significance. The function also validates its input: passing a matrix containing NA values will throw an error rather than silently producing nonsense.
Comparison with the mitch package
When the original mitch package is available we can verify that the two implementations produce identical results. Here we use the example data bundled with mitch:
data(myImportedData, genesetsExample, package = "mitch")
mitch_res <- suppressMessages(mitch::mitch_calc(
myImportedData,
genesetsExample,
priority = "significance",
minsetsize = 5,
cores = 2
))
bixverse_res <- calc_mitch(
contrast_mat = as.matrix(myImportedData),
gene_set_list = genesetsExample
)
# Pathway ordering and FDR values should match exactly
all(bixverse_res$pathway_names == mitch_res$enrichment_result$set) &&
all.equal(bixverse_res$manova_fdr, mitch_res$enrichment_result$p.adjustMANOVA)
#> [1] TRUEAnd let’s check the speed differences:
microbenchmark::microbenchmark(
mitch = suppressMessages(mitch::mitch_calc(
myImportedData,
genesetsExample,
priority = "significance",
minsetsize = 5,
cores = 2
)),
bixverse = calc_mitch(
contrast_mat = as.matrix(myImportedData),
gene_set_list = genesetsExample
),
times = 5L
)
#> Unit: milliseconds
#> expr min lq mean median uq max
#> mitch 106.560651 107.470870 110.704287 110.440091 114.384125 114.665696
#> bixverse 1.989689 2.071071 3.086893 2.176949 4.440206 4.756548
#> neval
#> 5
#> 5As with the GSVA and ssGSEA implementations, you should observe meaningful speed improvements from Rust, particularly as the number of gene sets or contrasts grows and the more cores/oomph your system has.