ATAC-Seq|ATAC-Seq Motif 富集分析 ATAC-SeqMotif富集分析

本文是《ATAC-Seq 分析流程》的 Motif 富集分析部分。
峰注释能够提供峰可能影响的基因和功能区域，但这不够直接解释开放区域的作用机制。开放区域可能被转录因子结合从而调控了基因的表达，转录因子结合的特定序列称之为 motif, 结合的位点称为 TFBS(TF binding sites)，转录因子可以允许 TFBS 有一定的可塑性（variations/flexible），所以 motif 序列不完全是固定死的。Motif 大部分不长，6-12 bp 左右。人类大约有 1600 个转录因子，其中超过 2/3 已鉴定了 Motif.

文章图片
TFs.png Motif 一般有 2 种模式，一是回文序列，比如 CACGTG 其反向互补序列也是 CACGTG. 二是 2 段保守序列被一段非保守序列分隔，这往往是因为结合的转录因子是二聚体，分别识别 2 段保守序列。
【ATAC-Seq|ATAC-Seq Motif 富集分析】Motif 的分析是容易受假阳性困扰的，这称为无效定理（Futility Theorem）。比如说往往在基因组序列中观察到大量的潜在转录因子结合位点，其中很少是真正起作用的，大部分预测的转录因子结合位点是无效的。
Position weight matrix
PWM 也被称为 PSWM(Position-specific weight matrix) 或 PSSM(Position-specific scoring matrix) 主要用于蛋白/核酸序列 Motif 的表征和预测。
构建 PWM 首先要有频数表 PFM(Position frequency matrix) 即一个 Motif 每个位置的碱基（以核酸为例）频数（这个数据可以来自于实验）。下面是一个例子。

文章图片
一个 PFM 例子
PWM 计算如下，公式 (1)

其中 i 是 1,2,3,4 分别表示 A,C,G,U 碱基，j 是位置。是该碱基的背景频率，是该碱基在 j 位置的频率，比如第一个位置碱基 A 的频率为

一个指定 Motif 的 PWMS(PWM score) 计算如下，公式 (2)

其中 L 为 Motif 长度，上面例子为 17.
另外，PWMS 是似然比的对数。假设一个 17bp 序列 S = ACGGTACCACGTAAGTT, 可以有 2 假设，第一该序列是一个 Motif，那么在特定位置就会有特定碱基（）；第二该序列是 17 随机碱基序列，那么特定位置不会是特定碱基模式（）。那么序列 S 在 2 种假设下的似然性（Likelihood）分别计算如下。公式 (3)

因此得到与公式 (2) 等同的公式 (4)

假设有了 PFM 了，缺少每个碱基背景频率也是不行的。最简单方法是假设背景频率相等，即都是 0.25，但这可能不符合许多生物学情景。第二种方法从输入序列计算，但这也有问题，比如说假设输入序列为 UUUUUUUU 那么会得到而 . 方法三是根据具体问题用相应的数据，比如说研究某物种的剪切位点模式，那么选该物种所有的内含子和外显子计算背景频率是适合的。
公式 (1) 可能会有些位点为 0，而这是不允许的，因为无法进行操作。解决这个问题，可以用下面公式 (5)

其中，是频数，这个方法在 N 很小时表现不好。
方法二是使用伪计数，下面公式 (6)

因此等于公式 (7)

这里的值应该取很小，如 0.0001.
HOMER
软件安装方法在Homer Software and Data Download. 安装后将 homer/bin 目录添加到 PATH 环境变量。
安装第三方软件，包括 Ghostscript, WebLogo WebLogo 3 - About, blat Index of /~kent/exe
安装需要的包。先查看有哪些可以下载的包。

perl ./configureHomer.pl -list # 只展示少部分结果作为举例 Version InstalledPackage Version Description SOFTWARE +homerv4.11.1 Code/Executables, ontologies, motifs for HOMER ORGANISMS +human-o v6.3Homo sapiens (human) accession and ontology information -hydrav6.3Hydra vulgaris (hydra) accession and ontology information -anemone v6.3Nematostella vectensis (anemone) accession and ontology information -yeast-o v6.3Saccharomyces cerevisiae (yeast) accession and ontology information

