单细胞去批次方法基础说明：Combat、Harmony、BBKNN

🔥 为什么要去批次？

在单细胞或空间转录组数据中，不同批次会产生技术噪声，让图上出现：

同一种细胞被分成多团（因为不同批次）
聚类受技术因素影响，而不是生物差异
UMAP “一团一个批次”

去批次就是：

去掉技术差异，保留真实生物差异。

🧩 一句话总结三者的作用

方法	一句话理解	在做什么？
Combat	把不同批次的表达值校正到同一水平	修正表达矩阵（每个基因的均值和方差）
Harmony	把同类细胞拉到一起，让不同批次对齐	修正 PCA 空间（细胞坐标）
BBKNN	强制每个细胞在所有批次中都找邻居	修正邻域图（KNN → 影响 UMAP/聚类）

三者的本质差异：

Combat → 修表达矩阵
Harmony → 修坐标（PCA embedding）
BBKNN → 修图结构（KNN graph）

1️⃣ Combat - 校准不同批次的表达值

想象你有两个温度计：

一个偏高 5 度
一个偏低 3 度

Combat 的作用就是：

校准温度计，让不同批次的基因表达读数一致。

它会对每个基因估计：

批次间的“整体偏高或偏低”（均值差）
批次间的“波动大小是否不同”（方差差）

然后全部统一到同一尺度。

✔ 擅长

简单、稳定
修正表达矩阵本身
适合批次差异不太“结构化”的情况（线性偏移）
批次数较少
样本间细胞类型分布相似

❌ 不擅长

复杂异质数据（多 donor、多细胞亚群）
空间转录组
希望整合 embedding（UMAP/聚类）

🌟 基础知识

Combat 基于 **贝叶斯分层模型 (Bayesian Hierarchical Model)**：

[
Y_{gij} = \alpha_g + X_{ij}\beta_g + \gamma_{gj} + \delta_{gj}\epsilon_{gij}
]

Combat 分别调整每个基因的：

均值偏移（location）
方差差异（scale）

📌 特点

特点	描述
线性模型	假设批次效应是偏移与缩放
基于 gene level 校正	对每个基因独立调整
不能对齐 embedding	不会调整 PCA/UMAP
简单稳定	常用于 bulk-like 数据

🔧 Scanpy 使用

sc.pp.combat(adata, key='batch')

🔧 Seurat 使用

library(sva)
combat_data <- ComBat(dat=as.matrix(expr), batch=batch, mod=NULL)

2️⃣ Harmony - 对齐细胞，使同类细胞靠在一起

假设：

batch A 的 T 细胞在左边
batch B 的 T 细胞在右边

但它们明明是同一种细胞！

Harmony 做的事情是：

调整 PCA 空间，让同类细胞（不同批次）对齐并重叠。

它不会改表达值，只会改“细胞在图上的坐标”。

✔ 擅长

不改变表达矩阵（DEG 安全）
细胞类型异质性大
能对齐不同批次、不同 donor 的同类细胞(多 batch、多个 donor )
需要统一 embedding 用于 UMAP/Louvain
效果强、应用最广

❌ 不擅长

极复杂的非线性结构（如强空间结构）
多 modality 整合不如 SCVI
极端非线性结构不如 BBKNN
空间坐标需要先 PCA 再整合

🌟 基础知识

Harmony 是 PCA embedding 层面的迭代校正算法：

利用 PCA 空间作为起点
聚类（soft clustering）
在 cluster 内校正 batch centroid
多次迭代直到 batch alignment 达到平衡

形式：

[
Z^{(t+1)} = Z^{(t)} - \lambda (Z^{(t)} - \mu_{b})
]

📌 特点

特点	描述
在 PCA 层面校正	不修改原始表达矩阵
对复杂多 batch 数据效果强	能自动对齐同类细胞
速度快，可扩展	大数据处理优秀
支持多个 batch 因子	donor + batch + chemistry

🔧 Scanpy 使用

sc.tl.pca(adata)
sc.external.pp.harmony_integrate(adata, key='batch')

🔧 Seurat 使用

obj <- RunHarmony(obj, group.by.vars="batch")

3️⃣ BBKNN (Batch Balanced KNN) - 让每个细胞都均衡地认识其他批次的邻居

普通 KNN：

一个 batch 特别大，就会支配邻居结构
小 batch 的细胞容易成为孤岛
UMAP 会按 batch 分成几坨

BBKNN 的做法是：

强制每个细胞在每个 batch 都找一定数量的邻居。

例如：

有 3 个 batch
你设置 k=30

BBKNN 会自动从每个 batch 各找 10 个邻居。

这样：

批次之间被“连接”在一起
细胞类型会更混合
图形结构更自然

✔ 擅长

Scanpy 主流程优先
多 batch 且分布不均匀
空间转录组（强烈推荐）
非线性能力强,处理非线性结构复杂数据
超快（适合几十万细胞）
空间组（Visium/Stereo）强烈推荐使用

