跟着Theis Lab学单细胞: 跨越批次的鸿沟（下）

· 10 min read · Tutorials

跟着Theis Lab学单细胞: 跨越批次的鸿沟（下）

跟着Theis Lab学单细胞: 跨越批次的鸿沟（下）

单细胞最佳实践

细胞整合（上）

本章导读 SECTION PREVIEW

1. BBKNN如何通过“修改邻接图”这一巧妙思路实现闪电般的快速整合？

2. 作为经典方法的Seurat MNN是如何通过寻找“锚点”来校正表达矩阵的？

3. 如何使用scib工具对整合质量进行定量打分？

4. 面对不同数据和需求，到底该如何选择最适合你的整合方法？

本节对应原书 12 Integration 12.7-12.10

在单细胞RNA测序的世界里，批次效应（Batch Effect） 是几乎所有分析者都会遇到的核心挑战。当我们希望整合来自不同实验批次、不同供体、甚至不同实验室的数据以获得更全面的生物学洞见时，如何有效去除批次效应，让数据“求同存异”，就成了整个分析的成败关键。本文将系统性地介绍单细胞数据整合的原理、方法与实战。

上篇详见跟着Theis Lab学单细胞: 跨越批次的鸿沟（上）

基于图的整合方法：BBKNN

BBKNN (Batch Balanced k-Nearest Neighbors) 采用了一种与众不同的整合思路。它并不去生成一个经过批次校正的低维嵌入或表达矩阵，而是直接在构建细胞近邻图（kNN graph） 的环节进行“干预”。

核心原理

标准kNN算法会在整个数据集中寻找与每个细胞最相似的邻居。而在存在强大批次效应时，这几乎不可避免地会导致细胞优先与同一批次的邻居连接。BBKNN的核心创新在于，它在为每个细胞寻找邻居时，会分别在每一个批次内部 寻找最近邻，然后将这些来自不同批次的邻居合并，从而强制在不同批次之间建立连接。

这个看似简单的修改，在应对强烈批次效应时却异常有效。

import bbknn

# 同样使用批次感知的HVG数据
adata_bbknn = adata_hvg.copy()

# BBKNN在PCA空间上运行
sc.tl.pca(adata_bbknn)

# 运行BBKNN整合
# 教程建议：细胞数 > 100,000时，neighbors_within_batch=25；否则设为3
neighbors_within_batch = 3 if adata_hvg.n_obs < 100000 else 25

bbknn.bbknn(
    adata_bbknn,
    batch_key=batch_key,
    neighbors_within_batch=neighbors_within_batch
)

# 直接基于BBKNN生成的整合邻接图计算UMAP
sc.tl.umap(adata_bbknn)
sc.pl.umap(
    adata_bbknn,
    color=[label_key, batch_key],
    wspace=0.5,
    title=["Cell Type (BBKNN integrated)", "Batch (BBKNN integrated)"]
)

BBKNN整合结果BBKNN整合结果

BBKNN最大的优势是速度极快 ，非常适合在分析初期进行快速的探索。但它的局限性也很明显：由于只输出一个整合后的图结构，它无法用于那些需要校正后表达矩阵的下游分析（如差异基因表达分析）。

线性嵌入整合方法：Seurat v3 (MNN)

Seurat v3的整合策略基于**“锚点”（Anchors）** 和互近邻（Mutual Nearest Neighbors, MNN） 的概念。简单来说，MNN是指一对分别来自不同批次的细胞，它们在降维后的空间中互为对方的最近邻。Seurat v3认为这些MNN对代表了跨批次的同一种细胞状态，并将它们作为“锚点”。

整合过程分为两步：

  1. 寻找锚点 ：在不同批次之间识别出足够多的MNN对。
  2. 整合数据 ：基于这些锚点计算出一个校正向量，将一个批次的数据“对齐”到另一个批次上，最终输出一个整合后的基因表达矩阵。

由于Seurat是R包，我们需要借助一些工具在Python和R之间传递数据。

1. Python端数据准备

 # 同样使用HVG数据
adata_seurat = adata_hvg.copy()

# 为与R兼容，确保批次和细胞类型列是字符串格式
adata_seurat.obs[batch_key] = adata_seurat.obs[batch_key].astype(str)
adata_seurat.obs[label_key] = adata_seurat.obs[label_key].astype(str)

# 删除uns槽中一些可能导致转换失败的复杂对象
del adata_seurat.uns

