scRNA-seq 多样本整合分析教程 (Seurat / Harmony)

作者: SeekGene

时长: 17 分钟

字数: 4.1k 字

更新: 2026-02-27

阅读: 0 次

3' 转录组 5' + 免疫组库 ATAC + RNA 双组学 FFPE 单细胞转录组 Notebooks 全序列转录组 分析指南 数据整合 甲基化 + RNA 双组学 空间转录组

环境准备

R 包加载

请选择 common_r 这个环境进行该整合教程的学习

#加载必要的R包
suppressPackageStartupMessages({
  library(Seurat)
  library(dplyr)
  library(ggplot2)
  library(patchwork)
  library(harmony)
})

# 设置随机种子
set.seed(1234)

# 设置Seurat选项（注意：8000 * 1024^2 实际上是8GB）
options(future.globals.maxSize = 8000 * 1024^2)  # 8GB

my36colors <-c( '#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87',
                '#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658',
                '#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398',
                '#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
                '#712820', '#DCC1DD', '#CCE0F5', '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
                '#968175', "#6495ED", "#FFC1C1",'#f1ac9d','#f06966','#dee2d1','#6abe83','#39BAE8','#B9EDF8','#221a12',
                '#b8d00a','#74828F','#96C0CE','#E95D22','#017890')

数据读取

本教程提供两种不同形式的数据输入方式，以满足不同用户的数据获取需求，选择适合自己的一种方式即可：

云平台 RDS 文件读取

数据特点：

• RDS 文件是标准的 Seurat 对象文件
• 数据已经过预处理和多个样本合并
• 可直接用于后续的下游分析，也可以提取表达矩阵，重新整合

适用场景：

• 当您无法获得标准的filtered_feature_bc_matrix表达矩阵时
• 希望整合云平台现有数据进行学习时
• 需要快速重新进行 scRNA-seq 数据整合分析时

注意事项：

• 具体挂载数据和 rds 文件的读取，请参照 jupyter 使用教程

例如下列项目数据/home/demo-SeekGene-com/workspace/data/AY1752565399550/

input <- readRDS("/home/demo-seekgene-com/workspace/data/AY1752565399550/input.rds")
meta <- read.table("/home/demo-seekgene-com/workspace/data/AY1752565399550/meta.tsv", header=TRUE, sep="\t", row.names = 1)
data <- AddMetaData(input,meta)
data <- CreateSeuratObject(counts = input@assay$RNA@counts, meta.data=meta)
seurat_list <- SplitObject(data, split.by = "Sample")
#删除多余的数据，清内存
rm(data,input)
gc()

标准 filtered_feature_bc_matrix 文件读取

适用场景：

• 当您拥有标准的基因表达矩阵文件时
• 希望自主完成 scRNA-seq 数据的多样本整合和批次矫正时
• 需要进行完整的从原始数据到整合分析的工作流程时

注意：

• 请保证样本与样本之间的文件结构如下：

scRNA-seq 数据（基因表达矩阵）
├── S127/
│ ├── filtered_feature_bc_matrix/
│ │ ├── barcodes.tsv.gz（细胞条形码）
│ │ ├── features.tsv.gz（peaks 列表）
│ │ └── matrix.mtx.gz（稀疏计数矩阵）
├── S44R/
│ ├── filtered_feature_bc_matrix/
│ │ ├── barcodes.tsv.gz
│ │ ├── features.tsv.gz
│ │ └── matrix.mtx.gz

# 定义样本名称
sample_names <- c('S127', 'S44R')
# 创建一个空的list来存储每个样本的Seurat对象
seurat_list <- list()

for (sample in sample_names) {
  # 构建文件路径
  data_path <- file.path(sample, 'filtered_feature_bc_matrix')

  # 读取表达矩阵文件
  counts <- Read10X(data.dir = data_path)
  
  # 创建Seurat对象
  rna_obj <- CreateSeuratObject(
    counts = counts,
    project = sample,
    min.cells = 3,      # 过滤在少于3个细胞中表达的基因
    min.features = 200   # 过滤表达基因少于200个的细胞
  )
  
  # 添加样本标识
  rna_obj$Sample <- sample
  rna_obj$orig.ident <- sample
  # 将Seurat对象添加到list中
  seurat_list[[sample]] <- rna_obj
  # 清除rna_obj对象
  rm(rna_obj)
  gc()
}

