单细胞分析（11）——scRNA-seq数据整合
单细胞 RNA-seq 数据整合：Seurat Integration and Harmony
1. 研究背景
在单细胞 RNA 测序（scRNA-seq）研究中，批次效应（batch effect） 是不可忽视的问题。不同样本来源（如多个实验室、不同测序平台、不同患者）可能会导致非生物学因素的影响，从而影响数据分析的准确性。
之前单独写过Harmony去除批次，为了更好地整合多个样本，这次使用以下两种方法进行批次校正：
- SCTransform + Seurat Integration
- SCTransform + Harmony
Luecken et al. (2022) 研究对 Seurat 方法的评价
在 Luecken et al. (2022)《Benchmarking atlas-level data integration in single-cell genomics》中，Seurat v3 被评估为适合小规模数据整合，尤其适用于细胞类型已知、批次效应较轻的数据 。然而，该研究指出，Seurat Integration 在高度异质性数据中可能会过度校正，从而丢失部分生物学变异。因此，对于跨多个实验室的大规模数据，可能需要结合其他方法（如 Harmony）。
SCTransform vs Harmony：核心对比
| 评估标准 | SCTransform + Seurat Integration | SCTransform + Harmony |
|---|---|---|
| 批次校正能力 | 强 ，适用于较小批次间的整合 | 适中 ，适用于大规模数据 |
| 生物学变异保留 | 高 ，适合免疫细胞、癌症数据等 | 较好 ，但可能保留部分批次效应 |
| 计算时间 | 较慢 （计算锚点消耗资源） | 较快 （基于 PCA 校正） |
| 适合的应用 | 细胞互作、DEG 分析 | 大规模数据整合、UMAP 计算 |
2. 数据预处理
数据筛选与清洗
- 合并不同样本的数据集
- 剔除低质量细胞 （如线粒体基因表达比例过高的细胞）
- 保留目标样本 （如某个特定处理组）
# 读取多个数据集
sample1 <- readRDS('sample1.rds')
sample2 <- readRDS('sample2.rds')
sample3 <- readRDS('sample3.rds')
# 合并样本
data_combined <- merge(sample1, y = c(sample2, sample3), add.cell.ids = c("S1", "S2", "S3"))
r
3. SCTransform + Seurat Integration
本方法借鉴于 A comprehensive single-cell breast tumor atlas defines epithelial and immune heterogeneity and interactions predicting anti-PD-1 therapy response 论文，应用于复杂肿瘤微环境的批次整合。
1️⃣ SCT 归一化
# 按来源拆分数据集
sample_list <- SplitObject(data_combined, split.by = "source")
# 设置多线程加速
library(future)
plan("multicore", workers = 4) # 根据硬件配置调整
options(future.globals.maxSize = 5 * 1024^3) # 允许更大数据处理（5GB）
# SCTransform 归一化
for (i in 1:length(sample_list)) {
sample_list[[i]] <- SCTransform(sample_list[[i]], verbose = TRUE, vars.to.regress = c("nCount_RNA", "percent.mt"))
}
r
2️⃣ 计算整合锚点
features <- SelectIntegrationFeatures(object.list = sample_list, nfeatures = 1500)
sample_list <- PrepSCTIntegration(object.list = sample_list, anchor.features = features)
integration_anchors <- FindIntegrationAnchors(object.list = sample_list, normalization.method = "SCT",
anchor.features = features, dims = 1:30)
data_integrated <- IntegrateData(anchorset = integration_anchors, normalization.method = "SCT")
r
为什么 FindIntegrationAnchors() 计算非常耗时？
- 计算多个数据集之间的共同特征基因匹配（anchor features）。
- 基于 PCA 计算最近邻（nearest neighbors）以建立整合参考框架。
- 大数据集（10W+ 细胞）可能需要几个小时计算，推荐使用多线程（future）加速。
3️⃣ 选择参考数据集
reference_samples <- c("S1", "S2")
integration_anchors <- FindIntegrationAnchors(object.list = sample_list, reference = reference_samples,
normalization.method = "SCT", anchor.features = features, dims = 1:30)
r
4️⃣ PCA + UMAP 降维
DefaultAssay(data_integrated) <- "integrated"
data_integrated <- RunPCA(data_integrated, npcs = 50, verbose = FALSE)
data_integrated <- RunUMAP(data_integrated, dims = 1:30, reduction = "pca")
r
什么时候用 “integrated” ，什么时候用 “SCT”？
| 使用数据 | 选择 Assay |
|---|---|
| 批次校正后的 UMAP、聚类分析 | integrated |
| 差异表达分析（DEG）、基因通路分析 | SCT |
4. SCTransform + Harmony
# SCTransform 归一化
data_combined <- SCTransform(data_combined, verbose = TRUE, vars.to.regress = c("nCount_RNA", "percent.mt"))
# PCA 降维
data_combined <- RunPCA(data_combined)
# Harmony 进行批次校正
data_combined <- RunHarmony(data_combined, group.by.vars = "source", max.iter.harmony = 50, theta = 3)
# UMAP
data_combined <- RunUMAP(data_combined, reduction = "harmony", dims = 1:30)
r
5. 结果对比
| 方法 | 批次校正效果 | 适合下游分析 | 计算速度 |
|---|---|---|---|
| SCTransform + Seurat Integration | 优秀，批次效应小，细胞类型分布清晰 | 适合 DEG、细胞互作分析 | 较慢 |
| SCTransform + Harmony | 一般，批次效应仍然存在 | 适合 UMAP 聚类分析，但不适合 DEG | 更快 |
6. 结论与推荐
✅ SCTransform + Seurat Integration 在批次校正方面更优，适用于 DEG 和细胞互作分析。
✅ SCTransform + Harmony 适用于大规模数据，但仍可能存在轻微批次效应，适用于 UMAP 聚类分析。
✅ 如果批次效应较小，建议 Seurat Integration；如果数据规模大且计算资源有限，Harmony 可能是更快的选择。
📌 对于单细胞 RNA-seq 研究，应该在比较去批次效果后，具体进行选择 🚀