包分为 4 大类 SOFTWARE, ORGANISM, PROMOTER, GENOME 其中 SOFTWARE 就是 homer 自己。有些包有 -p/-o/-g 后缀就是因为同名，增加一个后缀做区分。

-p - PROMOTER
-o - ORGANISM
-g - GENOME

用 -install 命令安装。

perl ./configureHomer.pl -install mm9

ATAC-seq 分析得到峰后就可以用 findMotifsGenome.pl 根据峰的位置进行 Motif 富集的分析。HOMER 自动分析已知 Motif 和新发现（de novo）Motif 的富集，还将新发现 Motif 和已知的进行对比。

findMotifsGenome.pl ${GenrichDir}/KO.narrowPeak hg38 ${MotifDir} -size 200

这是个封装的脚本，运行后自动进行许多步骤的分析。主要有前期处理，包含序列提取、背景提取、Normalization 等等，然后是进行新 Motif 的发现分析，然后是已知 Motif 的分析，最后输出结果。
输出主要是 2 个网页 homerResults.html, knownResults.html 分别是新发现的 motif 及已知 motif 的结果，同名目录下有相应的 Motif PFM 和 Sequence logo 图片；所有的 Motif PFM 在 homerMotifs 开头的文本文件。

$ ls homerMotifs.all.motifshomerMotifs.motifs8knownResultsmotifFindingParameters.txt homerMotifs.motifs10homerResultsknownResults.htmlseq.autonorm.tsv homerMotifs.motifs12homerResults.htmlknownResults.txt

HOMMER 的 Motif 主要分为 3 类。

ChIP-Seq Transcription Factor Motifs, 来自于 CHIP-Seq 数据的转录因子 Motif.
Promoter Motifs, 分析启动子区得到的 Motif.
General Factors X & Y (i.e. GFX & GFY), 尚不清楚与之结合的转录因子。

MEME-Chip
跟 HOMMER 一样，MEME-Chip 也是封装了 Motif 发现、Motif 富集分析、和已知 Motif 比较分析、输出结果等步骤的脚本。
MEME-Chip 用输入的序列最中间 100bp 进行 Motif 的发现分析，用输入序列的全长进行 Motif 的富集分析。
MEME-Chip 建议使用等长序列。可以选择去峰中心上下 250bp 共 500bp 进行分析。

awk -v FS="\t" -v OFS="\t" '{midpos=$2+$10; print $1,midpos-250,midpos+250; }' \ KO.narrowPeak > KO.MEME.bed

输入的是序列，所以从上面 bed 文件用 bedtools 从参考基因组提取序列。

bedtools getfasta -fo ${Genrich}/KO.MEME.fasta -fi ${GRCh38} -bed ${Genrich}/KO.MEME.bed meme-chip -meme-maxw 30 -oc ${MEMEDir}/KO -db ${MotifDb} ${GenrichDir}/KO.MEME.fasta

MEME-Chip 输出结果在 meme.html 查看，非常界面友好。
参考
Yan, Feng, et al. "From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis." Genome biology 21.1 (2020): 22.
Zia, Amin, and Alan M. Moses. "Towards a theoretical understanding of false positives in DNA motif finding." BMC bioinformatics 13.1 (2012): 1-9.
Van Loo, Peter, and Peter Marynen. "Computational methods for the detection of cis-regulatory modules." Briefings in bioinformatics 10.5 (2009): 509-524.
Wasserman, Wyeth W., and Albin Sandelin. "Applied bioinformatics for the identification of regulatory elements." Nature Reviews Genetics 5.4 (2004): 276-287.
Xia, Xuhua. "Position weight matrix, gibbs sampler, and the associated significance tests in motif characterization and prediction." Scientifica 2012 (2012).
Lambert, Samuel A., et al. "The human transcription factors." Cell 172.4 (2018): 650-665.
Das, Modan K., and Ho-Kwok Dai. "A survey of DNA motif finding algorithms." BMC bioinformatics. Vol. 8. No. S7. BioMed Central, 2007.
Position weight matrix - Dave Tang's blog
Homer Software and Data Download
MEME-ChIP - MEME Suite
Ma, Wenxiu, William S. Noble, and Timothy L. Bailey. "Motif-based analysis of large nucleotide data sets using MEME-ChIP." Nature protocols 9.6 (2014): 1428-1450.