HOMER | MEME | 转录因子的靶基因预测 | motif富集分析

转录因子motif是一些很短的模序（~10bp），在大基因组里很容易出现随机比对，而且是以position weight matrix (PWM)格式来呈现，说明它的可变性，因此研究motif有哪些binding区域是没有意义的，因为很难找到一个方法（阈值）来判断真正的比对和随机的比对。

换个思路，如果做富集分析，那就稳了，给定一个指定的区域（promoter或enhancer区域），根据统计学检验，我们很容易知道一个motif是否显著富集在这个区域（与背景区域相比），这就回答了一个很好的生物学问题：这个转录因子是否显著地结合到这片区域。

一个突出的矛盾：转录调控的稳定性和我们收集数据不确定性之间的矛盾。

为什么ChIP-seq和ATAC-seq能极大地助力motif研究：

• 真正的开放区域，得到的都是active的区域
• 过滤掉了大部分的无效或复杂区域，假阳性得到了极大的控制

如何根据这些信息来预测每个转录因子的binding区域以及靶基因？

• 不考虑远距离的调控作用，或者只考虑promoter的区域，我们就可以根据peak的注释信息找到最近的基因。
• 然后看这些promoter上分别有哪些富集的motif，然后与转录因子对应即可。
• 最后还需要基因表达来确认这些靶基因和转录因子确实是表达的！（如果这是一个抑制的转录因子，是否基因就不表达了）

转录因子的表达具有高度的组织特异性，而且已知的TF只有1000多个，基因有30000多个，所以一个TF的靶基因可能有几百个，具有高度的时空组织特异性。

实验的方法就暂且不说了，非常可靠，但成本高、耗费劳力。

最简单的预测就是基于基因表达，co-expressed就是可能的靶基因，预测软件一大把。

问题很多，首先理解假设：

1. TF的线性变化引起target gene的线性变化，他们线性相关；

2. TF的调控是sparse的，

问题：

1. 有人说这根本就不是线性的，TF的yes or no，决定target gene的表达；

2. 不是线性相关，存在shift，先后的shift是普遍存在的；

3. co-expression是无法得出target关系的；

所以，现在大家都开始结合motif enrichment了，TF的靶向作用是靠motif与基因组DNA结合来执行的。

但是我们不知道结合位点，所以大部分的富集都默认选择了10kb的flanking region，motif很短，随机比对会带来很多假阳性。

现在，大家有open chromatin的数据了，知道了候选的结合区域，我们就可以更有效的预测了，这就是HOMER的预测功能。

最终，open chromatin还是不准，因为DNA有三维结构，distal regulation是普遍存在的。【TAD很火，可以引入到模型中】

HOMER

HOMER Motif Analysis - 根据ChIP-seq和ATAC-seq的peak结果寻找可能binding的motif

Finding Enriched Motifs in Genomic Regions (findMotifsGenome.pl) - 核心的脚本

Finding Instance of Specific Motifs - 对人类和小鼠而言最有用的代码，因为大部分motif已知，没必要做denovo的motif预测。 

Motif Databases included in HOMER - HOMER最常使用的一些motif数据库

HOMER的motif格式

jaspar - 最全的转录因子数据库

下载所有human的TF对应的motif：链接

下载JASPER上的转录因子及其motif数据库，这部分很重要是因为我们不仅需要motif的信息，而且需要motif对应的人类转录因子的信息。【需要用homer的工具进行格式转换】

cd ~/softwares/miniconda3/share/homer-4.10-0/update
curl -O http://jaspar.genereg.net/download/CORE/JASPAR2020_CORE_vertebrates_non-redundant_pfms_jaspar.txt

perl ./motifs/parseJasparMatrix.pl JASPAR2020_CORE_vertebrates_non-redundant_pfms_jaspar.txt > jaspar.motifs

curl -O http://jaspar.genereg.net/download/CORE/JASPAR2020_CORE_vertebrates_redundant_pfms_jaspar.txt

