From 975d0c41f6cc12b5b75eb5d8f49bffb6e7c107f6 Mon Sep 17 00:00:00 2001 From: venkatk89 Date: Tue, 21 Jul 2026 18:07:23 +0930 Subject: [PATCH 1/3] removed transcriptopmics vignette --- ...paroscore_transcriptomics_technologies.Rmd | 768 ------------------ 1 file changed, 768 deletions(-) delete mode 100644 vignettes/sparoscore_transcriptomics_technologies.Rmd diff --git a/vignettes/sparoscore_transcriptomics_technologies.Rmd b/vignettes/sparoscore_transcriptomics_technologies.Rmd deleted file mode 100644 index d20017c..0000000 --- a/vignettes/sparoscore_transcriptomics_technologies.Rmd +++ /dev/null @@ -1,768 +0,0 @@ ---- -title: "SPAROscore: Transcriptomics Analysis Workflows" -package: "SPAROscore" -author: "Venkatesh Kamaraj" -date: "`r Sys.Date()`" -output: rmarkdown::html_vignette -vignette: > - %\VignetteIndexEntry{SPAROscore: Transcriptomics Analysis Workflows} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -```{r, include = FALSE} -knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>", - fig.width = 12, - fig.height = 6, - out.width = "100%" -) -``` - - -```{r setup, warning=FALSE, message=FALSE} -# load the necessary libraries -library(SPAROscore) -library(airway) -library(ggplot2) -library(dplyr) -library(edgeR) -library(Seurat) -library(SeuratData) -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 SeuratData 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](https://doi.org/10.1371/journal.pone.0099625) 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). - -```{r} -# 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")] - -# 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)) -``` -It 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 -```{r} -bulk_scores <- sparoscore( data = counts_data, signatures = gc_response_genes) - -bulk_scores -``` -The 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 -```{r} -metadata$SPAROscore <- bulk_scores[,1] - -metadata[, c("Run", "dex", "SPAROscore")] -``` -Appending 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. - -```{r} -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 - -```{r} -wilcox.test(SPAROscore ~ dex, data = metadata) -``` - -Although 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. - -```{r, warning=FALSE} -# 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) - -# 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](https://doi.org/10.1371/journal.pone.0099625). - -```{r} -# 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) - -# 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.](https://doi.org/10.1038/nbt.4042), which is available through the -SeuratData package will be used. The dataset consists of ~14,000 -peripheral blood mononuclear cells (PBMCs) collected from -lupus patients before and after stimulation with interferon-β (IFN-β). - -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 SeuratData - -```{r, warning=FALSE} -# install the dataset for first-time use -InstallData("ifnb") - -# load the ifnb dataset -ifnb <- LoadData("ifnb") - -# check the layer in the RNA assay -Layers(ifnb[["RNA"]]) - -# view the metadata of the ifnb seurat object -head(ifnb[[]]) -``` - -The 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 `seurat_annotations` metadata column. - -## Load the gene signature from MSigDB - -```{r} -ifn_signature <- msigdbr( species = "Homo sapiens", collection = "H" ) %>% - filter(gs_name == "HALLMARK_INTERFERON_ALPHA_RESPONSE") %>% - pull(gene_symbol) %>% - unique() - - -# check if the genes are present in the expression data -head(ifn_signature %in% rownames(ifnb)) -``` - -## Compute SPAROscore -The interferon-beta response for every cell can be quantified using SPAROscore - -```{r} -ifnb <- sparoscore(data = ifnb, - assay = "RNA", - layer = "counts", - signatures = ifn_signature, - prefix = "IFN_" ) - -# check the layer in the RNA assay after scoring -Layers(ifnb[["RNA"]]) - -# view the metadata after scoring -head(ifnb[[]]) -``` - -After 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. - -```{r} -VlnPlot(ifnb, - features = "IFN_SPAROscore", - group.by = "stim", - pt.size = 0) + - scale_fill_manual(values = c("#2166AC", "#B2182B"))+ - theme(legend.position = "none") - -``` - -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. - -```{r} -wilcox.test(IFN_SPAROscore ~ stim, ifnb[[]]) -``` - -The 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. - -```{r, warning=FALSE} -VlnPlot(ifnb, - features = "IFN_SPAROscore", - group.by = "seurat_annotations", - pt.size = 0, - flip = TRUE ) - -``` - -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. - -```{r, warning=FALSE, message=FALSE} -# 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) - -# 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. - -```{r} -# 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_" -) - -# view the metadata after scoring for the second signature -head(ifnb[[]]) -``` -This 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 (seurat_annotations) 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 -```{r} -pseudobulk_ifn <- AggregateExpression(ifnb, - group.by = c("seurat_annotations", "stim"), - assays = "RNA", - return.seurat = TRUE, - slot = "counts") - -pseudobulk_ifn[[]] -``` -Each 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 -```{r} -pseudobulk_ifn <- sparoscore(data = pseudobulk_ifn, - signatures = ifn_signature, - assay = "RNA", - layer = "counts") -``` -### Visualise pseudo-bulk level scores - -```{r} -pseudobulk_ifn[[]] %>% - ggplot(aes(y = seurat_annotations, 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. - -```{r} -# load the DLPFC Visium dataset -dlpfc_spe <- fetch_data(type = "spe") - -# check some domain-level metadata present in the object -head(colData(dlpfc_spe)[,c("sample_id", "spatialLIBD")]) - -``` -The object is a SpatialExperiment, containing spot-level gene expression counts, -spatial coordinates, tissue images, and manually annotated cortical layers. - -## Visualise the tissue morphology -```{r} -# 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) -``` - -## Defining a oligodendrocyte gene signature - -For demonstration purposes, we use a collection of well-established -oligodendrocyte marker genes. -```{r} -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)) -``` - - -## Compute SPAROscore - -Since SpatialExperiment objects are supported directly by SPAROscore, -pathway scores can be computed without extracting the expression matrix. - -```{r} -dlpfc_spe <- sparoscore(data = dlpfc_spe, - assay = "counts", - signatures = oligo_signature, - prefix = "Oligo_") - -# view the scores appended to colData -head(colData(dlpfc_spe)[,c("sample_id", "spatialLIBD", "Oligo_SPAROscore")]) -``` -The 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. - -```{r} -vis_gene(spe = dlpfc_spe, - geneid = "Oligo_SPAROscore", - spatial = TRUE, - point_size = 1.5) -``` -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. - -```{r} -# 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. - -```{r} -wilcox.test(Oligo_SPAROscore ~ Region, data = dlpfc_metadata) -``` -The 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](sparoscore_data_structures.html) - -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](sparoscore_advanced_options.html) - - -```{r} -sessionInfo() -``` From 5da470885639b9f236adf55262cb0592140c226e Mon Sep 17 00:00:00 2001 From: venkatk89 Date: Tue, 21 Jul 2026 18:11:00 +0930 Subject: [PATCH 2/3] merged gitignores --- .gitignore | 3 +++ vignettes/.gitignore | 2 -- 2 files changed, 3 insertions(+), 2 deletions(-) delete mode 100644 vignettes/.gitignore diff --git a/.gitignore b/.gitignore index cdf82c5..53d3b2c 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,6 @@ .Rbuildignore SPAROscore.Rproj inst/doc +/vignettes/*.html +/vignettes/*.R + diff --git a/vignettes/.gitignore b/vignettes/.gitignore deleted file mode 100644 index 097b241..0000000 --- a/vignettes/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -*.html -*.R From 615985934182c62f76052d644852afd44d0991f5 Mon Sep 17 00:00:00 2001 From: venkatk89 Date: Tue, 21 Jul 2026 19:32:54 +0930 Subject: [PATCH 3/3] moved from SeuratData to bioconductor --- DESCRIPTION | 49 +- R/get_scores_generic.R | 2 +- R/helper_functions.R | 2 +- R/sparoscore_generic.R | 6 +- vignettes/sparoscore_data_structures.Rmd | 33 +- vignettes/sparoscore_getting_started.Rmd | 10 +- ...paroscore_transcriptomics_technologies.Rmd | 774 ++++++++++++++++++ 7 files changed, 828 insertions(+), 48 deletions(-) create mode 100644 vignettes/sparoscore_transcriptomics_technologies.Rmd diff --git a/DESCRIPTION b/DESCRIPTION index e6f5c9f..323961c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,24 +1,26 @@ -Package: SPAROscore Type: Package -Title: A package to compute gene signature scores from transcriptomics data +Package: SPAROscore +Title: A package to compute gene signature scores from transcriptomics + data Version: 0.99.0 -Authors@R: c( - person( - "Venkatesh", "Kamaraj", - email = "venkatesh.kamaraj@adelaide.edu.au", - role = c("aut", "cre"), - comment = c(ORCID = "0009-0008-7309-9810") - ) - ) -Description: SPAROscore is a gene signature scoring method designed to be - robust across diverse gene expression datasets. SPAROscore adapts to - varying levels of sparsity, allowing signature scores to be efficiently - computed across bulk, single-cell, and spatial transcriptomic datasets. -biocViews: GeneExpression, SingleCell, Spatial, Transcriptomics, - GeneSetEnrichment, Pathways, Software +Authors@R: + person("Venkatesh", "Kamaraj", , "venkatesh.kamaraj@adelaide.edu.au", role = c("aut", "cre"), + comment = c(ORCID = "0009-0008-7309-9810")) +Description: SPAROscore is a gene signature scoring method designed to be + robust across diverse gene expression datasets. SPAROscore adapts to + varying levels of sparsity, allowing signature scores to be + efficiently computed across bulk, single-cell, and spatial + transcriptomic datasets. The resulting scores quantify the relative + expression of a gene signature compared with the background expression + of each sample/cell/domain, making the scores straightforward to + interpret biologically. Internally, SPAROscore ranks genes within + each column and computes the Spearman footrule distance between the + observed ranks of the signature genes and a background gene expression + rank estimated from the geometric mean expression of the sample, cell, + or spatial domain. License: GPL-3 -Encoding: UTF-8 -RoxygenNote: 7.3.3 +URL: https://github.com/MangiolaLaboratory/SPAROscore +BugReports: https://github.com/MangiolaLaboratory/SPAROscore/issues Imports: DelayedArray, DelayedMatrixStats, @@ -37,17 +39,20 @@ Suggests: GSEABase, knitr, msigdbr, + muscData, rmarkdown, scater, Seurat, - SeuratData, SingleCellExperiment, SpatialExperiment, spatialLIBD, SummarizedExperiment, testthat (>= 3.0.0), tidyr -VignetteBuilder: knitr +VignetteBuilder: + knitr +biocViews: GeneExpression, SingleCell, Spatial, Transcriptomics, + GeneSetEnrichment, Pathways, Software Config/testthat/edition: 3 -URL: https://github.com/MangiolaLaboratory/SPAROscore -BugReports: https://github.com/MangiolaLaboratory/SPAROscore/issues +Encoding: UTF-8 +RoxygenNote: 7.3.3 diff --git a/R/get_scores_generic.R b/R/get_scores_generic.R index f4abe61..a9ff6b2 100644 --- a/R/get_scores_generic.R +++ b/R/get_scores_generic.R @@ -59,7 +59,7 @@ #' Defaults to NULL. #' #' When down_signatures is not NULL, -#' Final Score = Score(signatures) - Score(down_signatures) +#' Final Score <- Score(signatures) - Score(down_signatures) #' #' #' @param prefix Character string to be appended before the headers of the diff --git a/R/helper_functions.R b/R/helper_functions.R index 328c45a..839e86f 100644 --- a/R/helper_functions.R +++ b/R/helper_functions.R @@ -229,7 +229,7 @@ append_to_matrix_like_object <- function(matrix_like_object, numeric_vector){ #' #' #' @examples -#' \dontrun{ +#' \donttest{ #' counts <- matrix( #' c(0, 5, 10, #' 2, 0, 20, diff --git a/R/sparoscore_generic.R b/R/sparoscore_generic.R index f5fc845..92e23c6 100644 --- a/R/sparoscore_generic.R +++ b/R/sparoscore_generic.R @@ -111,7 +111,7 @@ #' Defaults to NULL. #' #' When down_signatures is not NULL, -#' Final Score = Score(signatures) - Score(down_signatures) +#' Final Score <- Score(signatures) - Score(down_signatures) #' #' #' @param prefix Character string to be appended before the headers of the @@ -218,7 +218,7 @@ #' # Seurat object #' # ------------------------------------------------------------------ #' -#' \dontrun{ +#' \donttest{ #' seurat_object <- sparoscore( #' data = seurat_object, #' signatures = c("CCR7", "IL7R", "LTB") @@ -257,7 +257,7 @@ #' # SingleCellExperiment / SummarizedExperiment #' # ------------------------------------------------------------------ #' -#' \dontrun{ +#' \donttest{ #' sce <- sparoscore( #' data = sce, #' signatures = c("CCR7", "IL7R", "LTB") diff --git a/vignettes/sparoscore_data_structures.Rmd b/vignettes/sparoscore_data_structures.Rmd index 2ddb5df..4e656cb 100644 --- a/vignettes/sparoscore_data_structures.Rmd +++ b/vignettes/sparoscore_data_structures.Rmd @@ -10,7 +10,7 @@ vignette: > %\VignetteEncoding{UTF-8} --- -```{r, include = FALSE} +```{r global-options, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", @@ -97,7 +97,7 @@ as well as the hallmark inflammatory response genes taken from the Molecular Signatures Database (MSigDB). ### Load the gene expression dataset -```{r} +```{r load-data} # load the airway dataset data(airway) @@ -109,7 +109,7 @@ metadata <- as.data.frame(colData(airway)) ``` ### Load the gene signatures -```{r} +```{r load-signature} # set the glucocorticoid response markers from the original study as a signature response_genes <- c("ENSG00000112936", "ENSG00000170606", "ENSG00000120129", "ENSG00000152795", "ENSG00000211445", "ENSG00000163884", @@ -131,7 +131,7 @@ inflammation_genes <- msigdbr(species = "Homo sapiens", collection = "H") %>% # Gene expression as matrix-like data structures ## Dense matrix -```{r} +```{r matrix} # load gene expression data as dense matrix counts_matrix <- as.matrix(counts_data) @@ -147,7 +147,7 @@ Note that one signature gene is absent from the input expression data, and SPAROscore reports this as a warning during scoring. ## SparseMatrix -```{r} +```{r smatrix} # load gene expression data as sparseMatrix counts_sparse <- as(counts_matrix, "sparseMatrix") @@ -164,7 +164,7 @@ them to dense matrices, reducing memory usage and improving computational efficiency. ## DelayedArrays -```{r} +```{r dmatrix} # load gene expression data as delayedArray counts_delayed <- DelayedArray(counts_matrix) @@ -180,7 +180,7 @@ Note that the gene ranking is performed natively on the DelayedArray object, enabling improved computational efficiency. ## Data frame -```{r} +```{r dataframe} # load gene expression data as delayedArray counts_df <- as.data.frame(counts_matrix) @@ -197,7 +197,7 @@ because matrix operations are more efficient for this step. ### Verify if the output scores are the same with differing input formats -```{r} +```{r check-results} identical(scores_matrix, scores_sparse) identical(scores_matrix, scores_delayed) identical(scores_matrix, scores_df) @@ -206,7 +206,7 @@ identical(scores_matrix, scores_df) # Gene signature representations ## Single signature as a character vector -```{r} +```{r character} # compute scores from a character vector signature scores_for_vector <- sparoscore(data = counts_matrix, signatures = response_genes) @@ -217,20 +217,21 @@ scores_for_vector ## Single signature as a GeneSet object -```{r} +```{r geneset} # get a GeneSet object response_geneset <- GeneSet(response_genes) # compute scores from a GeneSet object scores_for_geneset <- sparoscore(data = counts_sparse, - signatures = response_geneset, prefix = "hello") + signatures = response_geneset, + prefix = "hello") # view scores scores_for_geneset ``` ## Multiple signatures as a named list -```{r} +```{r list} # create a named list of vectors to score for both signatures simultaneously gene_signatures <- list("glucocorticoid_response" = response_genes, "inflammation_response" = inflammation_genes) @@ -244,7 +245,7 @@ scores_for_list ``` ## Multiple signatures as a GeneSetCollection object -```{r} +```{r genesetcollection} # create a GeneSetCollection object to score for both signatures simultaneously gene_set_collection <- GeneSetCollection( mapply(function(g, n) GeneSet(unique(g), setName = n), @@ -268,7 +269,7 @@ named `ranks`. In this example, let us work with a SingleCellExperiment object -```{r} +```{r sce} # convert airways to a SingleCellExperiment object airway_sce <- as(airway, "SingleCellExperiment") @@ -309,7 +310,7 @@ For Seurat objects, after scoring, SPAROscore returns a Seurat object with the scores and the rank caps appended to the metadata. The ranks are stored in a layer named `ranks` in the same assay as the counts. -```{r} +```{r seurat} # convert airway to a seurat object airway_seurat <- CreateSeuratObject( counts = assay(airway), @@ -364,6 +365,6 @@ For examples of in-depth transcriptomics analyses using SPAROscore, refer here: For examples of additional options and advanced usage, see: [Advanced options](sparoscore_advanced_options.html) -```{r} +```{r session-info} sessionInfo() ``` diff --git a/vignettes/sparoscore_getting_started.Rmd b/vignettes/sparoscore_getting_started.Rmd index ffb3c3b..496b502 100644 --- a/vignettes/sparoscore_getting_started.Rmd +++ b/vignettes/sparoscore_getting_started.Rmd @@ -10,7 +10,7 @@ vignette: > %\VignetteEncoding{UTF-8} --- -```{r, include = FALSE} +```{r global-options, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", @@ -61,7 +61,7 @@ reported in the [original study](https://doi.org/10.1371/journal.pone.0099625). 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). -```{r} +```{r load-dataset} # load the airway dataset data(airway) @@ -87,7 +87,7 @@ head(signature_genes %in% rownames(counts_data)) ``` ## Compute signature scores using SPAROscore -```{r} +```{r compute-scores} scores <- sparoscore(data = counts_data, signatures = signature_genes) # view the scores @@ -99,7 +99,7 @@ scores Since dexamethasone activates glucocorticoid signalling, treated samples are expected to exhibit higher glucocorticoid response scores than untreated controls. -```{r} +```{r visualise-results} # add the signature scores to the metadata metadata$SPAROscore <- scores[, 1] @@ -159,6 +159,6 @@ For examples of additional options and more advanced usage, see: [Advanced options](sparoscore_advanced_options.html) -```{r} +```{r session-info} sessionInfo() ``` diff --git a/vignettes/sparoscore_transcriptomics_technologies.Rmd b/vignettes/sparoscore_transcriptomics_technologies.Rmd new file mode 100644 index 0000000..88ed0a0 --- /dev/null +++ b/vignettes/sparoscore_transcriptomics_technologies.Rmd @@ -0,0 +1,774 @@ +--- +title: "SPAROscore: Transcriptomics Analysis Workflows" +package: "SPAROscore" +author: "Venkatesh Kamaraj" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{SPAROscore: Transcriptomics Analysis Workflows} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r global-options, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 12, + fig.height = 6, + out.width = "100%" +) +``` + + +```{r setup, warning=FALSE, message=FALSE} +# 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 SeuratData 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](https://doi.org/10.1371/journal.pone.0099625) 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). + +```{r bulk-1} +# 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")] + +# 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)) +``` +It 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 +```{r bulk-2} +bulk_scores <- sparoscore( data = counts_data, signatures = gc_response_genes) + +bulk_scores +``` +The 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 +```{r bulk-3} +metadata$SPAROscore <- bulk_scores[,1] + +metadata[, c("Run", "dex", "SPAROscore")] +``` +Appending 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. + +```{r bulk-4} +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 + +```{r bulk-5} +wilcox.test(SPAROscore ~ dex, data = metadata) +``` + +Although 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. + +```{r bulk-6, warning=FALSE} +# 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) + +# 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](https://doi.org/10.1371/journal.pone.0099625). + +```{r bulk-7} +# 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) + +# 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.](https://doi.org/10.1038/nbt.4042), which is available through the +muscData package will be used. We will work with a seurat object of the same to +showcase the versatality of SPAROscore. The dataset consists of ~29,000 +peripheral blood mononuclear cells (PBMCs) collected from +lupus patients before and after stimulation with interferon-β (IFN-β). + +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 + +```{r sc-1, warning=FALSE} +# load the ifnb data from the muscData library +ifnb_sce_data <- Kang18_8vs8() + +# 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" "data" + +# view the metadata of the ifnb seurat object +head(ifnb[[]]) +``` + +The 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 `seurat_annotations` metadata column. + +## Load the gene signature from MSigDB + +```{r sc-2} +ifn_signature <- msigdbr( species = "Homo sapiens", collection = "H" ) %>% + filter(gs_name == "HALLMARK_INTERFERON_ALPHA_RESPONSE") %>% + pull(gene_symbol) %>% + unique() + + +# check if the genes are present in the expression data +head(ifn_signature %in% rownames(ifnb)) +``` + +## Compute SPAROscore +The interferon-beta response for every cell can be quantified using SPAROscore + +```{r sc-3} +ifnb <- sparoscore(data = ifnb, + assay = "RNA", + layer = "counts", + signatures = ifn_signature, + prefix = "IFN_" ) + +# check the layer in the RNA assay after scoring +Layers(ifnb[["RNA"]]) + +# view the metadata after scoring +head(ifnb[[]]) +``` + +After 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. + +```{r sc-4} +VlnPlot(ifnb, + features = "IFN_SPAROscore", + group.by = "stim", + pt.size = 0) + + scale_fill_manual(values = c("#2166AC", "#B2182B"))+ + theme(legend.position = "none") + +``` + +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. + +```{r sc-5} +wilcox.test(IFN_SPAROscore ~ stim, ifnb[[]]) +``` + +The 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. + +```{r sc-6, warning=FALSE} +VlnPlot(ifnb, + features = "IFN_SPAROscore", + group.by = "cell", + pt.size = 0, + flip = TRUE ) + +``` + +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. + +```{r sc-7, warning=FALSE, message=FALSE} +# 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) + +# 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. + +```{r sc-8} +# 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_" +) + +# view the metadata after scoring for the second signature +head(ifnb[[]]) +``` +This 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 (seurat_annotations) 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 +```{r pb-1} +pseudobulk_ifn <- AggregateExpression(ifnb, + group.by = c("cell", "stim"), + assays = "RNA", + return.seurat = TRUE, + slot = "counts") + +pseudobulk_ifn[[]] +``` +Each 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 +```{r pb-2} +pseudobulk_ifn <- sparoscore(data = pseudobulk_ifn, + signatures = ifn_signature, + assay = "RNA", + layer = "counts") +``` +### Visualise pseudo-bulk level scores + +```{r pb-3} +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. + +```{r sp-1} +# load the DLPFC Visium dataset +dlpfc_spe <- fetch_data(type = "spe") + +# check some domain-level metadata present in the object +head(colData(dlpfc_spe)[,c("sample_id", "spatialLIBD")]) + +``` +The object is a SpatialExperiment, containing spot-level gene expression counts, +spatial coordinates, tissue images, and manually annotated cortical layers. + +## Visualise the tissue morphology +```{r sp-2} +# 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) +``` + +## Defining a oligodendrocyte gene signature + +For demonstration purposes, we use a collection of well-established +oligodendrocyte marker genes. +```{r sp-3} +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)) +``` + + +## Compute SPAROscore + +Since SpatialExperiment objects are supported directly by SPAROscore, +pathway scores can be computed without extracting the expression matrix. + +```{r sp-4} +dlpfc_spe <- sparoscore(data = dlpfc_spe, + assay = "counts", + signatures = oligo_signature, + prefix = "Oligo_") + +# view the scores appended to colData +head(colData(dlpfc_spe)[,c("sample_id", "spatialLIBD", "Oligo_SPAROscore")]) +``` +The 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. + +```{r sp-5} +vis_gene(spe = dlpfc_spe, + geneid = "Oligo_SPAROscore", + spatial = TRUE, + point_size = 1.5) +``` +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. + +```{r sp-6} +# 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. + +```{r sp-7} +wilcox.test(Oligo_SPAROscore ~ Region, data = dlpfc_metadata) +``` +The 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](sparoscore_data_structures.html) + +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](sparoscore_advanced_options.html) + + +```{r session-info} +sessionInfo() +```