单细胞整合Harmony¶

假设你从3个医院收集了PBMC样本做单细胞测序：

不整合的问题：
  医院A的T细胞 ——和—— 医院B的T细胞
  虽然是同一种细胞，但在UMAP上可能分成两堆
  原因：不同的试剂批次、操作人员、测序平台等

整合的目的：
  把技术差异"抹掉"，让相同细胞类型聚在一起
  同时保留真正的生物差异（如疾病vs正常）

关键挑战：
  去太少 → 批次效应仍在
  去太多 → 生物差异也被消除了（过度校正）

Harmony的策略（迭代软聚类）：
1. 先在PCA空间中做软聚类
2. 计算每个cluster中不同批次的比例
3. 如果某batch在某cluster中占比异常 → 调整该batch细胞的坐标
4. 重复步骤1-3直到收敛
5. 输出校正后的低维嵌入（harmony坐标）

特点：只校正低维嵌入，不修改原始表达矩阵
优点：快速、保守、校准好（不太会过度校正）

# 安装Harmony
install.packages("harmony")  # 安装
library(harmony)               # 加载
library(Seurat)                # 加载Seurat

# 假设已有合并的多样本Seurat对象
# seurat_obj包含多个样本的数据，meta.data中有"batch"列

# 标准预处理
seurat_obj <- NormalizeData(seurat_obj)          # 归一化
seurat_obj <- FindVariableFeatures(seurat_obj)   # 找高变基因
seurat_obj <- ScaleData(seurat_obj)              # 缩放
seurat_obj <- RunPCA(seurat_obj, npcs = 30)      # PCA

# ===== 运行Harmony =====
seurat_obj <- RunHarmony(
  seurat_obj,
  group.by.vars = "batch",   # 批次变量名
  reduction = "pca",          # 基于PCA结果
  dims.use = 1:30,            # 使用的PC数
  max.iter.harmony = 20       # 最大迭代次数
)

# Harmony结果存储在reduction中
Reductions(seurat_obj)  # 查看可用的降维："pca", "harmony"

# 用Harmony嵌入做后续分析（不用PCA）
seurat_obj <- FindNeighbors(seurat_obj,
                             reduction = "harmony",  # 用harmony而非pca
                             dims = 1:30)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)
seurat_obj <- RunUMAP(seurat_obj,
                       reduction = "harmony",  # 基于harmony
                       dims = 1:30)

# 可视化：按批次着色（应该混合良好）
DimPlot(seurat_obj, group.by = "batch")  # 看批次是否混合
DimPlot(seurat_obj, group.by = "cell_type")  # 看细胞类型是否保留

# Seurat v5简化了整合接口
# 先按样本分层
seurat_obj[["RNA"]] <- split(seurat_obj[["RNA"]],
                              f = seurat_obj$batch)

# 标准预处理
seurat_obj <- NormalizeData(seurat_obj)
seurat_obj <- FindVariableFeatures(seurat_obj)
seurat_obj <- ScaleData(seurat_obj)
seurat_obj <- RunPCA(seurat_obj)

# 一行搞定Harmony整合
seurat_obj <- IntegrateLayers(
  seurat_obj,
  method = HarmonyIntegration,    # 使用Harmony方法
  orig.reduction = "pca",         # 基于PCA
  new.reduction = "harmony"       # 输出到"harmony"
)

# 合并层
seurat_obj <- JoinLayers(seurat_obj)

# 后续聚类和UMAP用harmony
seurat_obj <- FindNeighbors(seurat_obj, reduction = "harmony", dims = 1:30)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)
seurat_obj <- RunUMAP(seurat_obj, reduction = "harmony", dims = 1:30)

scVI（single-cell Variational Inference）：
1. 训练一个变分自编码器（VAE，深度学习模型）
2. 编码器把每个细胞的表达向量压缩成低维潜在变量
3. 在训练过程中，批次信息作为条件变量输入
4. 模型学会"什么是批次差异，什么是生物差异"
5. 输出的潜在空间就是去批次后的嵌入

特点：修改/重建表达矩阵（可以输出校正后的表达值）
优点：适合复杂整合任务、跨物种、多模态
缺点：需要GPU加速、参数较多、有时会过度校正

# scVI在Python中运行（通过scvi-tools）
# pip install scvi-tools

import scvi
import scanpy as sc

# 读取数据（AnnData格式）
adata = sc.read_h5ad("merged_data.h5ad")

# 设置scVI模型
scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",       # 使用原始计数
    batch_key="batch"     # 批次变量
)

# 训练模型
model = scvi.model.SCVI(adata, n_latent=30)  # 潜在维度30
model.train(max_epochs=400)                    # 训练400轮

# 获取潜在表示（去批次后的嵌入）
adata.obsm["X_scVI"] = model.get_latent_representation()

# 用scVI嵌入做聚类和UMAP
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"])

# 通过reticulate在R中调用scVI
library(reticulate)
library(Seurat)
library(SeuratDisk)

# 导出Seurat对象为h5ad
SaveH5Seurat(seurat_obj, filename = "data.h5seurat")
Convert("data.h5seurat", dest = "h5ad")

# 在Python中运行scVI（见上面的Python代码）
# 然后导入结果

# 或使用SeuratDisk导入scVI嵌入
adata <- ReadH5AD("scvi_result.h5ad")
scvi_embed <- adata@obsm$X_scVI

# 添加到Seurat对象
seurat_obj[["scvi"]] <- CreateDimReducObject(
  embeddings = scvi_embed,
  key = "scVI_",
  assay = "RNA"
)

1. 首先试Harmony（快速、保守、校准好）
2. 如果Harmony效果不好 → 试Seurat CCA
3. 如果涉及跨物种或很复杂的整合 → 用scVI
4. 如果有参考注释想迁移学习 → 用scANVI
5. 多方法比较取交集 → 最稳健的结果

# 方法1：目视检查
p1 <- DimPlot(seurat_obj, group.by = "batch") + ggtitle("By Batch")
p2 <- DimPlot(seurat_obj, group.by = "cell_type") + ggtitle("By Cell Type")
p1 | p2  # 并排显示

# 好的整合：batch颜色混合良好，cell_type颜色分离清晰

# 方法2：LISI分数（定量评估）
# install.packages("lisi")
library(lisi)

harmony_embed <- Embeddings(seurat_obj, "harmony")
meta <- seurat_obj@meta.data

lisi_scores <- compute_lisi(
  harmony_embed,
  meta,
  c("batch", "cell_type")  # 评估的变量
)
# iLISI（batch）: 越高越好（批次混合程度）
# cLISI（cell_type）: 越低越好（细胞类型保留程度）

# Harmony快速流程
NormalizeData() → FindVariableFeatures() → ScaleData()
→ RunPCA() → RunHarmony(group.by.vars="batch")
→ FindNeighbors(reduction="harmony")
→ FindClusters() → RunUMAP(reduction="harmony")

# Seurat v5整合流程
split() → NormalizeData() → FindVariableFeatures()
→ ScaleData() → RunPCA()
→ IntegrateLayers(method=HarmonyIntegration)
→ JoinLayers() → FindNeighbors() → FindClusters() → RunUMAP()

# 整合方法选择
简单批次: Harmony（首选）
中等复杂: Seurat CCA/RPCA
复杂/跨物种: scVI
有参考标签: scANVI

# 评估指标
iLISI: 批次混合程度（越高越好）
cLISI: 细胞类型保留（越低越好）
kBET: 批次混合统计检验
ASW: 平均轮廓系数