Introduction to tpSVG
19 April 2024
Contents
The objective of tpSVG is to detect spatially variable genes (SVG) when analyzing spatially-resolved transcriptomics data. This includes both unsupervised features where there’s not additional information is supplied besides the (normalized) gene counts and spatial coordinates, but also the spatial variation explained besides some covariates, such as tissue anatomy or possibly cell type composition.
Compared to previous SVG detection tools, tpSVG provides a scalable solution to model the gene expression as counts instead of logarithm-transformed counts. While log transformation provides convenience to model the spatial gene expression by mapping count data to the continuous domain, hence enabling well-understood Gaussian models, log transformation distorts low expressed genes counts and create bias populating high-expressed genes. For example, the rank of genes based on their effect size are commonly used for dimensional reduction, or its input. Hence, estimating gene ranking correctly is very important. Gaussian models, exemplified with nnSVG, often dissociates the mean-variance relationship which is commonly assumed for counts data, and hence often prioritizes the highly expressed genes over the lowly expressed genes. In the figure below, we saw that nnSVG is susceptible to such mean-rank relationship, meaning highly expressed genes are often ranked highly. In contrast, the proposed tpSVG with Poisson distribution is not susceptible to this mean-rank relationship when examining the DLPFC dataset.
Model counts data with tpSVG avoids mean-rank relationship.
1 Installation
1.1 GitHub
You can install the development version of tpSVG from GitHub with:
# install.packages("devtools")
devtools::install_github("boyiguo1/tpSVG")1.2 Bioconductor (pending)
The package is currently submitted to Bioconductor for review. Once the package
is accepted by Bioconductor, you can install the latest release
version of tpSVG from Bioconductor via the following code. Additional details
are shown on the Bioconductor page.
if (!require("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("tpSVG")The latest development version can also be installed from the devel version
of Bioconductor or from GitHub following
BiocManager::install(version='devel')2 Modeling spatially resolved gene expression using tpSVG
In the following section, we demonstrate the how to use tpSVG with a SpatialExperiment object.
2.1 Load Packages
To run the demonstration, there a couple of necessary packages to load. We will use the data set in STexampleData, which contains a 10x Visium dataset. Specifically, we will find one 10xVisium sample collected from the dorso lateral prefrontal region of a postmortem human brain from STexampleData package, indexed by brain number of “151673” from Maynard, Collado-Torres et al. (2021). For more information, please see the vignettes of STexampleData
library(tpSVG)
library(SpatialExperiment)
library(STexampleData) # Example data
library(scuttle) # Data preprocess2.2 Data processing
Before running the analysis using tpSVG, we will preprocess the data, which
includes 1) calculating normalization factor, or equivalently library size; 2)
down-sizing the number of genes to 6 to reduce running time. These preprocessing
step may not be necessary if real world data analysis.
spe <- Visium_humanDLPFC()
spe <- spe[, colData(spe)$in_tissue == 1]
spe <- logNormCounts(spe)
# Normalization factor
head(spe$sizeFactor)
#> AAACAAGTATCTCCCA-1 AAACAATCTACTAGCA-1 AAACACCAATAACTGC-1 AAACAGAGCGACTCCT-1
#> 1.8428941 0.3632188 0.8212187 1.1837838
#> AAACAGCTTTCAGAAG-1 AAACAGGGTCTATATT-1
#> 0.9321235 0.8724223
# Equivalently, library size
spe$total <- counts(spe) |> colSums()
# Down-sizing genes for faster computation
idx <- which(
rowData(spe)$gene_name %in% c("MOBP", "PCP4", "SNAP25",
"HBB", "IGKC", "NPY")
)
spe <- spe[idx, ]2.3 Modeling raw counts with Poisson model
The following example demonstrates how to model raw gene counts data with tpSVG.
The model fitting is simple, following stats::glm syntax to follow any
distributional assumption via the argument family. The only key point to
mention here is the model needs to account for any technical variation due to
the gene profiling procedure. To account for such techincal variation, we use
offset term in the model. In the following example, we use the commonly used
library size as the normalizing factor, and hence set offset = log(spe$total)
to account for the techinical variation in the data. Equivalently, it is also
possible/encouraged to use offset = log(spe$sizeFactor), where spe$sizeFactor
is calculated during logNormCounts and is a linear function of the library
size, i.e. spe$total. Note: it is very important to use the log function in
the offset making sure the data scale is conformable.
set.seed(1)
spe_poisson <- tpSVG(
spe,
family = poisson,
assay_name = "counts",
offset = log(spe$total) # Natural log library size
)
rowData(spe_poisson)
#> DataFrame with 6 rows and 6 columns
#> gene_id gene_name feature_type test_stat raw_p
#> <character> <character> <character> <numeric> <numeric>
#> ENSG00000211592 ENSG00000211592 IGKC Gene Expression 3447.12 0
#> ENSG00000168314 ENSG00000168314 MOBP Gene Expression 7143.34 0
#> ENSG00000122585 ENSG00000122585 NPY Gene Expression 3596.18 0
#> ENSG00000244734 ENSG00000244734 HBB Gene Expression 1536.56 0
#> ENSG00000132639 ENSG00000132639 SNAP25 Gene Expression 3529.46 0
#> ENSG00000183036 ENSG00000183036 PCP4 Gene Expression 1381.63 0
#> runtime
#> <numeric>
#> ENSG00000211592 1.224
#> ENSG00000168314 1.258
#> ENSG00000122585 1.297
#> ENSG00000244734 1.240
#> ENSG00000132639 1.289
#> ENSG00000183036 1.1692.4 Gaussian model for log-transformed normalized data
tpSVG provides flexibility regarding the distributional assumption. If interested,
it is possible to model log-transformed count data using Gaussian distribution.
spe_gauss <- tpSVG(
spe,
family = gaussian(),
assay_name = "logcounts",
offset = NULL
)2.5 Covariate-adjusted Model
It is scientifically interesting to understand if and how much additional spatial variation in gene expression when accounting some known biology. For example, if known anatomy is accounted for in the model, is there any additional spatial variation in the gene expression which can be informative to any unknown biology. Statistically, this question is known as covariate adjustment, where the known biology is quantified and accounted for in a model.
To address this question, tpSVG allows introducing covariates in the model via
the argument X, where X takes a vector of any kind, including categorical
variables.
The frist step is to remove any missing data in the dataset, specifically as
the covariate. This can be done via complete.cases.
# Check missing data
idx_complete_case <- complete.cases(spe$ground_truth)
# If multiple covariates
# idx_complete_case <- complete.cases(spe$ground_truth, spe$cell_count)
# Remove missing data
spe <- spe[, idx_complete_case]
# Create a design matrix
x <- spe$ground_truth
spe_poisson_cov <- tpSVG(
spe,
X = x,
family = poisson,
assay_name = "counts",
offset = log(spe$total) # Natural log library size
)2.5.1 image-based SRT in SpatialExperiment (e.g. SpatialFeatureExperiment)
tpSVG can be also used to model image-based SRT data. We use
the seqFISH data from Lohoff and Ghazanfar et al. (2020) to demonstrate tpSVG. Specifically, we use the curated example data in STexampleData package. For more information, please see the vignettes of STexampleData
library(STexampleData)
spe <- seqFISH_mouseEmbryo()
#> see ?STexampleData and browseVignettes('STexampleData') for documentation
#> loading from cache
#> Loading required package: BumpyMatrix
spe
#> class: SpatialExperiment
#> dim: 351 11026
#> metadata(0):
#> assays(2): counts molecules
#> rownames(351): Abcc4 Acp5 ... Zfp57 Zic3
#> rowData names(1): gene_name
#> colnames(11026): embryo1_Pos0_cell10_z2 embryo1_Pos0_cell100_z2 ...
#> embryo1_Pos28_cell97_z2 embryo1_Pos28_cell98_z2
#> colData names(14): cell_id embryo ... segmentation_vertices sample_id
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x y
#> imgData names(0):The example data set contains 351 genes for 11026 genes. To make the
demonstration computationally feasible, we down-size the number of genes to 1.
The average computation times for 11026 cells is roughly 2 minutes.
# Calculate "library size"
spe$total <- counts(spe) |> colSums()
# Down-size genes
idx_gene <- which(
rowData(spe)$gene_name %in%
c("Sox2")
)
library(tpSVG)
# Poisson model
tp_spe <- tpSVG(
input = spe[idx_gene,],
family = poisson(),
offset = log(spe$total),
assay_name = "counts")
rowData(tp_spe)
#> DataFrame with 1 row and 4 columns
#> gene_name test_stat raw_p runtime
#> <character> <numeric> <numeric> <numeric>
#> Sox2 Sox2 9588.51 0 3.6963 Session Info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.4.0 beta (2024-04-15 r86425)
#> os Ubuntu 22.04.4 LTS
#> system x86_64, linux-gnu
#> ui X11
#> language (EN)
#> collate C
#> ctype en_US.UTF-8
#> tz America/New_York
#> date 2024-04-19
#> pandoc 2.7.3 @ /usr/bin/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> abind 1.4-5 2016-07-21 [2] CRAN (R 4.4.0)
#> AnnotationDbi 1.65.2 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> AnnotationHub * 3.11.5 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> beachmat 2.19.4 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> Biobase * 2.63.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> BiocFileCache * 2.11.2 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> BiocGenerics * 0.49.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> BiocManager 1.30.22 2023-08-08 [2] CRAN (R 4.4.0)
#> BiocParallel 1.37.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> BiocStyle * 2.31.0 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> BiocVersion 3.19.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> Biostrings 2.71.5 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> bit 4.0.5 2022-11-15 [2] CRAN (R 4.4.0)
#> bit64 4.0.5 2020-08-30 [2] CRAN (R 4.4.0)
#> blob 1.2.4 2023-03-17 [2] CRAN (R 4.4.0)
#> bookdown 0.39 2024-04-15 [2] CRAN (R 4.4.0)
#> bslib 0.7.0 2024-03-29 [2] CRAN (R 4.4.0)
#> BumpyMatrix * 1.11.0 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> cachem 1.0.8 2023-05-01 [2] CRAN (R 4.4.0)
#> cli 3.6.2 2023-12-11 [2] CRAN (R 4.4.0)
#> codetools 0.2-20 2024-03-31 [3] CRAN (R 4.4.0)
#> crayon 1.5.2 2022-09-29 [2] CRAN (R 4.4.0)
#> curl 5.2.1 2024-03-01 [2] CRAN (R 4.4.0)
#> DBI 1.2.2 2024-02-16 [2] CRAN (R 4.4.0)
#> dbplyr * 2.5.0 2024-03-19 [2] CRAN (R 4.4.0)
#> DelayedArray 0.29.9 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> DelayedMatrixStats 1.25.3 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> digest 0.6.35 2024-03-11 [2] CRAN (R 4.4.0)
#> dplyr 1.1.4 2023-11-17 [2] CRAN (R 4.4.0)
#> evaluate 0.23 2023-11-01 [2] CRAN (R 4.4.0)
#> ExperimentHub * 2.11.3 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> fansi 1.0.6 2023-12-08 [2] CRAN (R 4.4.0)
#> fastmap 1.1.1 2023-02-24 [2] CRAN (R 4.4.0)
#> filelock 1.0.3 2023-12-11 [2] CRAN (R 4.4.0)
#> generics 0.1.3 2022-07-05 [2] CRAN (R 4.4.0)
#> GenomeInfoDb * 1.39.14 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> GenomeInfoDbData 1.2.12 2024-04-16 [2] Bioconductor
#> GenomicRanges * 1.55.4 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> glue 1.7.0 2024-01-09 [2] CRAN (R 4.4.0)
#> htmltools 0.5.8.1 2024-04-04 [2] CRAN (R 4.4.0)
#> httr 1.4.7 2023-08-15 [2] CRAN (R 4.4.0)
#> IRanges * 2.37.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> jquerylib 0.1.4 2021-04-26 [2] CRAN (R 4.4.0)
#> jsonlite 1.8.8 2023-12-04 [2] CRAN (R 4.4.0)
#> KEGGREST 1.43.0 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> knitr 1.46 2024-04-06 [2] CRAN (R 4.4.0)
#> lattice 0.22-6 2024-03-20 [3] CRAN (R 4.4.0)
#> lifecycle 1.0.4 2023-11-07 [2] CRAN (R 4.4.0)
#> magick 2.8.3 2024-02-18 [2] CRAN (R 4.4.0)
#> magrittr 2.0.3 2022-03-30 [2] CRAN (R 4.4.0)
#> Matrix 1.7-0 2024-03-22 [3] CRAN (R 4.4.0)
#> MatrixGenerics * 1.15.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> matrixStats * 1.3.0 2024-04-11 [2] CRAN (R 4.4.0)
#> memoise 2.0.1 2021-11-26 [2] CRAN (R 4.4.0)
#> mgcv * 1.9-1 2023-12-21 [3] CRAN (R 4.4.0)
#> mime 0.12 2021-09-28 [2] CRAN (R 4.4.0)
#> nlme * 3.1-164 2023-11-27 [3] CRAN (R 4.4.0)
#> pillar 1.9.0 2023-03-22 [2] CRAN (R 4.4.0)
#> pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.4.0)
#> png 0.1-8 2022-11-29 [2] CRAN (R 4.4.0)
#> purrr 1.0.2 2023-08-10 [2] CRAN (R 4.4.0)
#> R6 2.5.1 2021-08-19 [2] CRAN (R 4.4.0)
#> rappdirs 0.3.3 2021-01-31 [2] CRAN (R 4.4.0)
#> Rcpp 1.0.12 2024-01-09 [2] CRAN (R 4.4.0)
#> rjson 0.2.21 2022-01-09 [2] CRAN (R 4.4.0)
#> rlang 1.1.3 2024-01-10 [2] CRAN (R 4.4.0)
#> rmarkdown 2.26 2024-03-05 [2] CRAN (R 4.4.0)
#> RSQLite 2.3.6 2024-03-31 [2] CRAN (R 4.4.0)
#> S4Arrays 1.3.7 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> S4Vectors * 0.41.6 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> sass 0.4.9 2024-03-15 [2] CRAN (R 4.4.0)
#> scuttle * 1.13.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> sessioninfo 1.2.2 2021-12-06 [2] CRAN (R 4.4.0)
#> SingleCellExperiment * 1.25.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> SparseArray 1.3.5 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> sparseMatrixStats 1.15.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> SpatialExperiment * 1.13.2 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> STexampleData * 1.11.3 2024-04-18 [2] Bioconductor 3.19 (R 4.4.0)
#> SummarizedExperiment * 1.33.3 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> tibble 3.2.1 2023-03-20 [2] CRAN (R 4.4.0)
#> tidyselect 1.2.1 2024-03-11 [2] CRAN (R 4.4.0)
#> tpSVG * 0.99.7 2024-04-19 [1] Bioconductor 3.19 (R 4.4.0)
#> UCSC.utils 0.99.7 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> utf8 1.2.4 2023-10-22 [2] CRAN (R 4.4.0)
#> vctrs 0.6.5 2023-12-01 [2] CRAN (R 4.4.0)
#> withr 3.0.0 2024-01-16 [2] CRAN (R 4.4.0)
#> xfun 0.43 2024-03-25 [2] CRAN (R 4.4.0)
#> XVector 0.43.1 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#> yaml 2.3.8 2023-12-11 [2] CRAN (R 4.4.0)
#> zlibbioc 1.49.3 2024-04-19 [2] Bioconductor 3.19 (R 4.4.0)
#>
#> [1] /tmp/RtmplPMBQH/Rinst7ca2b71b210b0
#> [2] /home/biocbuild/bbs-3.19-bioc/R/site-library
#> [3] /home/biocbuild/bbs-3.19-bioc/R/library
#>
#> ──────────────────────────────────────────────────────────────────────────────