【原】基于Signac的motif分析流程学习

健明 2025-11-16 发布于广东

展开全文

在基因组学和生物信息学中，motif是指DNA、RNA或蛋白质序列中具有特定生物学功能的、保守的短序列模式。在scATAC-seq或染色质可及性分析中，motif通常指转录因子能够特异性结合的DNA序列模式。

由于不同细胞中基因的结合位点序列可能有差异，因此不会用一个固定序列来描述，而是用一个概率模型来表示每个位置上出现的A/T/C/G碱基的可能性。目前表示motif会使用两个参数，分别是PFM和PWM，PFM是指Position Frequency Matrix（位置频率矩阵），其用于记录在已知结合位点中，每个位置上A/T/C/G出现的实际次数，PWM是指Position Weight Matrix（位置权重矩阵），其由PFM转换而来，用对数似然比表示每个碱基在该位置的重要性，用于打分新序列是否可能是结合位点。

分析流程：

1.导入

rm(list=ls())
library(Signac)
library(Seurat)
library(JASPAR2020)
library(TFBSTools)
library(patchwork)
library(BSgenome.Hsapiens.UCSC.hg38)
library(qs2)

2.数据预处理

# 这个pbmc.qs2是基于Signac的scATAC-seq流程学习中得到的结果
pbmc <- qs2::qs_read("pbmc.qs2")
pbmc

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

3.向Seurat对象添加motif信息

#为了便于在 Signac 中进行 motif 分析，我们创建了 Motif 类，用于存储所有必需的信息，包括一组位置权重矩阵（PWMs）或位置频率矩阵（PFMs），以及一个 motif出现矩阵。
#在这里，AddMotifs() 函数会构建一个 Motif 对象，并将其添加到数据集中，同时还会附加其他信息，例如每个peak的碱基组成。

# 从JASPAR数据库中获取motif位置频率矩阵列表。
pfm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)

# 添加motif信息
pbmc <- AddMotifs(
  object = pbmc,
  genome = BSgenome.Hsapiens.UCSC.hg38,
  pfm = pfm
)

4.寻找过表达的motif

为了识别可能具有重要功能的细胞类型特异性调控序列，我们可以在不同细胞类型之间差异可及的染色质开放区域（peaks）中，搜索显著富集的DNAmotif。

在此，比较了pDC和NK dim两类细胞，寻找它们之间差异可及的peaks。由于scATAC-seq数据较为稀疏，发现通常需要将 FindMarkers() 函数中的 min.pct 阈值调低（默认值为 0.1，该参数最初是为 scRNA-seq 数据设计的），以提高检测灵敏度。

随后，采用超几何检验（hypergeometric test），评估某个motif在目标peak集中出现的频率是否显著高于随机预期。检验时，使用一组GC 含量匹配的背景peaks作为对照，以排除序列组成偏差带来的影响。

da_peaks <- FindMarkers(
  object = pbmc,
  ident.1 = 'pDC',
  ident.2 = 'NK dim',
  only.pos = TRUE,
  test.use = 'LR',
  min.pct = 0.05,
  latent.vars = 'nCount_peaks'
)

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

5.检验富集性

# test enrichment
enriched.motifs <- FindMotifs(
  object = pbmc,
  features = top.da.peak
)
head(enriched.motifs)
#             motif observed background percent.observed percent.background fold.enrichment        pvalue     motif.name      p.adjust
# MA0081.2 MA0081.2     4866      16020         62.58521            40.0500        1.562677  0.000000e+00           SPIB  0.000000e+00
# MA0080.5 MA0080.5     5070      18559         65.20900            46.3975        1.405442 2.032416e-302           SPI1 7.580912e-300
# MA0521.1 MA0521.1     2676       7852         34.41801            19.6300        1.753337 4.254708e-264          Tcf12 1.058004e-261
# MA0500.2 MA0500.2     1745       4353         22.44373            10.8825        2.062369 4.094314e-248           MYOG 7.635896e-246
# MA0816.1 MA0816.1     2260       6384         29.06752            15.9600        1.821273 4.588650e-240          Ascl2 6.846265e-238
# MA1635.1 MA1635.1     2145       5942         27.58842            14.8550        1.857181 1.295870e-238 BHLHE22(var.2) 1.611198e-236

6.绘制这些motif的位置权重矩阵（PWM），以便直观地查看不同motif的序列特征

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

7.计算motif活性

这里的MA0497.1是JASPAR数据库中Mef2c转录因子的DNA结合基序的唯一标识符，可以在JASPAR网站中查到相关信息：https://jaspar./matrix/MA0497.1/

pbmc <- RunChromVAR(
  object = pbmc,
  genome = BSgenome.Hsapiens.UCSC.hg38
)

DefaultAssay(pbmc) <- 'chromvar'

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

7.1 另一种方式

还可以直接对不同细胞类型之间的 chromVAR活性分数进行差异检验。这种方法通常会得到与前述方法——即在细胞类型间差异可及的peaks上进行motif富集分析——相似的结果。

在对chromVAR的z-score进行差异分析时，可以在FindMarkers函数中设置参数mean.fxn = rowMeans和fc.name = "avg_diff"，这样计算出的“倍数变化”实际上表示的是两组细胞之间z-score 的平均差值，而非传统的比值形式。

differential.activity <- FindMarkers(
  object = pbmc,
  ident.1 = 'pDC',
  ident.2 = 'NK dim',
  only.pos = TRUE,
  mean.fxn = rowMeans,
  fc.name = "avg_diff"
)

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