质量控制

在整合前，需要先把各样本数据分别"清洗干净"：

• QC（质量控制）：剔除低质量细胞（如线粒体占比过高、背景或者双细胞等）

质控指标计算

为每个样本计算关键的质控指标，用于后续的细胞过滤。

常用 QC 指标：

• nCount_ATAC：总 ATAC 计数，极低或极高都需警惕异常。
• percent.mt: 线粒体基因比例（通常<20%）
• nFeature_RNA: 总 ATAC 特征数（通常在 200-10000 之间）
• nCount_RNA: 总 ATAC 特征数（通常在 200-10000 之间）

# 对seurat_list中的每个对象进行质量控制
seurat_list <- lapply(seurat_list, function(x) {
  x[["percent.mt"]] <- PercentageFeatureSet(x, pattern = "^MT-")
  return(x)
})

质控指标可视化

用小提琴图查看各 QC 指标的分布/异常点，确定合适的阈值：

建议：

• 观察是否存在明显的长尾或双峰分布
• 尝试多组阈值并比较下游聚类/UMAP 是否更清晰

options(repr.plot.width = 15, repr.plot.height = 5)
for (sample_name in names(seurat_list)) {
  seurat_obj <- seurat_list[[sample_name]]

  # 小提琴图展示质控指标分布
  p1 <- VlnPlot(
    object = seurat_obj,
    features = c('nFeature_RNA', 'nCount_RNA', 'percent.mt'),
    pt.size = 0.1,
    ncol = 5
  )
  print(p1+plot_annotation(title = sample_name))
}

低质量细胞过滤

根据质控指标过滤低质量细胞。具体阈值应根据数据特征进行调整。具体过滤阈值参考上面的小提琴图分布情况。

# 使用lapply进行质量控制过滤
seurat_list <- lapply(seurat_list, function(x) {
  # 记录过滤前的细胞数量
  cells_before <- ncol(x)
  
  # 质量控制过滤
  x <- subset( x,subset = nFeature_RNA > 200 &     # 基因数量大于200
             nFeature_RNA < 8000 &    # 基因数量小于5000
             nCount_RNA > 500 &       # UMI数量大于500
             nCount_RNA < 30000 &     # UMI数量小于30000
             percent.mt < 20         # 线粒体基因比例小于20%
  )
  
  # 记录过滤后的细胞数量
  cells_after <- ncol(x)
  
  # 输出过滤信息
  cat('过滤完成: 过滤前', cells_before, '个细胞，过滤后', cells_after, '个细胞\n')
  
  return(x)
})

output

过滤完成: 过滤前 13719 个细胞，过滤后 12722 个细胞
过滤完成: 过滤前 11830 个细胞，过滤后 11122 个细胞

多样本合并

在 scRNA-seq 数据分析流程中，多样本合并是整合分析的重要前置步骤。

合并操作的目的：

• 将多个样本数据整合到同一个 Seurat 对象中
• 为后续的批次矫正和整合分析做准备
• 便于与批次矫正后的结果进行对比分析

重要说明：

• 此步骤仅进行简单的数据合并，尚未进行批次效应矫正
• 合并后的对象包含所有样本的原始数据
• 将合并的数据进行标准化、特征选择、PCA 降维和 UMAP 可视化等预处理步骤

# 合并所有样本
suppressWarnings({
  suppressMessages({
      obj_merge <- merge(seurat_list[[1]], seurat_list[-1], merge.data = FALSE)
      obj_merge <- NormalizeData(obj_merge)
      obj_merge <- FindVariableFeatures(obj_merge, nfeatures = 2000)
      obj_merge <- ScaleData(obj_merge)
      obj_merge <- RunPCA(obj_merge)
        })
})

数据整合

