scATAC-seq分析之Signac Motif富集分析

对scATAC-seq数据进行清洗质控【scATAC-seq数据分析之数据读入及质控】、注释以及进行差异peaks分析【scATAC-seq数据分析之标准分析流程】之后，后续的分析框架和方向与scRNA-seq分析（例如差异基因富集分析、拟时序分析和分化分析）及传统的bulk ATAC-seq分析（例如motif富集分析， motif footprinting）类似。本文基于Signac官方教程，介绍一下scATAC-seq数据的DNA 序列基序（motif）富集分析（https://stuartlab.org/signac/articles/motif_vignette）。

我们将探索两种互补的分析方法：

• 在一组差异可及（differentially accessible）peaks 中寻找富集的 motif；
• 在不同细胞群体之间进行差异 motif 活性分析。

Step1. 读入数据

在本次示例中，我们使用的是 成年小鼠脑部数据。有关用于生成该对象的代码和原始数据的链接，请参阅我们的 vignette 文档。首先，我们加载所需的 R 包和预先计算好的 Seurat 对象。

library(Signac)
library(Seurat)
library(JASPAR2020)
library(TFBSTools)
library(BSgenome.Mmusculus.UCSC.mm10)
library(patchwork)

mouse_brain <- readRDS("adult_mouse_brain.rds")
mouse_brain

p1 <- DimPlot(mouse_brain, label = TRUE, pt.size = 0.1) + NoLegend()
p1

image-20250516225703196

image-20250516225703196

Step2. 向 Seurat 对象添加基序（motif）信息

为了进行 motif 分析，我们需要将 DNA 序列的 motif 信息添加到 Seurat 对象中。可以通过运行 AddMotifs() 函数来实现：

# 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 = FALSE)
)

# add motif information
mouse_brain <- AddMotifs(
  object = mouse_brain,
  genome = BSgenome.Mmusculus.UCSC.mm10,
  pfm = pfm
)

为了方便在 Signac 中进行 motif 分析，我们创建了一个名为 Motif 的类，用于存储所有所需的信息，包括：

• 位置权重矩阵（PWMs） 或 位置频率矩阵（PFMs） 的列表
• 以及 motif 出现矩阵

在这里，AddMotifs() 函数会构建一个 Motif 对象，并将其添加到我们的 小鼠大脑数据集中，同时还会附加例如每个 peak 的碱基组成等信息。用户也可以使用 SetAssayData() 函数将 motif 对象添加到任何 Seurat 的 assay 中。

Step3. 查找富集的 motif

为了识别可能具有重要功能的 细胞类型特异性调控序列，我们可以在一组在不同细胞类型之间差异可及的 peaks 中，查找富集的 DNA motif。

在本示例中，我们比较 Pvalb 和 Sst 抑制性中间神经元之间的差异可及 peaks。

由于单细胞 ATAC-seq 数据较为稀疏，我们发现通常需要在 FindMarkers() 函数中 降低 min.pct 阈值（默认值为 0.1，适用于 scRNA-seq）。随后，我们使用 超几何检验（hypergeometric test），评估在一组 peaks 中观察到某个 motif 的频率是否高于随机背景（背景 peak 集经过 GC 含量匹配）。

da_peaks <- FindMarkers(
  object = mouse_brain,
  ident.1 = 'Pvalb',
  ident.2 = 'Sst',
  only.pos = TRUE,
  test.use = 'LR',
  min.pct = 0.05,
  latent.vars = 'nCount_peaks'
)

# get top differentially accessible peaks
top.da.peak <- rownames(da_peaks[da_peaks$p_val < 0.005 & da_peaks$pct.1 > 0.2, ])

我们还可以绘制 motif 的位置权重矩阵（Position Weight Matrix, PWM），以便可视化不同 motif 序列的图案（碱基偏好）：

# test enrichment
enriched.motifs <- FindMotifs(
  object = mouse_brain,
  features = top.da.peak
)

image-20250516225835365

image-20250516225835365

MotifPlot(
  object = mouse_brain,
  motifs = head(rownames(enriched.motifs))
)

image-20250508220747048

image-20250508220747048

Mef 家族的 motifs，尤其是 Mef2c，在 Pvalb 特异性的 scATAC-seq peaks 中表现出富集（参考文献：https://doi.org/10.1016/j.cell.2019.05.031；https://doi.org/10.1101/615179）。此外，Mef2c 对 Pvalb 抑制性中间神经元的发育是必需的（参考：https://www.nature.com/articles/nature25999）。

在本分析中，我们的结果与上述研究一致，在 FindMotifs() 的前几个显著结果中，观察到 Mef 家族 motifs 的强烈富集。

Step4. 计算 motif活性

我们还可以通过运行 chromVAR 来计算 每个细胞的 motif 活性评分。 这使我们能够在单细胞层面可视化 motif 活性，同时也提供了一种 识别不同细胞类型之间差异活性 motif 的替代方法。

chromVAR 可用于识别与 染色质可及性在细胞间变异相关的 motif。 有关该方法的完整描述，请参考 chromVAR 的原始论文。

mouse_brain <- RunChromVAR(
  object = mouse_brain,
  genome = BSgenome.Mmusculus.UCSC.mm10
)

DefaultAssay(mouse_brain) <- 'chromvar'

# look at the activity of Mef2c
p2 <- FeaturePlot(
  object = mouse_brain,
  features = "MA0497.1",
  min.cutoff = 'q10',
  max.cutoff = 'q90',
  pt.size = 0.1
)
p1 + p2

image-20250508220735768

image-20250508220735768

我们也可以直接在不同细胞类型之间比较 motif 活性评分是否存在差异。这种方法通常会得出与前面对差异可及 peaks 进行富集分析相似的结果。

在对 chromVAR 的 z-score 进行差异性检验时，我们可以在 FindMarkers() 函数中设置参数：

• mean.fxn = rowMeans
• fc.name = "avg_diff"

这样计算得到的 fold-change 就表示两组之间 z-score 的平均差值。

differential.activity <- FindMarkers(
  object = mouse_brain,
  ident.1 = 'Pvalb',
  ident.2 = 'Sst',
  only.pos = TRUE,
  mean.fxn = rowMeans,
  fc.name = "avg_diff"
)

MotifPlot(
  object = mouse_brain,
  motifs = head(rownames(differential.activity)),
  assay = 'peaks'
)

image-20250508220726376

image-20250508220726376