转录因子motif TSS区域富集分析 | motif enrichment | HOMER | FIMO | MEME

2021年09月15日 阅读数:7
这篇文章主要向大家介绍转录因子motif TSS区域富集分析 | motif enrichment | HOMER | FIMO | MEME,主要内容包括基础应用、实用技巧、原理机制等方面,希望对大家有所帮助。

【姐妹篇:HOMER | MEME | 转录因子的靶基因预测 | motif富集分析php

目的:假设你有一个感兴趣的TF,且这个TF的DNA binding motif已知,那么就能够预测这个TF的motif在其target gene的promoter(TSS附近)的富集信息。【需结合表观数据使用】html

 

主流的motif数据库node

  • JASPAR
  • dbcorrdb - SCENIC使用的
  • TRANSFAC® 7.0 Public 2005 and TRANSCompel 7.0 Public 2005 - 如今最新版要收费了

 

 HOMER很差用,且不更新,决定弃用。web

perl /home/lizhixin/softwares/miniconda3/bin/findMotifs.pl genelist.txt human motifResults/ -start -400 -end 100 -len 8,10,12 -find test.motif > result.txt数据库

/home/lizhixin/softwares/miniconda3/share/homer-4.10-0/update/motifs/vertebrates/jaspar.motifsapi

 

download
http://homer.ucsd.edu/homer/download.htmlbash

install
perl configureHomer.pl -installapp

dicard homer (stoped updateing)dom

 

如下为meme的教程:工具

 

motif格式问题

我关注的这个motif (ENCSR000ARI) 是Ezh2的,并无被收录在常规的TF数据库里,因此Homer里面没有。

dbcorrdb这个数据库里有,可是很老的JASPAR格式,能够转成meme格式。[dbcorrdb__EZH2__ENCSR000ARI_1__m5.png]

see format of motif:https://meme-suite.org/meme/doc/motif_conversion.html

prepare HDAC1.all.motif.sites file【须要手动修改格式】

>ENCSR000AQF-1-m1 HDAC1
A [1789 2376 3032 2941 1861 2873 1491 2709 2732 4134 5409 3443 1873 379 8456 2 1 576 1196 1808 ]
C [4094 5259 5024 3220 2666 1760 7390 4313 2199 2775 3023 3390 5021 8259 528 12442 4 2358 1081 3161 ]
G [4427 4089 4125 3207 4949 2036 2025 4972 7067 3844 2587 3396 1145 2814 3441 0 12439 215 8276 6114 ]
T [2134 720 263 3076 2968 5775 1538 450 446 1691 1425 2215 4405 992 19 0 0 9295 1891 1361 ]
jaspar2meme -pfm -bundle HDAC1.all.motif.sites > HDAC1.all.motif.meme
jaspar2meme -pfm -bundle ENCSR000ARI.all.sites > ENCSR000ARI.all.meme

 

能够用meme的Ceqlogo画出logo,plot motif

ceqlogo -i1 HDAC1.all.motif.meme -o logo.m1.png -f PNG
ceqlogo -i2 HDAC1.all.motif.meme -o logo.m2.png -f PNG
ceqlogo -i3 HDAC1.all.motif.meme -o logo.m3.png -f PNG
ceqlogo -i4 HDAC1.all.motif.meme -o logo.m4.png -f PNG
ceqlogo -i1 ENCSR000ARI.m5.meme -o logo.png -f PNG

 

check logo in website database
http://motifcollections.aertslab.org/v9/logos/

 

get TSS fasta
see project/scPipeline/motifEnrichment/get_TSS_fasta_R.ipynb

 

而后meme的FIMO能够直接把motif比对到指定的fasta

fimo --alpha 1 --max-strand -oc target HDAC1.all.motif.meme promoters_up500_down100.fasta
fimo --alpha 1 --max-strand -oc target HDAC1.all.motif.meme promoters_up1000_down500.fasta
fimo --alpha 1 --max-strand -oc target ENCSR000ARI.all.meme target.promt.TSS.fasta
fimo --alpha 1 --max-strand -oc random ENCSR000ARI.all.meme target.promt.TSS.random.fasta

  

提取promoter fasta:

这里就有个问题了,怎么从genome中提取出核心区域
由于promoter的定义不是很完善,这里我用的区域是TSS区 [-1000, 500]

其实Homer里面都写好了,但由于里面的perl过于复杂,最后仍是决定用R来作(核心就是防止位置溢出)

library(GenomicFeatures)

# library(TxDb.Hsapiens.UCSC.hg19.knownGene)
# library(BSgenome.Hsapiens.UCSC.hg19)

library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(BSgenome.Mmusculus.UCSC.mm10)

# PR <- promoters(TxDb.Mmusculus.UCSC.mm10.knownGene, upstream=2000, downstream=0)

# entrez geneID for a cell cycle control transcription
# factor, chr6 on the plus strand
# tmp.gene <- "429"  

# transcriptCoordsByGene.GRangesList <- transcriptsBy (TxDb.Mmusculus.UCSC.mm10.knownGene, by = "gene") [tmp.gene]
# a GrangesList of length one, describing three transcripts

# promoter.seqs <- getPromoterSeq (head(PR), Hsapiens, upstream=10, downstream=0)

trs <- transcriptsBy(TxDb.Mmusculus.UCSC.mm10.knownGene, by = "gene")

# trs <- keepStandardChromosomes(trs)