现在我们将多个样本的 scATAC-seq 数据进行整合,并对整合的数据进行降维聚类。scATAC-seq 数据整合主要有两种常用方法：

整合方法选择

CCA（Canonical Correlation Analysis）整合：

• 基于典型相关分析的整合方法
• 通过寻找样本间的共同变异模式进行整合
• 适用于样本间差异较大的情况
• Seurat 包的经典整合方法

Harmony 整合：

• 基于迭代聚类的快速整合方法
• 直接在降维空间中校正批次效应
• 计算效率高，适用于大规模数据
• 保留更多生物学变异信息

方法选择建议

• Harmony：同平台不同样本 scATAC-seq 数据；样本数量较多或数据量大，Harmony 方法更快。
• CCA：运行时间较久，通常批次效应严重的情况下推荐选 CCA，如跨平台的 scATAC-seq 数据

harmony 方法进行整合

注意：harmony 和 CCA 选择其中一种方法进行批次矫正

# 定义harmony整合函数
integrate_harmony <- function(obj) {
  DefaultAssay(obj) <- "RNA"
  obj <- RunHarmony(obj, "Sample")  
  obj <- RunUMAP(obj, reduction = "harmony",
                  dims = 1:30)
  #obj <- RunTSNE(obj, reduction = "harmony",dims = 1:30,check_duplicates = FALSE)
  # RNA数据聚类
  obj <- FindNeighbors(object = obj, reduction = 'harmony', dims = 1:30)
  obj <- FindClusters(object = obj, verbose = FALSE, algorithm = 3,resolution = 0.5)
  return(obj)
}

#执行harmony整合分析
suppressWarnings({
  suppressMessages({
      obj_integrated = integrate_harmony(obj_merge)
      })
    })

CCA 方法进行整合

注意：harmony 和 CCA 选择其中一种方法进行批次矫正

# 切换到RNA assay
suppressWarnings({
  suppressMessages({
      objs <- lapply(seurat_list, function(x) {
          DefaultAssay(x) <- "RNA"
          return(x)
      })

      cat('合并后细胞总数:', ncol(obj_merge), '\n')
      # 定义CCA整合函数
      integrate_cca <- function(objs) {
          objs <- lapply(objs, function(x) {
            x <- NormalizeData(x, verbose = FALSE)
            x <- FindVariableFeatures(x,nfeatures = 2000,selection.method = "vst")
            return(x)
          })
          features <- Seurat::SelectIntegrationFeatures(object.list = objs)
          anchors <- Seurat::FindIntegrationAnchors(
              object.list = objs,
              anchor.features = features
          ) 
          obj <- Seurat::IntegrateData(anchorset = anchors)
          DefaultAssay(obj) <- "integrated"
          obj <- Seurat::ScaleData(obj, verbose = FALSE) %>%
          RunPCA(verbose = FALSE)
          #对整合后的数据进行UMAP降维
          obj <- RunUMAP(obj, reduction = "pca", dims = 1:30)
          # RNA数据聚类
          obj <- FindNeighbors(object = obj, dims = 1:30)
          obj <- FindClusters(object = obj, verbose = FALSE, algorithm = 3,resolution = 0.5)
          #obj <- RunTSNE(obj, reduction = "pca", dims = 1:30, check_duplicates = FALSE)
          return(obj)
          }
          # 执行CCA整合
          obj_integrated <- integrate_cca(seurat_list)

  })
})

output

合并后细胞总数: 23844

整合效果评估

比较整合前后的结果，评估整合效果。

评估指标：

• 样本混合程度：不同样本的细胞在 UMAP 图中的分布
• 批次效应去除：技术重复样本的聚集情况
• 生物学信号保留：已知细胞类型的分离度

cat('开始可视化整合结果...', Sys.time(), '\n')
# 对比整合前后效果
suppressWarnings({
  suppressMessages({
      obj_merge <- RunUMAP(obj_merge, reduction = "pca", dims = 1:30)
      })
    })
p1 <- DimPlot(obj_merge, group.by = "Sample") +
  ggtitle("Before integrate")

