寫在前面
文獻(xiàn)中常用的有DREME和HOMER壤短，這次先搞定DREME缤剧，下次再寫HOMER客叉。
使用MEME套件中的DREME乃戈，用于鑒定meRIP-Seq數(shù)據(jù)中peak的motif屡立。motif是序列中反復(fù)出現(xiàn)的由幾個(gè)堿基組成的小片段直晨，可能有特定的生物學(xué)意義。不過看meRIP-seq的文章里，這塊的分析在結(jié)果中一般體現(xiàn)為一句話勇皇，交代找到的顯著motif與已有報(bào)道一致罩句，或者說多個(gè)發(fā)育階段鑒定的保守motif是什么。
DREME的原理（http://meme-suite.org/doc/dreme-tutorial.html?man_type=web）：
- 輸入文件：positive sequences和negative sequences(背景序列敛摘，如果沒有提供门烂，會(huì)根據(jù)positive sequences自動(dòng)生成）；
- 尋找positive sequences中的3-8bp(默認(rèn)長(zhǎng)度)片段, 計(jì)算每個(gè)片段在positive和negative中出現(xiàn)的次數(shù)着撩，通過Fisher's exact test 確定顯著性诅福，按照P-value排序成regular expression motifs.
- 選取前100個(gè)motif，將其中1個(gè)堿基替換成ambiguity code,重新計(jì)算p-value拖叙；重復(fù)該過程氓润，使得每個(gè)堿基都可以有ambiguity code，得到generalized RE motifs.（相當(dāng)于consensus sequences）
- 針對(duì)上一步得到的generalized RE motifs,計(jì)算每個(gè)motif在positive和negative中出現(xiàn)的次數(shù)薯鳍，通過Fisher's exact test 確定顯著性,按照設(shè)定的E-value值輸出motif及相應(yīng)的頻率矩陣咖气。
實(shí)戰(zhàn)
01 安裝
下載軟件：http://meme-suite.org/doc/download.html
#install
mkdir meme
tar zxf meme_5.0.1.tar.gz
cd meme_5.0.1
./configure --prefix=/root/cc/biosoft/bin/meme --with-url=http://meme-suite.org --enable-build-libxml2 --enable-build-libxslt
make
make test
make install
#加入PATH(~/.bash_profile)
PATH=/root/cc/biosoft/bin/meme/bin:$PATH
source ~/.bash_profile
#test
dreme
#根據(jù)提示需要安裝imagemagick
yum install ImageMagick
convert -version
02 實(shí)戰(zhàn)
輸入文件：fasta格式的序列文件
dreme -p input.fa -oc outDir -dna -png -eps -e 0.05 -t 18000 -m 20
#-dna, -rna : 默認(rèn)是dna
#-png, -eps ： 輸出圖片格式
# -e, -t, -m : 是設(shè)置軟件結(jié)束運(yùn)行的參數(shù)，不設(shè)置的話挖滤，默認(rèn)在e-value>e時(shí)停止崩溪，否則會(huì)一直運(yùn)行下去。
# -e: e-value斩松，0.05（default) 一直search伶唯，直到下一個(gè)motif 的e-value >e
# -t: 運(yùn)行時(shí)間 一直search，直到達(dá)到這個(gè)時(shí)間
# -m: 找到的motif數(shù)目 一直search惧盹，直到找到的motif數(shù)目達(dá)到這個(gè)值
#還可以設(shè)置Core Motif Width：--mink 最小值乳幸；--maxk 最大值
設(shè)置為-dna時(shí)，可以識(shí)別U(等同于T）,但輸出LOGO時(shí)使用T.
設(shè)置為-rna時(shí)钧椰，仍然能識(shí)別T(等同于U）粹断，但輸出LOGO時(shí)使用U. 同樣的序列和參數(shù)，設(shè)置為-dna和-rna嫡霞，找到的motif是不一樣的瓶埋。Why? 沒找到說明
是不是灰常簡(jiǎn)單？嘻嘻(#.#)
03 LOGO編輯
DREME生成的xml文件中诊沪，有每個(gè)motif的PWM matrix养筒，里面有motif的每個(gè)位置上每種堿基的probability《艘Γ可以用來直接畫LOGO晕粪，方便的工具有如下幾個(gè)。
seqLogo可以調(diào)整橫縱坐標(biāo)寄锐、堿基高度按比例/絕對(duì)大小展示等兵多。
source("https://bioconductor.org/biocLite.R")
biocLite("seqLogo")
require(seqLogo)
mFile <- read.table("inp")
mFile
# V1 V2 V3 V4 V5 V6
#1 0 0.000000 1 0 1 0.000000
#2 0 0.000000 0 1 0 0.000000
#3 1 0.907213 0 0 0 0.636825
#4 0 0.092787 0 0 0 0.363175
p<-makePWM(m)
seqLogo(p)
seqLogo(p,ic.scale=FALSE)
不足：目前seqLogo只能識(shí)別DNA alphabet尖啡，無法識(shí)別RNA，Protein剩膘。
motifStack用途更廣泛衅斩，這里僅用于把motif中的T替換成U。
需要安裝ghostscript怠褐，然后將安裝路徑設(shè)置為環(huán)境變量R_GSCMD.
比如我安裝的是64位的畏梆，安裝在C:\Program Files\gs\gs9.23，按如下設(shè)置：
Sys.setenv(R_GSCMD=file.path("C:", "Program Files", "gs",
"gs9.23", "bin", "gswin64c.exe"))
將矩陣中的T換成U奈懒。
source("https://bioconductor.org/biocLite.R")
biocLite("motifStack")
require(motifStack)
Sys.setenv(R_GSCMD=file.path("C:", "Program Files", "gs",
"gs9.23", "bin", "gswin64c.exe"))
pcm=read.table("inp")
#0 0.000000 1 0 1 0.000000
#0 0.000000 0 1 0 0.000000
#1 0.907213 0 0 0 0.636825
#0 0.092787 0 0 0 0.363175
rownames(pcm) <- c("A","C","G","U")
motif <- new("pcm", mat=as.matrix(pcm), name="KD")
plot(motif,ncex=2,xaxis=F) #ncex: fond size of name; xaxis 不顯示橫坐標(biāo)
04 Note
meRIP-seq的數(shù)據(jù)奠涌，分析得到的peak有時(shí)會(huì)跨區(qū)域。比如兩個(gè)相鄰的exon都是peak磷杏，但是得到的是一個(gè)包含這兩個(gè)exon的區(qū)域溜畅，那么，提取序列時(shí)极祸，應(yīng)該提取的是真實(shí)的peak區(qū)域慈格，即相應(yīng)的exon區(qū)域。如果是exomePeak的結(jié)果遥金，生成的BED12文件里包含block的信息浴捆，可以利用bedtools getfasta -split直接提取相應(yīng)block的序列，再“串聯(lián)”起來稿械；如果是普通的BED文件（chr start end)选泻，那么就需要結(jié)合exon的坐標(biāo)(比如GTF，GFF等注釋文件),自己生成相應(yīng)的BED12文件了美莫。
