主流的motif数据库
dbcorrdb - SCENIC使用的
TRANSFAC® 7.0 Public 2005 and TRANSCompel 7.0 Public 2005 - 现在最新版要收费了
motif格式问题
我关注的这个motif (ENCSR000ARI) 是Ezh2的,并没有被收录在常规的TF数据库里,所以Homer里面没有。
dbcorrdb这个数据库里有,但是很老的JASPAR格式,可以转成meme格式。[dbcorrdb__EZH2__ENCSR000ARI_1__m5.png]
jaspar2meme -pfm -bundle ENCSR000ARI.all.sites > ENCSR000ARI.all.meme
可以用meme的Ceqlogo画出logo
ceqlogo -i1 ENCSR000ARI.m5.meme -o logo.png -f PNG
然后meme的FIMO可以直接把motif比对到指定的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
这里就有个问题了,怎么从genome中提取出核心区域
因为promoter的定义不是很完善,这里我用的区域是TSS区 [-400, 100]
其实Homer里面都写好了,但因为里面的perl过于复杂,最后还是决定用R来做(核心就是防止位置溢出)
library(GenomicFeatures)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(BSgenome.Mmusculus.UCSC.mm10)
trs <- transcriptsBy(TxDb.Mmusculus.UCSC.mm10.knownGene, by = "gene")
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的对比
展望:
这里只考虑了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
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
jaspar2meme -pfm -bundle motifs_jaspar.sublime.txt > jaspar.all.meme
cat jaspar.all.meme | grep 'MOTIF' | wc -l
fimo --alpha 1 --max-strand -oc all.TF.targets jaspar.all.meme all.hs.hg19.TSS.up400.down100..fasta
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数据库
手机扫一扫
移动阅读更方便
你可能感兴趣的文章