Pathway Activity Analysis for Single Cells: Pathway Activity Analysis for spatial transcriptomics

loading packages

library(ggplot2)
library(ggrepel)
library(reshape2)
library(msigdbr)
library(Seurat)
library(SeuratData)
library(PaaSc)

Data preparation

We use the stxBrain data from SeuratData. Further details can be found in spatial_vignette.

# InstallData("stxBrain")
brain <- LoadData("stxBrain", type = "posterior1")

Preprocessing

We first preprocess the data using standard scRNA-seq pipeline.

# Data normalization by SCTransform
options(matrixStats.useNames.NA = "deprecated")
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)

# Dimensionality reduction and clustering
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30, verbose = FALSE)
brain <- FindClusters(brain, verbose = FALSE)
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30, verbose = FALSE)

# visualization the cluster labels in UMAP space
DimPlot(brain, reduction = "umap", label = TRUE)

# visualize the cluster label in image
SpatialDimPlot(brain, label = TRUE, label.size = 5)

Computing pathway activity scores

Here we tested all Gene Ontology (GO) temrs defined in MSigDB.

mm_c5 <- msigdbr(species = "Mus musculus", category = "C5")
mm_c5_gobp <- mm_c5[mm_c5$gs_subcat == "GO:BP", ]
background_geneset <- split(mm_c5_gobp$gene_symbol, mm_c5_gobp$gs_name)
length(background_geneset)

[1] 7656

brain <- computeMCA(brain, nmcs = 20, assay = "SCT")

2.94 sec elapsed

11.139 sec elapsed

2.153 sec elapsed

example <- c('GOBP_RESPONSE_TO_COCAINE',
             'GOBP_CEREBELLAR_CORTEX_FORMATION',
             'GOBP_CEREBELLAR_CORTEX_DEVELOPMENT',
             'GOBP_CEREBELLAR_CORTEX_MORPHOGENESIS',
             'GOBP_CEREBELLAR_PURKINJE_CELL_LAYER_DEVELOPMENT',
             'GOBP_CGMP_MEDIATED_SIGNALING',
             'GOBP_PROTEIN_KINASE_B_SIGNALING',
             'GOBP_MULTICELLULAR_ORGANISM_GROWTH')

pathway_geneset <- background_geneset[example]

gene_rate <- getGeneRate(background.geneset = background_geneset, pathway.geneset = pathway_geneset)
regression_data <- doRegression(brain, gene.rate = gene_rate)
score_data <- computeScore(brain, regression.data = regression_data, pvalue = 0.05, weight = FALSE, normalize = "z-score")
score_data[1:5, 1]

[1] -0.9497684  1.6081140 -0.8975759 -0.3700795 -0.5397234

Spatial mutual information

We use spatial mutual information (SMI) to measure the spatial pattern of pathway activity scores.

score_data <- score_data[, colSums(is.na(score_data)) == 0]
x <- GetTissueCoordinates(brain@images$posterior1)$x
y <- GetTissueCoordinates(brain@images$posterior1)$y
spatialMI <- calculateSpatialMI(x = x, y = y, z = score_data)
spatialMI <- data.frame(spatialMI, qvalue = p.adjust(spatialMI$pvalue))
spatialMI

                                                        MI pvalue
GOBP_RESPONSE_TO_COCAINE                        1.55902128      0
GOBP_CEREBELLAR_CORTEX_FORMATION                1.53493913      0
GOBP_CEREBELLAR_CORTEX_DEVELOPMENT              1.53493913      0
GOBP_CEREBELLAR_CORTEX_MORPHOGENESIS            1.53493913      0
GOBP_CEREBELLAR_PURKINJE_CELL_LAYER_DEVELOPMENT 1.53493913      0
GOBP_CGMP_MEDIATED_SIGNALING                    0.01831666      1
GOBP_PROTEIN_KINASE_B_SIGNALING                 0.01822519      1
GOBP_MULTICELLULAR_ORGANISM_GROWTH              0.01936508      1
                                                qvalue
GOBP_RESPONSE_TO_COCAINE                             0
GOBP_CEREBELLAR_CORTEX_FORMATION                     0
GOBP_CEREBELLAR_CORTEX_DEVELOPMENT                   0
GOBP_CEREBELLAR_CORTEX_MORPHOGENESIS                 0
GOBP_CEREBELLAR_PURKINJE_CELL_LAYER_DEVELOPMENT      0
GOBP_CGMP_MEDIATED_SIGNALING                         1
GOBP_PROTEIN_KINASE_B_SIGNALING                      1
GOBP_MULTICELLULAR_ORGANISM_GROWTH                   1

# visualize some interesting pathways
pathway <- "GOBP_RESPONSE_TO_COCAINE"
brain@meta.data[pathway] <- score_data[, pathway]
SpatialFeaturePlot(brain, features = pathway)

pathway <- c("GOBP_CEREBELLAR_CORTEX_FORMATION", "GOBP_CEREBELLAR_PURKINJE_CELL_LAYER_DEVELOPMENT")
brain@meta.data[pathway] <- score_data[, pathway]
SpatialFeaturePlot(brain, features = pathway)