基于Seurat的空转单样本数据分析流程学习(二)-SPOTlight

原创

凑齐六个字吧

发布于 2025-09-12 12:21:34

9400

代码可运行

文章被收录于专栏：单细胞单细胞

运行总次数：0

代码可运行

上一篇推文采用了RCTD工具，本次我们使用SPOTlight映射，SPOTlight 是一种用于解析（去卷积）细胞类型及其在每个捕获位置中所占比例的工具，这些位置通常由多种细胞混合组成。该方法最初是为 10X Visium 空间转录组学技术开发的，但它同样适用于所有能够返回细胞混合信号的技术。SPOTlight 的核心原理是：通过 NMFreg 模型学习每种细胞类型的主题特征（topic profile signatures），然后寻找能够最佳拟合目标点（spot）的细胞类型组合。

分析流程

aligned_fiducials.jpg：带有芯片定位点（fiducial markers）的图像。这些定位点是Visium芯片上的参照标记，用来帮助显微镜图像与芯片的spot array对齐；在图像配准过程中保证组织切片位置与空间条码坐标准确匹配。
detected_tissue_image.jpg：组织识别后的二值图像。用于标记哪些区域包含组织（组织蒙版 mask），帮助区分组织区域和背景区域；常用于筛选有效的 spots。
filtered_feature_bc_matrix.h5：核心的基因表达矩阵文件（HDF5 格式）。行表示基因，列表示空间位置（spot/barcode），数值为UMI counts；该矩阵已去除低质量基因和低质量spot，是Seurat、Scanpy等空间分析的主要输入。
scalefactors_json.json：存放图像与空间坐标之间的缩放因子。主要字段包括高/低分辨率图像的缩放比例和spot 在原始图像上的直径；用于将spot坐标准确映射到组织切片图像上，保证可视化时位置匹配。
tissue_hires_image.png：高分辨率 H&E 染色组织切片图像。用于在空间可视化中叠加基因表达信息，生成高质量的展示图像。
tissue_lowres_image.png：低分辨率 H&E 染色组织切片图像。常用于快速加载和浏览，在交互式分析中默认使用，节省内存和计算资源。
tissue_positions_list.csv：每个 spot 的空间位置信息文件。包括 spot 的条码、是否位于组织区域（1=组织内，0=背景）、芯片阵列上的行列位置以及在图像中的像素坐标；用于确定每个 spot 在组织图像上的具体位置。

在空间转录组分析中，必要的文件包括 filtered_feature_bc_matrix.h5（基因表达矩阵）、tissue_positions_list.csv（spot的空间位置信息）、scalefactors_json.json（坐标与图像的缩放比例）以及 tissue_lowres_image.png 或 tissue_hires_image.png（组织切片图像，用于可视化）。

1.导入

rm(list = ls())
library(Seurat)
library(future)
library(ggplot2)
library(hdf5r)
library(SingleCellExperiment)
library(SpatialExperiment)
library(scater)
library(scran)

dir.create("0-GSM8633891-2")   
setwd("0-GSM8633891-2")

把GSM8633891数据集整理成了这样的一个文件结构，并且把名称进行修改。由于空转数据基本上都是一个个处理，所以手动解压、改名、移动可能更快。

2.单样本读取数据

data <- Read10X_h5("./GSM8633891/filtered_feature_bc_matrix.h5")
dim(data)
object <- CreateSeuratObject(counts = data, 
                             assay = "Spatial", 
                             min.cells = 3, #过滤在少于3个细胞中表达的基因，以减少低表达基因的干扰
                             project = "GSM8633891")
object

# 再读取
image <- Read10X_Image(image.dir = "GSM8633891/", 
                       image.name = "tissue_hires_image.png",
                       filter.matrix = TRUE)
image
  
dim(image)
image <- image[Cells(x = object)]# 筛选图像对象中包含的SPOT
DefaultAssay(object = image) <- "Spatial" # 设置图像对象的默认assay为"Spatial"

# 添加图片到object中
object[["GSM8633891"]] <- image
object@images[["GSM8633891"]]@scale.factors$lowres <- object@images[["GSM8633891"]]@scale.factors$hires

sce.all <- object
## 简单探索一下数据结构
as.data.frame(sce.all[["Spatial"]]$counts[1:4,1:4])
as.data.frame(LayerData(sce.all, assay = "Spatial", layer = "counts")[1:5,1:5])
head(sce.all@meta.data)
table(sce.all$orig.ident)

Layers(sce.all)
Assays(sce.all)

library(qs)
qsave(sce.all, file="GSM8633891/sce.all.qs")

