This beta of SPACO is currently not yet available on Bioconductor or CRAN. Therefore it must be installed from github via devtools.
library(devtools)
devtools::install_github("IMSBCompBio/SpaCo")
SPACO uses the importer function of the Seurat library to import 10X Visium spatial data. The imported data is preprocessed using Variance Stabilizing Transformations by the sctransform package (URAL) The number of variable features to keep can be defined and any variables that should be regressed out. In out tutorial we correct for mitochondrial and hemoglobin genes detected in the spots
#Specify data paths
library(SPACO)
data_dir="~/path/to/directory"
slice = "~/path/to/directory/with/the/H&E/stain"
filename ="name of the feature matrix"
spatial_file = "name of the tissue positions list"
#Initialize the object from raw (non-normalized data and use the Seurat library for pre-processing)
SpaCoObject <- read_10x_for_spaco(data_dir = data_dir,slice = slice,filename = filename,
variable_features_n = 3000,
spatial_file = spatial_file,
vars_to_regress = c("^mt-","^Hbb-"))
library(Seurat)
library(SPACO)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
#install data from SeuratData
#InstallData(“stxBrain.SeuratData”)
brain <- LoadData("stxBrain", type = "anterior1")
brain <- PercentageFeatureSet(brain, pattern = "^mt-" ,col.name = "percent.mt")
brain <- PercentageFeatureSet(brain, pattern = "^Hbb-" ,col.name = "percent.hbb")
brain <- SCTransform(brain, assay = "Spatial", variable.features.n = 3000, verbose = F)
SpaCoObject <- seurat_to_spaco(Seurat = brain, assay = "SCT", n_image= 1, slot = "scale.data")
Run the Spatial component analysis as dimensionality reduction and computer the number of informative spatial components.
In this step we actually do the spatial component analysis. Further, we determine the number of relevant spatial components via sampling. These can be found in the @nSpacs slot in the object.
SpaCoObject <- RunSCA(SpaCoObject,compute_nSpacs = TRUE)
## computing number of releveant spacs
SpaCoObject@nSpacs
## [1] 40
We can plot the computed meta gene projections directly from the SPaCoObject or transfer the projections right back into an existing Seurat Object for downstream processing like clustering. In the object they will be treated as dimensionality reduction objects. It is also possible to plot visualize the normalized expression of specific genes as well. All plotting functions in the SPACO library can be modified with additional ggplot commands.
spacplot <- Spaco_plot(SpaCoObject, spac = 1:4, ncol = NULL, combine = T)+ggtitle("Meta genes 1 to 4")
spacplot
Ttr <- feature_plot(SpaCoObject, features = "Ttr", ncol = NULL, combine = TRUE)+ggtitle("original")
Ttr
First we remove all spots from the Seurat object which have no direct neighbors as they violate SPaCo assumptions. The computed SPaCo projections are stored in the object slot “object[[”spaco”]]”.
brain <- subset_non_neighbour_cells(SpaCoObject, brain)
brain <- spacs_to_seurat(SpaCoObject, brain)
## copying significan projections into reduction slot spaco
To compare the spatial component analysis to the analog process using PCA in the Seurat pipeline we do all the downstream processing right in the Seurat object.
cc <- SpatialFeaturePlot(brain,features = c("Spac_1","Spac_2","Spac_3","Spac_4"),combine = T)
brain_2 <- RunPCA(brain, verbose = F)
dd <- SpatialFeaturePlot(brain_2,features = c("PC_1","PC_2","PC_3","PC_4"),combine = T)
brain <- RunUMAP(brain,reduction = "spaco",dims = 1:SpaCoObject@nSpacs,n.neighbors = 45, verbose = F)
brain_2 <- RunUMAP(brain_2,reduction = "pca",dims = 1:30, verbose = F)
brain <- FindNeighbors(brain,reduction = "spaco" , dims = 1:SpaCoObject@nSpacs, verbose = F)
brain_2 <- FindNeighbors(brain_2,reduction = "pca" , dims = 1:30, verbose = F)
brain <- FindClusters(brain,resolution = 0.24, verbose = F)
brain_2 <- FindClusters(brain_2, verbose = F)
aa <- DimPlot(brain,group.by="seurat_clusters")+ggtitle("SpaCo")
bb <- DimPlot(brain_2,group.by="seurat_clusters")+ggtitle("Pca")
a <- DimPlot(brain,reduction = "spaco")+ggtitle("SpaCo")
b <- DimPlot(brain_2,reduction = "pca")+ggtitle("Pca")
rr <- SpatialDimPlot(brain)+ggtitle("SpaCo")
qq <- SpatialDimPlot(brain_2)+ggtitle("Pca")
We can plot all plots using standard ggplot and patchwork grammar.
patchwork::wrap_plots(a,b,aa,bb,rr,qq,ncol=2)
a+b
aa+bb
rr+qq
DE_genes<- SVGTest(SpaCoObject)
DE_genes_sort <- DE_genes[order(DE_genes$score, decreasing = TRUE),]
sigs <- rownames(DE_genes[DE_genes$p.adjust<0.05,])
We can use the significant spatially variable genes for functional enrichment methods like GO-terms or KEGG pathways. Here we use clusterProfiler for this.
library(org.Mm.eg.db)
library(clusterProfiler)
sigs_entres <- na.omit(AnnotationDbi::select(
x = org.Mm.eg.db,
keys = sigs,
keytype = "SYMBOL",
columns = c("ENTREZID"))[,"ENTREZID"])
enrichGO(sigs_entres,'org.Mm.eg.db',ont="BP") -> GOBP
dotplot_brain <- dotplot(GOBP, showCategory=10)
dotplot_brain
For the identification of functional areas on the tissue slide it is possible to use the relevant spatial components for denoising of genes. This should only be done with spatial variable genes as non-spatial genes will be highly distorted by artefacts
SpaCoObject <- denoise_profiles(SpaCoObject)
denoised <- denoised_projection_plot(SpaCoObject,features = "Ybx1")+ggtitle("denoised")
original <- feature_plot(SpaCoObject, features = "Ybx1", ncol = NULL, combine = TRUE)+ggtitle("original")
original + denoised