perl ./motifs/parseJasparMatrix.pl JASPAR2020_CORE_vertebrates_redundant_pfms_jaspar.txt > jaspar.motifs

　　

接下来就是全基因组的扫描了，找这些motif到底在哪binding【全基因组扫描过于费时，还是指定区域比较好】

perl ~/softwares/miniconda3/bin/scanMotifGenomeWide.pl ~/softwares/miniconda3/share/homer-4.10-0/update/motifs/vertebrates/jaspar.motifs hg38 -5p -bed -int -homer2 -p 10

选取promoter区域来扫描，看motif的结合区域

perl ~/softwares/miniconda3/bin/findMotifsGenome.pl encc-enhancer-atac.promt.Homer.bed hg38 promt.motif -p 10 -size 200 -find jaspar.motifs > promt.jaspar.txt

结果如何过滤？

For example: findMotifsGenome.pl ERalpha.peaks hg18 MotifOutputDirectory/ -find motif1.motif > outputfile.txt

The output file will contain the columns:

• Peak/Region ID
• Offset from the center of the region
• Sequence of the site
• Name of the Motif
• Strand
• Motif Score (log odds score of the motif matrix, higher scores are better matches)

根据Motif Score来过滤掉一些质量太低的比对。

这个结果仍然不是我想要的，我只想知道，某个motif是否在一个promoter或enhancer区域显著富集【相对于背景区域】

PositionID      Offset  Sequence        Motif Name      Strand  MotifScore
Peak_152271     -93     AGTAAG  Ahr::Arnt/MA0006.1/Jaspar       +       1.914416
Peak_152271     -85     CCCTTC  Ahr::Arnt/MA0006.1/Jaspar       +       2.895245
Peak_152271     -82     TTCAAG  Ahr::Arnt/MA0006.1/Jaspar       +       3.213699
Peak_152271     -75     GGCAGG  Ahr::Arnt/MA0006.1/Jaspar       +       4.644445
Peak_152271     -68     AGCTCC  Ahr::Arnt/MA0006.1/Jaspar       +       1.871856
Peak_152271     -65     TCCCTG  Ahr::Arnt/MA0006.1/Jaspar       +       6.391753
Peak_152271     -64     CCCTGG  Ahr::Arnt/MA0006.1/Jaspar       +       2.895245
Peak_152271     -63     CCTGGG  Ahr::Arnt/MA0006.1/Jaspar       +       2.895245
Peak_152271     -60     GGGATG  Ahr::Arnt/MA0006.1/Jaspar       +       4.687005
Peak_152271     -59     GGATGG  Ahr::Arnt/MA0006.1/Jaspar       +       1.508951
Peak_152271     -57     ATGGTG  Ahr::Arnt/MA0006.1/Jaspar       +       12.000225
Peak_152271     -55     GGTGAG  Ahr::Arnt/MA0006.1/Jaspar       +       11.552200

　　

HOMER的这个peak注释功能也是写得非常全面：

---

MEME

FIMO scans a set of sequences for individual matches to each of the motifs you provide - 看motif是否显著的富集在单独的序列里，需要把区域转换为fasta文件

CentriMo identifies known or user-provided motifs that show a significant preference for particular locations in your sequences - CentriMo可以做motif是否显著富集在一堆序列中，给出富集得分

看看这篇文章：Differential motif enrichment analysis of paired ChIP-seq experiments

FIMO

看一看sample output

FIMO Tutorial

首先提取出自己的fasta

可以用bedtools

bedtools flank -i encc-enhancer-atac.promt.bed -g chrom.sizes -b 300 > encc.enhc.atat.promt.f300.bed
bedtools sort -i tmp2.bed > tmp3.bed
bedtools merge -i tmp3.bed -c 4 -o collapse
bedtools getfasta -fo 

也可以用Homer的工具

---

参考：

Homer软件的介绍-最全面而详细的找motif教程