---
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 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](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. 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 

```{r sc-1, warning=FALSE}
# load the ifnb data from the muscData library
ifnb_sce_data <- Kang18_8vs8()

# view the colData
head(colData(ifnb_sce_data))

# check how many genes * cells are present in the original dataset
dim(ifnb_sce_data)

# 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)

# 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 `cell` 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 (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
```{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()
```
