空间转录组数据分析¶
一句话概述¶
利用10x Visium等空间转录组技术获取组织切片中基因表达的空间分布信息，通过Space Ranger处理原始数据，Seurat/Scanpy进行降维聚类，spot反卷积估计细胞组成，最终识别空间基因表达模式，实现基因表达与组织形态的精准映射。
核心知识点表格¶
| 知识点 | 关键内容 | 常用工具 |
|---|---|---|
| 数据生成 | 10x Visium/MERFISH/Slide-seq/Stereo-seq | - |
| 原始处理 | FASTQ→计数矩阵+空间坐标 | Space Ranger |
| 质控与标准化 | 过滤、SCTransform/LogNorm | Seurat, Scanpy |
| 降维聚类 | PCA→UMAP+空间信息图聚类 | BayesSpace, SpaGCN |
| Spot反卷积 | 从spot估计细胞类型比例 | cell2location, RCTD, SPOTlight |
| 空间变异基因 | 发现空间pattern的基因 | SpatialDE, SPARK, nnSVG |
| 细胞通讯 | 空间邻近细胞间信号 | CellChat, COMMOT, stLearn |
| 高分辨率方法 | 亚细胞分辨率空间组学 | Xenium, MERSCOPE |
各步骤详解¶
第一步：理解空间转录组技术¶
白话解释： 传统单细胞RNA-seq虽然能分析单个细胞的基因表达，但需要将组织消化成单个细胞，丢失了空间位置信息。空间转录组技术保留了基因表达在组织中的位置：就像在组织切片上的每个位置都测了一次RNA-seq，既知道表达了什么基因，又知道在哪里表达。
技术细节： - 10x Visium: 55μm spot，~1-10个细胞/spot，全转录组，约5000 spots/section - Slide-seq: 10μm beads，接近单细胞，全转录组 - MERFISH/Xenium: 亚细胞分辨率，但只检测预选基因（100-1000个） - Stereo-seq: 500nm分辨率，超高密度覆盖 - Visium HD: 2μm bin，接近单细胞分辨率
第二步：Space Ranger数据处理¶
白话解释： Space Ranger是10x官方的Visium数据处理软件。它做三件事：1）将FASTQ reads比对到参考基因组并生成计数矩阵；2）自动检测组织覆盖区域；3）将barcode坐标映射到H&E图像上的空间位置。
代码示例：
# 安装Space Ranger
# 从10x Genomics官网下载
# 运行Space Ranger
# Space Ranger v3.0+ 必须指定 --create-bam（true/false），否则命令报错
spaceranger count \
--id=sample_A \
--transcriptome=/ref/GRCh38-2024 \
--fastqs=/data/fastq/sample_A/ \
--sample=SampleA \
--image=/data/images/sample_A.tif \
--slide=V19L29-098 \
--area=B1 \
--create-bam false \
--localcores=16 \
--localmem=64
# 输出目录结构:
# sample_A/outs/
# filtered_feature_bc_matrix/ - 过滤后的计数矩阵
# spatial/ - 空间信息
# tissue_positions.csv - 每个spot的行列坐标和像素坐标
# scalefactors_json.json - 缩放因子
# tissue_hires_image.png - 高分辨率组织图像
# tissue_lowres_image.png - 低分辨率组织图像
# raw_feature_bc_matrix/ - 原始计数矩阵
# web_summary.html - QC报告
# 关键QC指标：
# - Reads Mapped Confidently to Transcriptome: >50%
# - Mean Reads per Spot: >20,000 (理想>50,000)
# - Median Genes per Spot: >2,000
# - Fraction of Spots Under Tissue: 组织覆盖率
第三步：Seurat空间数据分析¶
白话解释： 用Seurat加载Space Ranger的输出，进行标准的质控、标准化、降维、聚类，但特别加入空间维度——在UMAP聚类图和组织切片图上同时展示结果。
代码示例：
library(Seurat)
library(ggplot2)
library(patchwork)
# 加载Space Ranger输出
spatial_data <- Load10X_Spatial(
data.dir = "sample_A/outs/",
filename = "filtered_feature_bc_matrix.h5",
assay = "Spatial",
slice = "tissue_section")
# 质控
spatial_data[["percent.mt"]] <- PercentageFeatureSet(spatial_data, pattern="^MT-")
# 可视化QC指标在空间上的分布
SpatialFeaturePlot(spatial_data,
features = c("nCount_Spatial", "nFeature_Spatial", "percent.mt"))
# 过滤
spatial_data <- subset(spatial_data,
nFeature_Spatial > 200 & nFeature_Spatial < 8000 &
percent.mt < 25)
# SCTransform标准化（推荐用于空间数据）
spatial_data <- SCTransform(spatial_data,
assay = "Spatial",
verbose = FALSE)
# 降维和聚类
spatial_data <- RunPCA(spatial_data, assay = "SCT")
spatial_data <- FindNeighbors(spatial_data, reduction = "pca", dims = 1:30)
spatial_data <- FindClusters(spatial_data, resolution = 0.8)
spatial_data <- RunUMAP(spatial_data, reduction = "pca", dims = 1:30)
# 空间可视化
p1 <- DimPlot(spatial_data, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(spatial_data, label = TRUE, label.size = 3)
p1 + p2
# 单基因空间表达
SpatialFeaturePlot(spatial_data,
features = c("EPCAM", "COL1A1", "PTPRC", "CD3E"),
alpha = c(0.1, 1)) # 透明度渐变
# 保存
saveRDS(spatial_data, "spatial_seurat.rds")
第四步：空间变异基因检测¶
白话解释： 空间变异基因(Spatially Variable Genes, SVGs)是表达模式在空间上有特定规律的基因——比如只在肿瘤中心高表达、或沿着某个方向梯度变化。这些基因揭示了组织的空间功能分区。
代码示例：
# 方法1: Seurat内置的Moran's I
spatial_data <- FindSpatiallyVariableFeatures(spatial_data,
assay = "SCT",
features = VariableFeatures(spatial_data),
selection.method = "moransi", # 注意：markvariogram 在 Seurat 5.x 中已移除，仅支持 "moransi"
nfeatures = 1000)
# 查看top空间变异基因
top_svgs <- head(SpatiallyVariableFeatures(spatial_data), 20)
SpatialFeaturePlot(spatial_data, features = top_svgs[1:6])
# 方法2: nnSVG (推荐，基于高斯过程)
library(nnSVG)
# 从Seurat转为SpatialExperiment
library(SpatialExperiment)
spe <- as.SpatialExperiment(spatial_data)
# 运行nnSVG
spe <- nnSVG(spe, n_threads = 8)
# 结果
svg_results <- rowData(spe)
svg_results <- svg_results[order(svg_results$rank), ]
head(svg_results[, c("gene_name", "LR_stat", "padj", "rank")], 20)
# 方法3: SpatialDE (Python)
import scanpy as sc
import SpatialDE
# 读取数据
adata = sc.read_visium("sample_A/outs/")
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
# 准备坐标
coords = adata.obsm["spatial"]
# 运行SpatialDE
results = SpatialDE.run(coords, adata.to_df())
results = results.sort_values("qval")
print(results.head(20)[["g", "l", "qval", "FSV"]])
# 空间模式分组
patterns = SpatialDE.spatial_patterns(
coords, adata.to_df(), results[results["qval"] < 0.05],
n_patterns=5, length=1.0)
第五步：Spot反卷积¶
白话解释： Visium的每个spot通常覆盖1-10个细胞，得到的是混合表达谱。反卷积(deconvolution)利用单细胞RNA-seq作为参考，估计每个spot中各种细胞类型的比例——就像根据鸡尾酒的味道推断各种成分的比例。
代码示例：
# 方法1: RCTD (Robust Cell Type Decomposition)
library(spacexr)
# 准备单细胞参考
sc_data <- readRDS("scRNA_reference.rds")
sc_counts <- GetAssayData(sc_data, slot="counts")
sc_celltypes <- sc_data$celltype
reference <- Reference(sc_counts, sc_celltypes)
# 准备空间数据
sp_counts <- GetAssayData(spatial_data, slot="counts")
coords <- GetTissueCoordinates(spatial_data)
query <- SpatialRNA(coords, sp_counts)
# 运行RCTD
RCTD <- create.RCTD(query, reference, max_cores = 8)
RCTD <- run.RCTD(RCTD, doublet_mode = "full") # full=多细胞类型
# 提取结果
weights <- normalize_weights(RCTD@results$weights)
# weights: spots × cell_types 的比例矩阵
# 将结果加入Seurat
spatial_data@meta.data <- cbind(spatial_data@meta.data, as.data.frame(weights))
# 可视化各细胞类型空间分布
SpatialFeaturePlot(spatial_data,
features = c("T_cells", "B_cells", "Macrophages", "Epithelial"),
alpha = c(0.1, 1))
# 方法2: cell2location (Python, 推荐)
import cell2location
import scvi
# 准备单细胞参考（训练参考模型）
sc_adata = sc.read_h5ad("scRNA_reference.h5ad")
cell2location.models.RegressionModel.setup_anndata(
sc_adata, labels_key="celltype")
mod_ref = cell2location.models.RegressionModel(sc_adata)
mod_ref.train(max_epochs=250)
sc_adata = mod_ref.export_posterior(sc_adata)
# 提取细胞类型特征
inf_aver = sc_adata.varm["means_per_cluster_mu_fg"]
# 准备空间数据
sp_adata = sc.read_visium("sample_A/outs/")
intersect_genes = sp_adata.var_names.intersection(inf_aver.index)
sp_adata = sp_adata[:, intersect_genes]
inf_aver = inf_aver.loc[intersect_genes, :]
# 训练空间模型
cell2location.models.Cell2location.setup_anndata(sp_adata)
mod_spatial = cell2location.models.Cell2location(
sp_adata, cell_state_df=inf_aver,
N_cells_per_location=10, # 每个spot估计的细胞数
detection_alpha=200)
mod_spatial.train(max_epochs=30000)
# 导出结果
sp_adata = mod_spatial.export_posterior(sp_adata)
# sp_adata.obsm["q05_cell_abundance_w_sf"] 包含细胞类型丰度
# 可视化
cell2location.utils.plot_spatial(sp_adata,
color=["T_cells", "B_cells", "Macrophages"])
第六步：空间细胞通讯与模式分析¶
白话解释： 空间转录组的独特优势是能分析相邻细胞间的通讯。不像传统CellChat只能推断配体-受体关系是否可能存在，空间数据能直接验证配体和受体是否在空间上邻近——只有真正"碰得到"的细胞才能通讯。
代码示例：
# 使用squidpy进行空间分析
import squidpy as sq
# 构建空间邻域图
sq.gr.spatial_neighbors(adata, n_rings=2, coord_type="grid")
# 空间自相关 (Moran's I)
sq.gr.spatial_autocorr(adata, mode="moran",
genes=adata.var_names[:1000], n_jobs=8)
# 高Moran's I = 空间聚集性强
# 邻域富集分析
sq.gr.nhood_enrichment(adata, cluster_key="clusters")
sq.pl.nhood_enrichment(adata, cluster_key="clusters")
# 配体-受体空间共定位
sq.gr.ligrec(adata,
n_perms=100,
cluster_key="clusters",
interactions_params={"resources": "CellPhoneDB"})
sq.pl.ligrec(adata, cluster_key="clusters",
source_groups=["Tumor"], target_groups=["Macrophage"])
# 空间模式识别 (Ripley's statistics)
sq.gr.ripley(adata, cluster_key="clusters", mode="L")
sq.pl.ripley(adata, cluster_key="clusters", mode="L")
实战命令¶
#!/bin/bash
# === 10x Visium空间转录组完整分析流程 ===
# 1. Space Ranger处理
# Space Ranger v3.0+ 必须指定 --create-bam
spaceranger count --id=patient1 --transcriptome=/ref/GRCh38 \
--fastqs=/data/fastq/patient1/ --image=/data/images/patient1.tif \
--slide=V19L29-098 --area=A1 --create-bam false --localcores=16
# 2. R分析脚本
Rscript << 'EOF'
library(Seurat); library(spacexr)
# 加载和处理
sp <- Load10X_Spatial("patient1/outs/")
sp <- SCTransform(sp, assay="Spatial") %>%
RunPCA(assay="SCT") %>%
FindNeighbors(dims=1:30) %>%
FindClusters(resolution=0.8) %>%
RunUMAP(dims=1:30)
# 空间变异基因
sp <- FindSpatiallyVariableFeatures(sp, assay="SCT", selection.method="moransi")
# 反卷积
ref <- readRDS("sc_reference.rds")
reference <- Reference(GetAssayData(ref,"counts"), ref$celltype)
query <- SpatialRNA(GetTissueCoordinates(sp), GetAssayData(sp,"counts"))
rctd <- create.RCTD(query, reference, max_cores=8)
rctd <- run.RCTD(rctd, doublet_mode="full")
saveRDS(sp, "spatial_analyzed.rds")
EOF
面试常问点¶
Q1: Visium spot中有多个细胞，如何处理这个"混合"问题？ A: 两种策略：1）反卷积(deconvolution)：利用scRNA-seq参考估计每个spot的细胞类型比例（cell2location/RCTD）；2）亚细胞分辨率技术：使用Visium HD(2μm)、Stereo-seq或MERFISH直接达到单细胞甚至亚细胞分辨率。对于标准Visium，反卷积是必要的分析步骤。
Q2: 如何整合多个切片的空间数据？ A: 1）同一组织的连续切片：使用PASTE/STAligner进行3D重建和对齐；2）不同样本的对比：先独立分析再用Harmony等方法批次校正，或使用GraphST等工具直接整合。3）注意批次效应可能包含技术变异和真实生物学空间差异。
Q3: 空间变异基因检测和差异表达分析有什么不同？ A: 差异表达(DE)比较两组样本间的表达量差异（t-test/Wilcoxon）；空间变异(SV)检测基因表达在空间上是否有非随机分布模式（如梯度、岛状聚集），使用空间自相关统计量（Moran's I、高斯过程）。一个基因可以在所有spots中均匀表达但仍是两组间的DEG。
Q4: 如何验证反卷积结果的准确性？ A: 1）与H&E病理标注比对（如已知区域的细胞组成）；2）免疫组化/免疫荧光染色验证关键marker的空间分布；3）使用连续切片的单细胞数据验证；4）比较多种反卷积方法的一致性；5）模拟数据benchmarking。
易错点¶
忽略组织质量对数据的影响：组织切片的固定、染色、RNA质量直接影响数据。损伤区域的spot应标记/排除。
反卷积参考数据库不匹配：参考scRNA-seq必须来自相同或高度相似的组织/物种/疾病状态，否则反卷积结果不可靠。
空间聚类过度依赖表达信息：标准Seurat聚类只用表达信息，未利用空间邻近关系。应使用BayesSpace/SpaGCN等整合空间信息的方法。
spot过滤过于激进：组织边缘的spot质量低但可能有生物学意义（如肿瘤边界）。过滤需谨慎。
忽略切片方向和深度效应：2D切片是3D组织的投影，相邻spots的关系取决于切割角度。
补充知识¶
主流空间转录组技术比较¶
| 技术 | 分辨率 | 基因数 | 通量 | 适用场景 |
|---|---|---|---|---|
| Visium | 55μm(~5-10cells) | 全转录组 | ~5000 spots | 常规发现 |
| Visium HD | 2μm | 全转录组 | ~250K bins | 高分辨率 |
| Slide-seq V2 | 10μm | 全转录组 | ~50K beads | 近单细胞 |
| MERFISH | 亚细胞 | 100-10000基因 | 百万+细胞 | 单分子验证 |
| Stereo-seq | 500nm | 全转录组 | 超大面积 | 大组织/发育 |
| Xenium | 亚细胞 | 100-5000基因 | 数十万细胞 | 临床/验证 |