HOMER | MEME | 转录因子的靶基因预测 | motif富集分析-编程之家

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教程