prom <- getPromoterSeq(trs, Mmusculus, upstream = 500, downstream=100)

head(prom)

  

gene ID转换问题

这里的list的name就是NCBI的genename,即Entrez ID(长见识了,同一个基因在不一样物种中的entrez ID是不一样的,因此你找的genecard上确定是跟mouse的对不上)

如下是常见的geneID转换的脚本,ensemble ID转Entrez ID

mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
att <- listAttributes(mart)

transcripts <- getBM(attributes=c("chromosome_name", "transcript_start", "transcript_end", 
                                  "entrezgene_id", "ensembl_gene_id","ensembl_transcript_id"),
                      filters = "biotype",
                      values  = c("protein_coding"),
                      mart    = mart)

  

tmp.names <- select(TxDb.Mmusculus.UCSC.mm10.knownGene, keys = names(prom), 
                    columns=c("TXNAME"), keytype="GENEID")

tmp.names$transcriptName <- sapply(strsplit(tmp.names$TXNAME,"\\."), function(x) x[1])

tmp.names$entrezID <- transcripts[tmp.names$transcriptName,]$entrezgene_id

  

最终画出了这个图:

成图代码

options(repr.plot.width=6, repr.plot.height=6)
g <- ggplot(motif.bind, aes(x=pos, y=-log10(p.value))) + 
    facet_wrap( ~ gene, ncol=1) +
    geom_point() +
    geom_vline(xintercept=0, linetype="dashed", color = "blue") +
    geom_linerange(aes(x=pos, ymax=-log10(p.value), ymin=0),
                 position = position_jitter(height = 0L, seed = 1L), size=0.1) +
    theme_bw() +
    theme(strip.background = element_blank(), strip.text = element_text(face="italic", size = 15, color = "black")) +
    theme(axis.text.x  = element_text(face="plain", angle=0, size = 8, color = "black", vjust=-1),
            axis.text.y  = element_text(face="plain", size = 10, color = "black"),
            axis.title =element_text(size = 15)) +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank()) +
    theme(legend.title = element_blank()) +
    labs(x = "\nRelative distance to TSS (bp)",y = "- log10(motif binding P-value)\n", title = " ") +
    scale_y_continuous(limits = c(0, 7), expand = c(0, 0)) #+
    #scale_fill_manual(values = brewer.pal(8,"Set1")) 
g

  

其实要是Homer有对应的TF motif数据,能够作得更加高效。不用提fasta,不用id转换,一行代码搞定。

findMotifs.pl  genelist.txt mouse motifResults/ -start -400 -end 100 -len 8,10,12 -find test.motif > result.txt

  

Homer和meme的对比

  • Homer是用perl写的,一个脚本包含了不少功能,你想提取其中一个子功能很难
  • meme是C++的,执行效率更高,功能都拆分了 【一键化用Homer;精细化快速化用meme】
  • Homer彷佛不受待见,JASPAR里下载格式就没有Homer格式的,反而有meme的 【meme犹如瑞士军刀,功能很是全面,支持多种motif格式】

 

历史目录:

/home/lizhixin/project2/motif/homer - v1

/home/lizhixin/project/scPipeline/motifEnrichment - v2

 

发现的另外一个工具:https://www.regulatory-genomics.org/

 

展望:

这里只考虑了TSS flanking region

有ChIP-seq和ATAC-seq peak的数据能够confirm这个denovo的结果

再加入Hi-C等三维基因组的数据就能够整合promoter和enhancer了。  

 

下一专题:motif格式专题

 

了解一下数据分析的生物学基础:

『珍藏版』Nature综述 | 基因表达调控的新模型

Organization and regulation of gene transcription

 

 

全基因组(~20000个基因)批量扫描1541个motif

# download motif
wget http://www.jstacs.de/downloads/motifs_jaspar.txt

# change format in sublime, must have ',' at the end
# motif_(\d+)
# motif_$1 ,
# example
>CTCF_ENCSR000AKB-1_motif_1 ,
10474 6567 9974 12621 4732 79 32529 2514 7204 39703 37 18365 2206 427 1630 7094 19180 4523 6659 16417
11903 15167 9425 5070 36198 46507 1697 26396 19496 1126 22 46 1279 18 1817 33939 2120 23129 16325 11761
12421 7303 16734 21187 2705 14 5847 16524 3318 3237 46565 28151 26304 46054 39302 1153 23209 14621 4807 13486
11854 17615 10519 7774 3017 52 6579 1218 16634 2586 28 90 16863 153 3903 4466 2143 4379 18861 4988

# get wanted
# now run all motif: 1541

# to meme format
jaspar2meme -pfm -bundle motifs_jaspar.sublime.txt > jaspar.all.meme
# check
cat jaspar.all.meme  | grep 'MOTIF' | wc -l

# align to all promoter
fimo --alpha 1 --max-strand -oc all.TF.targets jaspar.all.meme all.hs.hg19.TSS.up400.down100..fasta

# submit job
echo "fimo --alpha 1 --max-strand -oc all.TF.targets.mm10 jaspar.all.meme all.Mm.mm10.TSS.up400.down100.fasta" | qsub -V -N Mm10_motif -q small_ext -l nodes=1:ppn=2,walltime=60:00:00,mem=10gb

 

数据库问题:

上面的DBcorrDB数据库里只有179个TF

想要更多能够查JASPAR数据库