p2 <- DimPlot(obj_integrated, group.by = "Sample") +
  ggtitle("After integrate")

# 保存对比图
pdf("integration_comparison.pdf", width = 16, height = 12)
print(p1 + p2)
dev.off()

options(repr.plot.width = 20, repr.plot.height = 8)
print(p1 + p2)

output

开始可视化整合结果... 1755238468

pdf: 2

细胞类型注释

基于 marker 基因表达进行细胞类型注释。先对各聚类做差异表达以获取候选 marker，再结合文献与数据库进行判读。

定义 marker 基因

整理常见细胞类型的参考 marker 集合（来源可包括文献、PanglaoDB、CellMarker 等），用于后续可视化与判读。

eye_marker_integrated <- list(
  # ========== 光感受器细胞 ==========
  "Rod_Photoreceptors" = c("RHO", "PDE6B", "CNGB1", "PDE6A", "NR2E3", "REEP6",
                          "RCVRN", "SAG", "NEB", "SLC24A3", "TRPM1"),
  "Cone_Photoreceptors" = c("ARR3", "GNAT2", "OPN1SW", "TRPM3"),

  # ========== 视网膜神经元 ==========
  "Bipolar_Cells" = c("VSX1", "OTX2", "GRM6", "PRKCA",
                     "DLG2", "NRXN3", "RBFOX3", "FSTL4"),
  "Amacrine_Cells" = c("GAD1", "GAD2", "C1QL2", "TFAP2B", "SLC32A1"),
  "Horizontal_Cells" = c("ONECUT1", "LHX1", "CALB1", "NFIA"),
  "Retinal_Ganglion_Cells" = c("RBPMS", "THY1", "NEFL", "POU4F2"),

  # ========== 胶质细胞 ==========
  "Muller_Glia" = c("SLC1A3", "RLBP1", "SOX9", "CRYAB"),
  "Astrocytes" = c("GFAP", "AQP4", "S100B"),
  "Microglia" = c("TMEM119", "C1QA", "AIF1", "CX3CR1", "CD74",
                 "PTPRC", "CSF1R", "CCL3L1", "SPP1", "P2RY12"),

  # ========== 血管与支持细胞 ==========
  "Endothelial" = c("FLT1", "PODXL", "PLVAP", "KDR", "EGFL7",
                   "CLDN5", "PECAM1", "CDH5"),
  "Pericytes" = c("PDGFRB", "RGS5", "CSPG4", "STAB1"),

  # ========== 上皮细胞 ==========
  "RPE" = c("RPE65", "BEST1", "PMEL", "TYRP1", "DCT", "SLC38A11"),
  #"Corneal_Epithelial" = c("KRT12", "KRT3"),
  "Corneal_Endothelial" = c("SLC4A11", "COL8A2", "ATP1A1"),
  "Conjunctival_Epithelial" = c("KRT13", "KRT19", "MUC5AC"),
  "Retinal_Progenitors" = c("VSX2", "PAX6","SOX2", "HES1", "NOTCH1"),  # 神经视网膜增殖期

  # ========== 晶状体细胞 ==========
  "Lens_Epithelial" = c("CRYAA", "BFSP1", "MIP"),
  "Lens_Fiber" = c("CRYGS", "LIM2", "FN1"),

  # ========== 其他细胞 ==========
  "Melanocytes" = c("TYR", "MLANA", "MITF"),
  #"Erythrocytes" = c("HBB", "HBA1", "HBA2"),
  "ECM" = c("COL1A1", "COL3A1", "COL12A1", "COL1A2", "COL6A3")#,
  #"Others" = c("TTN", "CLCN5", "DCC", "MIAT")
)

Marker 基因表达可视化

在不同聚类中展示候选 marker 的表达模式（小提琴图/点图/FeaturePlot），辅助确认类型归属。

cat('开始marker基因表达可视化...', Sys.time(), '\n')
# 创建DotPlot展示marker基因表达
DefaultAssay(obj_integrated)="RNA"
p_dot <- DotPlot(
  obj_integrated,
  features = eye_marker_integrated,
  cols = c("#f8f8f8","#ff3472"),
  dot.scale = 8
) +
  RotatedAxis() +
  ggtitle("Marker genes expression") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 保存DotPlot
