747. 单细胞批次整合Harmony/scVI¶

要去掉的：技术差异（批次效应）
  → 不同天实验导致的系统偏差
  → 不同平台的测序深度差异

要保留的：生物学差异（真实信号）
  → 疾病vs健康的基因表达差异
  → 不同组织的细胞类型差异

难点：这两种差异混在一起！
  过度校正 → 抹掉了真实的生物学差异
  校正不足 → 批次效应还在

Harmony的核心思想：
1. 先做PCA降维
2. 在PCA空间中，找到批次混合不好的地方
3. 迭代地"推"细胞，让每个聚类里各批次的比例均衡
4. 不改变原始表达矩阵，只修正PCA坐标

类比：
  不同班级（批次）的学生混坐在教室里
  Harmony就是"调座位"，让每个小组（聚类）里每班人数差不多
  但不改变学生的成绩（原始表达）

# ===== Python Harmony =====
import scanpy as sc  # 导入scanpy
import harmonypy  # 导入harmony
import matplotlib.pyplot as plt

# 读取多批次数据（已合并但未整合的）
adata = sc.read_h5ad("merged_data.h5ad")  # 读取合并数据
# adata.obs["batch"]列包含批次信息

# 标准预处理
sc.pp.normalize_total(adata, target_sum=1e4)  # 归一化
sc.pp.log1p(adata)  # log转换
sc.pp.highly_variable_genes(adata, n_top_genes=3000,
                             batch_key="batch")  # 按批次找高变基因
adata = adata[:, adata.var.highly_variable]  # 只保留高变基因
sc.pp.scale(adata, max_value=10)  # 标准化
sc.tl.pca(adata, n_comps=50)  # PCA降维

# ===== 运行Harmony =====
sc.external.pp.harmony_integrate(
    adata,
    key="batch",  # 批次变量名
    basis="X_pca",  # 在PCA空间上操作
    adjusted_basis="X_pca_harmony",  # 输出的校正后PCA
    max_iter_harmony=20  # 最大迭代次数
)

# 用校正后的PCA做下游分析
sc.pp.neighbors(adata, use_rep="X_pca_harmony")  # 用harmony后的PCA
sc.tl.umap(adata)  # UMAP降维
sc.tl.leiden(adata, resolution=0.5)  # Leiden聚类

# 对比整合前后
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# 整合前（用原始PCA）
sc.pp.neighbors(adata, use_rep="X_pca", key_added="neighbors_raw")
sc.tl.umap(adata, neighbors_key="neighbors_raw")
sc.pl.umap(adata, color="batch", title="Before Harmony",
           ax=axes[0], show=False)

# 整合后（用harmony PCA）
sc.pp.neighbors(adata, use_rep="X_pca_harmony", key_added="neighbors_harmony")
sc.tl.umap(adata, neighbors_key="neighbors_harmony")
sc.pl.umap(adata, color="batch", title="After Harmony",
           ax=axes[1], show=False)

plt.tight_layout()
plt.savefig("harmony_comparison.png", dpi=300)
plt.show()

# ===== R/Seurat中使用Harmony =====
library(Seurat)
library(harmony)

# 读取并合并多个样本的Seurat对象
seurat_list <- list(
  sample1 = readRDS("sample1.rds"),
  sample2 = readRDS("sample2.rds"),
  sample3 = readRDS("sample3.rds")
)

# 合并
merged <- merge(seurat_list[[1]], seurat_list[-1],
                add.cell.ids=names(seurat_list))

# 标准预处理
merged <- NormalizeData(merged)
merged <- FindVariableFeatures(merged)
merged <- ScaleData(merged)
merged <- RunPCA(merged, npcs=50)

# 运行Harmony
merged <- RunHarmony(
  merged,
  group.by.vars = "batch",  # 批次变量
  reduction = "pca",        # 输入降维
  assay.use = "RNA",        # 使用的assay
  plot_convergence = TRUE   # 显示收敛曲线
)

# 用Harmony结果做下游分析
merged <- FindNeighbors(merged, reduction="harmony", dims=1:30)
merged <- FindClusters(merged, resolution=0.5)
merged <- RunUMAP(merged, reduction="harmony", dims=1:30)