2. 通过rpy2接口接入R代码

 %%R -i adata_seurat -o integrated_expr

library(Seurat)
# 将AnnData对象转换为Seurat对象
seurat_obj <- as.Seurat(adata_seurat, counts = "X", data = NULL)

# 按批次拆分成列表
batch_list <- SplitObject(seurat_obj, split.by = batch_key)

# 寻找整合锚点
anchors <- FindIntegrationAnchors(
    object.list = batch_list,
    anchor.features = rownames(batch_list[[1]]) # 使用所有HVG作为特征
)

# 利用锚点进行数据整合
integrated <- IntegrateData(anchorset = anchors)

# 提取校正后的表达矩阵
integrated_expr <- GetAssayData(integrated, assay = "integrated", slot = "data")

# 转置以符合AnnData的格式 (cells x genes)
integrated_expr <- t(integrated_expr)

3. 返回Python并可视化

 # 将R返回的校正矩阵存入AnnData对象
adata_seurat.obsm["X_seurat"] = integrated_expr.copy()

# 重新计算PCA和UMAP
sc.pp.neighbors(adata_seurat, use_rep="X_seurat")
sc.tl.umap(adata_seurat)
sc.pl.umap(
    adata_seurat,
    color=[label_key, batch_key],
    wspace=0.5,
    title=["Cell Type (Seurat integrated)", "Batch (Seurat integrated)"]
)

Seurat MNN整合结果Seurat MNN整合结果

Seurat v3的优势在于能输出一个经过校正的、可解释的基因表达矩阵，这对于后续进行差异表达分析等任务至关重要。它在处理简单的批次效应时表现非常稳定。

整合质量的定量评估

只通过UMAP图来判断整合效果的好坏是远远不够的，因为视觉呈现可能具有误导性。一个真正好的整合，应该同时实现两个目标：

  1. 批次校正 ：最大程度地混合不同批次的细胞，消除技术差异。
  2. 生物学结构保留 ：保持细胞类型间的差异，不破坏真实的生物学异质性。

我们需要定量的指标来衡量这两个方面。scib  (single-cell integration benchmarking) 包为此提供了一个标准化的评估框架。

使用scib进行快速评估

 import scib

# 为之前运行的每种方法计算一系列评估指标
# 注意：需要传入原始数据(adata)和整合后的数据(adata_scvi等)进行比较

# 未整合的数据作为基线
metrics_hvg = scib.metrics.metrics_fast(
    adata, adata_hvg, batch_key, label_key, embed="X_pca"
)
# scVI
metrics_scvi = scib.metrics.metrics_fast(
    adata, adata_scvi, batch_key, label_key, embed="X_scVI"
)
# scANVI
metrics_scanvi = scib.metrics.metrics_fast(
    adata, adata_scanvi, batch_key, label_key, embed="X_scANVI"
)
# BBKNN (在PCA空间评估)
metrics_bbknn = scib.metrics.metrics_fast(
    adata, adata_bbknn, batch_key, label_key, embed='X_pca', type_='graph'
)
# Seurat
metrics_seurat = scib.metrics.metrics_fast(
    adata, adata_seurat, batch_key, label_key, embed="X_seurat"
)

单个方法的结果展示单个方法的结果展示

整合并可视化评估结果

 import pandas as pd

# 将所有结果合并到一个DataFrame中
metrics = pd.concat(
    [metrics_scvi, metrics_scanvi, metrics_bbknn,
     metrics_seurat, metrics_hvg],
    axis="columns"
)
metrics.columns = ["scVI", "scANVI", "BBKNN", "Seurat", "Unintegrated"]

# 选择几个关键指标进行展示
metrics_plot = metrics.loc[
    ["ASW_label", "NMI", "PCR_batch", "graph_conn", "kBET"], :
].T

# 将指标缩放到0-1范围，便于比较和计算总分
metrics_scaled = (metrics_plot - metrics_plot.min()) / (metrics_plot.max() - metrics_plot.min())

# 计算总分 (权重可根据需求调整)
# 生物保留指标 (60%权重)
bio_metrics = ["ASW_label", "NMI"]
# 批次校正指标 (40%权重)
batch_metrics = ["PCR_batch", "graph_conn", "kBET"]

metrics_scaled["Bio conservation"] = metrics_scaled[bio_metrics].mean(axis=1)
metrics_scaled["Batch correction"] = metrics_scaled[batch_metrics].mean(axis=1)

