ATAC + RNA 多组学：基于 Signac 的多样本整合分析

作者: SeekGene

时长: 28 分钟

字数: 6.9k 字

更新: 2026-07-07

阅读: 0 次

1. 教程简介

本教程基于 Seurat + Signac + Harmony 等单细胞多组学分析工具，用于对 SeekArc 单细胞多组学数据（RNA + ATAC） 进行 多样本整合分析、批次效应校正、联合聚类 和 细胞类型注释。

该教程面向单个细胞同时包含 转录组表达信息 和 染色质开放性信息 的单细胞数据，通过分别整合多个样本的 RNA 与 ATAC 模态，结合多模态分析策略，实现以下核心分析目标：

多样本整合：对来自不同样本的单细胞多组学数据进行统一处理，降低批次差异对下游分析的影响。
双模态联合分析：同时利用 RNA 表达特征和 ATAC 开放染色质特征，更全面地刻画细胞状态与生物学异质性。
批次效应校正：提供基于 CCA 和 Harmony 的整合方式，用于校正不同样本之间的技术偏差。
聚类与降维可视化：基于 RNA（PCA）、ATAC（LSI）及 WNN 多模态融合三种降维策略，分别构建细胞相似性网络，实现非线性降维（UMAP/t-SNE）与多分辨率聚类，支持结果交叉验证。
细胞类型注释：结合 marker 基因表达模式和基因组区域开放情况，对聚类结果进行生物学解释和细胞类型标注。

