在R语言中的 ATACseq 数据分析全流程实战（七）：Motif分析-阿南达文事网

在R语言中的 ATACseq 数据分析全流程实战（七）：Motif分析

本帖子学习资源：/

续上面第六期~

在ATAC实验中鉴定的Peak区域，为染色质开放的区域，含有转录因子结合位点（transcription factor binding site, TFBS），此位点对于转录因子介导的基因表达的调控起到重要的作用。因此，我们可以通过鉴定富集在Peak区域的代表性的Motif（over-represented）序列推测结合在此区域的潜在的转录因子。

motifs简单定义：

Motif 就是一段特定模式的DNA序列，可以理解为开放的染色质区域有着序列偏好性，而偏向的序列就是motif。识别Motif出发点是找出序列中或不同序列间的相似性片段，从而归结出序列片段中蕴涵的特征模式，进而推断出该特征模式与已知的结构和功能之间的内在联系。

ATACseq vs ChIPseq

与 ChIPseq 不同，ChIPseq 中我们知道免疫沉淀（IP）的目标是什么，而在 ATACseq 中，我们只知道某个区域是开放/可及的，但不知道可能存在哪些转录因子。

ChIPseq：适用于研究特定转录因子或组蛋白修饰的结合位点，适合验证性研究。
ATACseq：适用于发现性研究，例如在未知转录因子结合位点的情况下，寻找基因组中所有可能的开放区域。

接下来将学习如何使用 motifmatchr 来研究 ATACseq 中的motifs。

数据介绍

这里使用的数据为前面第一期：在R语言中的 ATACseq 数据分析全流程实战（一）中介绍的数据二，再重新温习一遍。

数据来自项目 ENCODE consortium 的一部分，为小鼠的组织，对应的文献为：ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature 2012 Sep 6;489(7414):57-74. PMID: 22955616

三种组织，每种组织两个重复：

Liver day 12：/
Kidney day 15：/
Hindbrain day 12：/

下载地址：

代码语言：javascript代码运行次数：0运行复制

## Liver day 12：/
# 样本1
r1：/@@download/ENCFF409BPW.fastq.gz
r2：/@@download/ENCFF016UWL.fastq.gz
# 样本2
r1：/@@download/ENCFF772GVP.fastq.gz
r2：/@@download/ENCFF987EVS.fastq.gz

## Kidney day 15：/
# 样本1
r1：/@@download/ENCFF340WOK.fastq.gz
r2：/@@download/ENCFF117TPW.fastq.gz
# 样本2
r1：/@@download/ENCFF395EDU.fastq.gz
r2：/@@download/ENCFF100YMT.fastq.gz

## Hindbrain day 12：/
# 样本1
r1：/@@download/ENCFF954MPR.fastq.gz
r2：/@@download/ENCFF824PVC.fastq.gz
# 样本2
r1：/@@download/ENCFF199WNL.fastq.gz
r2：/@@download/ENCFF968KUR.fastq.gz

数据下载完后，走之前的流程得到peaks。这里可以使用教程中准备好的多个样本的peaks结果count矩阵，链接为：myCounts.RData：/

已知Motif 数据库了解

已知的motif资源，Bioconductor提供了两个motif 数据库包：

MotifDb：.html
JASPAR databases, JASPAR2020是最新版：.html

1.MotifDb

简单看看 MotifDb 这个数据库包：相比作者笔记中2013年，这个包更新到了2022年。MotifDB 是一个专门用于收集和存储motif信息的 R 包，它整合了多种来源的数据，主要来源数据库：

FlyFactorSurvey:	614
hPDI:	437
JASPAR_CORE:	459
jolma2013:	843
ScerTF:	196
stamlab:	683
UniPROBE:	380
cisbp 1.02	874

主要支持的物种有：

Hsapiens:	2328
Dmelanogaster:	1008
Scerevisiae:	701
Mmusculus:	660
Athaliana:	160
Celegans:	44
other:	177

存储为 DB 对象：

代码语言：javascript代码运行次数：0运行复制

library(MotifDb)
MotifDb