3.数据预处理

plot1 <- VlnPlot(sce.all, 
                 features = c("nCount_Spatial","nFeature_Spatial"), 
                 pt.size = 0.1) + 
  NoLegend()
plot1
ggsave(filename = "1-QC/VlnPlot.pdf", 
       width = 12,height = 6, plot = plot1)

######
plot2 <- SpatialFeaturePlot(sce.all, 
                            features = "nCount_Spatial",
                            pt.size.factor = 2) 
plot2
ggsave(filename = "1-QC/FeaturePlot_nCount.pdf", 
       width = 18,height = 6, plot = plot2)

########
plot3 <- SpatialFeaturePlot(sce.all, 
                            features ="nFeature_Spatial",
                            pt.size.factor = 2)
plot3
ggsave(filename ="1-QC/FeaturePlot_nFeature.pdf", 
       width = 18,height = 6, plot = plot3)

4.过滤SPOT

## 过滤spot
sce.all <- subset(sce.all, subset = nFeature_Spatial > 200 )


### 计算线粒体/核糖体/特定基因集的比例
#mito_genes <- grep("^MT-", rownames(sce.all),ignore.case = T, value = T)
## 可能是13个线粒体基因
#print(mito_genes)
#sce.all <- PercentageFeatureSet(sce.all, 
#                                features = mito_genes, 
#                                col.name = "percent_mito")

5.标准化

sce.all <- SCTransform(sce.all, assay ="Spatial", verbose = T)
sce.all <- RunPCA(sce.all, assay ="SCT", verbose = FALSE)
sce.all <- FindNeighbors(sce.all, reduction ="pca", dims = 1:30)
sce.all <- FindClusters(sce.all, verbose = FALSE,resolution = 0.8)
sce.all <- RunUMAP(sce.all, reduction ="pca", dims = 1:30)

p1 <- DimPlot(sce.all, 
              reduction ="umap", 
              label = TRUE,
              label.size = 7);p1
ggsave(filename ="1-QC/DimPlot.pdf", width = 9,height = 6, plot = p1)

p2 <- SpatialDimPlot(sce.all, label = TRUE, 
                     label.size = 3,
                     pt.size.factor = 2,
                     image.alpha = 0.6);p2
ggsave(filename ="1-QC/SpatialDimPlot.pdf", width = 18,height = 6, plot = p2)

p1|p2

5.1 修改空转数据的数据样式，从Seurat变成SpatialExperiment

# 修改空转数据的数据样式，从Seurat变成SpatialExperiment
library(SpatialExperiment)
library(SingleCellExperiment)
sce <- as.SingleCellExperiment(sce.all)
# 检查 assay
assayNames(sce)       # 通常包含 counts 或 Spatial_counts
assay(sce, "counts")  # 确认 counts 是否存在
# 添加空间坐标（来自 Seurat object）
spatial_coords <- GetTissueCoordinates(sce.all)
# 保留数值列
spatial_coords_mat <- as.matrix(spatial_coords[, c("x","y")])
# 转换为SpatialExperiment
spe <- SpatialExperiment(
  assays = assays(sce),           # 保留 counts/logcounts
  colData = colData(sce),         # 保留 meta data
  spatialCoords = as.matrix(spatial_coords_mat)  # 添加空间坐标
)
class(spe)

6.基于现有单细胞数据集去映射

# 读取参考数据集
library(qs)
ref <- qread(".././0-reference/sce_celltype.qs")
# 建议大家都不要用特殊符号
# R语言不喜欢这些特殊符号
a <- ref@meta.data
a$celltype <- sub("/","_",a$celltype)
ref@meta.data <- a

# 该注释存储在对象元数据的“celltype”列中。
DimPlot(ref, group.by = "celltype", label = TRUE)

# 修改单细胞数据的数据格式
sce <- as.SingleCellExperiment(ref)

7.SPOTling分析前数据处理

# 研究者认为每个细胞身份的细胞数下采样至约100个并不会影响模型性能。
# 每个细胞身份包含超过100个细胞时，仅能带来边际性的提升，但会显著增加计算时间和资源消耗。
# 将基因集限制为每种细胞类型的标记基因，并结合最多3,000个高变基因，可以进一步优化性能和计算资源。

#特征选择
#第一步是获取每种细胞类型的标记基因。
#遵循OSCA中描述的标准化流程。
#首先执行文库大小标准化，以校正细胞特异性的偏差：
sce <- logNormCounts(sce)

8.获取重点基因

# 获取一个向量，用于指示哪些基因 既不是核糖体基因，也不是线粒体基因。
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的基因
hvg <- getTopHVGs(dec, n = 3000)

