企业官网建设流程全解析

单细胞数据整合实战：用Harmony消除SCTransform处理后的批次效应

当你在分析来自不同实验批次或供体的单细胞RNA测序数据时，是否遇到过这样的困扰：明明使用了Seurat的SCTransform进行标准化，但聚类结果仍然明显受到批次影响？这种现象在整合多源数据时尤为常见，而Harmony算法正是解决这一痛点的利器。本文将带你从实战角度，一步步掌握如何在使用SCTransform后正确调用Harmony进行数据整合。

1. 理解批次效应与整合原理

批次效应是单细胞数据分析中的"隐形杀手"。即使采用最严格的实验控制，不同时间点、不同操作者或不同试剂批次产生的数据仍会存在系统性差异。这种技术性变异会掩盖真实的生物学信号，导致聚类结果出现偏差。

SCTransform作为Seurat中的新一代标准化方法，通过正则化负二项回归模型（regularized negative binomial regression）对UMI计数数据进行建模，能有效处理测序深度差异和基因表达量的离散特性。但它主要解决的是样本内部的技术变异，对于样本间的批次效应仍需专门的处理。

Harmony算法的核心思想是通过迭代优化，在保留生物学变异的同时消除批次效应。它基于以下关键假设：

• 相同细胞类型在不同批次中应该具有相似的表达谱
• 批次效应可以通过线性模型进行校正
• 校正后的数据应该最大化保留生物学相关的聚类结构

提示：Harmony特别适合处理那些批次效应与生物学信号存在部分重叠的情况，这是许多线性校正方法难以解决的问题。

2. 数据准备与SCTransform标准化

在开始整合前，我们需要确保数据已经过适当的预处理。以下是一个典型的工作流程：

# 加载必要的R包 library(Seurat) library(dplyr) library(harmony) # 假设我们已经有了四个供体的PBMC数据 pbmc_combined <- merge(pbmc1, y = c(pbmc2, pbmc3, pbmc4), add.cell.ids = c("Donor1", "Donor2", "Donor3", "Donor4"), project = "pbmc_combined")

执行SCTransform标准化时，有几个关键参数需要注意：

# 执行SCTransform标准化 pbmc_combined <- SCTransform(pbmc_combined, vars.to.regress = "percent.mt", return.only.var.genes = FALSE, verbose = TRUE)

参数说明：

• vars.to.regress：指定需要回归掉的协变量（如线粒体基因百分比）
• return.only.var.genes：设为FALSE以确保保留所有基因用于后续分析
• verbose：显示详细运行信息

3. Harmony整合的关键步骤

完成SCTransform后，我们需要先进行PCA降维，然后再应用Harmony。以下是具体操作：

# 运行PCA pbmc_combined <- RunPCA(pbmc_combined, npcs = 50, verbose = FALSE) # 运行Harmony整合 pbmc_combined <- RunHarmony(pbmc_combined, group.by.vars = "donor.number", plot_convergence = TRUE, assay.use = "SCT")

关键参数解析：

注意：确保group.by.vars指定的批次变量已经存在于对象的元数据中。如果不存在，需要先添加到meta.data。

4. 结果评估与可视化

整合效果需要通过多种方式进行验证。以下是一些常用的评估方法：

4.1 UMAP可视化对比

# 整合前的UMAP pbmc_combined <- RunUMAP(pbmc_combined, reduction = "pca", dims = 1:30) DimPlot(pbmc_combined, group.by = "donor.number", reduction = "umap") # 整合后的UMAP pbmc_combined <- RunUMAP(pbmc_combined, reduction = "harmony", dims = 1:30) DimPlot(pbmc_combined, group.by = "donor.number", reduction = "umap")

4.2 混合度指标量化

除了视觉评估，我们还可以计算定量指标：

# 计算局部批次混合度 library(kBET) batch <- pbmc_combined$donor.number harmony_emb <- Embeddings(pbmc_combined, "harmony")[,1:30] batch_estimate <- kBET(harmony_emb, batch)

评估指标解读：

• kBET接受率：值越高表示批次混合越好（理想值>0.7）
• ASW（平均轮廓宽度）：衡量批次效应去除程度（值越小越好）
• LISI（局部逆辛普森指数）：评估细胞邻域中批次的多样性

5. 常见问题与解决方案

在实际应用中，我们可能会遇到各种挑战。以下是几个典型问题及其解决方法：

5.1 整合后生物学信号丢失

现象：细胞类型间的区分度降低，聚类质量下降
可能原因：

• 过度校正：Harmony将真实的生物学差异误认为批次效应
• PCA维度选择不当：使用了包含太多噪声的PCs

解决方案：

# 尝试调整theta参数（控制校正强度） pbmc_combined <- RunHarmony(pbmc_combined, group.by.vars = "donor.number", theta = 1, # 默认2，减小可降低校正强度 assay.use = "SCT") # 重新选择PCA维度 ElbowPlot(pbmc_combined, ndims = 50)

5.2 收敛速度慢或失败

现象：Harmony迭代次数过多或不收敛
可能原因：

• 批次间差异过大
• 数据质量存在问题

解决方案：

• 检查原始数据质量（如每个批次的细胞数、基因检出率）
• 尝试增加max.iter.harmony参数（默认10）
• 确保SCTransform参数设置合理

5.3 如何处理多个批次变量

当存在多个批次来源（如实验日期+测序平台）时：

# 合并多个批次变量为一个新变量 pbmc_combined$batch.combined <- paste(pbmc_combined$date, pbmc_combined$platform, sep = "_") # 使用合并后的变量进行整合 pbmc_combined <- RunHarmony(pbmc_combined, group.by.vars = "batch.combined", assay.use = "SCT")

6. 高级应用与优化技巧

对于更复杂的分析场景，我们可以考虑以下进阶方法：

6.1 与WNN整合框架结合

# 在Harmony整合后构建WNN图 pbmc_combined <- FindMultiModalNeighbors( pbmc_combined, reduction.list = list("harmony"), dims.list = list(1:30), modality.weight.name = "RNA.weight" ) # 基于WNN的UMAP pbmc_combined <- RunUMAP(pbmc_combined, nn.name = "weighted.nn", assay = "RNA", reduction.name = "wnn.umap")

6.2 处理超大单细胞数据集

当数据量极大时（>100万细胞），可以考虑：

# 启用Harmony的近似PCA模式 pbmc_combined <- RunHarmony(pbmc_combined, group.by.vars = "donor.number", approx_pca = TRUE, assay.use = "SCT") # 或者先进行子采样 set.seed(123) cells.use <- sample(colnames(pbmc_combined), size = 1e5) pbmc.sub <- subset(pbmc_combined, cells = cells.use)

6.3 整合多模态数据

对于CITE-seq等多模态数据，建议：

1. 对RNA数据单独进行SCTransform+Harmony
2. 对ADT数据使用DSB或CLR标准化
3. 使用Seurat的WNN框架进行跨模态整合

# 对ADT数据进行CLR标准化 DefaultAssay(pbmc_combined) <- "ADT" pbmc_combined <- NormalizeData(pbmc_combined, normalization.method = "CLR", margin = 2) # 切换回SCT assay进行Harmony整合 DefaultAssay(pbmc_combined) <- "SCT" pbmc_combined <- RunHarmony(pbmc_combined, group.by.vars = "donor.number", assay.use = "SCT")

在实际项目中，我发现当处理来自不同实验室的数据时，将theta参数设为1.5、同时使用前30个PCs通常能取得最佳平衡。对于特别敏感的细胞亚群（如过渡态细胞），建议在整合后手动检查这些群体的分布情况。