ATAC + RNA 多组学：基于 Signac 的单样本基础分析

作者: SeekGene

时长: 33 分钟

字数: 7.9k 字

更新: 2026-07-07

阅读: 0 次

1. 教程简介

本教程基于 Seurat + Signac 单细胞多组学分析工具，用于对 SeekArc 单细胞多组学数据（RNA + ATAC） 进行 单样本基础分析、联合降维聚类、细胞类型注释，并进一步开展 差异 peak 分析、LinkPeaks 分析、Motif 分析 和 Footprint 分析。

该流程面向一个细胞同时包含 转录组表达信息 和 染色质开放性信息 的单细胞数据，在单一样本内结合 RNA 与 ATAC 两种模态，完成从数据读取、质量控制到联合分析与调控机制解析的完整流程，实现以下核心分析目标：

单样本多组学联合分析：在同一样本内同时利用 RNA 表达和 ATAC 开放染色质信息，更全面地刻画细胞状态与生物学异质性。
数据质控与预处理：对 RNA 和 ATAC 两种模态分别进行质量控制、标准化和降维处理，为后续分析提供可靠的数据基础。
单模态与多模态聚类分析：分别开展 RNA 聚类、ATAC 聚类以及基于 WNN（Weighted Nearest Neighbor） 的联合聚类分析，用于比较不同模态下的细胞分群结果。
联合降维与可视化：融合 RNA 和 ATAC 两种模态信息，构建更稳定的低维表示，并通过 UMAP 展示细胞分布特征。
细胞类型注释：结合 marker 基因表达模式，对聚类结果进行生物学解释和细胞类型标注。
差异 peaks 分析：在不同细胞群或条件间鉴定差异开放染色质区域，筛选与细胞状态相关的调控元件。
差异 peaks 注释：将差异 peaks 注释最邻近的基因，辅助理解其潜在调控功能。
Linkpeaks 分析：关联染色质开放峰与目标基因表达，识别潜在调控关系，并结合基因区域可及性进行可视化验证。
Motif 富集分析：对差异开放 peaks 进行转录因子结合位点富集分析，筛选可能参与细胞状态调控的关键转录因子。
Footprint 分析：通过观察转录因子结合位点（motif）周围由于蛋白结合而形成的“足迹”信号，精准评估特定转录因子的实际结合活性，从而深入解析不同细胞类型特异的基因调控网络。

#加载必要的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(JASPAR2020)
  library(TFBSTools)
  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')

2. 输入文件准备

2.1 输入文件要求

单样本分析要求将 RNA 表达矩阵、ATAC peak 开放矩阵及对应的片段文件整理在同一个样本目录下。样本文件夹名称建议直接使用样本 ID，例如 S127。该目录下应包含以下内容：

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

S127/
├── 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
├── S127_A_fragments.tsv.gz
└── S127_A_fragments.tsv.gz.tbi

2.3 注意事项

输入文件需放置在同一样本目录下，保证路径结构正确，否则本教程可能无法正常读取。
ATAC 片段文件及其索引文件必须成对存在，否则无法正常构建 ChromatinAssay。
RNA 和 ATAC 数据应来自同一批细胞或同一多组学实验体系，以保证后续联合分析的准确性。
ATAC peak 坐标格式需与参考基因组注释的染色体命名风格保持一致。

3. 数据加载与预处理

在正式开展单样本多组学分析之前，需要先完成基因组注释信息加载，以及单个样本 RNA 和 ATAC 数据的读取与对象构建。

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 peak 开放矩阵以及片段文件，构建包含双组学信息的 Seurat 对象。

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

读取 RNA 数据：从 filtered_feature_bc_matrix 目录中读取 scRNA-seq 表达矩阵，并创建 RNA assay 的 Seurat 对象。
读取 ATAC 数据：从 filtered_peaks_bc_matrix 目录中读取 scATAC-seq peak 开放矩阵，同时指定片段文件路径。
构建 ATAC assay：利用 CreateChromatinAssay() 创建 ATAC assay，并引入前面生成的 annotation 注释信息，用于后续 TSS 富集分析、peak 注释及其他 ATAC 下游分析。
生成多组学对象：将 ATAC assay 添加到已有的 RNA Seurat 对象中，得到同时包含 RNA 和 ATAC 信息的单样本多组学对象。