# MotifDb object of length 12657
# | Created from downloaded public sources, last update: 2022-Mar-04
# | 12657 position frequency matrices from 22 sources:
#   |         cisbp_1.02:  874
# |    FlyFactorSurvey:  614
# |        HOCOMOCOv10: 1066
# | HOCOMOCOv11-core-A:  181
# | HOCOMOCOv11-core-B:   84
# | HOCOMOCOv11-core-C:  135
# | HOCOMOCOv11-secondary-A:   46
# | HOCOMOCOv11-secondary-B:   19
# | HOCOMOCOv11-secondary-C:   13
# | HOCOMOCOv11-secondary-D:  290
# |              HOMER:  332
# |               hPDI:  437
# |        JASPAR_2014:  592
# |        JASPAR_CORE:  459
# |         jaspar2016: 1209
# |         jaspar2018: 1564
# |         jaspar2022: 1956
# |          jolma2013:  843
# |             ScerTF:  196
# |            stamlab:  683
# |       SwissRegulon:  684
# |           UniPROBE:  380
# | 62 organism/s
# |           Hsapiens: 6075
# |          Mmusculus: 1554
# |      Dmelanogaster: 1437
# |          Athaliana: 1371
# |        Scerevisiae: 1221
# |                 NA:  184
# |              other:  815
# Scerevisiae-cisbp_1.02-M0001_1.02 
# Scerevisiae-cisbp_1.02-M0002_1.02 
# Scerevisiae-cisbp_1.02-M0003_1.02 
# Csativa-cisbp_1.02-M0004_1.02 
# Athaliana-cisbp_1.02-M0005_1.02 
# ...
# Mmusculus-UniPROBE-Zfp740.UP00022 
# Mmusculus-UniPROBE-Zic1.UP00102 
# Mmusculus-UniPROBE-Zic2.UP00057 
# Mmusculus-UniPROBE-Zic3.UP00006 
# Mmusculus-UniPROBE-Zscan4.UP00026 

class(MotifDb)
## [1] "MotifList"
## attr(,"package")
## [1] "MotifDb"

# 简单探索
length(MotifDb)
MotifNames <- names(MotifDb)
MotifNames[1:10]

# [1] "Scerevisiae-cisbp_1.02-M0001_1.02" "Scerevisiae-cisbp_1.02-M0002_1.02" "Scerevisiae-cisbp_1.02-M0003_1.02"
# [4] "Csativa-cisbp_1.02-M0004_1.02"     "Athaliana-cisbp_1.02-M0005_1.02"   "Athaliana-cisbp_1.02-M0006_1.02"  
# [7] "Athaliana-cisbp_1.02-M0007_1.02"   "Athaliana-cisbp_1.02-M0008_1.02"   "Athaliana-cisbp_1.02-M0009_1.02"  
# [10] "Athaliana-cisbp_1.02-M0010_1.02"

位置概率矩阵PPM：是一个矩阵，每一行代表一个核苷酸（A、C、G、T），每一列代表一个位置。矩阵中的每个值表示在该位置上特定核苷酸出现的概率。

代码语言：javascript代码运行次数：0运行复制

# extract the position probability matrix.
MotifDb[[1]]
colSums(MotifDb[[1]])

提取metadata信息：

代码语言：javascript代码运行次数：0运行复制

# metadata information
values(MotifDb)[1:2, ]

提取子集：以CTCF蛋白为例（/）

CTCF 是一种调控染色质三维结构的 DNA 结合因子，它的结合基序是什么呢？查询一下：

代码语言：javascript代码运行次数：0运行复制

# subset，检索所有被注释为“ctcf”（不区分大小写）的motif
CTCFMotifs <- query(MotifDb, "CTCF")
CTCFMotifs

# 限定一下检索条件
CTCFMotifs <- query(MotifDb, c("CTCF", "hsapiens", "jaspar2018"))
CTCFMotifs
length(CTCFMotifs)

# 比较：
# Biostrings 包中的 `consensusString` 函数
sapply(CTCFMotifs, consensusString)
CTCFMotifs[[1]]
CTCFMotifs[[2]]
# 可视化
library(seqLogo)
seqLogo(CTCFMotifs[[1]])  # Hsapiens-jaspar2018-CTCF-MA0139.1
seqLogo(CTCFMotifs[[2]])  # Hsapiens-jaspar2018-CTCFL-MA1102.1

来自不同来源的motifs有时会一致，有时则会有所不同，可以比较：Biostrings 包中的 consensusString 函数提供了一种快速且有时足够有效的比较方法。

对应的注释信息：

代码语言：javascript代码运行次数：0运行复制

meta <- values(CTCFMotifs)
meta <- as.data.frame(meta)

2.JASPAR2020

这个数据库本身的官网为：/，支持留个物种：

JASPAR2020相比MotifDb包存储信息更多：

JASPAR 包：主要存储motifs和位置概率矩阵（PPMs）的信息；
TFBSTools 包：提供了操作和交互这些数据的功能。

getMatrixByID 和 getMatrixByName 函数分别接受 JASPAR 数据库对象和 JASPAR ID 或转录因子名称作为参数。

以转录因子 GATA2为例，结果是一个新的对象类 PFMatrix。

代码语言：javascript代码运行次数：0运行复制

########################################
library(JASPAR2020)
JASPAR2020

library(TFBSTools)
`?`(getMatrixByID)

# 使用基因名字提取
GATA2mat <- getMatrixByName(JASPAR2020, "GATA2")
class(GATA2mat)