pdf("marker_genes_dotplot.pdf", width = 16, height = 8)
print(p_dot)
dev.off()
options(repr.plot.width = 26, repr.plot.height = 8)
print(p_dot)

cat('marker基因表达可视化完成!', Sys.time(), '\n')

output

开始marker基因表达可视化... 1755238806

pdf: 2

marker基因表达可视化完成! 1755238807 

p3 <- DimPlot(obj_integrated,group.by = "seurat_clusters",label = T, cols = my36colors)+ggtitle("Clusters")
options(repr.plot.width = 10, repr.plot.height = 8)
print(p3)

细胞类型注释

结合聚类层面差异基因与参考 marker，综合注释细胞类型；必要时细化到亚型或回退调整分辨率。

cat('开始细胞类型注释...', Sys.time(), '\n')

# 基于聚类结果进行细胞类型注释（需要根据实际的marker基因表达情况调整）
# 这里提供一个示例，实际使用时需要根据DotPlot结果进行调整
celltype_mapping <- c(
  "0" = "Retinal_Progenitors",
  "1" = "ECM",
  "2" = "ECM",
  "3" = "PRE",
  "4" = "Bipolar_Cells",
  "5" = "ECM",
  "6" = "Retinal_Progenitors",
  "7" = "Bipolar_Cells",
  "8" = "ECM",
  "9" = "Conjunctival_Epithelial",
  "10" = "ECM",
  "11" = "ECM",
  "12" = "Bipolar_Cells",
  "13" = "Doublets",
  "14" = "Endothelial",
  "15" = "Conjunctival_Epithelial",
  "16" = "Rod_Photoreceptors",
  "17" = "Astrocytes",
  "18" = "Microglia",
  "19" = "Muller_Glia"
)

# 应用细胞类型注释
obj_integrated$celltype <- recode(
  obj_integrated$seurat_clusters,
  !!!celltype_mapping
)

output

开始细胞类型注释... 1755239393

细胞类型注释可视化

以 UMAP 着色或按样本分面展示注释结果，并可统计各细胞类型的样本比例用于对比分析。

options(repr.plot.width = 10, repr.plot.height = 8)
# 细胞类型UMAP可视化
p1 <- DimPlot(
  obj_integrated,
  reduction = "umap",
  group.by = "celltype",
  label = TRUE,
  label.size = 3,
  cols = my36colors
) +
  theme(legend.position = "bottom")

# 按样本分组展示细胞类型分布
p2 <- DimPlot(
  obj_integrated,
  reduction = "umap",
  group.by = "celltype",
  split.by = "Sample",
  cols = my36colors,
  ncol = 2
) 

# 保存细胞类型注释图
pdf("celltype_annotation.pdf", width = 16, height = 12)
print(p1)
print(p2)
dev.off()

print(p1)
options(repr.plot.width = 20, repr.plot.height = 8)
print(p2)

pdf: 2

结果保存

保存整合后的对象与关键图表，建议同时保存中间对象以便复现与后续复查。

# 保存整合后的Seurat对象
saveRDS(obj_integrated, file = "scRNA_scRNA_integrated.rds")

总结

本教程演示了 scRNA-seq 多样本整合的完整流程，涵盖：

• 数据读取：导入 10x 格式表达矩阵，创建 Seurat 对象
• 质量控制：计算质控指标并筛除低质量细胞
• 数据预处理：归一化与高变基因识别
• 多样本合并：汇总样本并保留来源信息
• 数据整合：使用 CCA/RPCA 消除批次效应
• 降维与聚类：PCA、UMAP 及图聚类
• 细胞类型注释：基于 marker 的判读与可视化

最佳实践：

• 动态调整质控与整合参数（如阈值、dims、锚点数）
• 结合生物学先验与公共数据库进行注释
• 按阶段保存对象，确保分析可追溯与复现

后续分析建议：