❌ 不擅长

需要精确的 embedding（BBKNN不改PCA）
高维对齐任务（prefer Harmony/SCVI）

🌟 基础知识

BBKNN 基于 邻域图 (KNN graph) 去批次：

为每个细胞从 每个 batch 选固定数量的邻居
强制 batch 间“平衡连接”
修改邻域图，而不是修改表达或 PCA

例如 k=20，有两个 batch，则每批取10个最近邻。

📌 特点

特点	描述
图层面校正	修改邻域连接结构
强非线性能力	适合复杂数据分布
超快	可处理百万级单细胞
不影响原始表达矩阵	DEG 安全

🔧 Scanpy 使用

import bbknn
bbknn.bbknn(adata, batch_key='batch')
sc.tl.umap(adata)

🔥 三者横向对比总结

方法	校正层面	修改表达矩阵	非线性能力	空间组表现	适用场景
Combat	gene level	✔️	弱	较弱	简单批次效应
Harmony	PCA embedding	❌	中等	良好	多 donor & 多批次
BBKNN	邻域图	❌	强	强烈推荐	Scanpy/空间组

方法	类比
Combat	校准不同型号的秤，让重量读数一致
Harmony	拍全家福时，把同一类人（比如所有小孩）重新排到一起
BBKNN	班主任硬性规定：每个人必须和不同班级的人交朋友

📌 推荐方案（根据你的常用场景）

空间转录组（Visium/Stereo-seq）：
→ BBKNN > Harmony > Combat
常规 scRNA-seq（donor 多、异质性大）：
→ Harmony > BBKNN > Combat
简单线性批次效应：
→ Combat

场景	推荐方法
空间转录组（Visium/Stereo）	⭐⭐⭐ BBKNN 最优选
普通 scRNA-seq，多 donor、多批次	⭐⭐⭐ Harmony 最优选
批次效应很简单（只是线性偏移）	Combat 即可
数据很大（>30 万细胞）	BBKNN
需要做 DEG 分析	Harmony / BBKNN（不改变表达矩阵）

一个完整的 Scanpy 分析管线模板，在同一套预处理基础上，分别走三条分支：

Combat 去批次
Harmony 去批次
BBKNN 去批次

共同的基础预处理流程

import scanpy as sc
import pandas as pd
import numpy as np

# ============ 基础配置 ============
sc.set_figure_params(figsize=(5, 5))
sc.settings.verbosity = 3  # 日志等级

# ============ 读入数据 ============
# 例子：多个 batch 已经合并到一个 AnnData 中，并且在 adata.obs['batch'] 中标记
# 你可以根据自己实际情况改成 read_h5ad / read_10x_mtx / read_visium 等
adata = sc.read_h5ad("merged_raw.h5ad")

# adata.obs 应该至少有一列 batch 信息，例如：
# adata.obs['batch'] = ...
print(adata)

# ============ 基础 QC（可按需简化/替换） ============
sc.pp.calculate_qc_metrics(adata, inplace=True)

# 这里仅示意：按基因数和线粒体比例过滤细胞
adata = adata[adata.obs['n_genes_by_counts'] > 200, :].copy()
adata = adata[adata.obs['pct_counts_mt'] < 10, :].copy()

# 过滤低表达基因
sc.pp.filter_genes(adata, min_cells=3)

# ============ 归一化 & log ============
# 按每个细胞总 counts 归一化到 1e4
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# ============ 高变基因 ============
sc.pp.highly_variable_genes(
    adata,
    flavor="seurat_v3",
    n_top_genes=3000,
    subset=True,    # 只保留高变基因
)

# ============ 保存原始 log 表达矩阵（用于 DEG 等） ============
adata.raw = adata

# 从这里开始，kaishi走三条不同的去批次分支

分支一：Combat 去批次（sc.pp.combat）

思路：在 log 表达矩阵上做 Combat 校正，然后再 scale + PCA + 邻域 + UMAP

# ========= Combat 分支 =========
adata_combat = adata.copy()

# 1. Combat 直接作用于 adata.X（当前是 log-normalized 数据）
#    adata_combat.obs['batch'] 必须存在
sc.pp.combat(adata_combat, key='batch')

# 2. 标准化 / 缩放
sc.pp.scale(adata_combat, max_value=10)

# 3. PCA
sc.tl.pca(adata_combat, svd_solver='arpack')

# 4. 邻域图
sc.pp.neighbors(adata_combat, n_neighbors=15, n_pcs=30)

# 5. UMAP / 聚类
sc.tl.umap(adata_combat)
sc.tl.leiden(adata_combat, resolution=0.5, key_added='leiden_combat')

# 6. 可视化
sc.pl.umap(
    adata_combat,
    color=['batch', 'leiden_combat'],
    wspace=0.4,
    title=['Combat: batch', 'Combat: clusters']
)