# 最终总分
metrics_scaled["Overall score"] = 0.6 * metrics_scaled["Bio conservation"] + 0.4 * metrics_scaled["Batch correction"]

# 绘制条形图比较总分
metrics_scaled.sort_values("Overall score", ascending=False).plot.bar(
    y="Overall score", legend=None
)

★

小编注：  scib评估体系中的权重分配（例如，60%生物保留 vs 40%批次校正）是一个经验性的选择。在你的实际项目中，可以根据核心目标进行调整。如果你更关心的是构建一个无缝拼接的细胞图谱，可以适当提高“批次校正”的权重；反之，如果你的目标是发现新的、稀有的细胞亚群，那么“生物保留”的权重就应该更高。

关键评估指标解读（越高越好，kBET除外）：

  • 生物保留类
    • NMI: 比较整合前后的细胞类型聚类结果，衡量细胞类型标签的保留程度。
    • ASW_label: 轮廓宽度系数，衡量相同细胞类型的细胞是否紧密聚集。
  • 批次校正类
    • PCR_batch: 主成分回归，衡量整合后的主成分中，有多少方差仍然能被批次解释（该指标越高，说明批次解释力越弱，整合越好）。
    • graph_conn: 图连通性，衡量整合后，不同批次的同种细胞是否在图上连通。
    • kBET: （该指标越低越好 ）衡量局部区域细胞的批次混合程度。

从结果来看，scANVI和scVI在本次任务中取得了最高的综合得分，展现了强大的平衡能力。

方法选择的终极建议

没有一种整合方法是万能的。选择哪种方法，取决于你的数据规模、复杂程度和分析目标。

场景/需求 首选方法 备选/特点
大数据集 ( >100k细胞) scVI / scANVI 扩展性极佳，能有效处理海量数据。
简单批次校正  (如同技术平台、不同批次) Harmony, Seurat v3 速度快，稳定可靠，易于使用。
复杂数据整合  (跨技术、跨物种) scVI / scANVI, Scanorama 能够处理非线性的、复杂的批次效应。
需要校正后的表达矩阵 Seurat v3, scVI (可重建) 对于差异表达等下游分析至关重要。
快速探索性分析 BBKNN 速度最快，适合初步查看整合潜力。

一些实战经验：

  1. 先观察，再动手 ：在整合前，务必先通过PCA或UMAP评估原始数据的批次效应强度。微弱的效应可能无需校正。
  2. 多方比较 ：对于重要的分析项目，建议尝试2-3种不同类型的整合方法，并使用scib进行定量评估。
  3. 参数并非一成不变 ：所有方法都有可调参数，默认值是一个好的起点，但针对你的特定数据进行微调可能会获得更好的结果。
  4. 回归生物学 ：整合后，务必检查已知的细胞类型Marker基因的表达模式是否正常，确保整合没有引入或消除关键的生物学信号。

核心要点总结

  • 批次效应普遍存在，但在整合前必须想清楚：哪些是需要移除的技术噪音，哪些是需要保留的生物学差异。
  • 方法各有所长 ：
    • 深度学习 (scVI/scANVI)  ：处理复杂、大规模数据的王者。
    • 线性嵌入 (Seurat/Harmony)  ：应对简单任务的稳定之选。
    • 基于图 (BBKNN)  ：追求速度的探索性分析利器。
  • 评估是金标准 ：永远不要只相信你的眼睛（UMAP）。使用scib等工具进行定量评估，从“批次校正”和“生物保留”两个维度综合判断。
  • 没有万能的方法 ：最好的方法永远是最适合你当前数据特点、分析目标和计算资源 的那一个。

后记

scib是Theis Lab开发的一个评估整合效果的Python包，包含了多种常用的标准（比如图连接性、标签一致性等）。很妙的一点是它可以自定义/评估批次效应和样本本身的生物学多样性的占比，从而更好的平衡批次效应的去除程度。有兴趣的读者可以去GitHub（theislab/scib）和文章Benchmarking atlas-level data integration in single-cell genomics了解更多细节。

Workflow

来源：scib GitHub repository

大家怎么看待去批次？什么工具效果比较好，什么工具容易过拟合？欢迎在评论区分享你的经验！本章英文原文及完整代码，可以在后台回复scbp12 获取。

Bioinfo自习室交流群欢迎各位生信爱好者加入：

----------

单细胞领域大牛编写的最强硬核单细胞教程来啦！想要追更学习笔记，请给我们加个✨，不错过最新推送：

Single Cell Tutorial