Seurat中使用Harmony进行批次效应校正的最佳实践

概述

在单细胞RNA测序数据分析中，批次效应是一个常见的技术变异来源。Seurat作为单细胞分析的主流工具，提供了多种批次效应校正方法，其中Harmony因其高效和易用性而广受欢迎。本文将详细介绍在Seurat工作流中整合使用Harmony进行数据整合的方法和注意事项。

Harmony集成的基本流程

1. 数据准备与预处理

首先需要将原始数据按照批次变量进行分割。例如，如果批次信息存储在"Method"列中：

obj[["RNA"]] <- split(obj[["RNA"]], f = obj$Method)

然后进行标准的预处理步骤：

• 数据标准化(NormalizeData)
• 寻找高变基因(FindVariableFeatures)
• 数据缩放(ScaleData)
• PCA降维(RunPCA)

2. 使用IntegrateLayers进行Harmony整合

预处理完成后，可以直接调用IntegrateLayers函数进行Harmony整合：

obj <- IntegrateLayers(
  object = obj, 
  method = HarmonyIntegration,
  orig.reduction = "pca", 
  new.reduction = "harmony",
  verbose = FALSE
)

3. 下游分析

整合完成后，可以进行UMAP降维和聚类分析：

obj <- RunUMAP(obj, reduction = "harmony", dims = 1:30, reduction.name = "umap.harmony")

4. 数据清理

最后建议进行数据清理以释放内存：

obj <- JoinLayers(obj)
obj[["RNA"]]$scale.data <- NULL

直接使用harmony包的替代方案

除了通过Seurat的IntegrateLayers接口，也可以直接调用harmony包：

# 合并数据集
obj_1 = merge(x=o1, y=o2)

# 标准预处理
obj_1 = NormalizeData(obj_1)
obj_1 = FindVariableFeatures(obj_1)
obj_1 = ScaleData(obj_1)
obj_1 = RunPCA(obj_1)

# 直接调用RunHarmony，指定批次变量
obj_1 = RunHarmony(obj_1, c("dataset","flowcell"))

# 下游分析
obj_1[["RNA"]] = JoinLayers(obj_1[["RNA"]])
obj_1 = FindNeighbors(obj_1, dims = 1:20, reduction = "harmony")
obj_1 = FindClusters(obj_1)
obj_1 = RunUMAP(obj_1, dims = 1:20, reduction = "harmony")

方法比较与选择

1. IntegrateLayers方法：

   • 更符合Seurat的工作流
   • 自动处理批次信息
   • 适合Seurat v5的层级数据结构

2. 直接RunHarmony方法：

   • 更灵活，可以指定多个批次变量
   • 适合需要精细控制的情况
   • 适合熟悉harmony包的用户

注意事项

1. 批次变量的选择至关重要，应该基于实验设计和技术因素确定。

2. 在整合前后都应该检查数据质量，可以通过批次混合指标或可视化评估整合效果。

3. 对于大型数据集，Harmony通常比CCA等方法更高效，但仍需注意内存使用。

4. 整合后的数据通常用于降维和聚类，但不建议直接用于差异表达分析。

通过合理应用Harmony整合方法，可以有效地消除单细胞数据中的批次效应，使不同实验条件下的细胞能够被准确比较和分析。