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

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

目的：假设你有一个感兴趣的TF，且这个TF的DNA binding motif已知，那么就可以预测这个TF的motif在其target gene的promoter（TSS附近）的富集信息。【需结合表观数据使用】

主流的motif数据库

dbcorrdb - SCENIC使用的

TRANSFAC ® 7.0 Public 2005 and TRANSCompel 7.0 Public 2005 - 现在最新版要收费了

HOMER不好用，且不更新，决定弃用。

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.motifs

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

install
perl configureHomer.pl -install

dicard homer (stoped updateing)

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")) 
其实要是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数据库