# 使用具有唯一性的JASPAR IDs提取想要的motif，MA0036.2为motif ID
GATA2mat <- getMatrixByID(JASPAR2020, "MA0036.2")
GATA2mat

# 提取名字
ID(GATA2mat)

# 提取the Position Frequency Matrix (PFM)
myMatrix <- Matrix(GATA2mat)
myMatrixToo <- as.matrix(myMatrix)
myMatrix

还可以使用 getMatrixSet() 限定搜索提取motif

代码语言：javascript代码运行次数：0运行复制

# 使用 getMatrixSet() 提取motif
opts <- list()
opts[["collection"]] <- "CORE"
opts[["tax_group"]] <- "vertebrates"
opts
motifList <- getMatrixSet(JASPAR2020, opts)
motifList
motifList[1]
motifList[[1]]

提取的一个motif的信息：

代码语言：javascript代码运行次数：0运行复制

An object of class PFMatrix
ID: MA0004.1
Name: Arnt
Matrix Class: Basic helix-loop-helix factors (bHLH)
strand: +
Tags: 
$alias
[1] "HIF-1beta,bHLHe2"

$description
[1] "aryl hydrocarbon receptor nuclear translocator"

$family
[1] "PAS domain factors"

$medline
[1] "7592839"

$remap_tf_name
[1] "ARNT"

$symbol
[1] "ARNT"

$tax_group
[1] "vertebrates"

$tfbs_shape_id
[1] "11"

$type
[1] "SELEX"

$unibind
[1] "1"

$collection
[1] "CORE"

$species
         10090 
"Mus musculus"

$acc
[1] "P53762"

Background: 
   A    C    G    T 
0.25 0.25 0.25 0.25 
Matrix: 
  [,1] [,2] [,3] [,4] [,5] [,6]
A    4   19    0    0    0    0
C   16    0   20    0    0    0
G    0    1    0   20    0   20
T    0    0    0    0   20    0

一个motif的相关信息都在这张图里面了，包括ID，ID命名方式，name，关联的转录因子，碱基位置频率矩阵等：

motifs可视化

SeqLogo 包通过 seqLogo 函数提供一种简单直观的方式来可视化motifs中的碱基频率，SeqLogo 图通过比较同一位置上其他碱基的相对大小，展示了每个motif位置上碱基的相对频率。

可视化一个motif看看：

代码语言：javascript代码运行次数：0运行复制

library(seqLogo)
# 可视化 来自 MotifDb数据库的motif
CTCFMotifs <- query(MotifDb, "CTCF")
seqLogo::seqLogo(CTCFMotifs[[1]], ic.scale = FALSE)

生成的 SeqLogo 图展示每个位置上碱基的相对频率，频率越高的碱基在图中显示得越大：

将概率信息转换为 information content，取值为0-2，每个位置上概率相等的碱基转换后打分为0，只有一个碱基的打分为2，这样可以突出重要的碱基：

代码语言：javascript代码运行次数：0运行复制

seqLogo::seqLogo(CTCFMotifs[[1]])

对JASPAR2020的motif可视化：

代码语言：javascript代码运行次数：0运行复制

# JASPAR2020 提供的PFM矩阵需要转换一下，然后可视化
myMatrix
ppm <- myMatrix/colSums(myMatrix)
ppm
seqLogo::seqLogo(ppm)

# 或者 TFBSTools工具自己进行转换后可视化
GATA2mat
GATA2_IC <- toICM(GATA2mat)
TFBSTools::seqLogo(GATA2_IC)

还可以使用 ggseqlogo 包：

代码语言：javascript代码运行次数：0运行复制

## ggplot风格
library(ggseqlogo)
library(ggplot2)
ggseqlogo(myMatrix) + theme_minimal()

ATAC中Motifs鉴定

ATACseq 是一种用于检测染色质开放区域的技术，它只能告诉我们哪些区域是开放的，而无法直接告诉我们哪些转录因子可能结合在这些区域。MotifDb 和 JASPAR2020 这些数据库提供了大量已知的motif，这些motif可以用于在 ATACseq 峰区域中识别潜在的转录因子结合位点。

使用的方法为：motifmatchr 包。

motifmatchr 的主要方法是 matchMotifs函数。这个函数有两个必需的参数：

位置权重矩阵（PWM）或位置频率矩阵（PFM）：存储在 TFBSTools 包中的 PWMatrix、PFMatrix、PWMatrixList 或 PFMatrixList 对象中。
一组基因组范围（GenomicRanges 或 RangedSummarizedExperiment 对象）或一组序列（可以是 DNAStringSet、DNAString 或简单的字符向量）。

1.第一个参数

首先，我们可以检索一组适合在小鼠组织的 ATACseq 数据中扫描的motifs，得到 matchMotifs 函数的第一个参数。

