---
title: "Use of SpatialDecon in a Spatial Transcriptomics dataset"
output:
rmarkdown::html_vignette:
toc: true
vignette: >
%\VignetteIndexEntry{Use of SpatialDecon in a Spatial Transcriptomics dataset}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
<style>
p.caption {
font-size: 1.5em;
}
</style>
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
### Installation
```{r installation, eval=FALSE}
if(!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SpatialDecon")
```
### Overview
This vignette demonstrates the use of the SpatialDecon package to estimate cell
abundance in spatial gene expression studies.
The workflow demonstrated here focuses on analysis of Seurat objects, which
are commonly used to store Visium and Spatial Transcriptomics data. See our other
vignettes for examples in GeoMx data.
We'll analyze a Spatial Transcriptomics dataset from a HER2+ breast tumor, looking for abundance of
different immune cell types.
### Data preparation
First, we load the package:
```{r setup}
library(SpatialDecon)
library(SeuratObject)
```
Let's load the example ST Seurat object and examine it:
```{r loaddata}
# download data:
con <- gzcon(url("https://github.com/almaan/her2st/raw/master/data/ST-cnts/G1.tsv.gz"))
txt <- readLines(con)
temp <- read.table(textConnection(txt), sep = "\t", header = TRUE, row.names = 1)
# parse data
raw = t(as.matrix(temp))
norm = sweep(raw, 2, colSums(raw), "/") * mean(colSums(raw))
x = as.numeric(substr(rownames(temp), 1, unlist(gregexpr("x", rownames(temp))) - 1))
y = -as.numeric(substr(rownames(temp), unlist(gregexpr("x", rownames(temp))) + 1, nchar(rownames(temp))))
# put into a seurat object:
andersson_g1 = CreateSeuratObject(counts = raw, assay="Spatial")
andersson_g1@meta.data$x = x
andersson_g1@meta.data$y = y
```
### Cell profile matrices
A "cell profile matrix" is a pre-defined matrix that specifies the expected
expression profiles of each cell type in the experiment.
The SpatialDecon library comes with one such matrix pre-loaded, the "SafeTME"
matrix, designed for estimation of immune and stroma cells in the tumor
microenvironment.
(This matrix was designed to avoid genes commonly expressed by cancer cells;
see the SpatialDecon manuscript for details.)
Let's take a glance at the safeTME matrix:
```{r showsafetme, fig.height=5, fig.width=10, fig.cap = "The safeTME cell profile matrix"}
data("safeTME")
data("safeTME.matches")
signif(safeTME[seq_len(3), seq_len(3)], 2)
heatmap(sweep(safeTME, 1, apply(safeTME, 1, max), "/"),
labRow = NA, margins = c(10, 5))
```
For studies of other tissue types, we have provided a library of cell profile
matrices, available on Github and downloadable with the "download_profile_matrix" function.
For a complete list of matrices, see [CellProfileLibrary GitHub Page](https://github.com/Nanostring-Biostats/CellProfileLibrary/tree/NewProfileMatrices).
Below we download a matrix of cell profiles derived from scRNA-seq of a mouse
spleen.
```{r downloadmatrix, fig.height=7, fig.width=10, fig.cap = "The Mouse Spleen profile matrix", eval=T}
mousespleen <- download_profile_matrix(species = "Mouse",
age_group = "Adult",
matrixname = "Spleen_MCA")
dim(mousespleen)
mousespleen[1:4,1:4]
head(cellGroups)
metadata
heatmap(sweep(mousespleen, 1, apply(mousespleen, 1, max), "/"),
labRow = NA, margins = c(10, 5), cexCol = 0.7)
```
For studies where the provided cell profile matrices aren't sufficient or if a specific single cell dataset is wanted, we can make a custom profile matrix using the function create_profile_matrix().
This mini single cell dataset is a fraction of the data from Kinchen, J. et al. Structural Remodeling of the Human Colonic Mesenchyme in Inflammatory Bowel Disease. Cell 175, 372-386.e17 (2018).
```{r single cell data}
data("mini_singleCell_dataset")
mini_singleCell_dataset$mtx@Dim # genes x cells
as.matrix(mini_singleCell_dataset$mtx)[1:4,1:4]
head(mini_singleCell_dataset$annots)
table(mini_singleCell_dataset$annots$LabeledCellType)
```
**Pericyte cell** and **smooth muscle cell of colon** will be dropped from this matrix due to low cell count. The average expression across all cells of one type is returned so the more cells of one type, the better reflection of the true gene expression. The confidence in these averages can be changed using the minCellNum filter.
```{r creatematrix, fig.height=7, fig.width=10, fig.cap = "Custom profile matrix"}
custom_mtx <- create_profile_matrix(mtx = mini_singleCell_dataset$mtx, # cell x gene count matrix
cellAnnots = mini_singleCell_dataset$annots, # cell annotations with cell type and cell name as columns
cellTypeCol = "LabeledCellType", # column containing cell type
cellNameCol = "CellID", # column containing cell ID/name
matrixName = "custom_mini_colon", # name of final profile matrix
outDir = NULL, # path to desired output directory, set to NULL if matrix should not be written
normalize = FALSE, # Should data be normalized?
minCellNum = 5, # minimum number of cells of one type needed to create profile, exclusive
minGenes = 10, # minimum number of genes expressed in a cell, exclusive
scalingFactor = 5, # what should all values be multiplied by for final matrix
discardCellTypes = TRUE) # should cell types be filtered for types like mitotic, doublet, low quality, unknown, etc.
head(custom_mtx)
heatmap(sweep(custom_mtx, 1, apply(custom_mtx, 1, max), "/"),
labRow = NA, margins = c(10, 5), cexCol = 0.7)
```
Custom matrices can be created from all single cell data classes as long as a counts matrix and cell annotations can be passed to the function. Here is an example of creating a matrix using a Seurat object.
```{r createSeuratmatrix}
library(SeuratObject)
data("mini_singleCell_dataset")
rownames(mini_singleCell_dataset$annots) <- mini_singleCell_dataset$annots$CellID
seuratObject <- CreateSeuratObject(counts = mini_singleCell_dataset$mtx, meta.data = mini_singleCell_dataset$annots)
Idents(seuratObject) <- seuratObject$LabeledCellType
rm(mini_singleCell_dataset)
annots <- data.frame(cbind(cellType=as.character(Idents(seuratObject)),
cellID=names(Idents(seuratObject))))
custom_mtx_seurat <- create_profile_matrix(mtx = seuratObject@assays$RNA@counts,
cellAnnots = annots,
cellTypeCol = "cellType",
cellNameCol = "cellID",
matrixName = "custom_mini_colon",
outDir = NULL,
normalize = FALSE,
minCellNum = 5,
minGenes = 10)
head(custom_mtx_seurat)
paste("custom_mtx and custom_mtx_seurat are identical", all(custom_mtx == custom_mtx_seurat))
```
### Deconvolving a Seurat object with the runspatialdecon function
Now our data is ready for deconvolution.
First we'll show how to use spatialdecon under the basic settings, omitting
optional bells and whistles.
```{r runiss}
res = runspatialdecon(object = andersson_g1,
bg = 0.01,
X = safeTME,
align_genes = TRUE)
str(res)
```
We're most interested in "beta", the matrix of estimated cell abundances.
```{r plotissres, fig.height = 5, fig.width = 8, fig.cap = "Cell abundance estimates"}
heatmap(res$beta, cexCol = 0.5, cexRow = 0.7, margins = c(10,7))
```
### Using the advanced settings of spatialdecon
spatialdecon has several abilities beyond basic deconvolution:
1. If given the nuclei counts for each region/observation, it returns results on
the scale of total cell counts. This option is generally not available for ST/Visium data.
2. If given the identities of pure tumor regions/observations, it infers a
handful of tumor-specific expression profiles and appends them to the cell
profile matrix. Doing this accounts for cancer cell-derived expression from any
genes in the cell profile matrix, removing contaminating signal from cancer
cells. This operation is complicated in ST/Visium studies, since they generally
do not contain regions previously marked as pure tumor. As a heuristic for identifying pure
tumor regions, we identify a small subset of regions with relatively low counts from
genes in the safeTME matrix vs. from the rest of the transcriptome.
3. If given raw count data, it derives per-data-point weights. For Visium/ST data,
we assume Poisson error.
4. If given a "cellmatches" argument, it sums multiple closely-related cell
types into a single score. E.g. if the safeTME matrix is used with the
cell-matching data object "safeTME.matches", it e.g. sums the "T.CD8.naive" and
"T.CD8.memory" scores into a single "CD8.T.cells" score.
Let's take a look at an example cell matching object:
```{r showmatches}
str(safeTME.matches)
```
Now let's mark some regions as pure tumor:
```{r puretumor}
sharedgenes = intersect(rownames(raw), rownames(safeTME))
plot(colSums(raw), colSums(raw[sharedgenes, ]), log = "xy")
hist(colSums(raw[sharedgenes, ]) / colSums(raw), breaks = 20)
alltumor = colSums(raw[sharedgenes, ]) / colSums(raw) < 0.03 # for alma data
table(alltumor)
```
Calculate weights:
```{r wts}
sd_from_noise <- runErrorModel(counts = raw, platform = "st")
wts <- 1 / sd_from_noise
```
Now let's run spatialdecon:
```{r runisstils}
# run spatialdecon with all the bells and whistles:
restils = runspatialdecon(object = andersson_g1, # Seurat object
X = safeTME, # safeTME matrix, used by default
bg = 0.01, # Recommended value of 0.01 in Visium/ST data
wts = wts, # weight
cellmerges = safeTME.matches, # safeTME.matches object, used by default
is_pure_tumor = alltumor, # identities of the Tumor segments/observations
n_tumor_clusters = 5) # how many distinct tumor profiles to append to safeTME
str(restils)
```
There are quite a few readouts here. Let's review the important ones:
* beta: the cell abundance scores of the rolled-up/major cell types
* beta.granular: the cell abundance scores of the granular cell types,
corresponding to the columns of the cell profile matrix
* yhat, resids: the fitted values and log2-scale residuals from the
deconvolution fit. Can be used to measure each observation's goodness-of-fit, a
possible QC metric.
* prop_of_nontumor: the beta matrix rescaled to give the proportions of
non-tumor cells in each observation.
* X: the cell profile matrix used, including newly-derived tumor-specific
columns.
To illustrate the derivation of tumor profiles, let's look at the cell profile
matrix output by spatialdecon:
```{r shownewX, fig.height=5, fig.width=8, fig.cap = "safeTME merged with newly-derived tumor profiles"}
heatmap(sweep(restils$X, 1, apply(restils$X, 1, max), "/"),
labRow = NA, margins = c(10, 5))
```
Note the new tumor-specific columns.
### Plotting deconvolution results
The "florets" function plotting cell abundances atop some
2-D projection.
Here, we'll plot cell abundances atop the first 2 principal components of the
data:
```{r florets, fig.width=8, fig.height=6, fig.cap = "TIL abundance plotted on PC space"}
# PCA of the normalized data:
pc = prcomp(t(log2(pmax(norm, 1))))$x[, c(1, 2)]
# run florets function:
par(mar = c(5,5,1,1))
layout(mat = (matrix(c(1, 2), 1)), widths = c(6, 2))
florets(x = pc[, 1], y = pc[, 2],
b = restils$beta, cex = .5,
xlab = "PC1", ylab = "PC2")
par(mar = c(0,0,0,0))
frame()
legend("center", fill = cellcols[rownames(restils$beta)],
legend = rownames(restils$beta), cex = 0.7)
```
```{r floretsxy, fig.width=8, fig.height=6, fig.cap = "TIL abundance plotted on physical space"}
# and plot florets in space:
par(mar = c(5,5,1,1))
layout(mat = (matrix(c(1, 2), 1)), widths = c(6, 2))
florets(x = x, y = y,
b = restils$beta, cex = .5,
xlab = "", ylab = "")
par(mar = c(0,0,0,0))
frame()
legend("center", fill = cellcols[rownames(restils$beta)],
legend = rownames(restils$beta), cex = 0.7)
```
So we can see that PC1 roughly tracks many vs. few immune cells, and PC2 tracks
the relative abundance of lymphoid/myeloid populations.
### Other functions
The SpatialDecon library includes several helpful functions for further
analysis/fine-tuning of deconvolution results.
#### Combining cell types:
When two cell types are too similar, the estimation of their abundances becomes
unstable. However, their sum can still be estimated easily.
The function "collapseCellTypes" takes a deconvolution results object and
collapses any colsely-related cell types you tell it to:
```{r collapse, fig.width=5, fig.height=5, fig.cap="Cell abundance estimates with related cell types collapsed"}
matching = list()
matching$myeloid = c( "macrophages", "monocytes", "mDCs")
matching$T.NK = c("CD4.T.cells","CD8.T.cells", "Treg", "NK")
matching$B = c("B")
matching$mast = c("mast")
matching$neutrophils = c("neutrophils")
matching$stroma = c("endothelial.cells", "fibroblasts")
collapsed = collapseCellTypes(fit = restils,
matching = matching)
heatmap(collapsed$beta, cexRow = 0.85, cexCol = 0.75)
```
### Session Info
```{r sessioninfo}
sessionInfo()
```