完成该步骤后，即可获得一个标准的单样本 SeekArc 多组学 Seurat 对象，为后续的数据质量控制、标准化、降维、聚类和细胞注释分析做好准备。

# --- 输入参数配置 ---
file_path = "/path/to/sampleDir/"
sample_name = "pbmc_1"

# load the RNA and ATAC data
atac_path <- file.path(file_path, sample_name, 'filtered_peaks_bc_matrix')
rna_path <- file.path(file_path, sample_name, 'filtered_feature_bc_matrix')
frag_path <- file.path(file_path, sample_name, paste0(sample_name, '_A_fragments.tsv.gz'))
# 读取RNA数据
rna_counts <- Read10X(data.dir = rna_path)
# create a Seurat object containing the RNA adata
seu <- CreateSeuratObject(
  counts = rna_counts,
  assay = "RNA"
)
# 读取ATAC数据
atac_counts <- Read10X(data.dir = atac_path)
atac_counts <- atac_counts[Matrix::rowSums(atac_counts > 0) >= 3, ]  
# create ATAC assay and add it to the object
seu[["ATAC"]] <- CreateChromatinAssay(
  counts = atac_counts,
  sep = c(":", "-"),
  fragments = frag_path,
  annotation = annotation
)

4. 数据质量控制

4.1 质量指标计算

开始分析前，需要对数据进行质量控制，以去除低质量细胞，保证后续分析结果的可靠性。本步骤会计算 RNA 和 ATAC 两个模态的常用质控指标。

RNA 质控指标：

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

ATAC 质控指标：

TSS.enrichment：TSS 富集分数，用于评估开放染色质信号在转录起始位点附近的富集程度。
nucleosome_signal：核小体信号，用于反映片段分布特征，通常越低越好。
nCount_ATAC：每个细胞的 ATAC 总计数，用于衡量染色质开放信号强度。
nFeature_RNA：每个细胞检测到的 RNA 特征数，可辅助判断细胞质量。

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

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

4.2 质量指标可视化

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

建议重点关注以下内容：

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

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

4.3 低质量细胞过滤

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

# 质控过滤
 
cells_before <- ncol(seu)
  