motif提取筛选：the vertebrate, JASPAR CORE motifs，并设置 all_versions=F 保证提取的为最新的版本：

提取到了 746个motifs

代码语言：javascript代码运行次数：0运行复制

## retrieve the vertebrate, JASPAR CORE motifs
opts <- list()
opts[["tax_group"]] <- "vertebrates"
opts[["collection"]] <- "CORE"
opts[["all_versions"]] <- FALSE
opts
motifsToScan <- getMatrixSet(JASPAR2020, opts)
motifsToScan
motifsToScan[[1]]

# PFMatrixList of length 746
# names(746): MA0004.1 MA0006.1 MA0019.1 MA0029.1 MA0030.1 MA0031.1 ... MA1123.2 MA0093.3 MA0526.3 MA0748.2 MA0528.2 MA0609.2

2.第二个参数

接着从自己的peaks中提取序列得到 matchMotifs 函数的第二个参数：

导入前面鉴定到的the ENCODE mouse tissue ATAC data：summarizedExperiment对象

myCounts.RData：/

代码语言：javascript代码运行次数：0运行复制

# 导入前面教程中ATACseq的数据
load("RU_ATACseq-master/r_course/data/myCounts.RData")
myCounts

# class: RangedSummarizedExperiment 
# dim: 36320 6 
# metadata(0):
#   assays(1): counts
# rownames: NULL
# rowData names(6): Sorted_Hindbrain_day_12_1_Small_Paired_peaks.narrowPeak
# Sorted_Hindbrain_day_12_2_Small_Paired_peaks.narrowPeak ... Sorted_Liver_day_12_1_Small_Paired_peaks.narrowPeak
# Sorted_Liver_day_12_2_Small_Paired_peaks.narrowPeak
# colnames(6): Sorted_Hindbrain_day_12_1_Small_Paired.bam Sorted_Hindbrain_day_12_2_Small_Paired.bam ...
# Sorted_Liver_day_12_1_Small_Paired.bam Sorted_Liver_day_12_2_Small_Paired.bam
# colData names(0):

提取peaks的范围，将peaks重新归一化到100bp,并提取需要的序列：

代码语言：javascript代码运行次数：0运行复制

# retrieve the ranges of our peaks 
library(SummarizedExperiment)
peakRanges <- rowRanges(myCounts)
peakRanges[1, ]

# re-center our peaks to 100bp
# and extract the required sequences 
library(BSgenome.Mmusculus.UCSC.mm10)
peakRangesCentered <- resize(peakRanges, fix = "center", width = 100)
peakSeqs <- getSeq(BSgenome.Mmusculus.UCSC.mm10, peakRangesCentered)
names(peakSeqs) <- as.character(peakRangesCentered)
peakSeqs

准备好的序列如下：

3.鉴定已知motif

matchMotifs 函数可以提供motif匹配的输出，形式可以是匹配结果、分数或位置。

使用 ATACseq peaks前100个peaks匹配来自 JASPAR2020 的motif，结果输出类型为 positions：

代码语言：javascript代码运行次数：0运行复制

# 匹配 motif
library(motifmatchr)
motif_positions <- matchMotifs(motifsToScan[1:4], peakSeqs[1:100], out = "positions")
class(motif_positions)
length(motif_positions)
motif_positions$MA0029.1

鉴定到了四个motif：

挑一个出来看看：

代码语言：javascript代码运行次数：0运行复制

MA0029hits <- motif_positions$MA0029.1
names(MA0029hits) <- names(peakSeqs[1:100])
unlist(MA0029hits, use.names = TRUE)

指定输出类型为输出类型为 matches 看看：

代码语言：javascript代码运行次数：0运行复制

# map motifs to their ATACseq peaks.
motifHits <- matchMotifs(motifsToScan, peakSeqs, out = "matches")
class(motifHits)
motifHits

# retrieve a matrix of matches by motif and peak
mmMatrix <- motifMatches(motifHits)
dim(mmMatrix)
mmMatrix[1:10, 1:8]
as.matrix(mmMatrix[1:10, 1:8])

# identify the total occurrence of motifs in our peak sequences
totalMotifOccurence <- colSums(mmMatrix)
totalMotifOccurence[1:4]

# identify peaks which contain a hit for a selected motif.
peaksWithMA0912 <- peakRangesCentered[mmMatrix[, "MA0912.2"] == 1]
peaksWithMA0912

mmMatrix匹配结果：|表示匹配到了，每一行为一个序列，每一列为一个已知的motif ID

总的motif数：

提取某个motif的peaks：

完~

下次就上具体的文献应用实战啦~

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-05-01，如有侵权请联系 cloudcommunity@tencent 删除数据库数据分析函数可视化数据

在R语言中的 ATACseq 数据分析全流程实战（七）：Motif分析