9.获取每个细胞身份的标志基因

# 可以使用任何方法，只要它能返回一个权重，表示该基因对于该细胞类型的重要性。
# 常见的指标包括 avgLogFC、AUC、pct.expressed、p-value 等。
colLabels(sce) <- colData(sce)$celltype # 是否使用celltype要根据自己的数据集进行修改

# 计算标记基因
mgs <- scoreMarkers(sce, subset.row = genes)

10.保留那些与每个细胞身份相关的基因

# 接下来，想保留那些与每个细胞身份相关的基因。
mgs_fil <- lapply(names(mgs), function(i) {
    x <- mgs[[i]]
    # 筛选并保留相关的标志基因（marker genes），即 AUC > 0.8 的基因。
    x <- x[x$mean.AUC > 0.8, ]
    # 将基因按权重从高到低排序。
    x <- x[order(x$mean.AUC, decreasing = TRUE), ]
    # 向数据框中添加基因和簇（cluster）ID。
    x$gene <- rownames(x)
    x$cluster <- i
    data.frame(x)
})
mgs_df <- do.call(rbind, mgs_fil)

11.提取细胞

#随机从每个细胞身份中选择最多100个细胞。如果某个细胞类型少于 100 个细胞，则会使用该类型的全部细胞。对于生物学差异非常大的细胞身份，可以使用更少的细胞，因为它们的转录谱差异很大。而当细胞身份在转录水平上较为相似时，需要增加细胞数量以捕捉它们之间的生物学异质性。

# 按细胞身份划分细胞索引。
idx <- split(seq(ncol(sce)), ref$celltype)
# 对每个细胞身份最多采样100个细胞并进行子集选择。

# 分析
n_cells <- 100
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)]

12.SPOTlight解卷积

res <- SPOTlight(
    x = sce,
    y = spe,
    groups = as.character(sce$celltype),
    mgs = mgs_df,
    hvg = hvg,
    weight_id = "mean.AUC",
    group_id = "cluster",
    gene_id = "gene")

13 可视化

13.1 查看Topic profiles

#查看Topic profiles。通过设置 facet = FALSE，希望评估每个topic 签名对每个细胞身份的特异性。理想情况下，每个细胞身份都会对应一个独特的topic profile
plotTopicProfiles(
    x = mod,
    y = sce$celltype,
    facet = FALSE,
    min_prop = 0.01,
    ncol = 1) +
    theme(aspect.ratio = 1)

13.2 Topic profiles分面权重值

plotTopicProfiles(
    x = mod,
    y = sce$celltype,
    facet = TRUE,
    min_prop = 0.01,
    ncol = 3)

13.3 相关性分析

plotCorrelationMatrix(mat)

13.4 共定位

现在已经知道每个spot中存在的细胞类型，就可以绘制一个表示空间相互作用的图。在图中，细胞类型之间的权越大，表示它们在同一个spot中出现的频率越高。

# 共定位（Co-localization）
plotInteractions(mat, which = "heatmap", metric = "prop")
plotInteractions(mat, which = "heatmap", metric = "jaccard")
plotInteractions(mat, which = "network")

13.5 散点饼图

绘制每个SPOT中不同细胞的百分比饼图

#还可以将每个 spot 中的细胞类型比例可视化为饼图的不同扇区。
ct <- colnames(mat)
mat[mat < 0.1] <- 0

# 定义多种颜色
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 = spatialCoords(spe),
    y = mat,
    cell_types = colnames(mat),
    img = F,
    scatterpie_alpha = 1,
    pie_scale = 0.4) +
    scale_fill_manual(
        values = pal,
        breaks = names(pal))

对比了一下RCTD的结果，稍有差别。

13.6 残差法评估准确度

还可以查看模型对每个spot预测的比例效果如何。通过观察每个spot平方和残差来评估，数值越高说明不准确度越高。

spe$res_ss <- res[[2]][colnames(spe)]
xy <- spatialCoords(spe)
spe$x <- xy[, 1]
spe$y <- xy[, 2]
ggcells(spe, aes(y, x, color = res_ss)) +
    geom_point() +
    scale_color_viridis_c() +
    coord_fixed() +
    theme_bw()

参考资料：

SPOTlight: seeded NMF regression to deconvolute spatial transcriptomics spots with single-cell transcriptomes. Nucleic Acids Res 49(9):e50.
SPOTight: https://github.com/MarcElosua/SPOTlight

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据挖掘

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

数据挖掘

#单细胞

登录后参与评论

0 条评论

热度