seu <- subset(seu,
    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(seu)
cat('样本', sample_name, '过滤完成: 过滤前', cells_before, '个细胞，过滤后', cells_after, '个细胞\n')

output

样本 pbmc_1 过滤完成: 过滤前 8433 个细胞，过滤后 7515 个细胞

5. 数据标准化处理

5.1 RNA 数据标准化和线性降维

在完成质控和低质量细胞过滤后，需要对 RNA 数据进行标准化和降维处理，为后续聚类分析提供基础。包括：

数据标准化：使用 Seurat 默认的 LogNormalize 方法，对每个细胞的表达量进行归一化处理，并进行对数转换，以减少测序深度差异带来的影响。
高变基因筛选：筛选在细胞间具有较高变异性的基因，用于后续降维和聚类分析。
数据缩放：对基因表达矩阵进行线性变换，使不同基因具有可比性，为 PCA 分析做好准备。
PCA 降维：基于标准化后的 RNA 数据进行主成分分析（PCA），提取主要变异信息，用于后续邻居搜索、聚类和可视化分析。

suppressWarnings({
  suppressMessages({
      DefaultAssay(seu) <- "RNA"
      seu <- NormalizeData(seu, assay = "RNA")
      seu <- FindVariableFeatures(seu, assay = "RNA", selection.method = "vst", nfeatures = 2000)
      seu <- ScaleData(seu, assay = "RNA")
      seu <- RunPCA(seu, assay = "RNA", npcs = 50)
      })
    })

5.2 ATAC 数据标准化和线性降维

在完成 RNA 数据预处理后，还需要对 ATAC 数据进行标准化和降维处理，以提取染色质开放性特征，为后续聚类和联合分析提供基础。包括：

数据标准化：使用 Signac 提供的 TF-IDF 方法对 ATAC 数据进行标准化，以校正不同细胞测序深度差异，并提升稀有 peak 的识别能力。
特征筛选：通过 FindTopFeatures() 选择后续分析使用的 peak 特征，通常保留信号较强或在一定数量细胞中出现的 peaks。
SVD 降维：基于标准化后的 ATAC 矩阵进行奇异值分解（SVD），提取主要变异信息，用于后续邻居搜索、聚类和联合分析。

suppressWarnings({
  suppressMessages({
      seu <- RunTFIDF(seu, assay = "ATAC")
      seu <- FindTopFeatures(seu, assay = "ATAC", min.cutoff = 'q0')
      seu <- RunSVD(seu, assay = "ATAC")
  })
})

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

完成 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(seu, reduction = "lsi",n = 50)

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

ElbowPlot(seu, ndims =30, reduction = "pca")

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

6.1 WNN 联合降维聚类

确定好 PCA 和 LSI 的有效维数后，可开始进行降维聚类分析。

suppressWarnings({
  suppressMessages({
      seu <- FindMultiModalNeighbors(seu, reduction.list = list("pca", "lsi"), dims.list = list(1:30, 3:30))  #这里即前面确定的pca和lsi的实际可用于降维聚类的维数
      # 基于WNN进行聚类
      seu <- FindClusters(seu, graph.name = "wknn", resolution = 0.5)
      # WNN UMAP
      seu <- RunUMAP(seu, nn.name = "weighted.nn", reduction.name = "wnn.umap")
        })
})

output

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 7515
Number of edges: 108882

Running Louvain algorithm...n Maximum modularity in 10 random starts: 0.9158
Number of communities: 13
Elapsed time: 0 seconds

options(repr.plot.width = 9, repr.plot.height = 7)
DimPlot(seu, reduction = "wnn.umap", group.by = "wknn_res.0.5",label=T, cols = my36colors) + ggtitle("WNN")

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

suppressWarnings({
  suppressMessages({
      seu <- RunUMAP(seu, reduction = "pca", dims = 1:30, assay = "RNA",reduction.name="rnaumap")
      seu <- FindNeighbors(seu, reduction = "pca", dims = 1:30, assay = "RNA",graph.name = "rnaneigobr")
      seu <- FindClusters(seu, resolution = 0.5, algorithm = 1,graph.name = "rnaneigobr")
        })
})

output

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 7515
Number of edges: 71144

Running Louvain algorithm...n Maximum modularity in 10 random starts: 0.9152
Number of communities: 23
Elapsed time: 0 seconds

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

6.3 基于 ATAC 数据进行降维与聚类

suppressWarnings({
  suppressMessages({
      seu <- RunUMAP(seu, reduction = "lsi", dims = 1:30, assay = "ATAC",reduction.name="atacumap")
      seu <- FindNeighbors(seu, reduction = "lsi", dims = 1:30, assay = "ATAC",graph.name = "atacneigobr")
      seu <- FindClusters(seu, resolution = 0.5, algorithm = 1,graph.name = "atacneigobr")
        })
})

output

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 7515
Number of edges: 70618

Running Louvain algorithm...n Maximum modularity in 10 random starts: 0.9129
Number of communities: 15
Elapsed time: 0 seconds

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

💡 Note
注意：不同方法可能产生不同的聚类结果，WNN 聚类通常能发现更细致的细胞亚群，建议优先使用 WNN 结果进行细胞注释及其下游分析。

7. 细胞类型注释

7.1 差异表达基因筛选

识别各细胞簇中特异性高表达的基因，基于表达比例和差异倍数筛选代表性标记基因，为细胞类型鉴定提供初步依据。可通过气泡图查看不同 cluster 中 top 差异基因的表达情况，从而辅助判断各个 cluster 对应的细胞类型。

Idents(seu) <- seu@meta.data$wknn_res.0.5
pbmc.markers <- FindAllMarkers(seu, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.58)

head(pbmc.markers)

A data.frame: 6 × 7
	p_val	avg_log2FC	pct.1	pct.2	p_val_adj	cluster	gene
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<fct>	<chr>
ACSM3	0	2.909427	0.861	0.272	0	0	ACSM3
RUBCNL	0	2.567519	0.937	0.379	0	0	RUBCNL
ROR1	0	2.547254	0.779	0.193	0	0	ROR1
RAPGEF5	0	2.422696	0.874	0.290	0	0	RAPGEF5
PCDH9	0	2.418539	0.883	0.416	0	0	PCDH9
DTX1	0	2.367083	0.879	0.254	0	0	DTX1

top.pbmc.markers <- pbmc.markers %>% 
group_by(cluster) %>% 
  top_n(n = 3, wt = avg_log2FC)

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

# 绘制DotPlot
DefaultAssay(seu)="RNA"
pbmc.markers <- pbmc.markers[!duplicated(pbmc.markers$gene),]
DotPlot(seu, 
        group.by = "wknn_res.0.5", 
        features = unique(top.pbmc.markers$gene),
        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("Top DEG Genes Expression") +
  labs(color = "Expression\nLevel")  # 修改图例标题

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

通过已知细胞类型的经典标记基因，在 RNA 表达层面验证各细胞簇的生物学属性，实现初步的细胞类型匹配。本示例数据为 PBMC，因此整理了常见免疫细胞类型及其对应的 marker 基因集。实际分析中，应根据样本来源和组织类型，自定义相应的 marker 基因。

pbmc_marker_classic <- 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")
)

单细胞双组学数据细胞注释一般基于 wnn 多模态降维聚类情况，依据 marker 基因在 wknn 聚类结果中表达情况进行注释。

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

# 绘制DotPlot
DefaultAssay(seu)="RNA"
DotPlot(seu, 
        group.by = "wknn_res.0.5", 
        features = pbmc_marker_classic,
        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
注意：①气泡颜色反映基因在该亚群中的平均表达水平（浅色→低表达，深色→高表达）；②气泡大小表示表达该基因的细胞比例。横坐标标签表示marker基因；顶部表示对应的marker基因属于什么细胞类型

7.3 基于染色质可及性评估 marker 基因活性

利用 ATAC-seq 数据计算相同标记基因的染色质可及性活性，从表观遗传层面验证 RNA 表达注释结果，提高细胞类型鉴定的准确性。

DefaultAssay(seu) <- "ATAC"
pbmc_marker_genes <- unlist(pbmc_marker_classic, use.names = FALSE)
gene.activities <- GeneActivity(seu,features = pbmc_marker_genes)

# add the gene activity matrix to the Seurat object as a new assay and normalize it
seu[['GeneActivity']] <- CreateAssayObject(counts = gene.activities)
seu <- NormalizeData(
  object = seu,
  assay = 'GeneActivity',
  normalization.method = 'LogNormalize',
  scale.factor = median(seu$nCount_GeneActivity)
)

DefaultAssay(seu)="GeneActivity"
existing_genes <- rownames(seu)

# 过滤pbmc_marker_classic，只保留在seu中存在的基因
pbmc_marker_classic <- lapply(pbmc_marker_classic, function(gene_vec) {
  gene_vec[gene_vec %in% existing_genes]
})

# 删除空列表元素（即该细胞类型的所有标记基因都不在seu中）
pbmc_marker_classic <- pbmc_marker_classic[sapply(pbmc_marker_classic, length) > 0]

# 设置图形大小
options(repr.plot.width=12, repr.plot.height=6)
# 绘制DotPlot
DotPlot(seu, 
        group.by = "wknn_res.0.5", 
        features = pbmc_marker_classic,
        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:
"Removed 50 rows containing missing values or values outside the scale range
(\`geom_point()\`)."

7.4 marker 基因组区域开放情况

通过 CoveragePlot 同步展示特定标记基因的染色质可及性模式与 RNA 表达水平，在基因组坐标上直观验证两种组学数据的一致性，为细胞类型特异性 marker 表达和开放提供可视化证据。

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

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

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

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

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

P6 <- CoveragePlot(
  object = seu,
  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()\`)."

8. 细胞注释

8.1 细胞类型注释

综合差异表达基因筛选、经典标记基因表达模式及染色质可及性活性的多组学分析结果，对 wknn 聚类得到的各细胞簇进行系统性的细胞类型鉴定。通过建立聚类编号与已知细胞类型（如 B 细胞、T 细胞、单核细胞、NK 细胞等）之间的映射关系，完成对全部细胞的生物学分类。

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

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

# 应用细胞类型注释
seu$celltype <- recode(
  seu$wknn_res.0.5,
  !!!celltype_mapping
)

output

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

8.2 注释结果可视化

将细胞类型注释结果在 WNN-UMAP 降维空间中进行可视化，直观展示不同细胞类型在二维空间的分布特征，验证了注释的合理性和细胞群体的分离效果。

cat('细胞类型注释可视化...', '\n')

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

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

output

细胞类型注释可视化...

9. 进阶分析

9.1 差异 peak 分析

为识别细胞类型特异的调控序列，在细胞类型间差异可及的染色质区域中进行富集分析。对于稀疏数据（scATAC-seq），通常需要在 FindMarkers()中将 min.pct阈值适当降低。

Idents(seu) <- "celltype"
DefaultAssay(seu) <- "ATAC"
da_peaks <- FindMarkers(
  object = seu,
  ident.1 = 'B cells',
  ident.2 = 'T cells',
  only.pos = TRUE,  # 只保留正差异峰
  test.use = 'LR',  # 使用逻辑回归检验
  min.pct = 0.1,  # 提高最小细胞比例阈值
  logfc.threshold = 0.58  # 添加logFC
)

# 获取显著差异可及性峰
top.da.peak <- rownames(da_peaks[da_peaks$p_val < 0.005 & da_peaks$pct.1 > 0.2, ])

head(top.da.peak)

'chr2-231671510-231673282'
'chr10-1575356-1577892'
'chr15-90190821-90192528'
'chr12-113074961-113076104'
'chr7-74100205-74102036'
'chr2-136206351-136207902'

DefaultAssay(seu) <- "ATAC"
P1 = CoveragePlot(
    object = seu,
    region = top.da.peak[1],
    extend.upstream = 8000,
    extend.downstream = 5000
)
P2 = CoveragePlot(
    object = seu,
    region = top.da.peak[2],
    extend.upstream = 8000,
    extend.downstream = 5000
)

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

output

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

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

9.2 差异 peak 的 motif 富集分析

对筛选出的差异可及性峰集合，检测其中是否存在过表达的 DNA 基序。通过超几何检验，评估基序在差异峰中相对于 GC 含量匹配的背景峰是否显著富集。

# Get a list of motif position frequency matrices from the JASPAR database
pfm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = TRUE)
)

# add motif information
seu <- AddMotifs(
  object = seu,
  genome = BSgenome.Hsapiens.UCSC.hg38,
  pfm = pfm
)

head(rownames(da_peaks))

'chr2-231671510-231673282'
'chr10-1575356-1577892'
'chr15-90190821-90192528'
'chr12-113074961-113076104'
'chr7-74100205-74102036'
'chr2-136206351-136207902'

# 只使用有完整元数据的峰
peak_names <- rownames(da_peaks)
peak_metadata <- seu@assays$ATAC@meta.features[peak_names, ]

# 检查哪些峰有NA值
has_na <- apply(peak_metadata, 1, function(x) any(is.na(x)))
# 移除有NA值的峰
valid_peaks <- peak_names[!has_na]
# 进行motif富集分析
if(length(valid_peaks) > 0) {
  enriched.motifs <- FindMotifs(
    object = seu,
    features = valid_peaks
  )
} else {
  message("没有可用的峰进行motif分析")
}

options(repr.plot.width = 15, repr.plot.height = 5)
MotifPlot(
  object = seu,
  motifs = head(rownames(enriched.motifs))
)

output

Warning message:
"The \`\` argument of \`guides()\` cannot be \`FALSE\`. Use "

none" instead as
of ggplot2 3.3.4.
ℹ The deprecated feature was likely used in the ggseqlogo package.
Please report the issue at ."

9.3 差异 peak 注释

通过将差异可及性峰注释到其邻近的基因，建立染色质开放性与基因表达之间的关联。该分析有助于解释差异可及性峰在转录调控中的作用，识别受其调控的潜在靶基因，并为后续的生物学功能解析提供依据。

closest_genes <- ClosestFeature(seu, rownames(da_peaks))

head(closest_genes)

A data.frame: 6 × 8
	tx_id	gene_name	gene_id	gene_biotype	type	closest_region	query_region	distance
	<chr>	<chr>	<chr>	<chr>	<fct>	<chr>	<chr>	<int>
ENSE00001940067	ENST00000466801	PTMA	ENSG00000187514	protein_coding	exon	chr2-231706895-231707228	chr2-231671510-231673282	33612
ENST00000381312	ENST00000381312	ADARB2	ENSG00000185736	protein_coding	gap	chr10-1379161-1737050	chr10-1575356-1577892	0
ENST00000559792	ENST00000559792	SEMA4B	ENSG00000185033	protein_coding	utr	chr15-90191928-90192026	chr15-90190821-90192528	0
ENST00000257600	ENST00000257600	DTX1	ENSG00000135144	protein_coding	gap	chr12-113058452-113077423	chr12-113074961-113076104	0
ENST00000538333	ENST00000538333	LIMK1	ENSG00000106683	protein_coding	gap	chr7-74099239-74105874	chr7-74100205-74102036	0
ENSE00001587001	ENST00000241393	CXCR4	ENSG00000121966	protein_coding	exon	chr2-136118046-136118165	chr2-136206351-136207902	88185

9.4 峰-基因关联分析

通过计算基因表达量与邻近染色质开放区域（Peaks）可及性的相关性，识别潜在的顺式调控元件。分析过程中已针对 GC 含量、整体可及性背景及峰大小等潜在混杂因素进行了校正。

DefaultAssay(seu) <- "ATAC"

# first compute the GC content for each peak
seu <- RegionStats(seu, genome = BSgenome.Hsapiens.UCSC.hg38)

# link peaks to genes
seu <- LinkPeaks(
  object = seu,
  peak.assay = "ATAC",
  expression.assay = "RNA",
  genes.use = c("MS4A1", "BACH2", "PAX5")
)

p3 <- CoveragePlot(
  object = seu,
  region = "MS4A1",
  features = "MS4A1",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

p4 <- CoveragePlot(
  object = seu,
  region = "BACH2",
  features = "BACH2",
  expression.assay = "RNA",
  extend.upstream = 2000,
  extend.downstream = 2000
)

patchwork::wrap_plots(p3, p4, ncol = 2)

output

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

9.5 Footprint足迹分析

差异peak的motif富集分析可提示潜在的调控转录因子，而足迹分析则更进一步：通过检测转录因子结合时产生的DNA酶切保护信号，直接鉴定那些在细胞中实际发挥调控功能的活跃转录因子结合位点，从"可能结合"升级到"正在作用"的证据层面。

# gather the footprinting information for sets of motifs
seu <- Footprint(
  object = seu,
  motif.name = "CEBPA"),
  genome = BSgenome.Hsapiens.UCSC.hg38
)

# plot the footprint data for each group of cells
p2 <- PlotFootprint(seu, features = "CEBPA")

options(repr.plot.width = 12, repr.plot.height = 7)
p2 + patchwork::plot_layout(ncol = 1)

output

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

10. 保存结果

# 保存整合后的Seurat对象
saveRDS(seu, file = "processed.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. 数据标准化处理 ​

5.1 RNA 数据标准化和线性降维 ​

5.2 ATAC 数据标准化和线性降维 ​

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

三种聚类策略 ​

注意： ​

6.1 WNN 联合降维聚类 ​

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

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

7. 细胞类型注释 ​

7.1 差异表达基因筛选 ​

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

7.3 基于染色质可及性评估 marker 基因活性 ​

7.4 marker 基因组区域开放情况 ​

8. 细胞注释 ​

8.1 细胞类型注释 ​

8.2 注释结果可视化 ​

9. 进阶分析 ​

9.1 差异 peak 分析 ​

9.2 差异 peak 的 motif 富集分析 ​

9.3 差异 peak 注释 ​

9.4 峰-基因关联分析 ​

9.5 Footprint足迹分析 ​

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. 数据标准化处理

5.1 RNA 数据标准化和线性降维

5.2 ATAC 数据标准化和线性降维

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

三种聚类策略

注意：

6.1 WNN 联合降维聚类

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

6.3 基于 ATAC 数据进行降维与聚类

7. 细胞类型注释

7.1 差异表达基因筛选

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

7.3 基于染色质可及性评估 marker 基因活性

7.4 marker 基因组区域开放情况

8. 细胞注释

8.1 细胞类型注释

8.2 注释结果可视化

9. 进阶分析

9.1 差异 peak 分析

9.2 差异 peak 的 motif 富集分析

9.3 差异 peak 注释

9.4 峰-基因关联分析

9.5 Footprint足迹分析

10. 保存结果