Harmony 批次整合 — 快速高效的单细胞批次效应校正工具¶

# Python 版本（推荐，和 Scanpy 无缝整合）
conda activate scanpy_env
pip install harmonypy                  # 安装 harmonypy（Python 版 Harmony）

# 验证安装
python -c "import harmonypy; print(harmonypy.__version__)"

# Scanpy 内置接口（自动调用 harmonypy）
# 已包含在 Scanpy 的 sc.external.pp 模块中，无需额外安装

# R 版本（原版 Harmony，性能更优）
# install.packages('harmony')   # R 版本

import scanpy as sc                    # Scanpy 主包
import harmonypy as hm                # harmonypy 直接接口

# ── 方法1：Scanpy 内置 Harmony 接口（推荐）─────────────
adata = sc.read_h5ad('multi_batch.h5ad')

# 先做基础预处理
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=3000)
sc.pp.scale(adata)
sc.tl.pca(adata, n_comps=50)           # 必须先做 PCA

# 用 Harmony 校正批次（在 PCA 空间进行）
sc.external.pp.harmony_integrate(
    adata,
    key='batch',                       # 批次信息列名（adata.obs 中的列）
    basis='X_pca',                     # 输入的 PCA 嵌入（默认值）
    adjusted_basis='X_pca_harmony'     # 输出的校正后嵌入（新键名）
)

# 使用校正后的嵌入构建邻居图
sc.pp.neighbors(
    adata,
    use_rep='X_pca_harmony'            # 用 Harmony 校正后的空间
)
sc.tl.umap(adata)                      # UMAP 可视化
sc.tl.leiden(adata, resolution=0.5)    # 聚类

# 对比批次校正效果
sc.pl.umap(
    adata,
    color=['batch', 'leiden'],         # 查看批次分布是否混合良好
    ncols=2
)

import scanpy as sc
import harmonypy as hm
import pandas as pd
import matplotlib.pyplot as plt

# ── 完整多批次整合流程 ────────────────────────────────

# 1. 读取多个样本数据（假设已合并到一个 h5ad 文件）
# 每个样本有 'sample' 列标注批次信息
adata = sc.read_h5ad('merged_samples.h5ad')
print(f"样本信息：\n{adata.obs['sample'].value_counts()}")

# 2. 统一预处理（所有样本用相同流程）
adata.layers['counts'] = adata.X.copy()       # 保存原始计数
sc.pp.normalize_total(adata, target_sum=1e4)  # 标准化
sc.pp.log1p(adata)                             # log 变换

# 3. 选高变基因（按批次分别计算，再取并集）
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=3000,
    batch_key='sample',                        # 考虑批次的高变基因选择
    subset=False                               # 不直接过滤，标记即可
)
# 只使用高变基因做 PCA
adata_hvg = adata[:, adata.var.highly_variable].copy()
sc.pp.scale(adata_hvg, max_value=10)
sc.tl.pca(adata_hvg, n_comps=50)

# 将 PCA 结果复制回原始 adata
adata.obsm['X_pca'] = adata_hvg.obsm['X_pca']
adata.uns['pca'] = adata_hvg.uns['pca']

# 4. 校正前的 UMAP（查看批次效应）
sc.pp.neighbors(adata, use_rep='X_pca')
sc.tl.umap(adata)
sc.pl.umap(adata, color='sample', title='校正前（批次效应明显）')

# 5. Harmony 批次校正
sc.external.pp.harmony_integrate(
    adata,
    key='sample',                              # 按样本批次校正
    max_iter_harmony=20,                       # 最大迭代次数
    random_state=42                            # 固定随机种子
)

# 6. 校正后的 UMAP（应该批次混合在一起）
sc.pp.neighbors(adata, use_rep='X_pca_harmony', n_neighbors=20)
sc.tl.umap(adata, random_state=42)
sc.tl.leiden(adata, resolution=0.5)

# 7. 对比校正前后效果
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sc.pl.umap(
    adata, color='sample',
    title='Harmony 校正后', ax=axes[0], show=False
)
sc.pl.umap(
    adata, color='leiden',
    title='细胞聚类（校正后）', ax=axes[1], show=False
)
plt.tight_layout()
plt.savefig('harmony_integration.pdf', dpi=150, bbox_inches='tight')

# 8. 评估批次混合质量（LISI 指数，越高越好）
# 使用 scib 包评估（可选）
try:
    import scib
    lisi = scib.me.lisi(
        adata,
        batch_key='sample',
        label_key='leiden'
    )
    print(f"批次 LISI 分数：{lisi['iLISI'].mean():.3f}（越接近批次数越好）")
    print(f"细胞类型 LISI：{lisi['cLISI'].mean():.3f}（越接近1越好）")
except ImportError:
    print("安装 scib 包可评估整合质量：pip install scib")

# 9. 保存结果
adata.write_h5ad('harmony_integrated.h5ad')
print("Harmony 批次整合完成！")

# 安装
# pip install harmonypy scanpy

# 核心三行代码
sc.tl.pca(adata, n_comps=50)                        # 先 PCA
sc.external.pp.harmony_integrate(adata, key='batch') # Harmony 校正
sc.pp.neighbors(adata, use_rep='X_pca_harmony')      # 用校正后嵌入

# 多协变量校正（同时校正批次和供体）
sc.external.pp.harmony_integrate(adata, key=['batch', 'donor'])

# 批次校正方法对比
# Harmony     → 快，轻量，推荐首选
# scVI        → 深度学习，适合大数据集（见323教程）
# BBKNN       → 批次感知近邻图，保守
# scArches    → 参考映射（见331教程）