In tsinghua-ZTX/SPER: Spatial Paired Expression Ratio

library(SPER)
##  Other necessary libraries
library(Seurat)
if (!require("ggpubr")) install.packages("ggpubr")
library(ggpubr)

Example: adult mouse brain spatial transciptomics dataset

Dataset

We have provided a small example where you have 37 genes across 2695 spots in the spatial transcriptomics (ST) data 'ST_expression_data'. We also included the physical coordinates of all these spots as 'ST_coordinates'.

####  Spatial transcriptomics data
head(ST_expression_data)

####  Coordinates & Distance matrix
head(ST_coordinates)
spot_dist <- as.matrix(dist(ST_coordinates))
round_dist_mat <- ceiling(spot_dist / 120) * 120
round_dist_mat[which(round_dist_mat > 1800)] <- 1800

If wanted, one can plot the gene's spatial distribution among the spots using 'gg_gene' function

gg_gene("Grp",
        ST_expression_data,
        coord = ST_coordinates,
        title.anno = "Gene")

Another important component of SPER is the cell-type spatial compositional data. Usually, one can get them by using cell type decomposition algorithms on a ST dataset. Here we just provide the results which is obtained from RCTD algorithms (now as a part of spacexr).

####  spatial compositional data
head(spatial_comp_data)

Additionally, on can also generate the cell-type spatial distribution through the function used above

p1 <- gg_gene("L2.3.IT",
              spatial_comp_data,
              coord = ST_coordinates,
              title.anno = "Compositional")
p2 <- gg_gene("L4",
              spatial_comp_data,
              coord = ST_coordinates,
              title.anno = "Compositional")
p3 <- gg_gene("L5.IT",
              spatial_comp_data,
              coord = ST_coordinates,
              title.anno = "Compositional")
p4 <- gg_gene("L6.IT",
              spatial_comp_data,
              coord = ST_coordinates,
              title.anno = "Compositional")
ggarrange(p1, p2, p3, p4)

SPER:

With given data, we can run SPER and get the outcome as a list. The list has each gene as the name of an element where it actually contains the spatial paired expression matrix between the gene and all cell types.

res <- SPER(dist.mat = round_dist_mat,
            dist.list = (0:15)*120,
            ST.mat = ST_expression_data,
            CoDa.data = spatial_comp_data)
res[["L2.3.IT"]][1:10, 1:10]

Above results shows the SPER between L2/3 IT cell type and all genes. For example, for the 1st column, the values represents the spatial association between Grp and L2/3 IT cells at the given distance (0, 120 microns, 240 microns, etc.). Higher value means the gene is spatially more associated with the cells.
However, it will be impractical to compare the whole column between pairs. Therefore, we can apply a customized weight on the column so that we can integrate the values into one score, which is SPER score.

# custom_weights <- c(rep(1/5, 5), rep(0, 11))
custom_weights <- c(0, 0, rep(1/3, 3), rep(0, 11))
SPER_res <- weightSPER(res,
                       custom_weights)
head(SPER_res)

After applying the weight, we can finally get the matrix of SPER scores for each gene/cell-type pair. In the matrix, the values represent, depending the weights you have applied, the strength of spatial relationship between the gene and cells. In our example, we used the 'paracrine' weight which has its focus on the distance of 240, 360, and 480 microns.

Find potential paracrine signals by SPER scores

Once we get the SPER scores, the next step is to find genes that may potentially be paracrine signals. Several criteria for such a gene can be:
Theoretically, a ligand should not be a marker genes (cell-type-specific gene) for the target cells, although they both have a high spatial dependence to the target cells. Therefore, we want to remove the markers first when finding our markers.
A gene should be located at extracellular space to serve as a signals. Furthermore, if it is a well-known gene, we should expect it appearing in a ligand-receptor pair database as a ligand gene.

##  Provided cell-type-specific gene sets: P values are unnecessary, as long as it includes 'Gene' and 'Type'
marker_gene <- read.table("~/Desktop/Dropbox_NYU/Dropbox/SPER/SPER_public_data/Data/Marker_gene.txt", 
                          sep = "\t",
                          check.names = F)

##  Ligand-receptor pair dataset: gene1 is the ligand and gene2 is the receptor
LRP_data <- read.table("~/Desktop/Dropbox_NYU/Dropbox/SPER/SPER_public_data/Data/Mouse_LRP_data.txt", 
                       sep = "\t",
                       check.names = F)
##  Expression fraction data: showing genes' expression prevalence in each cell types
example_expr_frac <- read.table("~/Desktop/Dropbox_NYU/Dropbox/SPER/SPER_public_data/Data/Expr_frac_matrix.txt",
                                sep = "\t",
                                check.names = F)
# example_expr_frac <- example_expr_frac[rownames(SPER_res),]
head(example_expr_frac)

L5_IT_signals <- findSPERsignals(target.cell = "L5.IT",
                                 score.mat = SPER_res,
                                 expr.frac.mat = example_expr_frac,
                                 marker.gene = marker_gene,
                                 LRP.data = LRP_data)
L5_IT_signals

Visualization:

To visualize our SPER results, we may need more data than what have been provided in the package. Specifically, we will need: 1) an external single-cell RNA-seq reference, which can be used as a validation to tell you gene's expression profile across cell types; 2) cell-type-specific gene list (marker gene list). As SPER detects spatial association, there is no wonder that cell-type-specific genes will be highly associated with their own cells and thus having a high SPER score. In order to remove them from our result, it would be better if we can have such a list of genes and find the difference between it and our SPER results. 3) ligand-receptor pair database (or other specified background). Since we are trying to find potential paracrine signals, they should be expected to be enriched in the ligand gene set or at least extracellular gene set. Furthermore, with the help of ligand-receptor pair information, we can evaluate the receptor's expression patterns in the target cells as a validation.
All of the above datasets are too large to be included in a package, so one can download them from the Dropbox link and run the below codes by themselves.
As a reminder, one should change the file directory to where you download the data.

##  The Allen Mouse Brain Atlas: WHOLE CORTEX & HIPPOCAMPUS - 10X GENOMICS
allen_reference <- readRDS("~/Desktop/Dropbox_NYU/Dropbox/SPER/SPER_public_data/allen_reference_SCT.rds")
allen_reference <- SetIdent(allen_reference,
                            value = "subclass")

plotSPER(cell.type = "L5.IT",
         gene.id = "Grp",
         coord = ST_coordinates,
         ST.data = ST_expression_data,
         CoDa.data = spatial_comp_data,
         reference.data = allen_reference,
         SPER.data = res,
         dist.list = (0:15)*120,
         save.plots = F)