#加载必要的R包
suppressPackageStartupMessages({
  library(Seurat)
  library(Signac)
  library(EnsDb.Hsapiens.v86)
  library(BSgenome.Hsapiens.UCSC.hg38)
  library(biovizBase)
  #library(BSgenome.Mmusculus.UCSC.mm10) #小鼠
  #library(EnsDb.Mmusculus.v79)
  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')

# --- 输入参数配置 ---

## file_path：输入的 Seurat 对象路径
file_path = "/path/to/sampleDir"

# integration_method：用于整合的算法，本教程提供 CCA 和 harmony 两种数据整合的方式，只能二选一，整合的目的是去除样本间批次效应
integration_method = "harmony" 

## sample_name：需要分析的样本名称。
##   - 多样本整合模式：填入多个样本名，用逗号分隔（例如："sample1,sample2,sample3"）。系统会自动触发多样本整合与批次校正。
sample_names = c('pbmc_1', 'pbmc_2')

2. 输入文件准备

2.1 输入文件要求

请保证每个样本均按照统一的目录结构进行整理。每个样本使用一个独立文件夹，文件夹名称建议直接使用样本 ID，例如 pbmc_1、pbmc_2 等。每个样本目录下应包含以下内容：

filtered_feature_bc_matrix：scRNA-seq 表达矩阵目录，包含 barcodes.tsv.gz、features.tsv.gz 和 matrix.mtx.gz 文件。
filtered_peaks_bc_matrix：scATAC-seq peak 开放矩阵目录，包含 barcodes.tsv.gz、features.tsv.gz 和 matrix.mtx.gz 文件。
{样本ID}_A_fragments.tsv.gz：ATAC 片段文件，用于后续构建 ChromatinAssay、计算 TSS 富集和核小体信号等分析。
{样本ID}_A_fragments.tsv.gz.tbi：ATAC 片段文件对应的索引文件。

2.2 目录结构示例

text

pbmc_1/
├── filtered_feature_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_peaks_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── pbmc_1_A_fragments.tsv.gz
└── pbmc_1_A_fragments.tsv.gz.tbi

pbmc_2/
├── filtered_feature_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_peaks_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── pbmc_2_A_fragments.tsv.gz
└── pbmc_2_A_fragments.tsv.gz.tbi

2.3 注意事项

所有样本的目录结构必须保持一致，否则程序在批量读取时可能报错。
sample_names 中填写的样本名称应与实际文件夹名称完全一致。
ATAC 片段文件及其索引文件必须成对存在，否则无法正常构建 ChromatinAssay。
RNA 和 ATAC 数据应来自同一批细胞或同一多组学实验体系，以保证后续联合分析的准确性。

3. 数据加载与预处理

在正式开展多样本整合分析之前，需要先完成基因组注释信息加载，以及每个样本 RNA 和 ATAC 数据的读取与对象构建。该步骤的目标是：为每个样本分别建立包含 RNA 和 ATAC 两种模态的 Seurat 对象，并统一存入列表中，供后续质量控制、整合和聚类分析使用。

3.1 获取基因注释信息

首先，通过 EnsDb 数据库读取参考基因组注释信息，并生成 annotation 对象。该注释信息包含基因位置、转录起始位点（TSS）等内容，可用于后续的 TSS 富集分析、基因活性分析和 peak 注释。

本示例使用人类注释数据库 EnsDb.Hsapiens.v86，并将染色体名称统一为带 chr 前缀的格式，同时指定参考基因组版本为 `hg38``。这样可以保证注释信息与 ATAC 数据的染色体命名保持一致，避免后续分析出错。

# 获取基因注释信息（静默处理警告和消息）
suppressWarnings({
  suppressMessages({
    annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
    seqlevels(annotation) <- paste0('chr', seqlevels(annotation))
    genome(annotation) <- 'hg38'
  })
})

3.2 数据读取与预处理

本步骤主要完成多样本 RNA 和 ATAC 数据的读取、基因组注释信息加载，以及多组学 Seurat 对象的构建，为后续质控、整合分析和聚类分析做好准备，包含以下环节：

读取多样本 RNA 数据：依次读取每个样本的 filtered_feature_bc_matrix 表达矩阵，并创建 RNA assay 的 Seurat 对象。
读取多样本 ATAC 数据：读取每个样本的 filtered_peaks_bc_matrix 开放矩阵及片段文件，构建 ATAC assay，并添加到对应的 Seurat 对象中。其中，构建 ATAC assay 时引入上述 annotation 基因组注释信息，可用于ATAC 相关下游分析。
保存多组学对象列表：将每个样本构建完成的 RNA + ATAC Seurat 对象保存到 seurat_list 中，供后续多样本整合分析使用。

# 创建空列表存储Seurat对象
seurat_list <- list()
# 循环读取每个样本的数据
for(sample in sample_names){
  cat('正在处理样本:', sample, '\n')
  
  # 构建文件路径 - 修复：将中文逗号改为英文逗号
  atac_path <- file.path(file_path, sample, 'filtered_peaks_bc_matrix')
  rna_path <- file.path(file_path, sample, 'filtered_feature_bc_matrix')
  frag_path <- file.path(file_path, sample, paste0(sample, '_A_fragments.tsv.gz'))
  
  # 读取RNA数据
  rna_counts <- Read10X(data.dir = rna_path)
  
  # 创建Seurat对象
  obj <- CreateSeuratObject(
    counts = rna_counts,
    assay = "RNA")
  
  obj$Sample <- sample
  rm(rna_counts)
  gc()
  
  # 读取ATAC数据
  atac_counts <- Read10X(data.dir = atac_path)
  atac_counts <- atac_counts[Matrix::rowSums(atac_counts > 0) >= 3, ]   
  
  # 注意：需要先定义 annotation
  # 例如：annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
  ChromatinAssay <- CreateChromatinAssay(
    counts = atac_counts,
    sep = c(":", "-"),
    fragments = frag_path,
    annotation = annotation  # 这里需要提前定义annotation
  )
  
  obj[["ATAC"]] <- ChromatinAssay
  rm(atac_counts, ChromatinAssay)
  gc()
  
  # 将对象添加到列表
  seurat_list[[sample]] <- obj
  cat('样本', sample, '处理完成\n')
}

4. 数据质量控制

4.1 质量指标计算

在完成多样本数据读取后，需要对每个样本分别进行质量控制，以去除低质量细胞，保证后续整合分析结果的可靠性。本步骤会计算 RNA 和 ATAC 两个模态的常用质控指标。

RNA 质控指标：

percent.mt：线粒体基因比例，通常用于评估细胞状态，过高可能提示低质量或受损细胞。

ATAC 质控指标：

TSS.enrichment：TSS 富集分数，用于评估开放染色质信号在转录起始位点附近的富集程度。
nucleosome_signal：核小体信号，用于反映片段分布特征，通常越低越好。

# 对每个样本进行质量控制
suppressWarnings({
  suppressMessages({
      for(i in names(seurat_list)){
          # RNA质控指标
          seurat_list[[i]][["percent.mt"]] <- PercentageFeatureSet(seurat_list[[i]], pattern = "^MT-")
          # ATAC质控指标
          DefaultAssay(seurat_list[[i]]) <- "ATAC"
          # 计算TSS富集分数
          seurat_list[[i]] <- TSSEnrichment(object = seurat_list[[i]], fast = FALSE) 
          # 计算核小体信号
          seurat_list[[i]] <- NucleosomeSignal(object = seurat_list[[i]])
      }
  })
})

说明： seekARC 双组学数据质控指标还包含 nCount_RNA、nCount_ATAC、nFeature_RNA、nFeature_ATAC 等指标，这些指标在前面构建 Seurat 对象的时候已经自动生成。

4.2 质量指标可视化

在完成质控指标计算后，可以通过小提琴图和散点图查看各项指标的分布情况，从而判断数据质量并选择合适的过滤阈值。

建议重点关注以下内容：

是否存在明显的异常值、长尾分布或双峰分布；
不同样本之间质控指标分布是否一致；
根据可视化结果合理调整过滤阈值，以获得更清晰的下游聚类和 UMAP 结果。

# 可视化各质控指标，该cell内容选择性执行，非必要
options(repr.plot.width = 16, repr.plot.height = 10)
suppressWarnings({
for (sample_name in names(seurat_list)) {
  seurat_obj <- seurat_list[[sample_name]]
  p1=DensityScatter(seurat_obj, x = 'nCount_ATAC', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)
  p2=VlnPlot(
      object = seurat_obj,
      features = c('nCount_ATAC', 'TSS.enrichment',  'nucleosome_signal',"nFeature_RNA", "nCount_RNA", "percent.mt"),
      pt.size = 0.1,
      ncol = 6
  )
 print(p1 / p2 + plot_annotation(title = sample_name))
}
})

4.3 低质量细胞过滤

根据前面计算得到的质控指标，对低质量细胞进行过滤，以去除异常细胞并提高后续分析结果的可靠性。具体过滤阈值需要结合样本实际情况进行设置，通常可参考质控指标的小提琴图和散点图分布结果进行调整。

# 质控过滤
for(i in names(seurat_list)){
  cat('正在对样本', i, '进行质控过滤...\n')
  
  cells_before <- ncol(seurat_list[[i]])
  
  seurat_list[[i]] <- subset(
    seurat_list[[i]],
    subset = nFeature_RNA > 200 &
             nFeature_RNA < 8000 &
             nCount_RNA > 500 &
             nCount_RNA < 30000 &
             percent.mt < 20 &
             nCount_ATAC > 500 &
             nCount_ATAC < 50000 &
             TSS.enrichment > 1 &
             nucleosome_signal < 1
  )
  
  cells_after <- ncol(seurat_list[[i]])
  cat('样本', i, '过滤完成: 过滤前', cells_before, '个细胞，过滤后', cells_after, '个细胞\n')
}

output

正在对样本 pbmc_1 进行质控过滤...n 样本 pbmc_1 过滤完成: 过滤前 8433 个细胞，过滤后 7515 个细胞
正在对样本 pbmc_2 进行质控过滤...n 样本 pbmc_2 过滤完成: 过滤前 7703 个细胞，过滤后 6805 个细胞

5. 计算样本间的共有 peaks

由于不同样本通常是分别进行 peak calling 的，因此各样本的 ATAC peak 集合并不完全一致。为了开展多样本整合分析，需要先将所有样本的 peaks 进行合并，构建统一的共有 peak 集合，并基于该集合重新生成每个样本的 ATAC 矩阵。

本步骤主要包括以下内容：

合并所有样本的 peaks：提取每个样本中的 peak 区间信息，并合并为一个统一的 peak 集合。
过滤 peak 长度：去除过短或过长的 peaks，通常保留长度在 20 bp 到 10 kb 之间的区间。
重新量化共有 peaks：以统一的 peak 集合为基础，对每个样本重新计算 ATAC counts，保证不同样本使用相同的特征空间。

# 提取所有样本的peaks
all_peaks <- lapply(seurat_list, function(x) {
  DefaultAssay(x) <- "ATAC"
  granges(x@assays$ATAC)
})

# 合并所有peaks
gr_list <- GRangesList(all_peaks)
all_granges <- unlist(gr_list, use.names = FALSE)
combined.peaks <- reduce(x = all_granges)

# 过滤peak长度
peakwidths <- width(combined.peaks)
common_peaks <- combined.peaks[peakwidths < 10000 & peakwidths > 20]

cat('合并后的peak数量:', length(common_peaks), '\n')

# 为每个样本创建统一的peak矩阵
seurat_list <- lapply(seurat_list, function(x) {
  combined_counts <- FeatureMatrix(
    fragments = Fragments(x@assays$ATAC),
    features = common_peaks,
    cells = colnames(x)
  )
  
  combined_peaks_assay <- CreateChromatinAssay(
    counts = combined_counts,
    fragments = Fragments(x@assays$ATAC),
    annotation = Annotation(x@assays$ATAC)
  )
  
  # 添加到Seurat对象中
  x[["combinedpeaks"]] <- combined_peaks_assay
  DefaultAssay(x) <- "combinedpeaks"
  x[["ATAC"]] <- NULL
  return(x)
})

6. 多样本合并

在完成各样本数据预处理后，需要先将多个样本合并到同一个 Seurat 对象中，作为后续批次校正和多组学整合分析的输入。

本步骤的主要目的包括：

将多个样本的数据统一整合到同一个对象中；
为后续的批次效应校正、降维和聚类分析做准备；
便于比较整合前后数据的差异。

需要注意的是，此步骤仅进行样本合并，并不涉及批次效应校正。合并后的对象仍保留各样本的原始特征信息，随后通常还需要进行标准化、特征选择、降维等基础预处理。

# 合并所有样本
suppressWarnings({
  suppressMessages({
      obj_merge <- merge(seurat_list[[1]], seurat_list[-1], merge.data = FALSE)
      DefaultAssay(obj_merge) <- "RNA" 
      obj_merge <- NormalizeData(obj_merge)
      obj_merge <- FindVariableFeatures(obj_merge, nfeatures = 2000)
      obj_merge <- ScaleData(obj_merge)
      obj_merge <- RunPCA(obj_merge)
      obj_merge <- RunUMAP(obj_merge, reduction = "pca", dims = 2:30,reduction.name="rnaumap")
      DefaultAssay(obj_merge) <- "combinedpeaks"
      obj_merge <- FindTopFeatures(obj_merge, min.cutoff = 10)
      obj_merge <- RunTFIDF(obj_merge)
      obj_merge <- RunSVD(obj_merge)
      obj_merge <- RunUMAP(obj_merge, reduction = "lsi", dims = 2:30,reduction.name="atacumap")
        })
})

P1=DimPlot(obj_merge, reduction = "rnaumap", group.by = "Sample")
P2=DimPlot(obj_merge, reduction = "atacumap", group.by="Sample")

options(repr.plot.width=17, repr.plot.height=8)
patchwork::wrap_plots(P1, P2, ncol = 2)

💡 Note
说明：当前阶段仅完成初步数据合并（scRNA 与 scATAC 数据分别按样本合并），尚未进行批次矫正，左侧是RNA数据合并后的降维图，右侧是ATAC数据合并后的降维图。图中展示的样本间分布差异可初步反映批次效应强弱：若同类型细胞因样本来源不同而明显分离，则表明存在批次效应，需进行后续矫正；若样本间细胞分布高度重叠，则可跳过批次矫正步骤。实际分析中，多数多样本数据存在不同程度的批次效应，通常需要进行矫正处理的。

7. 数据整合

在完成多样本合并后，需要进一步对不同样本之间的批次差异进行校正。对于 SeekArc 单细胞多组学数据，需要分别对 RNA 和 ATAC 模态各自进行整合，再用于后续的联合分析。

7.1 整合方法选择

本教程提供两种常用的批次校正方法：CCA 和 Harmony。

CCA（Canonical Correlation Analysis）整合：

基于 Seurat 的锚点整合思路，通过寻找不同样本之间的共同特征实现数据对齐；
适合样本间差异较明显、批次效应较强的情况；
计算开销相对较大，运行时间通常较长。

Harmony 整合：

在已有降维空间中直接进行批次校正，例如 RNA 的 PCA 空间或 ATAC 的 LSI 空间；
运行速度较快，适合样本数量较多或数据规模较大的场景；通常能够较好地保留生物学差异，同时降低批次效应；
更适合同平台、同类型数据的多样本整合分析。

注意事项

RNA 和 ATAC 需要分别进行整合，再用于后续的多模态联合分析。
实际分析中建议结合 UMAP 混样效果、聚类结果和已知生物学信息综合判断整合效果。

7.2 开始对 RNA 数据整合

if(integration_method == "CCA"){
    # 切换到RNA assay
    rnaintegratedumap="rnaintegratedumap"
    rnaintegratedtsne="rnaintegratedtsne"
    suppressWarnings({  
        suppressMessages({
            objs <- lapply(seurat_list, function(x) {
                DefaultAssay(x) <- "RNA"
                return(x)
            })
            # 定义CCA整合函数
            integrate_cca <- function(objs) {
                # 使用lapply替代for循环
                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 <- ScaleData(obj, verbose = FALSE)
                obj <- RunPCA(obj, verbose = FALSE)
                obj <- RunUMAP(obj, reduction = "pca",reduction.name=rnaintegratedumap, dims = 1:30)
                #obj <- RunTSNE(obj, reduction = "pca", reduction.name=rnaintegratedtsne,dims = 1:30, check_duplicates = FALSE)
                return(obj)
            }
            # 执行CCA整合
            obj_integrated <- integrate_cca(seurat_list)
        })
    })      
} else {
    # RNA数据整合函数
    rnaintegratedumap="rnaharmonyumap"
    rnaintegratedtsne="rnaharmonytsne"
    integrate_harmony <- function(obj) {
        DefaultAssay(obj) <- "RNA"
        obj <- RunHarmony(obj, "Sample")
        obj <- RunUMAP(obj, reduction = "harmony",reduction.name=rnaintegratedumap,dims = 1:30)
        #obj <- RunTSNE(obj, reduction = "harmony",reduction.name=rnaharmonytsne,dims = 1:30,check_duplicates = FALSE)
        return(obj)
    }
    suppressWarnings({
        suppressMessages({
            obj_integrated = integrate_harmony(obj_merge)
        })
    }) 
}

##可视化RNA数据去批次效果
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(obj_integrated, reduction = rnaintegratedumap, group.by = "Sample")
options(repr.plot.width=16, repr.plot.height=8)
DimPlot(obj_integrated, reduction = rnaintegratedumap, split.by = "Sample")

💡 Note
说明：此图展示经批次矫正后 RNA 数据的分布情况。评估批次效应去除效果时，应在相同细胞类型内观察：若同类型细胞已实现样本间混合分布，表明批次矫正有效；若同类型细胞仍按样本分离成簇，则批次去除不完全。请注意，评估需基于细胞类型层面，而非整体 UMAP 分布，因为样本间可能存在真实的生物学差异（如细胞类型组成或比例不同）。

7.3 ATAC 数据整合

if(integration_method == "CCA"){
    # 定义ATAC整合函数
    atacintegratedumap="atacintegratedumap"
    atacintegratedtsne="atacintegratedtsne"
    integrated_lsi="integrated_lsi"
    options(future.globals.maxSize = 20 * 1024^3)
    integrate_atac_cca <- function(object.list){
        object.list <- lapply(object.list, function(x) {
            DefaultAssay(x)="combinedpeaks"
            x=RunTFIDF(x)
            x=FindTopFeatures(x, min.cutoff = 5)
            x=RunSVD(x)
            return(x)
        })
        anchors <- FindIntegrationAnchors(
            object.list = object.list,
            anchor.features = rownames(obj_merge),
            reduction = "rlsi",
            dims = 2:30
        )
        integrated <- IntegrateEmbeddings(
            anchorset = anchors,
            reductions = obj_merge[["lsi"]],
            new.reduction.name = "integrated_lsi",
            dims.to.integrate = 1:30)
        integrated<- RunUMAP(integrated, reduction = integrated_lsi, dims = 2:30,reduction.name=atacintegratedumap)
        integrated<- RunTSNE(integrated, reduction = integrated_lsi, dims = 2:30,reduction.name=atacintegratedtsne)
        return(integrated)
    }
    # 执行ATAC整合
    suppressWarnings({
        suppressMessages({
            integrated_atac <- integrate_atac_cca(seurat_list)
        })
    })  
    # 将整合后的降维结果添加到主对象
    obj_integrated[[integrated_lsi]] <- integrated_atac[[integrated_lsi]]
    obj_integrated[[atacintegratedumap]] <- integrated_atac[[atacintegratedumap]]
    obj_integrated[[atacintegratedtsne]] <- integrated_atac[[atacintegratedtsne]]
    DefaultAssay(obj_integrated) <- "integrated"
} else {
    atacintegratedumap="atacharmonyumap"
    atacintegratedtsne="atacharmonytsne"
    integrated_lsi="harmonylsi"
    atac_integrate_harmony <- function(obj) {
        library(harmony)
        DefaultAssay(obj) <- "combinedpeaks"
        obj<- RunHarmony(
            object = obj,
            group.by.vars = 'Sample',
            reduction = 'lsi',
            assay.use = 'combinepeaks',
            reduction.save = integrated_lsi,
            project.dim = FALSE
        )
        obj <- RunUMAP(obj,reduction = integrated_lsi,reduction.name=atacintegratedumap, dims = 2:50)
        #obj <- RunTSNE(obj, reduction = integrated_lsi,reduction.name=atacintegratedtsne,dims = 2:50,check_duplicates = FALSE)
        return(obj)
    }
    suppressWarnings({
        suppressMessages({
            obj_integrated=atac_integrate_harmony(obj_integrated)
        })
    })
}
#可视化ATAC去批次效果
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(obj_integrated, reduction = atacintegratedumap, group.by = "Sample")
options(repr.plot.width=16, repr.plot.height=8)
DimPlot(obj_integrated, reduction = atacintegratedumap, split.by = "Sample")

💡 Note
说明：此图展示经批次矫正后 ATAC 数据的分布情况。评估批次矫正效果的方法同上所述。

8. 非线性降维与聚类分析

在完成 RNA 和 ATAC 数据的整合和批次矫正后，我们可以通过非线性降维方法进一步揭示细胞在高维空间中的复杂关系。常用的非线性降维技术包括 t-SNE 和 UMAP，它们能够将高维数据映射到二维或三维空间，便于可视化观察细胞亚群结构。

三种聚类策略

RNA 聚类：基于基因表达相似性，使用 RNA-seq 数据的 PCA 降维结果构建最近邻图。
ATAC 聚类：基于染色质可及性相似性，使用 ATAC-seq 数据的 LSI 降维结果构建最近邻图。
WNN 聚类：整合两种模态信息（推荐），基于加权最近邻算法，综合考虑 RNA 和 ATAC 两种模态对细胞相似性的贡献。

注意：

在做降维聚类分析前，需分别确定 RNA 和 ATAC 数据的最佳降维维度，以避免噪音干扰或信息丢失。

options(repr.plot.width = 13, repr.plot.height = 7)
DepthCor(obj_integrated, reduction = integrated_lsi,n = 50)
#DepthCor(obj_integrated，reduction = integrated_lsi, n = 50) #如果前面用多个CCA整合数据，reduction则选择"integrated_lsi"

output

Warning message in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
"font width unknown for character 0x9"
Warning message in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
"font width unknown for character 0x9"

💡 Note
LSI 维度选择说明：在 ATAC 数据分析中，第一维 LSI 成分有时会更强地反映测序深度等技术因素，而不是生物学差异。因此，建议使用 DepthCor() 评估各个 LSI 维度与测序深度之间的相关性。如果发现第一维与测序深度存在明显相关性，则在后续邻居搜索、UMAP 降维和聚类分析中通常不使用该维度，有时候第二个 LSI 也可能是测序深度等影响的，实际选择哪些维数，根据改图来确定。

if(integration_method == "CCA"){
    tmp_reduction="pca"
}else{
    tmp_reduction="harmony"
}
ElbowPlot(obj_integrated, ndims =30, reduction = tmp_reduction)

💡 Note
主成分维度选择说明：在 RNA 数据分析中，可通过 ElbowPlot() 查看各主成分对总体变异的解释情况，从而辅助判断后续分析应保留的维度数。通常在曲线出现明显“拐点”后，新增主成分对变异的贡献会逐渐减弱，因此可将该位置附近作为维度截断的参考范围。需要注意的是，数据真实维度往往并不存在绝对固定的阈值，建议结合拐点图、已知生物学信息以及下游结果综合判断。实际分析中通常可以在相邻范围内尝试不同维度数，例如 10、15 或更高，并优先选择略偏大的维度范围，以避免遗漏潜在的生物学信号。

8.1 WNN 联合降维聚类

# 寻找多模态邻居
obj_integrated <- FindMultiModalNeighbors(
          object = obj_integrated,
          reduction.list = list(tmp_reduction, integrated_lsi),
          dims.list = list(1:30, 3:50),
          modality.weight.name = "RNA.weight",
          verbose = TRUE)
obj_integrated <- RunUMAP(
          object = obj_integrated,
          nn.name = "weighted.nn",
          assay = "RNA",  #assay = "RNA" 主要是一个 关联/上下文设置 ，不是说 UMAP 只基于 RNA。表示这个新生成的降维结果挂在RNA这个 assay 上
          verbose = TRUE,
          reduction.name = "wnnintergratedumap"
      )
obj_integrated <- FindClusters(obj_integrated, 
                                     graph.name = "wknn", 
                                     algorithm = 3, resolution = 0.5, 
                                     verbose = TRUE)
options(repr.plot.width = 9, repr.plot.height = 7)

DimPlot(obj_integrated, reduction = "wnnintergratedumap", group.by = "wknn_res.0.5",label=T, cols = my36colors) + ggtitle("WNN")

8.2 基于 RNA 数据进行降维聚类

obj_integrated <- FindNeighbors(object = obj_integrated, 
                                      reduction = 'pca',graph.name = "rnaneigobr", 
                                      dims = 1:30)
obj_integrated <- FindClusters(object = obj_integrated, 
                                     verbose = FALSE,graph.name = "rnaneigobr", 
                                     algorithm = 3, resolution = 0.5)

DimPlot(obj_integrated, reduction = rnaintegratedumap, group.by = "rnaneigobr_res.0.5",label=T, cols = my36colors) + ggtitle("RNA")

8.3 基于 ATAC 数据进行降维聚类

obj_integrated <- FindNeighbors(object = obj_integrated, 
                                reduction = integrated_lsi, graph.name = "atacneigobr", 
                                dims = 2:30)
obj_integrated <- FindClusters(object = obj_integrated, 
                               verbose = FALSE, graph.name = "atacneigobr",
                               resolution = 0.5, algorithm = 3)

DimPlot(obj_integrated, reduction = atacintegratedumap, group.by = "atacneigobr_res.0.5",label=T, cols = my36colors) + ggtitle("ATAC")

9. 细胞类型注释

9.1 在 scRNA-seq 数据层面查看 marker 基因的表达情况

本示例数据为 PBMC，因此整理了常见免疫细胞类型及其对应的 marker 基因集。实际分析中，应根据样本来源和组织类型，自定义相应的 marker 基因。后续可通过气泡图查看不同 cluster 中 marker 基因的表达情况，从而辅助判断各个 cluster 对应的细胞类型。

pbmc_marker_integrated <- list(
  "Plasma Cells" = c("IGHG1", "IGKC"),
  "Monocytes" = c("VCAN", "FCN1", "TREM1"),
  "B cells" = c("MS4A1", "BACH2", "PAX5"),
  "CMP" = c("MPO"),
  "Pro B cells" = c("DNTT"),
  "Erythroblast" = c("ALAS2", "NFIA", "SOX6"),
  "T cells" = c("PLXNA4", "THEMIS", "ITK", "BCL11B"),
  "NK cells" = c("NCAM1", "KLRC3"),
  "Dividing B cells" = c("MKI67", "TOP2A"),
  "pDC" = c("RUNX2")
)

# 设置图形大小
options(repr.plot.width=15, repr.plot.height=7)

# 绘制DotPlot
DefaultAssay(obj_integrated)="RNA"
DotPlot(obj_integrated, 
        group.by = "wknn_res.0.5", 
        features = pbmc_marker_integrated,
        cols = c("#f8f8f8","#ff3472"),
       #dot.min = 0.05,
       dot.scale = 8)+  # 应用自定义配色
  RotatedAxis() + 
  scale_x_discrete("") + 
  scale_y_discrete("") +
  theme(
    axis.text.x = element_text(size = 12, face = "bold", 
                              angle = 45, hjust = 1, vjust = 1),
    axis.text.y = element_text(size = 12, face = "bold"),
    plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
    legend.title = element_text(size = 10, face = "bold")
  ) +
  ggtitle("Marker Genes Expression") +
  labs(color = "Expression\nLevel")  # 修改图例标题

output

Warning message:
"The \`facets\` argument of \`facet_grid()\` is deprecated as of ggplot2 2.2.0.
ℹ Please use the \`rows\` argument instead.
ℹ The deprecated feature was likely used in the Seurat package.
Please report the issue at ."

💡 Note
细胞注释说明：在 SeekArc 单细胞双组学分析中，细胞类型注释通常基于 WNN 多模态降维结果进行，并结合 marker 基因在 wknn 聚类结果中的表达模式，对各个 cluster 进行细胞类型判定。

9.2 在 scATAC-seq 数据层面查看 marker 基因组区域开放情况

除 RNA 表达外，还可以从 scATAC-seq 数据层面查看 marker 基因附近调控区域的开放情况，以进一步辅助细胞类型注释。通常可结合基因附近的 peak 信号、染色质开放程度以及不同 cluster 间的可及性差异，判断特定 marker 基因是否在对应细胞群中具有活跃的调控状态。

DefaultAssay(obj_integrated) <- "combinedpeaks"
P1 <- CoveragePlot(
  object = obj_integrated,
  region = "VCAN",
  features = "MS4A1",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

P2 <- CoveragePlot(
  object = obj_integrated,
  region = "RUNX2",
  features = "RUNX2",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

P3 <- CoveragePlot(
  object = obj_integrated,
  region = "MPO",
  features = "MPO",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

P4 <- CoveragePlot(
  object = obj_integrated,
  region = "ALAS2",
  features = "ALAS2",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

P5 <- CoveragePlot(
  object = obj_integrated,
  region = "THEMIS",
  features = "THEMIS",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

P6 <- CoveragePlot(
  object = obj_integrated,
  region = "KLRC3",
  features = "KLRC3",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

options(repr.plot.width=15, repr.plot.height=4)
patchwork::wrap_plots(P1, P2, P3, ncol = 3)
patchwork::wrap_plots(P4, P5, P6, ncol = 3)

output

Warning message:
"Removed 54 rows containing missing values or values outside the scale range
(\`geom_segment()\`)."

Warning message:
"Removed 90 rows containing missing values or values outside the scale range
(\`geom_segment()\`)."

Warning message:
"Removed 11 rows containing missing values or values outside the scale range
(\`geom_segment()\`)."

Warning message:
"Removed 3 rows containing missing values or values outside the scale range
(\`geom_segment()\`)."

9.3 细胞类型标注及其可视化

# 基于聚类结果进行细胞类型注释（需要根据实际的marker基因表达情况调整）
# 这里提供一个示例，实际使用时需要根据DotPlot结果进行调整
cat('开始细胞类型注释...', Sys.time(), '\n')

# 基于聚类结果进行细胞类型注释（需要根据实际的marker基因表达情况调整）
# 这里提供一个示例，实际使用时需要根据DotPlot结果进行调整
celltype_mapping <- c(
  "15" = "Plasma Cells",
  "1" = "Monocytes",
  "14" = "Monocytes",
  "4" = "B cells",
  "6" = "B cells",
  "13" = "B cells",
  "9" = "CMP",
  "8" = "Pro B cells",
  "3" = "Erythroblast",
  "0" = "T cells",
  "5" = "T cells",
  "10" = "T cells",
  "12" = "T cells",    
  "2" = "NK cells",
  "16" = "NK cells",
  "7" = "Dividing B cells",
  "11" = "pDC"
)
# 应用细胞类型注释
obj_integrated$celltype <- recode(
  obj_integrated$wknn_res.0.5,
  !!!celltype_mapping
)

output

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

# 细胞类型UMAP可视化
p1 <- DimPlot(
  obj_integrated,
  reduction = "wnnintergratedumap",
  group.by = "celltype",
  label = TRUE,
  label.size = 3,
  cols = my36colors
) +
  ggtitle("celltype") +
  theme(legend.position = "bottom")

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

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

cat('细胞类型注释结果分样本展开!', '\n')

细胞类型注释结果分样本展开!

10. 保存结果

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

ATAC + RNA 多组学：基于 Signac 的多样本整合分析 ​

1. 教程简介 ​

2. 输入文件准备 ​

2.1 输入文件要求 ​

2.2 目录结构示例 ​

2.3 注意事项 ​

3. 数据加载与预处理 ​

3.1 获取基因注释信息 ​

3.2 数据读取与预处理 ​

4. 数据质量控制 ​

4.1 质量指标计算 ​

4.2 质量指标可视化 ​

4.3 低质量细胞过滤 ​

5. 计算样本间的共有 peaks ​

6. 多样本合并 ​

7. 数据整合 ​

7.1 整合方法选择 ​

注意事项 ​

7.2 开始对 RNA 数据整合 ​

7.3 ATAC 数据整合 ​

8. 非线性降维与聚类分析 ​

三种聚类策略 ​

注意： ​

8.1 WNN 联合降维聚类 ​

8.2 基于 RNA 数据进行降维聚类 ​

8.3 基于 ATAC 数据进行降维聚类 ​

9. 细胞类型注释 ​

9.1 在 scRNA-seq 数据层面查看 marker 基因的表达情况 ​

9.2 在 scATAC-seq 数据层面查看 marker 基因组区域开放情况 ​

9.3 细胞类型标注及其可视化 ​

10. 保存结果 ​

ATAC + RNA 多组学：基于 Signac 的多样本整合分析

1. 教程简介

2. 输入文件准备

2.1 输入文件要求

2.2 目录结构示例

2.3 注意事项

3. 数据加载与预处理

3.1 获取基因注释信息

3.2 数据读取与预处理

4. 数据质量控制

4.1 质量指标计算

4.2 质量指标可视化

4.3 低质量细胞过滤

5. 计算样本间的共有 peaks

6. 多样本合并

7. 数据整合

7.1 整合方法选择

注意事项

7.2 开始对 RNA 数据整合

7.3 ATAC 数据整合

8. 非线性降维与聚类分析

三种聚类策略

注意：

8.1 WNN 联合降维聚类

8.2 基于 RNA 数据进行降维聚类

8.3 基于 ATAC 数据进行降维聚类

9. 细胞类型注释

9.1 在 scRNA-seq 数据层面查看 marker 基因的表达情况

9.2 在 scATAC-seq 数据层面查看 marker 基因组区域开放情况

9.3 细胞类型标注及其可视化

10. 保存结果