空间转录组+单细胞联合分析---（SPOTlight最新版教程）
作者在2022年更新了SPOTlight脚本代码，或许有伙伴在运行以往的版本教程（https://marcelosua.github.io/）spotlight_deconvolution函数的时候会出现报错：could not find function "spotlight_deconvolution"，这就说明版本更新了，以往的函数代码作者也已经更新。
什么是SPOTlight？
由于10x Visium数据初衷是想将每个spot看作一个细胞，但是现实实验中单个spot中包含不仅是一个细胞。如何确定每个spot中包含的细胞，有利于空间转录组的分析。
SPOTlight 是一个工具，能够解读每个捕获位置细胞混合物内存在的细胞类型和比例，最初是为10X的Visium-空间转录组学技术而开发。
SPOTlight可以结合单细胞转录组RNA测序信息反卷积deconvolute空间数据，识别每个spot中的细胞类型和比例。SPOTlight利用这两种数据类型的优势，能够将ST与scRNA-seq数据集成，从而推断出复杂组织中细胞类型和状态的位置。SPOTlight基于一个种子的非负矩阵因子分解回归（Seeded NMF regression ），使用细胞类型标记基因和非负最小二乘(NNLS)初始化，随后去卷积ST捕获位置(spot)。作者文章：Elosua-Bayes M, Nieto P, Mereu E, Gut I, Heyn H (2021): SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes. Nucleic Acids Res 49(9):e50. doi:10.1093/nar/gkab043.
研究流程图如下：
image.png
研究代码如下：
首先安装SPOTlight包
install.packages("BiocManager")
BiocManager::install("SPOTlight")
# Or the devel version
BiocManager::install("SPOTlight", version = "devel")
#或者从Github中安装
install.packages("devtools")
library(devtools)
install_github("https://github.com/MarcElosua/SPOTlight")
加载包开始分析
library(ggplot2)
library(SPOTlight)
library(SingleCellExperiment)
library(SpatialExperiment)
library(scater)
library(scran)
library(scuttle)
library(Matrix)
library(data.table)
library(Seurat)
library(SeuratData)
library(dplyr)
library(gt)
library(igraph)
library(RColorBrewer)
## Loading the spatial and single cell data:
spe <- readRDS("path:/scRNA")
sce <- readRDS("path:/scRNA.rds") ##单细胞数据经scater处理的,对接logNormCounts
##代码空转和单细胞演示数据是作者从Biocondcutor packages下载
library(TENxVisiumData) ##空转数据
spe <- MouseKidneyCoronal()
# Use symbols instead of Ensembl IDs as feature names
rownames(spe) <- rowData(spe)$symbol
##单细胞数据集来自：[Tabula Muris Sensis](https://www.nature.com/articles/s41586-020-2496-1)
library(TabulaMurisSenisData) ##单细胞数据
sce <- TabulaMurisSenisDroplet(tissues = "Kidney")$Kidney
##快速查看数据:
table(sce$free_annotation, sce$age)
# 作者为了去除批次影响挑选特定年龄段的小鼠细胞
sce <- sce[, sce$age == "18m"]
# 挑选具有明确注释的细胞类别,这个可以根据自身情况调整
sce <- sce[, !sce$free_annotation %in% c("nan", "CD45")]
## Preprocessing
### Feature selection(第一步获得每种细胞类型的标记基因,归一化处理)
sce <- logNormCounts(sce)
##获得既不是核糖体也不是线粒体的基因向量
genes <- !grepl(pattern = "^Rp[l|s]|Mt", x = rownames(sce))
dec <- modelGeneVar(sce, subset.row = genes)
plot(dec$mean, dec$total, xlab = "Mean log-expression", ylab = "Variance")
curve(metadata(dec)$trend(x), col = "blue", add = TRUE)
# 获得3000个的高变基因HVGs.
hvg <- getTopHVGs(dec, n = 3000)
colLabels(sce) <- colData(sce)$celltype
# 计算标记基因，表明该基因对该细胞类型的重要性
mgs <- scoreMarkers(sce, subset.row = genes)
mgs_fil <- lapply(names(mgs), function(i) {
x <- mgs[[i]]
# 筛选并保留相关标记基因，即AUC>0.8的基因
x <- x[x$mean.AUC > 0.8, ]
# 将基因从高到低的权重排序
x <- x[order(x$mean.AUC, decreasing = TRUE), ]
# 将基因和聚类ID添加到数据框中
x$gene <- rownames(x)
x$cluster <- i
data.frame(x)
})
mgs_df <- do.call(rbind, mgs_fil)
##Cell Downsampling(每个细胞类别最多随机选择100个细胞,<100个细胞都将被使用。转录谱明显不同的(如B和T),可使用较少N的细胞;转录谱相似的细胞增加N)
# 按细胞类别分割细胞索引
idx <- split(seq(ncol(sce)), sce$celltype)
# 每个类别和亚群最多降样到20个
# 我们在这里使用5个来加快进程，但对于你的实际情况，可以设置为75-100个。
# 生活模式分析
n_cells <- 5
cs_keep <- lapply(idx, function(i) {
n <- length(i)
if (n < n_cells)
n_cells <- n
sample(i, n_cells)
})
sce <- sce[, unlist(cs_keep)]
## Deconvolution(运行 "SPOTlight 函数"来解构spot)
res <- SPOTlight(
x = sce,
y = spe,
groups = as.character(sce$free_annotation),
mgs = mgs_df,
hvg = hvg,
weight_id = "mean.AUC",
group_id = "cluster",
gene_id = "gene")
##或者分两步运行`SPOTlight'得到训练好的模型,在其他数据集上可重复使用。
mod_ls <- trainNMF(
x = sce,
y = spe,
groups = sce$type,
mgs = mgs,
weight_id = "weight",
group_id = "type",
gene_id = "gene")
# Run deconvolution(运行去卷积)
res <- runDeconvolution(
x = spe,
mod = mod_ls[["mod"]],
ref = mod_ls[["topic"]])
##从`SPOTlight`中提取数据
# 提取去卷积矩阵
head(mat <- res$mat)[, seq_len(3)]
# 提取NMF模型拟合
mod <- res$NMF
## 可视化过程
## Topic profiles(设置`facet = FALSE`,topic profile表征细胞类别,每个细胞类别都有一个独特的主题特征)
plotTopicProfiles(
x = mod,
y = sce$free_annotation,
facet = FALSE,
min_prop = 0.01,
ncol = 1) +
theme(aspect.ratio = 1)
##确保所有来自同一细胞类别的细胞都有类似的topic profile
plotTopicProfiles(
x = mod,
y = sce$free_annotation,
facet = TRUE,
min_prop = 0.01,
ncol = 6)
##查看模型为每个主题学习了哪些基因,更高的数值表明该基因与该主题更相关
library(NMF)
sign <- basis(mod)
colnames(sign) <- paste0("Topic", seq_len(ncol(sign)))
head(sign)
# 这可以用DT动态可视化，如下所示
# DT::datatable(sign, fillContainer = TRUE, filter = "top")
## Spatial Correlation Matrix(空间相关矩阵)
#参照(http://www.sthda.com/english/wiki/ggcorrplot-visualization-of-a-correlation-matrix-using-ggplot2)获得参数
plotCorrelationMatrix(mat)
## Co-localization(空间共定位及相互作用)
#参照(https://www.r-graph-gallery.com/network.html) 获取参数.
plotInteractions(mat, which = "heatmap", metric = "prop")
plotInteractions(mat, which = "heatmap", metric = "jaccard")
plotInteractions(mat, which = "network")
## Scatterpie(将细胞类型的比例可视化为饼状图)
ct <- colnames(mat)
mat[mat < 0.1] <- 0
# 定义调色板
# 使用'colorBlindness'软件包中的'paletteMartin'
paletteMartin <- c(
"#000000", "#004949", "#009292", "#ff6db6", "#ffb6db",
"#490092", "#006ddb", "#b66dff", "#6db6ff", "#b6dbff",
"#920000", "#924900", "#db6d00", "#24ff24", "#ffff6d")
pal <- colorRampPalette(paletteMartin)(length(ct))
names(pal) <- ct
plotSpatialScatterpie(
x = spe,
y = mat,
cell_types = colnames(mat),
img = FALSE,
scatterpie_alpha = 1,
pie_scale = 0.4) +
scale_fill_manual(
values = pal,
breaks = names(pal))
##在图像下方--将其逆时针旋转90度,并在横轴上做镜像,以显示如果spot不在图像上覆盖该如何对齐
plotSpatialScatterpie(
x = spe,
y = mat,
cell_types = colnames(mat),
img = FALSE,
scatterpie_alpha = 1,
pie_scale = 0.4,
# 将图像逆时针旋转90度
degrees = -90,
# 将图像在其X轴上进行旋转
axis = "h") +
scale_fill_manual(
values = pal,
breaks = names(pal))
## Residuals
##最后观察模型对每个点的比例的预测情况,通过观察每个斑点的平方之和的残差来实现,这表明模型不能解释的生物信号的数量
spe$res_ss <- res[[2]][colnames(spe)]
xy <- spatialCoords(spe)
spe$x <- xy[, 1]
spe$y <- xy[, 2]
ggcells(spe, aes(x, y, color = res_ss)) +
geom_point() +
scale_color_viridis_c() +
coord_fixed() +
theme_bw()
# Session information(查看目前相关包的版本信息)
sessionInfo()
