SPAROscore: Transcriptomics Analysis Workflows
# load the necessary libraries
library(SPAROscore)
library(airway)
library(ggplot2)
library(dplyr)
library(edgeR)
library(Seurat)
library(muscData)
library(msigdbr)
library(scater)
library(spatialLIBD)
library(ggspavis)Introduction
Gene signature scoring has become an integral component of transcriptomics analyses, enabling the activity of biological pathways, cellular processes, and molecular phenotypes to be quantified at the level of individual samples, cells, or spatial domains.
Rather than focusing on the expression of individual genes, signature scoring summarizes the coordinated behaviour of set of genes, providing a measure of biological activity that complements conventional differential expression and clustering analyses.
SPAROscore is a rank-based gene signature scoring method designed to provide robust signature scores across a wide range of transcriptomics technologies. By ranking genes within each column and comparing the observed ordering of signature genes against the background gene expression profile, SPAROscore produces scores that are resilient to differences in sequencing depth, expression sparsity, and technical variation.
Consequently, the same scoring framework can be applied to bulk RNA sequencing, single-cell RNA sequencing, and spatial transcriptomics datasets without requiring technology-specific modifications. This vignette demonstrates how SPAROscore can be incorporated into typical transcriptomics analysis workflows across these three experimental modalities.
Example datasets
The following publicly available datasets are used in this vignette:
Bulk RNA-seq: the airway Bioconductor
dataset, containing RNA-seq data from human airway smooth muscle cell
lines treated with dexamethasone alongside untreated controls. We
quantify the activity of a glucocorticoid response signature derived
from the original publication.
Single-cell RNA-seq: the interferon-β stimulation dataset from Kang et al., available through the muscData package. This dataset consists of peripheral blood mononuclear cells (PBMCs) collected before and after IFN-β stimulation, allowing the activation of interferon signalling pathways to be examined at single-cell resolution.
Spatial transcriptomics: the human dorsolateral prefrontal cortex (DLPFC) Visium dataset from the SpatialLIBD package. This dataset combines spatially resolved gene expression with histological tissue images, cell-types presence to be visualized directly within its anatomical context.
Bulk Transcriptomics Analysis
For this analysis, the gene expression matrix from the
airway package is used. The data contains RNA-seq read
counts from eight human airway smooth muscle (ASM) cell line samples:
four treated with dexamethasone (a synthetic glucocorticoid), and four
untreated controls.
Rather than performing differential expression analysis, the objective is to quantify the activity of a glucocorticoid response gene signature within each sample. The set of top glucocorticoid-responsive genes reported in the original study is used as the signature to compute a SPAROscore for each sample.
Load the dataset and the gene signature
The expression matrix contains genes in rows and samples in columns. The accompanying metadata indicates whether each sample received dexamethasone treatment (trt) or served as an untreated control (untrt).
# load the airway dataset
data(airway)
# extract the counts matrix
counts_data <- assay(airway)
# extract the sample metadata
metadata <- as.data.frame(colData(airway))
# inspect the experimental design
metadata[, c("Run", "dex")]
#> Run dex
#> SRR1039508 SRR1039508 untrt
#> SRR1039509 SRR1039509 trt
#> SRR1039512 SRR1039512 untrt
#> SRR1039513 SRR1039513 trt
#> SRR1039516 SRR1039516 untrt
#> SRR1039517 SRR1039517 trt
#> SRR1039520 SRR1039520 untrt
#> SRR1039521 SRR1039521 trt
# set the glucocorticoid response markers as the gene signature
gc_response_genes <- c("ENSG00000112936", "ENSG00000170606", "ENSG00000120129",
"ENSG00000152795", "ENSG00000211445", "ENSG00000163884",
"ENSG00000189221", "ENSG00000101347", "ENSG00000196136",
"ENSG00000152583", "ENSG00000120658", "ENSG00000157514",
"ENSG00000103196", "ENSG00000179094", "ENSG00000152763",
"ENSG00000103310", "ENSG00000127954")
# check if the gene ids in the signature are also in the expression data
mean(gc_response_genes %in% rownames(counts_data))
#> [1] 1It is better to have most signature genes to be present in the expression data. The number of signature missing from the dataset will also be reported by SPAROscore as a warning.
Compute SPAROscore
bulk_scores <- sparoscore( data = counts_data, signatures = gc_response_genes)
#> SPAROscore says: Calculating column-wise geometric averages
#> SPAROscore says: Ranking a matrix object
bulk_scores
#> SPAROscore
#> SRR1039508 0.5238153
#> SRR1039509 0.7654385
#> SRR1039512 0.5495932
#> SRR1039513 0.7473585
#> SRR1039516 0.5337068
#> SRR1039517 0.7162212
#> SRR1039520 0.5222398
#> SRR1039521 0.7593574The resulting score represents the relative enrichment of the glucocorticoid response signature within each sample. Larger values indicate that the signature genes have higher expression relative to the background expression of the sample, consistent with stronger pathway activation.
Combine scores with the sample metadata
metadata$SPAROscore <- bulk_scores[,1]
metadata[, c("Run", "dex", "SPAROscore")]
#> Run dex SPAROscore
#> SRR1039508 SRR1039508 untrt 0.5238153
#> SRR1039509 SRR1039509 trt 0.7654385
#> SRR1039512 SRR1039512 untrt 0.5495932
#> SRR1039513 SRR1039513 trt 0.7473585
#> SRR1039516 SRR1039516 untrt 0.5337068
#> SRR1039517 SRR1039517 trt 0.7162212
#> SRR1039520 SRR1039520 untrt 0.5222398
#> SRR1039521 SRR1039521 trt 0.7593574Appending the scores to the sample metadata makes them readily available for downstream statistical analyses and visualization.
Compare treated and untreated samples
Since dexamethasone activates glucocorticoid signalling, treated samples are expected to exhibit higher glucocorticoid response scores than untreated controls.
metadata %>%
ggplot(aes(x = dex, y = SPAROscore, fill = dex)) +
geom_boxplot(width = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.12, size = 3) +
xlab("Treatment") +
ylab("SPAROscore")+
ggtitle("Glucocorticoid response signature activity")+
scale_fill_manual(values = c("#CC79A7", "#009E73"))+
theme_classic()
As expected, dexamethasone-treated samples display substantially higher signature scores than untreated controls, demonstrating that SPAROscore successfully captures the biological response induced by the treatment.
Statistical comparison
The observed difference can be formally evaluated using the non-parametric Wilcoxon ranksum test
wilcox.test(SPAROscore ~ dex, data = metadata)
#>
#> Wilcoxon rank sum exact test
#>
#> data: SPAROscore by dex
#> W = 16, p-value = 0.02857
#> alternative hypothesis: true location shift is not equal to 0Although the dataset contains only eight samples, the treatment effect is sufficiently strong that the glucocorticoid response signature clearly separates the two experimental groups.
Visualising pathway activity on a PCA
Principal component analysis (PCA) is commonly used to visualise sample-to-sample variation in transcriptomic datasets. By colouring samples according to their SPAROscore rather than their treatment group, we can determine whether pathway activity explains the dominant source of variation.
# Create the DGEList object
bulk_dge <- DGEList(counts = counts_data)
# filter genes that are sufficiently informative
bulk_keep <- filterByExpr(bulk_dge)
bulk_dge <- bulk_dge[bulk_keep, keep.lib.sizes = FALSE]
# Calculate normalization factors
bulk_dge <- calcNormFactors(bulk_dge)
#> calcNormFactors has been renamed to normLibSizes
# Obtain log2 counts per million (logCPM)
bulk_logCPM <- edgeR::cpm(bulk_dge, log = TRUE, prior.count=1)
# perform PCA
bulk_pca_res <- prcomp(t(bulk_logCPM), scale. = TRUE)
# create a dataframe with PCs
bulk_pca_df <- data.frame(PC1 = bulk_pca_res$x[,1],
PC2 = bulk_pca_res$x[,2],
SPAROscore = metadata$SPAROscore,
Treatment = metadata$dex)
# calculate the variance explained by PC1 and PC2
pc1_exp <- round(summary(bulk_pca_res)$importance[2,1]*100,1)
pc2_exp <- round(summary(bulk_pca_res)$importance[2,2]*100,1)
# visualise the PCA
bulk_pca_df %>%
ggplot(aes(x = PC1, y = PC2, color = SPAROscore, shape = Treatment)) +
geom_point(size = 4) +
labs(
title = "Principal component analysis",
x=paste0("PC1 (", pc1_exp, "%)"),
y=paste0("PC2 (", pc2_exp, "%)"))+
theme_classic()
The continuous colour gradient reveals that samples with similar glucocorticoid response scores cluster together in PCA space. This illustrates how pathway-level activity can explain the major sources of biological variation without relying solely on experimental labels.
Relationship between SPAROscore and individual marker genes
Gene signature scores summarize the coordinated behaviour of multiple genes. To illustrate this, we compare the SPAROscore with the expression of most differentially expressed glucocorticoid response gene (C7) from the original study.
# assign the top gene from the original study
top_gene <- "ENSG00000112936"
bulk_expr_df <- data.frame(C7_log_expression = log2(counts_data[top_gene, ]+1),
SPAROscore = metadata$SPAROscore,
Treatment = metadata$dex)
# compute correlation of top gene expression vs signature scores
cor(bulk_expr_df$C7_log_expression, bulk_expr_df$SPAROscore)
#> [1] 0.8088692
# plot top gene expression and SPAROscore along with treatment
bulk_expr_df %>%
tidyr::pivot_longer(cols = c("C7_log_expression", "SPAROscore")) %>%
ggplot(aes(x = Treatment, y = value, fill = Treatment)) +
geom_boxplot(width = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.12, size = 3) +
facet_wrap(~name, scales = "free") +
ggtitle("Glucocorticoid response signature activity")+
scale_fill_manual(values = c("#2297E6", "#E69F00"))+
theme_classic()
It can be seen that the values overlap for the marker gene expression, but not for the signature scores.
Although expression of individual marker genes often correlates with pathway activity, SPAROscore integrates information across the entire gene signature, making it less sensitive to noise or variability affecting any single gene. This provides a more stable and biologically interpretable measure of pathway activation.
Summary
In this example, SPAROscore quantified glucocorticoid pathway activity from bulk RNA-seq data, clearly distinguishing treated from untreated samples. The scores aligned with the main transcriptomic variation and provided an interpretable metric for downstream analyses and visualisation. This illustrates how SPAROscore complements conventional bulk RNA-seq analyses by summarising coordinated biological responses into a single interpretable metric.
Single-cell RNA-seq analysis
Single-cell RNA sequencing (scRNA-seq) enables gene expression to be measured at the resolution of individual cells, revealing cellular heterogeneity that is often masked in bulk RNA-seq experiments. In addition to identifying distinct cell populations, scRNA-seq provides an opportunity to quantify biological pathway activity within each cell, allowing functional differences between cell types and experimental conditions to be explored.
In this example, the IFN-β stimulation dataset from Kang et al., which is
available through the muscData package will be used. The dataset
consists of ~29,000 peripheral blood mononuclear cells (PBMCs) collected
from lupus patients before and after stimulation with interferon-β
(IFN-β). The colData contains the cell type annotations in the
cell column, treatment condition in the stim
column, donor_id in the ind column and multiplet status in
the multiplets column. For the simplicity of this vignette,
cells from a a single sample will be analysed. We will work with a
seurat object of the same data to demonstrate the versatility of
SPAROscore.
To quantify interferon signalling, the Hallmark Interferon Alpha Response gene set from the Molecular Signatures Database (MSigDB) will be used. Although the cells were stimulated with IFN-β, the Hallmark collection contains a single canonical type I interferon response signature, which captures the conserved transcriptional response induced by both IFN-α and IFN-β signalling.
Load the dataset from muscData
# load the ifnb data from the muscData library
ifnb_sce_data <- Kang18_8vs8()
#> see ?muscData and browseVignettes('muscData') for documentation
#> downloading 1 resources
#> retrieving 1 resource
#> loading from cache
# view the colData
head(colData(ifnb_sce_data))
#> DataFrame with 6 rows and 5 columns
#> ind stim cluster cell multiplets
#> <integer> <factor> <integer> <factor> <factor>
#> AAACATACAATGCC-1 107 ctrl 5 CD4 T cells doublet
#> AAACATACATTTCC-1 1016 ctrl 9 CD14+ Monocytes singlet
#> AAACATACCAGAAA-1 1256 ctrl 9 CD14+ Monocytes singlet
#> AAACATACCAGCTA-1 1256 ctrl 9 CD14+ Monocytes doublet
#> AAACATACCATGCA-1 1488 ctrl 3 CD4 T cells singlet
#> AAACATACCTCGCT-1 1256 ctrl 9 CD14+ Monocytes singlet
# check how many genes * cells are present in the original dataset
dim(ifnb_sce_data)
#> [1] 35635 29065
# subset the original data to work with only one sample and singlets
ifnb_sce_data <- ifnb_sce_data[,
(colData(ifnb_sce_data)$ind == "107" &
colData(ifnb_sce_data)$multiplets == "singlet")
]
# check how many genes * cells are present in the filtered dataset
dim(ifnb_sce_data)
#> [1] 35635 1042
# working with a seurat object
ifnb <- CreateSeuratObject(
counts = counts(ifnb_sce_data),
meta.data = as.data.frame(colData(ifnb_sce_data)),
assay = "RNA"
)
# check the layer in the RNA assay
Layers(ifnb[["RNA"]])
#> [1] "counts"
#> [1] "counts" "data"
# view the metadata of the ifnb seurat object
head(ifnb[[]])
#> orig.ident nCount_RNA nFeature_RNA ind stim cluster
#> AAACATACTGCGTA-1 SeuratProject 2747 980 107 ctrl 3
#> AAACGGCTTTCGTT-1 SeuratProject 1271 464 107 ctrl 5
#> AAATCCCTAAGAGT-1 SeuratProject 1333 482 107 ctrl 5
#> AAATTGACCTGAGT-1 SeuratProject 3892 1049 107 ctrl 9
#> AACAATACGGATTC-1 SeuratProject 932 505 107 ctrl 10
#> AACACTCTGGAGGT-1 SeuratProject 2771 810 107 ctrl 9
#> cell multiplets
#> AAACATACTGCGTA-1 CD4 T cells singlet
#> AAACGGCTTTCGTT-1 CD4 T cells singlet
#> AAATCCCTAAGAGT-1 CD4 T cells singlet
#> AAATTGACCTGAGT-1 CD14+ Monocytes singlet
#> AACAATACGGATTC-1 CD4 T cells singlet
#> AACACTCTGGAGGT-1 CD14+ Monocytes singletThe Seurat object contains gene expression of PBMCs under two
experimental conditions: unstimulated cells (CTRL), and celle stimulated
with interferron-beta (STIM). Along with the experimental condition, the
cell types are also given under the cell metadata
column.
Load the gene signature from MSigDB
Compute SPAROscore
The interferon-beta response for every cell can be quantified using SPAROscore
ifnb <- sparoscore(data = ifnb,
assay = "RNA",
layer = "counts",
signatures = ifn_signature,
prefix = "IFN_" )
#> SPAROscore says: Calculating column-wise geometric averages
#> SPAROscore says: Ranking a sparseMatrix object
#> Warning: SPAROscore says: The following 11 signature genes are missing in the
#> input dataset: C1S, CNP, DHX58, HLA-C, PSMB8, PSMB9, TAP1, TENT5A, TRIM26,
#> TXNIP, WARS1
# check the layer in the RNA assay after scoring
Layers(ifnb[["RNA"]])
#> [1] "counts" "ranks"
# view the metadata after scoring
head(ifnb[[]])
#> orig.ident nCount_RNA nFeature_RNA ind stim cluster
#> AAACATACTGCGTA-1 SeuratProject 2747 980 107 ctrl 3
#> AAACGGCTTTCGTT-1 SeuratProject 1271 464 107 ctrl 5
#> AAATCCCTAAGAGT-1 SeuratProject 1333 482 107 ctrl 5
#> AAATTGACCTGAGT-1 SeuratProject 3892 1049 107 ctrl 9
#> AACAATACGGATTC-1 SeuratProject 932 505 107 ctrl 10
#> AACACTCTGGAGGT-1 SeuratProject 2771 810 107 ctrl 9
#> cell multiplets rank_caps IFN_SPAROscore
#> AAACATACTGCGTA-1 CD4 T cells singlet 361 0.01311115
#> AAACGGCTTTCGTT-1 CD4 T cells singlet 148 0.02659397
#> AAATCCCTAAGAGT-1 CD4 T cells singlet 143 0.02419072
#> AAATTGACCTGAGT-1 CD14+ Monocytes singlet 355 0.03814999
#> AACAATACGGATTC-1 CD4 T cells singlet 107 0.03918696
#> AACACTCTGGAGGT-1 CD14+ Monocytes singlet 276 0.04831208After scoring, SPAROscore returns the original Seurat object with the calculated pathway scores added to the cell metadata. In addition, the computed gene ranks are stored as a new layer within the RNA assay, allowing subsequent signatures to be scored without repeating the ranking step.
The newly created IFN_SPAROscore column can now be used in the same way as any other cell-level metadata throughout the Seurat analysis workflow.
Comparing stimulated and control cells
Since interferon signalling is expected to increase following stimulation, the distribution of SPAROscore between the two experimental conditions should reflect that.
VlnPlot(ifnb,
features = "IFN_SPAROscore",
group.by = "stim",
pt.size = 0) +
scale_fill_manual(values = c("#2166AC", "#B2182B"))+
theme(legend.position = "none")
#> Warning: Default search for "data" layer in "RNA" assay yielded no results;
#> utilizing "counts" layer instead.
As expected, stimulated cells exhibit higher interferon response scores than unstimulated controls.
Statistical comparison
The difference between the two experimental groups can also be assessed statistically using the non-parametric Wilcoxon ranksum test.
wilcox.test(IFN_SPAROscore ~ stim, ifnb[[]])
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: IFN_SPAROscore by stim
#> W = 6611.5, p-value < 2.2e-16
#> alternative hypothesis: true location shift is not equal to 0The highly significant difference reflects the robust activation of interferon-responsive genes following stimulation, and the ability of SPAROscore to quantify that biological phenomenon.
Interferon activity across cell types
Because pathway scores are calculated for every individual cell, they can be compared across annotated immune cell populations.
This visualisation highlights differences in interferon signalling across immune cell types. It summarises the pathway activity within each cell population under each experimental condition, making it straightforward to use SPAROscore and identify populations that exhibit the greater transcriptional response following interferon stimulation.
Visualising pathway activity on the UMAP
One of the most useful ways to explore pathway activity is to visualise the scores on the UMAP embeddings.
# get umap embedding from the gene expression data
# normalise the gene expression
ifnb <- NormalizeData(ifnb)
# identify the variable genes
ifnb <- FindVariableFeatures(ifnb, selection.method = "vst", nfeatures = 2000)
# scale the data
ifnb <- ScaleData(ifnb, features = rownames(ifnb))
# perform PCA
ifnb <- RunPCA(ifnb, features = VariableFeatures(object = ifnb))
# Run network-based clustering algorithm
ifnb <- FindNeighbors(ifnb, dims = 1:10)
ifnb <- FindClusters(ifnb, resolution = 0.5)
#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
#>
#> Number of nodes: 1042
#> Number of edges: 34050
#>
#> Running Louvain algorithm...
#> Maximum modularity in 10 random starts: 0.9049
#> Number of communities: 11
#> Elapsed time: 0 seconds
# run umap
ifnb <- RunUMAP(ifnb, dims = 1:10)
# visualise the experimental condition in umap
DimPlot(ifnb, reduction = "umap", group.by = "stim")
# visualise the SPAROscores in umap
FeaturePlot( ifnb, features = "IFN_SPAROscore",
reduction = "umap") 
It can be seen that the cells exhibiting strong activation of the
interferon pathway are readily identified by SPAROscore, while
neighbouring cells with weaker responses remain distinguishable.
Efficiently scoring additional pathways
Since SPAROscore stores the calculated gene ranks within the Seurat object, subsequent pathway analyses can reuse these ranks without repeating the computationally intensive ranking step.
For example, an inflammatory response signature could be scored as follows.
# load the inflammatory response hallmark genes from msigdb
inflammation_signature <- msigdbr(species = "Homo sapiens",collection = "H") %>%
filter(gs_name == "HALLMARK_INFLAMMATORY_RESPONSE") %>%
pull(gene_symbol) %>%
unique()
# compute sparoscores for the new signature with pre-computed ranks
ifnb <- sparoscore(
data = ifnb,
data_has_ranks = TRUE,
assay = "RNA",
layer = "ranks",
signatures = inflammation_signature,
prefix = "Inflammation_"
)
#> Warning: SPAROscore says: The following 12 signature genes are missing in the
#> input dataset: ADGRE1, CCL24, CCL5, CLEC5A, CXCL8, CYBB, GABBR1, IFNGR2, ITGB3,
#> LTA, SELENOS, TAPBP
# view the metadata after scoring for the second signature
head(ifnb[[]])
#> orig.ident nCount_RNA nFeature_RNA ind stim cluster
#> AAACATACTGCGTA-1 SeuratProject 2747 980 107 ctrl 3
#> AAACGGCTTTCGTT-1 SeuratProject 1271 464 107 ctrl 5
#> AAATCCCTAAGAGT-1 SeuratProject 1333 482 107 ctrl 5
#> AAATTGACCTGAGT-1 SeuratProject 3892 1049 107 ctrl 9
#> AACAATACGGATTC-1 SeuratProject 932 505 107 ctrl 10
#> AACACTCTGGAGGT-1 SeuratProject 2771 810 107 ctrl 9
#> cell multiplets rank_caps IFN_SPAROscore
#> AAACATACTGCGTA-1 CD4 T cells singlet 361 0.01311115
#> AAACGGCTTTCGTT-1 CD4 T cells singlet 148 0.02659397
#> AAATCCCTAAGAGT-1 CD4 T cells singlet 143 0.02419072
#> AAATTGACCTGAGT-1 CD14+ Monocytes singlet 355 0.03814999
#> AACAATACGGATTC-1 CD4 T cells singlet 107 0.03918696
#> AACACTCTGGAGGT-1 CD14+ Monocytes singlet 276 0.04831208
#> RNA_snn_res.0.5 seurat_clusters Inflammation_SPAROscore
#> AAACATACTGCGTA-1 8 8 0.016885553
#> AAACGGCTTTCGTT-1 2 2 0.010942041
#> AAATCCCTAAGAGT-1 2 2 0.018591348
#> AAATTGACCTGAGT-1 0 0 0.071711520
#> AACAATACGGATTC-1 2 2 0.008340747
#> AACACTCTGGAGGT-1 0 0 0.073530274This approach is particularly advantageous when evaluating many gene signatures, as the gene ranking only needs to be performed once.
Pseudo-bulk analyses
While single-cell analyses reveal pathway activity at the level of individual cells, it is often informative to summarise expression profiles across biologically similar cells. The pseudo-bulk approach aggregates gene expression within each cell population, reducing cell-to-cell variability while retaining cell type specificity. Pseudo-bulk analyses are commonly used for differential expression and pathway analyses, as they provide more stable estimates of average expression for each biological group.
Here, we aggregate the IFN-β dataset by both the cell type annotation (cell) and stimulation condition (stim). We then compute SPAROscore on the resulting pseudo-bulk expression matrix to quantify interferon pathway activity across each cell population.
Aggregate cells into pseudo-bulk profiles
pseudobulk_ifn <- AggregateExpression(ifnb,
group.by = c("cell", "stim"),
assays = "RNA",
return.seurat = TRUE,
slot = "counts")
#> Centering and scaling data matrix
pseudobulk_ifn[[]]
#> orig.ident cell stim
#> B cells_ctrl B cells_ctrl B cells ctrl
#> B cells_stim B cells_stim B cells stim
#> CD14+ Monocytes_ctrl CD14+ Monocytes_ctrl CD14+ Monocytes ctrl
#> CD14+ Monocytes_stim CD14+ Monocytes_stim CD14+ Monocytes stim
#> CD4 T cells_ctrl CD4 T cells_ctrl CD4 T cells ctrl
#> CD4 T cells_stim CD4 T cells_stim CD4 T cells stim
#> CD8 T cells_ctrl CD8 T cells_ctrl CD8 T cells ctrl
#> CD8 T cells_stim CD8 T cells_stim CD8 T cells stim
#> Dendritic cells_ctrl Dendritic cells_ctrl Dendritic cells ctrl
#> Dendritic cells_stim Dendritic cells_stim Dendritic cells stim
#> FCGR3A+ Monocytes_ctrl FCGR3A+ Monocytes_ctrl FCGR3A+ Monocytes ctrl
#> FCGR3A+ Monocytes_stim FCGR3A+ Monocytes_stim FCGR3A+ Monocytes stim
#> Megakaryocytes_ctrl Megakaryocytes_ctrl Megakaryocytes ctrl
#> Megakaryocytes_stim Megakaryocytes_stim Megakaryocytes stim
#> NK cells_ctrl NK cells_ctrl NK cells ctrl
#> NK cells_stim NK cells_stim NK cells stimEach column in the expression data now corresponds to a unique combination of cell type and treatment (for example, “B_CTRL” or “B_STIM”), while rows represent genes.
Compute scores at pseudo-bulk level
pseudobulk_ifn <- sparoscore(data = pseudobulk_ifn,
signatures = ifn_signature,
assay = "RNA",
layer = "counts")
#> SPAROscore says: Calculating column-wise geometric averages
#> SPAROscore says: Ranking a sparseMatrix object
#> Warning: SPAROscore says: The following 11 signature genes are missing in the
#> input dataset: C1S, CNP, DHX58, HLA-C, PSMB8, PSMB9, TAP1, TENT5A, TRIM26,
#> TXNIP, WARS1Visualise pseudo-bulk level scores
pseudobulk_ifn[[]] %>%
ggplot(aes(y = cell, x = SPAROscore, fill = stim)) +
geom_col(position = position_dodge(width = 0.8)) +
theme_classic() +
scale_fill_manual(values = c("#2166AC", "#B2182B"))+
theme()
Compared with the single-cell visualisations, the pseudo-bulk analysis
emphasises the average interferon response within each immune cell
population. Cell types exhibiting the strongest pathway activation
following IFN-β stimulation can be readily identified, while differences
between stimulated and control samples remain easy to interpret.
Overall, pseudo-bulk analyses with SPAROscore complement cell-level investigations by providing a robust summary of pathway activity within biologically meaningful groups.
Since SPAROscore computes scores in a robust, sparsity-adaptive manner (and operates on any matrix-like expression object), it can be seamlessly applied to bulk RNA-seq data, pseudo-bulk profiles, and individual cells without any modification.
In such cross-modal analyses, SPAROscore computes biologically interprettable scores that are feasible to be compared and analysed together.
Summary
This example demonstrates how SPAROscore integrates naturally into a standard single-cell RNA-seq workflow. The computed pathway scores provide an interpretable measure of interferon signalling for every individual cell, enabling biological responses to be visualised on low-dimensional embeddings, compared across experimental conditions, and summarised within annotated cell populations using pseudo-bulk aggregations.
By reducing the expression of an entire gene signature to a single quantitative metric, SPAROscore complements conventional scRNA-seq gene analyses and facilitates the exploration of functional heterogeneity in single-cell datasets.
Spatial transcriptomics analysis
Spatial transcriptomics extends conventional RNA sequencing by preserving the spatial locations of transcriptional measurements within intact tissue sections. In addition to quantifying gene expression, spatial transcriptomics enables biological processes to be investigated within their anatomical context, allowing molecular activity to be related directly to tissue morphology and organization.
In this example, the Visium human dorsolateral prefrontal cortex (DLPFC) dataset available through the SpatialLIBD package is used. This dataset was generated using the 10x Genomics Visium platform and comprises multiple human brain sections with manually annotated cortical layers.
Rather than scoring a signalling pathway, the oligodendrocyte marker gene signature will be quantified. Oligodendrocyte are abundant in the white matter, but not abundant in the cortical layers. Consequently, oligodendrocyte signature scores are expected to be high within the white matter and substantially lower in the cortex.
Load the DLPFC dataset from SpatialLIBD
For computational efficiency, we analyse one representative tissue section.
# load the DLPFC Visium dataset
dlpfc_spe <- fetch_data(type = "spe")
#> adding rname 'https://www.dropbox.com/s/f4wcvtdq428y73p/Human_DLPFC_Visium_processedData_sce_scran_spatialLIBD.Rdata?dl=1'
#> 2026-07-22 08:07:34.754444 loading file /github/home/.cache/R/BiocFileCache/275d2e8a03b5_Human_DLPFC_Visium_processedData_sce_scran_spatialLIBD.Rdata%3Fdl%3D1
# check some domain-level metadata present in the object
head(colData(dlpfc_spe)[,c("sample_id", "spatialLIBD")])
#> DataFrame with 6 rows and 2 columns
#> sample_id spatialLIBD
#> <character> <factor>
#> AAACAACGAATAGTTC-1 151507 L1
#> AAACAAGTATCTCCCA-1 151507 L3
#> AAACAATCTACTAGCA-1 151507 L1
#> AAACACCAATAACTGC-1 151507 WM
#> AAACAGCTTTCAGAAG-1 151507 L6
#> AAACAGGGTCTATATT-1 151507 L6The object is a SpatialExperiment, containing spot-level gene expression counts, spatial coordinates, tissue images, and manually annotated cortical layers.
Visualise the tissue morphology
# for simplicity of the vignette, work with one sample
dlpfc_spe <- dlpfc_spe[ , dlpfc_spe$sample_id == "151675"]
# visualise the cortical annotations
vis_clus(spe = dlpfc_spe,
colors = libd_layer_colors,
clustervar = "spatialLIBD",
point_size = 1.5)
#> adding rname 'https://spatial-dlpfc.s3.us-east-2.amazonaws.com/images/151675_tissue_lowres_image.png'
Defining a oligodendrocyte gene signature
For demonstration purposes, we use a collection of well-established oligodendrocyte marker genes.
oligo_signature <- c(
"ENSG00000167748", "ENSG00000166923", "ENSG00000174607", "ENSG00000135678",
"ENSG00000089123", "ENSG00000183018", "ENSG00000129538", "ENSG00000136937",
"ENSG00000105695", "ENSG00000168314", "ENSG00000184113", "ENSG00000124920",
"ENSG00000054274", "ENSG00000122507", "ENSG00000086205", "ENSG00000114867",
"ENSG00000169532", "ENSG00000166292", "ENSG00000141338"
)
# verify if the gene representationas are consistent in signature and data
head(oligo_signature %in% rownames(dlpfc_spe))
#> [1] TRUE TRUE TRUE TRUE TRUE TRUECompute SPAROscore
Since SpatialExperiment objects are supported directly by SPAROscore, pathway scores can be computed without extracting the expression matrix.
dlpfc_spe <- sparoscore(data = dlpfc_spe,
assay = "counts",
signatures = oligo_signature,
prefix = "Oligo_")
#> SPAROscore says: Calculating column-wise geometric averages
#> SPAROscore says: Ranking a sparseMatrix object
#> Warning: SPAROscore says: The following 2 signature genes are missing in the
#> input dataset: ENSG00000054274, ENSG00000169532
# view the scores appended to colData
head(colData(dlpfc_spe)[,c("sample_id", "spatialLIBD", "Oligo_SPAROscore")])
#> DataFrame with 6 rows and 3 columns
#> sample_id spatialLIBD Oligo_SPAROscore
#> <character> <factor> <numeric>
#> AAACAACGAATAGTTC-1 151675 L1 0.0000000
#> AAACAAGTATCTCCCA-1 151675 L3 0.0593458
#> AAACACCAATAACTGC-1 151675 WM 0.1759591
#> AAACAGAGCGACTCCT-1 151675 L3 0.0000000
#> AAACAGGGTCTATATT-1 151675 L5 0.0000000
#> AAACAGTGTTCCTGGG-1 151675 WM 0.0881556The resulting pathway scores are appended to the colData of the SpatialExperiment object, while the computed gene ranks are stored as an additional assay for reuse in subsequent analyses.
Visualising pathway activity across the tissue
One of the major advantages of spatial transcriptomics is the ability to visualise biological processes within their anatomical context.
Using the spot coordinates stored in the SpatialExperiment object, the SPAROscore can be projected directly onto the tissue section.
Higher oligodendrocyte signature scores are observed throughout the
white matter, whereas substantially lower scores (mostly zero) are
evident within the cortex layers as illustrated by the colours being
empty for most spatial spots outside the white matter.
This agrees with the known cellular composition of the tissue and demonstrates that SPAROscore successfully captures the spatial distribution of oligodendrocyte gene expression.
Pathway activity across anatomical regions
The manual annotations for each cortical layer enables quantitative comparison of oligodendrocyte pathway activity.
# differentiate white matter from cortex based on annotations
dlpfc_metadata <- as.data.frame(colData(dlpfc_spe)) %>%
filter(layer_guess != "NA") %>% # remove missing annotations
mutate(Region = ifelse(layer_guess == "WM",
"White_Matter", "Cortex"))
dlpfc_metadata %>%
ggplot(aes(y = reorder(Region, Oligo_SPAROscore, median),
x = Oligo_SPAROscore,
fill = Region)) +
geom_boxplot() +
labs(x = "Oligo SPAROscore",
y = "DLPFC Regions",
title = "Pathway activity in cortical layers") +
scale_fill_manual(values = c("#56B4E9", "#D55E00"))+
theme_classic() +
theme(legend.position = "none") 
The oligodendrocyte signature exhibits high scores in the white matter
but
lower scores within the cortex layers. This observation is consistent
with the known neuroanatomy of the cerebral cortex, where
oligodendrocyte cell bodies are concentrated within the white matter
while the cortical grey matter predominantly contains neuronal
cells.
Statistical comparison
To further illustrate the difference, oligodendrocyte signature scores between cortical spots and white matter can be compared statistically.
wilcox.test(Oligo_SPAROscore ~ Region, data = dlpfc_metadata)
#>
#> Wilcoxon rank sum test with continuity correction
#>
#> data: Oligo_SPAROscore by Region
#> W = 139971, p-value < 2.2e-16
#> alternative hypothesis: true location shift is not equal to 0The significant reduction in oligodendrocyte signature scores within grey matter provides quantitative evidence that SPAROscore accurately recovers known anatomical organization from spatial transcriptomics data.
Summary
This example demonstrates how SPAROscore can be integrated directly into spatial transcriptomics analyses. By quantifying a oligodendrocyte gene signature, the method successfully distinguishes cortical grey matter from white matter and recovers known tissue architecture using only gene expression information.
Such analyses extend naturally to signalling pathways, cell-type signatures and disease-associated gene programmes, providing a flexible framework for investigating biological processes in their spatial context.
Conclusion
This vignette demonstrates how SPAROscore can be integrated into analysis workflows spanning the three major transcriptomics technologies: bulk RNA sequencing, single-cell RNA sequencing, and spatial transcriptomics.
Although these technologies differ substantially in their data structures and biological resolution, SPAROscore provides a unified framework for quantifying gene signature activity by ranking genes within each sample, cell, or spatial domain.
Together, these analyses demonstrate that SPAROscore provides a flexible and scalable approach for pathway-level analysis across diverse transcriptomics applications. By producing interpretable signature scores that are robust to differences in sequencing depth, sparsity, and data modality, SPAROscore enables biologically meaningful comparisons across samples, cells, and spatial domains while integrating naturally into existing Bioconductor and Seurat workflows.
For detailed information on the supported input data structures, gene signature representations, and object-specific behaviour, see the Data Structures vignette: Data structures
A subsequent vignette, on advanced usage, describes additional functionality including usage of custom caps to fine tune the adaptive rank capping strategy of SPAROscore, performance optimisation, tertiary parameter tuning, and advanced scoring options for specialised analysis workflows. The advance usage vignette can be found here: Advanced options
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] future_1.75.0 ggspavis_1.19.1
#> [3] spatialLIBD_1.25.1 SpatialExperiment_1.23.0
#> [5] muscData_1.27.1 ExperimentHub_3.3.1
#> [7] AnnotationHub_4.3.2 BiocFileCache_3.3.0
#> [9] dbplyr_2.6.0 edgeR_4.11.4
#> [11] limma_3.69.2 scater_1.41.2
#> [13] ggplot2_4.0.3 scuttle_1.23.1
#> [15] SingleCellExperiment_1.35.2 Seurat_5.5.1
#> [17] SeuratObject_5.4.0 sp_2.2-3
#> [19] GSEABase_1.75.0 graph_1.91.0
#> [21] annotate_1.91.0 XML_3.99-0.23
#> [23] AnnotationDbi_1.75.2 DelayedArray_0.39.3
#> [25] SparseArray_1.13.2 S4Arrays_1.13.0
#> [27] abind_1.4-8 Matrix_1.7-5
#> [29] dplyr_1.2.1 msigdbr_26.1.0
#> [31] airway_1.33.2 SummarizedExperiment_1.43.0
#> [33] Biobase_2.73.1 GenomicRanges_1.65.1
#> [35] Seqinfo_1.3.0 IRanges_2.47.2
#> [37] S4Vectors_0.51.5 BiocGenerics_0.59.10
#> [39] generics_0.1.4 MatrixGenerics_1.25.0
#> [41] matrixStats_1.5.0 SPAROscore_0.99.3
#> [43] rmarkdown_2.31
#>
#> loaded via a namespace (and not attached):
#> [1] bitops_1.0-9 spatstat.sparse_3.2-0
#> [3] httr_1.4.8 RColorBrewer_1.1-3
#> [5] doParallel_1.0.17 tools_4.6.1
#> [7] sctransform_0.4.3 DT_0.34.0
#> [9] R6_2.6.1 lazyeval_0.2.3
#> [11] uwot_0.2.4 GetoptLong_1.1.1
#> [13] withr_3.0.3 gridExtra_2.3.1
#> [15] progressr_1.0.0 cli_3.6.6
#> [17] spatstat.explore_3.8-1 fastDummies_1.7.6
#> [19] labeling_0.4.3 sass_0.4.10
#> [21] S7_0.2.2 spatstat.data_3.1-9
#> [23] ggridges_0.5.7 pbapply_1.7-4
#> [25] Rsamtools_2.29.0 sessioninfo_1.2.4
#> [27] parallelly_1.48.0 attempt_0.3.1
#> [29] RSQLite_3.53.3 BiocIO_1.23.3
#> [31] shape_1.4.6.1 ica_1.0-3
#> [33] spatstat.random_3.5-0 ggbeeswarm_0.7.3
#> [35] lifecycle_1.0.5 yaml_2.3.12
#> [37] Rtsne_0.17 paletteer_1.7.0
#> [39] grid_4.6.1 blob_1.3.0
#> [41] promises_1.5.0 crayon_1.5.3
#> [43] miniUI_0.1.2 lattice_0.22-9
#> [45] beachmat_2.29.0 cowplot_1.2.0
#> [47] cigarillo_1.3.1 KEGGREST_1.53.5
#> [49] magick_2.9.1 sys_3.4.3
#> [51] maketools_1.3.2 pillar_1.11.1
#> [53] knitr_1.51 ComplexHeatmap_2.29.0
#> [55] rjson_0.2.23 future.apply_1.20.2
#> [57] codetools_0.2-20 glue_1.8.1
#> [59] spatstat.univar_3.2-0 data.table_1.18.4
#> [61] vctrs_0.7.3 png_0.1-9
#> [63] spam_2.11-4 gtable_0.3.6
#> [65] assertthat_0.2.1 rematch2_2.1.2
#> [67] cachem_1.1.0 xfun_0.60
#> [69] mime_0.13 ggside_0.4.1
#> [71] survival_3.8-9 iterators_1.0.14
#> [73] statmod_1.5.2 fitdistrplus_1.2-6
#> [75] ROCR_1.0-12 nlme_3.1-170
#> [77] bit64_4.8.2 filelock_1.0.3
#> [79] RcppAnnoy_0.0.23 bslib_0.11.0
#> [81] irlba_2.3.7 vipor_0.4.7
#> [83] KernSmooth_2.23-26 otel_0.2.0
#> [85] colorspace_2.1-3 DBI_1.3.0
#> [87] tidyselect_1.2.1 bit_4.6.0
#> [89] compiler_4.6.1 curl_7.1.0
#> [91] httr2_1.3.0 BiocNeighbors_2.7.2
#> [93] plotly_4.12.0 rtracklayer_1.73.0
#> [95] scales_1.4.0 lmtest_0.9-40
#> [97] rappdirs_0.3.4 stringr_1.6.0
#> [99] digest_0.6.39 goftest_1.2-3
#> [101] spatstat.utils_3.2-4 benchmarkmeData_2.0.0
#> [103] XVector_0.53.0 htmltools_0.5.9
#> [105] pkgconfig_2.0.3 sparseMatrixStats_1.25.0
#> [107] fastmap_1.2.0 GlobalOptions_0.1.4
#> [109] rlang_1.3.0 htmlwidgets_1.6.4
#> [111] shiny_1.14.0 DelayedMatrixStats_1.35.0
#> [113] farver_2.1.2 jquerylib_0.1.4
#> [115] zoo_1.8-15 jsonlite_2.0.0
#> [117] BiocParallel_1.47.0 config_0.3.2
#> [119] RCurl_1.98-1.19 BiocSingular_1.29.0
#> [121] magrittr_2.0.5 dotCall64_1.2
#> [123] patchwork_1.3.2 Rcpp_1.1.2
#> [125] babelgene_22.9 viridis_0.6.5
#> [127] reticulate_1.46.0 stringi_1.8.7
#> [129] MASS_7.3-66 plyr_1.8.9
#> [131] parallel_4.6.1 listenv_1.0.0
#> [133] ggrepel_0.9.8 deldir_2.0-4
#> [135] Biostrings_2.81.5 splines_4.6.1
#> [137] tensor_1.5.1 circlize_0.4.18
#> [139] locfit_1.5-9.12 igraph_2.3.3
#> [141] spatstat.geom_3.8-1 RcppHNSW_0.7.0
#> [143] buildtools_1.0.0 reshape2_1.4.5
#> [145] ScaledMatrix_1.21.0 BiocVersion_3.24.0
#> [147] evaluate_1.0.5 golem_1.0.1
#> [149] BiocManager_1.30.27 foreach_1.5.2
#> [151] httpuv_1.6.17 RANN_2.6.2
#> [153] tidyr_1.3.2 purrr_1.2.2
#> [155] polyclip_1.10-7 benchmarkme_1.0.8
#> [157] clue_0.3-68 scattermore_1.2
#> [159] BiocBaseUtils_1.15.1 rsvd_1.0.5
#> [161] xtable_1.8-8 restfulr_0.0.17
#> [163] RSpectra_0.16-2 later_1.4.8
#> [165] viridisLite_0.4.3 tibble_3.3.1
#> [167] GenomicAlignments_1.49.1 memoise_2.0.1
#> [169] beeswarm_0.4.0 cluster_2.1.8.2
#> [171] shinyWidgets_0.9.1 globals_0.19.1