# 可视化
DimPlot(merged, group.by="batch", label=FALSE)
ggsave("harmony_umap.png", width=10, height=8, dpi=300)

scVI使用变分自编码器（VAE）：
1. 编码器：将高维表达数据压缩到低维潜在空间
2. 在潜在空间中，显式建模"批次效应"
3. 解码器：从潜在空间重建表达数据
4. 批次作为条件变量，在潜在空间中"消除"批次效应

优势：
  - 直接建模count数据（负二项分布）
  - 可以处理复杂的批次效应
  - 更好地保留生物学差异
  - 可以做差异表达、去卷积等下游分析

# ===== scVI分析流程 =====
import scvi  # 导入scvi-tools
import scanpy as sc

# 读取数据
adata = sc.read_h5ad("merged_data.h5ad")

# 预处理（scVI需要原始counts）
# 确保adata.X是原始counts或adata.layers["counts"]有counts
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=3000,
    batch_key="batch",  # 按批次计算高变基因
    flavor="seurat_v3",  # 使用Seurat v3方法
    layer="counts"  # 从counts层计算
)

# 设置scVI模型
scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",  # 使用原始counts
    batch_key="batch",  # 批次变量
    categorical_covariate_keys=["sample"],  # 额外的分类协变量
    continuous_covariate_keys=["percent_mt"]  # 连续协变量（如线粒体比例）
)

# 创建并训练模型
model = scvi.model.SCVI(
    adata,
    n_latent=30,  # 潜在空间维度
    n_layers=2,   # 神经网络层数
    gene_likelihood="nb"  # 负二项分布（适合count数据）
)

# 训练
model.train(
    max_epochs=200,  # 最大训练轮数
    early_stopping=True,  # 早停
    train_size=0.9  # 90%训练，10%验证
)

# 查看训练历史
model.history["elbo_train"].plot()  # 绘制训练损失曲线

# 获取整合后的潜在表示
adata.obsm["X_scVI"] = model.get_latent_representation()

# 用scVI表示做下游分析
sc.pp.neighbors(adata, use_rep="X_scVI")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

# 可视化
sc.pl.umap(adata, color=["batch", "cell_type"], ncols=2)

# 保存模型（以后可以加载继续分析）
model.save("scvi_model/")

# ===== 定量评估整合效果 =====
import scib  # 导入scib（单细胞整合基准测试工具）

# 安装：pip install scib

# 计算整合指标
# 1. 批次混合指标（高=整合好）
ari = scib.metrics.ari(adata, cluster_key="leiden",
                        label_key="cell_type")
nmi = scib.metrics.nmi(adata, cluster_key="leiden",
                        label_key="cell_type")
asw_batch = scib.metrics.silhouette_batch(
    adata, batch_key="batch", label_key="cell_type",
    embed="X_scVI"
)

# 2. 生物学保留指标（高=保留好）
asw_cell = scib.metrics.silhouette(
    adata, label_key="cell_type", embed="X_scVI"
)

print(f"ARI (聚类与真实标签的一致性): {ari:.3f}")
print(f"NMI (归一化互信息): {nmi:.3f}")
print(f"ASW_batch (批次混合): {asw_batch:.3f}")
print(f"ASW_cell (细胞类型分离): {asw_cell:.3f}")

# 理想结果：ASW_batch高（批次混合好）+ ASW_cell高（类型分离好）

# ===== 批次整合速查 =====

# Harmony（Python）
sc.external.pp.harmony_integrate(adata, key="batch")
sc.pp.neighbors(adata, use_rep="X_pca_harmony")

# Harmony（R/Seurat）
merged <- RunHarmony(merged, group.by.vars="batch")
merged <- FindNeighbors(merged, reduction="harmony")

# scVI（Python）
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="batch")
model = scvi.model.SCVI(adata, n_latent=30)
model.train(max_epochs=200)
adata.obsm["X_scVI"] = model.get_latent_representation()

# 整合方法选择指南：
# 批次效应轻 → Harmony/BBKNN
# 批次效应重 → scVI/Scanorama
# 速度优先 → Harmony
# 精度优先 → scVI
# Seurat生态 → CCA/rpca

# 评估指标：
# ASW_batch: 高=批次混合好（去批次效应成功）
# ASW_cell: 高=细胞类型分离好（保留了生物学）