• 组织类型偏好性分析
• 差异富集分析
• 拟时序分析
• 细胞互作分析

sessionInfo()

output

R version 4.3.3 (2024-02-29)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Debian GNU/Linux 12 (bookworm)

Matrix products: default
BLAS/LAPACK: /jp_envs/envs/common/lib/libopenblasp-r0.3.29.so; LAPACK version 3.12.0

locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C

time zone: Asia/Shanghai
tzcode source: system (glibc)

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] future_1.40.0 harmony_1.2.3 Rcpp_1.0.14 patchwork_1.3.0
[5] ggplot2_3.5.2 dplyr_1.1.4 SeuratObject_4.1.4 Seurat_4.4.0
[9] repr_1.1.7

loaded via a namespace (and not attached):
[1] deldir_2.0-4 pbapply_1.7-2 gridExtra_2.3
[4] rlang_1.1.5 magrittr_2.0.3 RcppAnnoy_0.0.22
[7] spatstat.geom_3.3-6 matrixStats_1.5.0 ggridges_0.5.6
[10] compiler_4.3.3 png_0.1-8 vctrs_0.6.5
[13] reshape2_1.4.4 stringr_1.5.1 pkgconfig_2.0.3
[16] crayon_1.5.3 fastmap_1.2.0 labeling_0.4.3
[19] promises_1.3.2 ggbeeswarm_0.7.2 purrr_1.0.4
[22] jsonlite_2.0.0 goftest_1.2-3 later_1.4.2
[25] uuid_1.2-1 spatstat.utils_3.1-3 irlba_2.3.5.1
[28] parallel_4.3.3 cluster_2.1.8.1 R6_2.6.1
[31] ica_1.0-3 stringi_1.8.7 RColorBrewer_1.1-3
[34] spatstat.data_3.1-6 reticulate_1.42.0 parallelly_1.43.0
[37] spatstat.univar_3.1-2 lmtest_0.9-40 scattermore_1.2
[40] IRkernel_1.3.2 tensor_1.5 future.apply_1.11.3
[43] zoo_1.8-14 R.utils_2.13.0 base64enc_0.1-3
[46] sctransform_0.4.1 httpuv_1.6.15 Matrix_1.6-5
[49] splines_4.3.3 igraph_2.0.3 tidyselect_1.2.1
[52] abind_1.4-5 spatstat.random_3.3-3 codetools_0.2-20
[55] miniUI_0.1.1.1 spatstat.explore_3.4-2 listenv_0.9.1
[58] lattice_0.22-7 tibble_3.2.1 plyr_1.8.9
[61] withr_3.0.2 shiny_1.10.0 ROCR_1.0-11
[64] ggrastr_1.0.2 evaluate_1.0.3 Rtsne_0.17
[67] survival_3.8-3 polyclip_1.10-7 fitdistrplus_1.2-2
[70] pillar_1.10.2 KernSmooth_2.23-26 plotly_4.10.4
[73] generics_0.1.3 sp_2.2-0 IRdisplay_1.1
[76] munsell_0.5.1 scales_1.3.0 globals_0.16.3
[79] xtable_1.8-4 glue_1.8.0 lazyeval_0.2.2
[82] tools_4.3.3 data.table_1.17.0 pbdZMQ_0.3-13
[85] RANN_2.6.2 leiden_0.4.3.1 Cairo_1.6-2
[88] cowplot_1.1.3 grid_4.3.3 tidyr_1.3.1
[91] colorspace_2.1-1 nlme_3.1-168 beeswarm_0.4.0
[94] vipor_0.4.7 cli_3.6.4 spatstat.sparse_3.1-0
[97] viridisLite_0.4.2 uwot_0.2.3 gtable_0.3.6
[100] R.methodsS3_1.8.2 digest_0.6.37 progressr_0.15.1
[103] ggrepel_0.9.6 htmlwidgets_1.6.4 farver_2.1.2
[106] R.oo_1.27.0 htmltools_0.5.8.1 lifecycle_1.0.4
[109] httr_1.4.7 mime_0.13 MASS_7.3-60.0.1