Recovering Spatial Ecotypes from Single-Cell Spatial Transcriptomics Data
Source:vignettes/Recovery_scST.Rmd
Recovery_scST.RmdOverview
In this tutorial, we will illustrate how to recover spatial ecotypes (SEs) from single-cell-scale spatial transcriptomics (ST) data, generated by platforms such as MERSCOPE, Xenium, CosMx, and Visium HD.
First load required packages for this vignette
Data preparation
SE recovery for single-cell ST data requires two inputs:
- Gene expression matrix: a numeric matrix with genes as rows and cells as columns.
- Metadata: a data frame containing at least three columns: “X” (x-coordinate), “Y” (y-coordinate), and “CellType” (cell type annotation). Row names must match the column names (cell IDs) of the expression matrix.
Text files as input
# Load metadata
scmeta <- read.table("https://spatialecotyper.stanford.edu/inc/inc.public.vignettes.php?file=CRC2_subset_scmeta.tsv",
sep = "\t", header = TRUE, row.names = 1)
head(scmeta[, c("X", "Y", "CellType")])## X Y CellType
## HumanColonCancerPatient2__cell_1085 5057.912 -1721.871 Macrophage
## HumanColonCancerPatient2__cell_1205 5035.339 -1783.352 Macrophage
## HumanColonCancerPatient2__cell_1312 5012.093 -1908.390 Macrophage
## HumanColonCancerPatient2__cell_1625 5073.743 -1861.740 Macrophage
## HumanColonCancerPatient2__cell_1629 5030.235 -1867.441 Macrophage
## HumanColonCancerPatient2__cell_1639 5020.448 -1882.751 Macrophage
# Load expression matrix
scdata <- fread("https://spatialecotyper.stanford.edu/inc/inc.public.vignettes.php?file=CRC2_subset_counts.tsv.gz",
sep = "\t",header = TRUE, data.table = FALSE)
rownames(scdata) <- scdata[, 1]
scdata <- as.matrix(scdata[, -1])
head(scdata[,1:5])## HumanColonCancerPatient2__cell_1085 HumanColonCancerPatient2__cell_1205
## PDK4 0 0
## CCL26 0 0
## CX3CL1 0 0
## PGLYRP1 0 0
## CD4 0 0
## SNAI2 0 0
## HumanColonCancerPatient2__cell_1312 HumanColonCancerPatient2__cell_1625
## PDK4 1 0
## CCL26 0 0
## CX3CL1 0 0
## PGLYRP1 0 0
## CD4 0 0
## SNAI2 0 0
## HumanColonCancerPatient2__cell_1629
## PDK4 0
## CCL26 0
## CX3CL1 0
## PGLYRP1 0
## CD4 0
## SNAI2 0Data normalization
The data should be normalized to a log-transformed CPM (counts per
million) or a similar scale. The NormalizeData function is
applied at this step to perform this normalization.
normdata = NormalizeData(scdata)SE recovery
The RecoverSE function will be used to assign single cells into SEs.
Using pan-cancer SE recovery model
For recovery of the nine pan-cancer spatial ecotypes defined in our study,
the pre-trained model file (SE_recovery_model.rds) can be
downloaded from Spatial
EcoTyper web platform. Alternatively, you can upload your
single-cell expression data directly to the web platform for the
recovery.
The model was trained on discovery MERSCOPE data, encompassing five cancer types: melanoma, and four carcinomas, and nine distinct cell types: B cells, CD4+ T cells, CD8+ T cells, NK cells, plasma cells, macrophages, dendritic cells, fibroblasts, and endothelial cells.
Therefore, cells in the query data should be grouped into one of “B”, “CD4T”, “CD8T”, “NK”, “Plasma”, “Macrophage”, “DC”, “Fibroblast”, and “Endothelial”, case sensitive. All the other cell types will be considered non-SE compartments.
model_Ws = readRDS("SE_recovery_model.rds")
sepreds <- RecoverSE(normdata, celltypes = scmeta$CellType, Ws = model_Ws,
custom = FALSE, min.score = 0.6, ncores = 8)
head(sepreds)## CID
## HumanColonCancerPatient2__cell_1085 HumanColonCancerPatient2__cell_1085
## HumanColonCancerPatient2__cell_1205 HumanColonCancerPatient2__cell_1205
## HumanColonCancerPatient2__cell_1312 HumanColonCancerPatient2__cell_1312
## HumanColonCancerPatient2__cell_1625 HumanColonCancerPatient2__cell_1625
## HumanColonCancerPatient2__cell_1629 HumanColonCancerPatient2__cell_1629
## HumanColonCancerPatient2__cell_1639 HumanColonCancerPatient2__cell_1639
## CellType InitSE SE PredScore
## HumanColonCancerPatient2__cell_1085 Macrophage SE11 SE9 0.7484101
## HumanColonCancerPatient2__cell_1205 Macrophage SE09 NonSE 0.8450224
## HumanColonCancerPatient2__cell_1312 Macrophage SE11 SE9 0.7659342
## HumanColonCancerPatient2__cell_1625 Macrophage SE07 SE6 0.6217793
## HumanColonCancerPatient2__cell_1629 Macrophage SE08 NonSE 0.5823970
## HumanColonCancerPatient2__cell_1639 Macrophage SE11 SE9 0.8723833The SE recovery output contains five columns: cell ID (CID), cell type (CellType), initial SE assignment (InitSE), final SE assignment (SE), and prediction score (PredScore).
Using custom model
To use custom model, users should first develop a model following the tutorial NMF Model Development for Spatial Ecotype Recovery from Single-Cell and Spatial Transcriptomics Data. The resulting model can be used for SE recovery.
Ws <- readRDS("path/to/SE_recovery_model.rds")
sepreds <- RecoverSE(normdata, celltypes = scmeta$CellType,
Ws = Ws, custom = TRUE, min.score = 0.5)After prediction, we recommend retaining only assignments corresponding to SE-specific cell states identified in Tutorial 4. Cells not assigned to these SE-specific states should be grouped as “NonSE,” as they cannot be robustly mapped to the corresponding SE.
Combine SE recovery results into metadata
Visualization of SE recovery results in the tissue
To visualize the spatial SE landscape, we recommend computing SE abundances within spatial neighborhoods and displaying the continuum of SE levels across the tissue, rather than directly plotting SE assignments. Direct visualization of SE assignments may bias interpretation due to variation in cell density and the presence of a large number of non-SE cells, including cells with low detection of SE marker genes and cells belonging to cell states that are not specific to any individual SE.
# Compute SE abundances by SNs
seabunds = ComputeSEAbundanceBySN(scmeta)
colors = rev(RColorBrewer::brewer.pal(11, "Spectral"))
plist = lapply(colnames(seabunds)[-(1:3)], function(se){
p = SpatialView(seabunds, by = se, X = "X", Y = "Y", pt.size = 0.5) +
coord_fixed() + scale_color_gradientn(colors = colors) +
theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) +
labs(title = se)
p
})
p = patchwork::wrap_plots(plist, nrow = 3)
plot(p)
Validation: spatial colocalization of SE cell states
After recovering SE-specific cell states from single-cell ST data, users can validate the SEs by testing whether cell states belonging to the same SE are more spatially colocalized than expected by random chance using the Colocalization function. The spatial coordinates (column X and Y), SE and cell type (CellType) are required for the analysis. The function also evaluates the significance of cell state colocalization within each SE and returns a vector of p-values.
# Verify the four columns required for the analysis.
scmeta$CellState = paste0(scmeta$SE, "_", scmeta$CellType)
head(scmeta[, c("X", "Y", "CellType", "CellState")])## X Y CellType
## HumanColonCancerPatient2__cell_1085 5057.912 -1721.871 Macrophage
## HumanColonCancerPatient2__cell_1205 5035.339 -1783.352 Macrophage
## HumanColonCancerPatient2__cell_1312 5012.093 -1908.390 Macrophage
## HumanColonCancerPatient2__cell_1625 5073.743 -1861.740 Macrophage
## HumanColonCancerPatient2__cell_1629 5030.235 -1867.441 Macrophage
## HumanColonCancerPatient2__cell_1639 5020.448 -1882.751 Macrophage
## CellState
## HumanColonCancerPatient2__cell_1085 SE9_Macrophage
## HumanColonCancerPatient2__cell_1205 NonSE_Macrophage
## HumanColonCancerPatient2__cell_1312 SE9_Macrophage
## HumanColonCancerPatient2__cell_1625 SE6_Macrophage
## HumanColonCancerPatient2__cell_1629 NonSE_Macrophage
## HumanColonCancerPatient2__cell_1639 SE9_Macrophage
## For quick test run, 100 permutations were performed to compute the colocalization scores.
coloc_res = Colocalization(scmeta, coords = c("X", "Y"),
CellType = "CellType",
CellState = "CellState",
nperm = 100, ## 1000 or higher is recommended.
test = TRUE,
ncores = 8)
head(coloc_res$ColocIndex[, 1:5]) ## Matrix of colocalization indices between all cell states## NonSE_B NonSE_CD4T NonSE_CD8T NonSE_DC NonSE_Endothelial
## NonSE_B 3.0273914 -0.3871890 -2.4146753 -2.7740600 -1.7139008
## NonSE_CD4T -0.2230683 1.4884956 -0.5378198 -3.1516158 1.8466398
## NonSE_CD8T -2.1822585 -0.5866724 1.2915895 -3.4119072 -0.5020747
## NonSE_DC -2.5797764 -2.6971647 -3.2207830 -1.1411956 1.3444704
## NonSE_Endothelial -1.6912665 1.7783736 -0.2552831 0.7904064 1.1589863
## NonSE_Fibroblast -8.8121407 3.6152242 -2.1043831 3.9147400 3.7981699
coloc_res$Pval ## Vector of p-values for colocalization significance## NonSE SE1 SE2 SE3 SE4 SE5
## 1.433468e-01 8.368894e-07 1.385676e-02 1.068561e-03 1.498422e-13 1.887150e-02
## SE6 SE7 SE8 SE9
## 3.670351e-04 7.344488e-05 1.462465e-05 5.689382e-02
## attr(,"Zscore")
## NonSE SE1 SE2 SE3 SE4 SE5 SE6 SE7
## 1.463441 4.926562 2.460955 3.271824 7.387298 2.348060 3.562718 3.964841
## SE8 SE9
## 4.334267 1.904126The ColocalizationMetaAnalysis function can be used to perform meta-analysis of cell-state colocalization results across samples. It applies Stouffer’s method to combine colocalization indices across ST samples and to aggregate SE-specific z-scores, generating an integrated colocalization index matrix and a meta P value for each SE.
# Generate a second colocalization result for demonstration of integrative analysis
coloc_res2 = Colocalization(scmeta, nperm = 100, ncores = 8)
coloc_res_list = list(Sample1 = coloc_res, Sample2 = coloc_res2)
# Combine colocalization results from multiple samples using meta-analysis
meta_coloc_res = ColocalizationMetaAnalysis(coloc_res_list)
# Integrated matrix of cell-state colocalization indices
head(meta_coloc_res$MetaColocIndex[, 1:5])## NonSE_B NonSE_CD4T NonSE_CD8T NonSE_DC NonSE_Endothelial
## NonSE_B 4.4622035 -0.3652733 -3.1267282 -4.017636 -2.3871742
## NonSE_CD4T -0.1283318 1.8418853 -0.8880536 -4.478315 2.5915767
## NonSE_CD8T -2.8516299 -0.9126206 1.9696086 -4.691937 -0.6623149
## NonSE_DC -3.7810804 -3.7719211 -4.2842524 -1.601143 1.9710879
## NonSE_Endothelial -2.3630828 2.4818032 -0.3262626 1.268528 1.3948120
## NonSE_Fibroblast -7.0710673 5.3684188 -3.2468506 5.626446 5.2800988
# Meta-analysis P values
meta_coloc_res$MetaPval## NonSE SE1 SE2 SE3 SE4 SE5
## 4.386579e-02 7.013388e-10 1.156336e-05 9.443509e-06 3.080129e-20 8.225041e-04
## SE6 SE7 SE8 SE9
## 2.592443e-06 6.409643e-09 4.933816e-09 5.817645e-03You can visualize the results using the CooccurrenceHeatmapView function.
## Visualization of single-sample result
idx = grepl("NonSE", rownames(coloc_res$ColocIndex))
p = CooccurrenceHeatmapView(coloc_res$ColocIndex[!idx, !idx],
coloc_res$Pval)
p = ComplexHeatmap::draw(p)
ses = gsub("_.*", "", rownames(coloc_res$ColocIndex)[!idx])
drawRectangleAnnotation(p, ses, ses)
## Visualization of integrated result
idx = grepl("NonSE", rownames(meta_coloc_res$MetaColocIndex))
p = CooccurrenceHeatmapView(meta_coloc_res$MetaColocIndex[!idx, !idx],
meta_coloc_res$MetaPval)
p = ComplexHeatmap::draw(p)
ses = gsub("_.*", "", rownames(meta_coloc_res$MetaColocIndex)[!idx])
drawRectangleAnnotation(p, ses, ses)
Validation: spatial autocorrelation
Users can also validate the SEs by assessing the spatial autocorrelation with Moran’s I. To control for bias, permutation experiments will be performed to normalize Moran’s I into z-scores using the ComputeNormalizedMoranI function.
## For quick test run, test with 10 permutations were performed.
morani_zscore = ComputeNormalizedMoranI(scmeta, coords = c("X", "Y"),
SE = "SE", CellType = "CellType",
nperm = 10, ncores = 8)
morani_zscore## NonSE SE1 SE2 SE3 SE4 SE5 SE6 SE7
## 14.70323 105.02635 40.29611 12.14710 141.02512 12.12789 13.99375 12.05357
## SE8 SE9
## 50.89906 11.39071Validation: average expression of SE markers
Users can also validate the SEs by visualizing the average expression of SE consensus markers across recovered SEs using the AverageMarkerExpression function.
## The function is built on Seurat functions, so first create a Seurat object.
obj = CreateSeuratObject(scdata, meta.data = scmeta)
obj = NormalizeData(obj)By default, the AverageMarkerExpression function compute the average expression of SE consensus markers from Zhang et al., Nature 2026 paper.
## Visualize average expression of SE consensus markers across recovered SEs
res = AverageMarkerExpression(obj, group.by = "SE")
plot(res$p)
# Underlying data
res$AvgExp## SE1 SE2 SE3 SE4 SE5 SE6
## SE1 1.64070080 -0.7060161 -0.4663656 0.8807441 -1.0590788 1.21249206
## SE2 -0.84354729 1.7105596 -0.6344529 -1.6443932 -0.1966020 0.97375468
## SE3 1.25589643 -1.2176720 1.1718960 -0.5889244 -0.8418603 -0.84898871
## SE4 0.05998728 -2.1657621 0.3102000 1.6572110 -0.0312405 0.22618885
## SE5 -0.04530647 1.3943294 -0.2890947 -2.3422362 0.1779251 0.28174520
## SE6 0.01167333 0.3181417 -0.8555862 0.8529600 -0.8583647 1.64279726
## SE7 0.55327809 -0.1473391 0.6339091 -2.0968811 -0.7236552 -0.03412487
## SE8 -1.22614497 0.8803489 -0.6618139 -0.6590871 -0.1434522 0.07370003
## SE9 -1.01769519 0.9995106 -0.3650149 -0.6427398 -0.9947289 -0.38007387
## SE7 SE8 SE9
## SE1 -0.06865693 -0.4070571 -1.0267623
## SE2 -0.07751705 0.4245898 0.2876084
## SE3 0.73128916 -0.6048618 0.9432257
## SE4 -0.36954015 0.4923614 -0.1794057
## SE5 -0.03893458 0.4159551 0.4456171
## SE6 0.03605022 0.5086691 -1.6563407
## SE7 1.37233632 0.6498305 -0.2073538
## SE8 0.56940084 1.9651776 -0.7981293
## SE9 -0.41772201 1.1459626 1.6725014The AverageMarkerExpression
function can also be used to visualize the average expression of
user-defined gene sets by specifying the genesets parameter
as a list of gene lists.
genesets = list(set1 = c("FOS", "EGR1", "DUSP1", "JUNB", "JUN", "FOSB"),
set2 = c("SPP1", "CCL2", "MMP9"),
set3 = c("CXCL12", "TGFBR2", "PDK4", "HLA-DMA"),
set4 = c("ACTA2", "CAV1", "CLDN5"),
set5 = c("LRP1", "MAFB"),
set6 = c("CXCL2", "NFKBIA", "CDKN1A"),
set7 = c("STAT1", "TAP1", "HLA-A", "HLA-B", "HLA-C", "B2M"),
set8 = c("PKM", "PCNA", "BST2"),
set9 = c("NRP1", "ITGB1", "CD276"))
res = AverageMarkerExpression(obj, group.by = "SE", genesets = genesets)
plot(res$p)
Session info
The session info allows users to replicate the exact environment and identify potential discrepancies in package versions or configurations that might be causing problems.
## R version 4.4.1 (2024-06-14)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.6.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] SpatialEcoTyper_1.0.5 pals_1.9 NMF_0.28
## [4] Biobase_2.64.0 BiocGenerics_0.50.0 cluster_2.1.6
## [7] rngtools_1.5.2 registry_0.5-1 RANN_2.6.2
## [10] Matrix_1.7-0 data.table_1.18.4 Seurat_5.1.0
## [13] SeuratObject_5.0.2 sp_2.1-4 ggplot2_4.0.3
## [16] dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] RcppAnnoy_0.0.22 splines_4.4.1 later_1.3.2
## [4] R.oo_1.26.0 tibble_3.3.1 polyclip_1.10-7
## [7] fastDummies_1.7.4 lifecycle_1.0.5 sf_1.1-0
## [10] doParallel_1.0.17 globals_0.16.3 lattice_0.22-6
## [13] MASS_7.3-60.2 magrittr_2.0.5 plotly_4.10.4
## [16] sass_0.4.9 rmarkdown_2.28 jquerylib_0.1.4
## [19] yaml_2.3.10 httpuv_1.6.15 sctransform_0.4.1
## [22] spam_2.10-0 spatstat.sparse_3.1-0 reticulate_1.39.0
## [25] cowplot_1.1.3 mapproj_1.2.11 pbapply_1.7-2
## [28] DBI_1.3.0 RColorBrewer_1.1-3 maps_3.4.2
## [31] abind_1.4-8 Rtsne_0.17 R.utils_2.12.3
## [34] purrr_1.0.2 circlize_0.4.18 IRanges_2.38.1
## [37] S4Vectors_0.42.1 ggrepel_0.9.6 irlba_2.3.5.1
## [40] listenv_0.9.1 spatstat.utils_3.1-0 units_1.0-1
## [43] goftest_1.2-3 RSpectra_0.16-2 spatstat.random_3.3-1
## [46] fitdistrplus_1.2-1 parallelly_1.38.0 pkgdown_2.1.0
## [49] leiden_0.4.3.1 codetools_0.2-20 tidyselect_1.2.1
## [52] shape_1.4.6.1 farver_2.1.2 matrixStats_1.5.0
## [55] stats4_4.4.1 spatstat.explore_3.3-2 jsonlite_1.8.8
## [58] GetoptLong_1.1.1 e1071_1.7-16 progressr_0.14.0
## [61] ggridges_0.5.6 survival_3.6-4 iterators_1.0.14
## [64] systemfonts_1.1.0 foreach_1.5.2 tools_4.4.1
## [67] ragg_1.3.2 ica_1.0-3 Rcpp_1.1.2
## [70] glue_1.8.1 gridExtra_2.3.1 xfun_0.52
## [73] withr_3.0.3 BiocManager_1.30.25 fastmap_1.2.0
## [76] boot_1.3-30 spData_2.3.4 digest_0.6.39
## [79] R6_2.6.1 mime_0.12 wk_0.9.5
## [82] textshaping_0.4.0 colorspace_2.1-2 scattermore_1.2
## [85] tensor_1.5 dichromat_2.0-0.1 spatstat.data_3.1-2
## [88] R.methodsS3_1.8.2 tidyr_1.3.1 generics_0.1.4
## [91] class_7.3-22 httr_1.4.7 htmlwidgets_1.6.4
## [94] spdep_1.4-2 uwot_0.2.2 pkgconfig_2.0.3
## [97] gtable_0.3.6 ComplexHeatmap_2.20.0 lmtest_0.9-40
## [100] S7_0.2.2 htmltools_0.5.8.1 dotCall64_1.1-1
## [103] clue_0.3-68 scales_1.4.0 png_0.1-9
## [106] spatstat.univar_3.0-1 knitr_1.48 rstudioapi_0.16.0
## [109] reshape2_1.4.4 rjson_0.2.23 nlme_3.1-164
## [112] proxy_0.4-27 cachem_1.1.0 zoo_1.8-12
## [115] GlobalOptions_0.1.4 stringr_1.5.1 KernSmooth_2.23-24
## [118] miniUI_0.1.1.1 s2_1.1.9 desc_1.4.3
## [121] pillar_1.11.1 grid_4.4.1 vctrs_0.7.3
## [124] promises_1.3.0 xtable_1.8-4 evaluate_0.24.0
## [127] magick_2.8.5 cli_3.6.6 compiler_4.4.1
## [130] rlang_1.3.0 crayon_1.5.3 future.apply_1.11.2
## [133] labeling_0.4.3 classInt_0.4-11 plyr_1.8.9
## [136] fs_1.6.4 stringi_1.8.4 viridisLite_0.4.3
## [139] deldir_2.0-4 gridBase_0.4-7 lazyeval_0.2.2
## [142] spatstat.geom_3.3-2 RcppHNSW_0.6.0 patchwork_1.2.0
## [145] future_1.34.0 shiny_1.9.1 highr_0.11
## [148] ROCR_1.0-11 igraph_2.0.3 bslib_0.8.0