# 如需保存
adata_combat.write_h5ad("adata_combat_corrected.h5ad")

要点：

Combat 修改了 adata_combat.X（表达矩阵）；DEG 分析建议用校正后的矩阵或者用 adata.raw 原始矩阵对比两种情况。
Combat 适合批次效应比较“线性”的情况。

分支二：Harmony 去批次（sc.external.pp.harmony_integrate）

思路：在 PCA 空间上用 Harmony 调整 embedding，不改表达矩阵。

需要额外安装：pip install harmonypy

# ========= Harmony 分支 =========
import scanpy.external as sce  # Harmony 的接口在 scanpy.external

adata_harmony = adata.copy()

# 1. scale + PCA
sc.pp.scale(adata_harmony, max_value=10)
sc.tl.pca(adata_harmony, svd_solver='arpack')

# 2. Harmony 整合
#    结果会写入 adata_harmony.obsm['X_pca_harmony']
sce.pp.harmony_integrate(
    adata_harmony,
    key='batch'  # 使用哪一列 obs 作为批次因子
)

# 3. 用 Harmony PCA 结果建邻域图
sc.pp.neighbors(
    adata_harmony,
    use_rep='X_pca_harmony',  # 使用 Harmony 校正过的 PCA
    n_neighbors=15,
    n_pcs=30
)

# 4. UMAP / 聚类
sc.tl.umap(adata_harmony)
sc.tl.leiden(adata_harmony, resolution=0.5, key_added='leiden_harmony')

# 5. 可视化
sc.pl.umap(
    adata_harmony,
    color=['batch', 'leiden_harmony'],
    wspace=0.4,
    title=['Harmony: batch', 'Harmony: clusters']
)

# 如需保存
adata_harmony.write_h5ad("adata_harmony_integrated.h5ad")

要点：

Harmony 只改 obsm['X_pca_harmony']，不改 .X，所以 DEG 可以安心用 adata_harmony.raw（和原始 adata.raw 一样）。
key='batch' 也可以换成 ['batch', 'donor'] 等更复杂组合（scanpy/harmony 有支持）。

分支三：BBKNN 去批次（bbknn.bbknn）

思路：用 BBKNN 替换 neighbors 步骤，只改邻域图，不改 PCA 和表达矩阵。

需要安装：pip install bbknn

# ========= BBKNN 分支 =========
import bbknn

adata_bbknn = adata.copy()

# 1. scale + PCA
sc.pp.scale(adata_bbknn, max_value=10)
sc.tl.pca(adata_bbknn, svd_solver='arpack')

# 2. 使用 BBKNN 构建 batch-balanced 邻域图
#    会覆盖 adata_bbknn.obsp['connectivities'] 和 ['distances']
bbknn.bbknn(
    adata_bbknn,
    batch_key='batch',   # 批次信息所在列
    n_pcs=30,
    neighbors_within_batch=10,  # 每个 batch 内部取多少邻居
    n_pcs_resid=0               # 一般默认即可
)

# 3. UMAP / 聚类
sc.tl.umap(adata_bbknn)
sc.tl.leiden(adata_bbknn, resolution=0.5, key_added='leiden_bbknn')

# 4. 可视化
sc.pl.umap(
    adata_bbknn,
    color=['batch', 'leiden_bbknn'],
    wspace=0.4,
    title=['BBKNN: batch', 'BBKNN: clusters']
)

# 如需保存
adata_bbknn.write_h5ad("adata_bbknn_integrated.h5ad")

要点：

BBKNN 不改 .X，也不改 PCA，只重建邻域图。
特别适合空间组 / 细胞数很大的单细胞数据。

bajiu

去批次的三种方法Combat、Harmony、BBKNN介绍

🔥 为什么要去批次？

🧩 一句话总结三者的作用

1️⃣ Combat - 校准不同批次的表达值

✔ 擅长

❌ 不擅长

🌟 基础知识

📌 特点

🔧 Scanpy 使用

🔧 Seurat 使用

2️⃣ Harmony - 对齐细胞，使同类细胞靠在一起

✔ 擅长

❌ 不擅长

🌟 基础知识

📌 特点

🔧 Scanpy 使用

🔧 Seurat 使用

3️⃣ BBKNN (Batch Balanced KNN) - 让每个细胞都均衡地认识其他批次的邻居

✔ 擅长

❌ 不擅长

🌟 基础知识

📌 特点

🔧 Scanpy 使用

🔥 三者横向对比总结

📌 推荐方案（根据你的常用场景）

一个完整的 Scanpy 分析管线模板，在同一套预处理基础上，分别走三条分支：

共同的基础预处理流程

分支一：Combat 去批次（sc.pp.combat）

分支二：Harmony 去批次（sc.external.pp.harmony_integrate）

分支三：BBKNN 去批次（bbknn.bbknn）