HOMER | MEME | 转录因子的靶基因预测 | motif富集分析
2024年04月30日
根据fasta来找motif enrichment
findMotifs.pl prom1_e3.fasta fasta motifResults/ -fasta background.fa
2023年03月20日
想找到一个peak文件里富集的motif，MEME可以用网页版，HOMER可以用本地版，速度也很快，最多4行命令。
conda install -c bioconda homer perl /home/zz950/softwares/miniconda3/envs/cutrun/share/homer/configureHomer.pl -list perl /home/zz950/softwares/miniconda3/envs/cutrun/share/homer/configureHomer.pl -install mm10 findMotifsGenome.pl HDAC1.narrowPeak mm10 motif -e 0.01
BETA也能做，还能整合RNA-seq，但不确定其motif数据库是否全面，还是用“金标准”homer吧。
http://localhost:17435/notebooks/projects/cut_run/HDAC-b2/pipeline/DiffBind-merge-integrate.ipynb#venn-ATAC-and-AC-loss
findMotifsGenome.pl HT115_M60a_ATAC.HDAC.gain.bed hg38 motif -e 0.01 findMotifsGenome.pl HT115_M60a_ATAC.HDAC.loss.bed hg38 motif -e 0.01
findMotifsGenome.pl HT115_M60a_ATAC.HDAC.gain.bed hg38 HT115_M60a_ATAC.HDAC.gain_homer_motif -e 0.01 findMotifsGenome.pl HT115_M60a_ATAC.HDAC.loss.bed hg38 HT115_M60a_ATAC.HDAC.loss_homer_motif -e 0.01
【姐妹篇，更新：转录因子motif TSS区域富集分析 | motif enrichment | HOMER | FIMO | MEME】
转录因子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
首先提取出自己的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的工具
参考：
