单细胞数据整合方法¶
一句话概述¶
利用Harmony、scVI、SCANORAMA和MNN等批次校正算法，将来自不同实验批次、技术平台或样本的单细胞RNA-seq数据整合到统一的低维空间中，消除技术变异同时保留生物学差异，是大规模单细胞图谱构建的关键步骤。
核心知识点表格¶
| 知识点 | 关键内容 | 常用工具 |
|---|---|---|
| 批次效应来源 | 试剂批次/测序run/文库制备差异 | - |
| 线性校正 | 回归去除批次变量 | Seurat CCA, limma |
| 非线性嵌入校正 | 在低维空间对齐 | Harmony, SCANORAMA |
| 深度学习整合 | VAE学习批次不变表示 | scVI, scANVI |
| 锚点方法 | 找跨批次对应细胞对 | Seurat CCA/RPCA, MNN |
| 评估指标 | kBET/LISI/ASW/ARI | scIB benchmark |
| 过校正风险 | 去除真实生物学差异 | - |
| 参考映射 | 将查询数据映射到参考图谱 | Seurat reference mapping |
各步骤详解¶
第一步：理解批次效应与整合需求¶
白话解释： 不同批次产生的单细胞数据就像用不同相机拍的照片——颜色和亮度不一样，但拍的是同样的场景。整合就是给这些照片做"色彩校正"，使它们可以放在一起比较。关键挑战是区分技术差异（需要去除）和真实的生物学差异（需要保留）。
技术细节： - 技术批次效应：文库制备、测序深度、酶批次、实验人员 - 生物学差异：不同样本/条件的真实细胞组成差异 - 过校正：把生物学差异当技术差异去除了 - 欠校正：技术差异未完全消除，仍在UMAP上分开
代码示例：
library(Seurat)
library(ggplot2)
# 加载多批次数据
batch1 <- Read10X("batch1/filtered_feature_bc_matrix/")
batch2 <- Read10X("batch2/filtered_feature_bc_matrix/")
batch3 <- Read10X("batch3/filtered_feature_bc_matrix/")
# 分别创建Seurat对象
s1 <- CreateSeuratObject(batch1, project="batch1", min.cells=3, min.features=200)
s2 <- CreateSeuratObject(batch2, project="batch2", min.cells=3, min.features=200)
s3 <- CreateSeuratObject(batch3, project="batch3", min.cells=3, min.features=200)
# 标记批次
s1$batch <- "batch1"; s2$batch <- "batch2"; s3$batch <- "batch3"
# 简单合并（不做批次校正）
merged <- merge(s1, y=c(s2, s3))
merged <- NormalizeData(merged) %>% FindVariableFeatures() %>%
ScaleData() %>% RunPCA() %>% RunUMAP(dims=1:30)
# 查看批次效应
DimPlot(merged, group.by="batch")
# 如果不同批次明显分开→需要整合
第二步：Harmony整合（推荐首选）¶
白话解释： Harmony是最快、最简单且效果好的整合方法。它在PCA空间中迭代地将不同批次的细胞"混合"到一起。核心思想：找到PCA空间的线性校正，使得每个聚类中不同批次的细胞比例接近。运行速度快（几分钟处理百万细胞），内存占用小。
技术细节： - 在PCA embedding上操作（不改变原始表达矩阵） - 迭代算法：聚类→校正→聚类→校正 - 软聚类：每个细胞可属于多个聚类（模糊分配） - 保守校正：不会强制所有批次完全混合 - 参数：theta控制校正强度，sigma控制软聚类宽度
代码示例：
# 安装: install.packages("harmony")
library(harmony)
# 在已有Seurat对象上运行Harmony
# 先做标准预处理
merged <- NormalizeData(merged) %>%
FindVariableFeatures(nfeatures=3000) %>%
ScaleData() %>%
RunPCA(npcs=50)
# 运行Harmony
merged <- RunHarmony(merged,
group.by.vars = "batch", # 批次变量
theta = 2, # 校正强度(默认2)
lambda = 1, # ridge回归正则化
sigma = 0.1, # 软聚类宽度
nclust = 50, # 聚类数
max.iter.harmony = 20, # 最大迭代次数
plot_convergence = TRUE)
# 使用校正后的embedding进行下游分析
merged <- RunUMAP(merged, reduction="harmony", dims=1:30) %>%
FindNeighbors(reduction="harmony", dims=1:30) %>%
FindClusters(resolution=0.5)
# 可视化比较
p1 <- DimPlot(merged, group.by="batch", reduction="umap") + ggtitle("After Harmony")
p2 <- DimPlot(merged, group.by="cell_type", reduction="umap") + ggtitle("Cell Types")
p1 + p2
# 多批次变量同时校正
merged <- RunHarmony(merged,
group.by.vars = c("batch", "donor"), # 同时校正批次和个体
theta = c(2, 1)) # 不同变量不同校正强度
第三步：scVI深度学习整合¶
白话解释： scVI使用变分自编码器(VAE)学习一个"批次不变"的低维表示。它同时建模基因表达的技术噪声（如测序深度、批次）和生物学信号，在潜在空间中自动分离两者。特别适合大规模数据和复杂批次结构。
代码示例：
import scanpy as sc
import scvi
import numpy as np
# 1. 准备数据
adata = sc.read_h5ad("merged_adata.h5ad")
# 标准预处理
sc.pp.filter_genes(adata, min_cells=10)
sc.pp.highly_variable_genes(adata, n_top_genes=3000,
batch_key="batch", flavor="seurat_v3")
adata_hvg = adata[:, adata.var["highly_variable"]].copy()
# 2. 设置scVI模型
scvi.model.SCVI.setup_anndata(
adata_hvg,
layer="counts", # 原始计数
batch_key="batch", # 批次变量
categorical_covariate_keys=["donor"], # 其他分类协变量
continuous_covariate_keys=["percent_mito"]) # 连续协变量
# 3. 训练模型
model = scvi.model.SCVI(
adata_hvg,
n_latent=30, # 潜在维度数
n_layers=2, # 神经网络层数
n_hidden=128, # 隐藏层节点数
dropout_rate=0.1,
dispersion="gene", # 离散参数按基因估计
gene_likelihood="zinb") # 零膨胀负二项分布
model.train(
max_epochs=400,
early_stopping=True,
batch_size=128,
plan_kwargs={"lr": 1e-3})
# 查看训练曲线
model.history["elbo_train"].plot()
# 4. 获取校正后的低维表示
adata_hvg.obsm["X_scVI"] = model.get_latent_representation()
# 5. 下游分析
sc.pp.neighbors(adata_hvg, use_rep="X_scVI", n_neighbors=30)
sc.tl.umap(adata_hvg)
sc.tl.leiden(adata_hvg, resolution=0.8)
# 可视化
sc.pl.umap(adata_hvg, color=["batch", "cell_type", "leiden"])
# 6. 获取校正后的表达矩阵（可选）
# 用于差异表达分析
denoised = model.get_normalized_expression(
library_size=1e4, # 标准化到10000
transform_batch="batch1") # 校正到batch1的水平
# 7. scANVI（半监督整合，使用已知标签）
# 注意：使用from_scvi_model时不需要单独调用SCANVI.setup_anndata
# scANVI会继承SCVI的setup，并自动注册labels_key
scanvi_model = scvi.model.SCANVI.from_scvi_model(
model,
labels_key="cell_type", # 细胞类型标签列
unlabeled_category="Unknown") # 未标记细胞的标签值
scanvi_model.train(max_epochs=20)
adata_hvg.obsm["X_scANVI"] = scanvi_model.get_latent_representation()
# scANVI通常比scVI整合效果更好（有标签指导）
第四步：Seurat CCA/RPCA整合¶
白话解释： Seurat的整合流程先在不同批次间找"锚点"——对应到同种细胞的细胞对，然后根据锚点的关系校正数据。CCA (Canonical Correlation Analysis) 找最相关的跨批次维度；RPCA (Reciprocal PCA) 更快但稍保守。Seurat v5 推荐使用 IntegrateLayers() 新API（见上方第二步代码），以下是 v4 旧API写法供参考。
代码示例：
library(Seurat)
# === Seurat v4 旧API（仍可用但不推荐，v5 请用 IntegrateLayers） ===
# 分别预处理每个批次
obj_list <- list(s1, s2, s3)
obj_list <- lapply(obj_list, function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, nfeatures=3000)
x
})
# 找整合特征
features <- SelectIntegrationFeatures(obj_list, nfeatures=3000)
# 方法1: CCA整合 (标准，适合<200K细胞)
anchors <- FindIntegrationAnchors(obj_list,
anchor.features=features,
reduction="cca", # CCA方法
dims=1:30)
integrated <- IntegrateData(anchors, dims=1:30)
# 在integrated assay上做下游分析
DefaultAssay(integrated) <- "integrated"
integrated <- ScaleData(integrated) %>% RunPCA() %>%
RunUMAP(dims=1:30) %>% FindNeighbors(dims=1:30) %>%
FindClusters(resolution=0.5)
# 方法2: RPCA整合 (更快，适合大数据)
obj_list <- lapply(obj_list, function(x) {
x <- ScaleData(x, features=features)
x <- RunPCA(x, features=features)
x
})
anchors_rpca <- FindIntegrationAnchors(obj_list,
anchor.features=features,
reduction="rpca", # RPCA方法
dims=1:30)
integrated_rpca <- IntegrateData(anchors_rpca, dims=1:30)
# 方法3: Seurat V5 整合（推荐，更高效）
# 注意：V5的关键是先split再整合，整合后再JoinLayers
merged_v5 <- merge(s1, y=c(s2, s3))
# 按批次拆分layers（V5核心概念：每个批次是一个layer）
merged_v5[["RNA"]] <- split(merged_v5[["RNA"]], f = merged_v5$batch)
# 标准预处理（NormalizeData等会自动在每个layer上分别运行）
merged_v5 <- NormalizeData(merged_v5) %>%
FindVariableFeatures() %>% ScaleData() %>% RunPCA()
# IntegrateLayers支持多种整合方法
merged_v5 <- IntegrateLayers(merged_v5,
method = CCAIntegration, # 或RPCAIntegration/HarmonyIntegration
orig.reduction = "pca",
new.reduction = "integrated.cca")
# 整合完成后再JoinLayers（差异表达分析前必须rejoin）
merged_v5[["RNA"]] <- JoinLayers(merged_v5[["RNA"]])
# 使用整合后的reduction进行下游分析
merged_v5 <- FindNeighbors(merged_v5, reduction = "integrated.cca", dims = 1:30) %>%
FindClusters(resolution = 0.5) %>%
RunUMAP(reduction = "integrated.cca", dims = 1:30)
第五步：SCANORAMA和MNN整合¶
代码示例：
# === SCANORAMA ===
import scanorama
import scanpy as sc
# 按批次分割数据
batches = []
genes_list = []
for batch in adata.obs["batch"].unique():
batch_data = adata[adata.obs["batch"] == batch]
batches.append(batch_data.X.toarray())
genes_list.append(batch_data.var_names.tolist())
# 运行SCANORAMA
corrected, genes = scanorama.correct(batches, genes_list,
return_dimred=True, dimred=100)
# 合并结果
corrected_concat = np.concatenate(corrected)
# 添加到adata
adata.obsm["X_scanorama"] = np.array(corrected_concat)
# === MNN (Mutual Nearest Neighbors) ===
# 使用mnnpy或Scanpy内置
sc.external.pp.mnn_correct(adata, batch_key="batch")
# 或
import mnnpy
corrected_data = mnnpy.mnn_correct(
*[adata[adata.obs["batch"]==b].X for b in batches],
var_index=adata.var_names)
第六步：整合质量评估¶
白话解释： 整合后需要评估两个方面：1）批次混合程度（不同批次的细胞是否在同一聚类中混合良好？）；2）生物学保守性（相同细胞类型是否仍然聚在一起？是否丢失了真实的生物学差异？）。好的整合应该两者兼顾。
代码示例：
# 使用scib评估整合质量
# pip install scib # 经典版本（仍可用）
# pip install scib-metrics # 新版本（GPU加速，推荐）
# === 方法A: 使用scib（经典API） ===
import scib
metrics = scib.metrics.metrics(
adata,
adata_int=adata_integrated,
batch_key="batch",
label_key="cell_type",
embed="X_scVI", # 整合后的embedding
cluster_key="leiden",
type_="embed",
organism="human",
n_jobs=8)
# 关键指标:
# 批次校正:
# - kBET: 0-1, 越高批次混合越好
# - Graph connectivity: 图连通性
# - Batch ASW: 批次轮廓系数(越接近0越好)
# - iLISI: 整合局部逆辛普森指数(越高越好)
# 生物保守:
# - ARI: 调整兰德指数（聚类与真实标签一致性）
# - NMI: 标准化互信息
# - Cell type ASW: 细胞类型轮廓系数(越接近1越好)
# - cLISI: 细胞类型局部逆辛普森指数(越接近1越好)
print(metrics)
# === 方法B: 使用scib-metrics（新版，推荐） ===
# import scib_metrics
# 提供与scib兼容的API，但运行更快
# 详见: https://scib-metrics.readthedocs.io/
# 简单的LISI评估（scib经典API）
from scib.metrics import lisi
ilisi = lisi.ilisi_graph(adata_integrated, batch_key="batch",
type_="embed", use_rep="X_scVI")
clisi = lisi.clisi_graph(adata_integrated, label_key="cell_type",
type_="embed", use_rep="X_scVI")
print(f"iLISI (批次混合): {ilisi:.3f}")
print(f"cLISI (生物保守): {clisi:.3f}")
# R中的评估
library(kBET)
# kBET检验
data_embed <- Embeddings(integrated, "harmony")
batch_labels <- integrated$batch
kbet_result <- kBET(data_embed, batch_labels,
k0=30, testSize=500, do.pca=FALSE)
cat("kBET rejection rate:",
mean(kbet_result$summary$kBET.observed), "\n")
# 越低越好（接近0=批次混合完美）
实战命令¶
#!/usr/bin/env Rscript
# === 单细胞数据整合完整流程 ===
library(Seurat); library(harmony)
# 1. 加载多批次数据
obj_list <- lapply(list.files("data/", full.names=TRUE), function(f) {
obj <- Read10X(f) %>% CreateSeuratObject(min.cells=3, min.features=200)
obj$batch <- basename(f)
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern="^MT-")
subset(obj, nFeature_RNA>200 & nFeature_RNA<6000 & percent.mt<20)
})
# 2. 合并和预处理
merged <- Reduce(merge, obj_list)
merged <- NormalizeData(merged) %>% FindVariableFeatures(nfeatures=3000) %>%
ScaleData() %>% RunPCA(npcs=50)
# 3. Harmony整合
merged <- RunHarmony(merged, "batch", theta=2)
merged <- RunUMAP(merged, reduction="harmony", dims=1:30) %>%
FindNeighbors(reduction="harmony", dims=1:30) %>%
FindClusters(resolution=0.5)
# 4. 可视化和评估
DimPlot(merged, group.by=c("batch","seurat_clusters"), ncol=2)
saveRDS(merged, "integrated_seurat.rds")
面试常问点¶
Q1: Harmony和scVI的选择依据？ A: Harmony：快速（分钟级）、简单、对中等批次效应效果好、不改变表达矩阵（只校正embedding）。scVI：处理复杂/强批次效应更好、提供概率框架（可以做差异表达）、适合大数据、支持多协变量。小数据集/简单批次用Harmony；大数据集/复杂设计用scVI。
Q2: 如何判断是否过校正？ A: 1）已知不同条件/疾病状态的细胞被混在一起→过校正；2）特异性标记基因在细胞类型间区分度降低→过校正；3）cLISI过低（细胞类型特异性被破坏）；4）比较整合前后的已知marker表达模式。最佳实践：多种theta/lambda参数对比。
Q3: 整合是否改变了原始表达值？ A: Harmony和SCANORAMA只改变低维embedding(PCA空间)，不改变原始表达矩阵——差异表达分析应使用原始数据。scVI提供校正后的表达估计但也保留原始数据。Seurat v4 CCA的integrated assay创建了校正后的表达矩阵（只用于聚类，DE应用原始RNA assay）。Seurat v5的IntegrateLayers同样只返回校正后的低维reduction，原始counts/data layers不变，DE前需要JoinLayers合并回原始RNA assay。
Q4: 哪些情况不适合做数据整合？ A: 1）数据间不共享细胞类型时（如完全不同的组织）；2）感兴趣的信号就是"批次"差异时（如疾病vs对照的差异就是研究目标）；3）技术差异极端（不同物种、不同建库方法差异太大）。这些情况下整合可能消除感兴趣的信号。
Q5: scIB benchmark中哪些指标最重要？ A: 整体分数是batch correction (40%) + bio conservation (60%)的加权平均。关键指标：iLISI（批次混合）、cLISI（生物保守）、ARI（聚类准确性）、Graph connectivity（图连通性）。实际中最关心：已知细胞类型是否正确聚类 + UMAP上不同批次是否混合。
易错点¶
不做整合就直接合并分析：不同批次的技术差异会主导聚类结果，导致按批次而非细胞类型分群。
过度校正：theta/lambda参数设置过大会消除真实生物学差异。建议从默认参数开始。
在整合后的assay上做差异表达：Seurat v4的integrated assay校正了表达值，不适合做DE。Seurat v5使用IntegrateLayers只产生新reduction，DE前必须先调用JoinLayers()合并原始RNA layers。DE应该使用原始RNA assay/layer。
高变基因选择不当：应该在每个批次内分别选高变基因再取并集，而非合并后选择（会被批次效应主导）。
参考与查询不匹配：reference mapping时，参考图谱必须包含查询数据中的所有细胞类型，否则查询中的新类型会被错误分配。
补充知识¶
整合方法选择指南¶
| 场景 | 推荐方法 | 原因 |
|---|---|---|
| 快速整合(<100K cells) | Harmony | 快，简单，效果好 |
| 大规模(>500K cells) | scVI/Harmony | 可扩展 |
| 强批次效应 | scVI/scANVI | 深度学习更强力 |
| 有细胞标签 | scANVI | 半监督利用标签 |
| 保守校正 | Seurat RPCA | 不过度校正 |
| 跨平台(10X vs Smart-seq) | scVI/Harmony | 处理平台差异 |
| 参考图谱映射 | Seurat reference | 已有参考 |
scIB benchmark总结(Luecken et al., 2022, Nature Methods)¶
整体排名(综合得分 = 40%批次校正 + 60%生物保守):
1. scVI/scANVI - 深度学习, 综合最优
2. Harmony - 快速, 接近最优
3. Scanorama - 稳定
4. Seurat CCA/RPCA - 经典
5. BBKNN - 简单
6. Combat - 过于保守
7. MNN - 较旧
关键发现:
- 没有一种方法在所有场景都最好
- scVI在复杂设计中优势明显
- Harmony在简单设计中最佳性价比
- 过校正比欠校正危害更大
注意：scIB工具现改名为scib-metrics（pip install scib-metrics），
提供更高效的GPU加速实现，API有所变化。