Motif預(yù)測.docx

Motif預(yù)測 在許多生物學(xué)研究過程中，我們需要對真核和原核生物的轉(zhuǎn)錄調(diào)控，蛋白質(zhì)結(jié)構(gòu)活性位點，以及DNA、RNA的酶切位點進(jìn)行識別。在這些實際過程中，我們通常會碰到下面兩個問題： （1）我們通常需要從未知的數(shù)據(jù)盡可能多的發(fā)掘有用的信息，也就是說，在探索實驗過程中，由于科研工作者一開始并不知道哪些信息是我們真正“有用”的，在沒有其他太多的信息可以利用的時候，一個很自然的想法就是我們是否能知道哪些因素是這些數(shù)據(jù)共有部分，也就是我們通常所謂的數(shù)據(jù)的公共特征提取，或者叫motif的預(yù)測或模式識別，這里的“模式”（motif）可以簡單的理解為特定數(shù)據(jù)的共同特征。 （2）另一方面，假設(shè)當(dāng)我們已經(jīng)知道了某個特征，我們需要把具有這些特征的數(shù)據(jù)都收集起來，那么，我們需要在大量數(shù)據(jù)中挑選符合我們要求地數(shù)據(jù)，這個過程就是所謂的“模式匹配”?！澳Ｊ阶R別”和“模式匹配”是生物信息學(xué)輔助實驗生物學(xué)的一個重要手段。關(guān)于這些模式的具體形式可以是相當(dāng)廣泛的。從廣義上講，即使是多序列比對尋找保守區(qū)、蛋白質(zhì)結(jié)構(gòu)預(yù)測等問題，也可以納入模式識別的范疇。但下面我們要介紹是以圍繞轉(zhuǎn)錄因子結(jié)合位點（TFBS）為代表的一類模式識別的程序，是大家傳統(tǒng)意義上的模式識別問題。這一類模式（motif）的最大的特點就是大家共有特征較短，一般一個TFBS位點的長度在5-20bp左右，而且信號比較靈活多變，所以這些調(diào)控元件在通常情況下不適合直接使用多序列聯(lián)配的方式來尋找，而是需要一些專門的算法來解決這個問題。 下面我們介紹幾個較著名的模式（motif）預(yù)測軟件。7.3.1 MEME/MAST系統(tǒng) MEME和MAST是由T.L.Bailey、Charles Elkan和Bill Noble合作開發(fā)的一套搜索motif的程序組合套件。是目前生物信息學(xué)領(lǐng)域，對motif預(yù)測方面最著名的程序之一。 這個組件的兩個程序分別執(zhí)行的是motif的預(yù)測和搜索兩個不同的功能。其中MEME的全稱是Multi EM for motif elicitation，是一個基于EM算法的一個motif預(yù)測程序。而MAST的全稱是Motif Alignment & search Tool，是一個motif搜索組件。下載： 這個軟件包可以從MEME的官方網(wǎng)站免費下載，我們這里使用的最新版本是version 3.5.4，其源代碼可以從如下鏈接下載：/downloads/meme_3.5.4.tar.gz安裝： 這里介紹在Linux/unix下的安裝步驟： 1、把源代碼文件上載本地的Linux/Unix指定目錄: /mnt/disk1/motif_workplace/ 2、解壓： $ gunzip -c meme_3.5.4.tar.gz | tar xvf - $ cd meme_3.5.4 3、配置：通常情況下，可以使用默認(rèn)參數(shù)。如果需要指定安裝路徑，可以使用參數(shù)-prefix=your-install-path$ ./configure或$ ./configure -prefix=your-install-path 如果配置成功，會出現(xiàn)以下畫面：=Configuration parameters=Install path: your-install-pathInstall UID: usridVersion: 3.5.4C compiler: gccC compiler flags: -Wall -DUNIX -D_USE_FIXED_PROTOTYPES_ -O3Linker: /usr/bin/ldSpecial Libs: -lmMPICCMPIINCMPILIBDIRMPIFLAGS -DPARALLEL -OMPIRUNURL: http:/usr_home/memePERL: /usr/bin/perlMEME_LOGS: $prefix/LOGSRun the following commands to compile, test and install meme: make make test make install 4、編譯和安裝： 在命令行鍵入make開始編譯， $ make 當(dāng)編譯完成后，鍵入命令：make test，測試， $ make test 當(dāng)返回提示： Meme test for-Dataset Model result-crp0 oops OKcrp0 zoops OKcrp0 tcm OK.PASS: runcheck=All 1 tests passed=.說明測試成功，可以鍵入命令：make install進(jìn)行安裝。$ make install 其中測試畫面中出現(xiàn)的三個模型oops、zoops和tcm是MEMmotif預(yù)測中最重要的選擇參數(shù)，分別表示對任意給定的motif來講，每條序列當(dāng)且僅當(dāng)出現(xiàn)一次、每條序列至多有一個以及每條序列出現(xiàn)兩個或兩個以上三種情況。對于一個給定的序列，我們可以根據(jù)我們需要研究的問題，選擇這三個模型中的一個。關(guān)于meme參數(shù)的選擇將在后面參數(shù)一部分詳細(xì)討論。使用： 當(dāng)我們在命令行鍵入meme命令，可以看到meme的參數(shù)列表和幫助文檔。$ memeUSAGE:meme optional arguments file containing sequences in FASTA format-h print this message-dna sequences use DNA alphabet-protein sequences use protein alphabet-mod oops|zoops|anr distribution of motifs-nmotifs maximum number of motifs to find-evt stop if motif E-value greater than -nsites number of sites for each motif-minsites minimum number of sites for each motif-maxsites maximum number of sites for each motif-wnsites weight on expected number of sites-w motif width-minw minumum motif width-maxw maximum motif width-nomatrim do not adjust motif width using multiple alignment-wg gap opening cost for multiple alignments-ws gap extension cost for multiple alignments-noendgaps do not count end gaps in multiple alignments-bfile name of background Markov model file-revcomp allow sites on + or - DNA strands-pal force palindromes (requires -dna)-maxiter maximum EM iterations to run-distance EM convergence criterion-prior dirichlet|dmix|mega|megap|addone type of prior to use-b strength of the prior-plib name of Dirichlet prior file-spfuzz fuzziness of sequence to theta mapping-spmap uni|pam starting point seq to theta mapping type-cons consensus sequence to start EM from-text output in text format (default is HTML)-maxsize maximum dataset size in characters-nostatus do not print progress reports to terminal-p use parallel version with processors-time quit before CPU seconds consumed-sf print as name of sequence file. .當(dāng)我們在命令行鍵入mast命令，可以看到mast的參數(shù)列表和幫助文檔。$ mastUSAGE:mast optional arguments . file containing motifs to use; may be a MEME output file or a file with the format given below or-d database to search with motifs or-stdin read database from standard input; Default: reads database specified inside -c only use the first motifs-a is assumed to contain motifs in the format output by make_logodds and is their alphabet; -d or -stdin must be specified when this option is used-stdout print output to standard output instead of file-text output in text (ASCII) format; (default: hypertext (HTML) format)-sep score reverse complement DNA strand as a separate sequence-norc do not score reverse complement DNA strand-dna translate DNA sequences to protein-comp adjust p-values and E-values for sequence composition-rank print results starting with best (default: 1)-smax print results for no more than sequences (default: all)-ev print results for sequences with E-value (default: 10)-mt show motif matches with p-value mt (default: 0.0001)-w show weak matches (mtp-valuemt*10) in angle brackets-bfile read background frequencies from -seqp use SEQUENCE p-values for motif thresholds (default: use POSITION p-values)-mf print as motif file name-df print as database name-minseqs lower bound on number of sequences in db-mev + use only motifs with E-values less than -m + use only motif(s) number (overrides -mev)-diag nominal order and spacing of motifs-best include only the best motif in diagrams-remcorr remove highly correlated motifs from query-brief brief output-do not print documentation-b print only sections I and II-nostatus do not print progress report-hit_list print machine-readable list of all hits only; implies-text. . 下面我們開始逐一解釋這些基本參數(shù)的含義以及在實際應(yīng)用中需要注意哪些問題。 MEME參數(shù)的說明： 1）基本參數(shù)篇 我們可以看到，meme程序的參數(shù)很多，為了更好的說明問題，meme把參數(shù)大致地分幾類，我們按照meme文檔里的分類模式，簡要說明： 1. ：meme輸入文件格式是fasta格式，主要形式如下： seq1 GDIFYPGYCPDVKPVNDFDLSAFAGAWHEIAK seq2 GDMFCPGYCPDVKPVGDFDLSAFAGAWHELAK seq3 QKVAGTWYSLAMAASDISLLDAQSAPLRVYVEELKPTPEGDLEILLQKW 可以是DNA序列的fasta文件，也可以蛋白質(zhì)序列的fasta格式文件。 2.可選的幾個參數(shù)類：字母表類；分布類；搜索類；系統(tǒng)類 字母表參數(shù)類： -dna/-protein，說明輸入序列類型是DNA還是蛋白質(zhì)。 分布模型類： -mod有3個可選的值，就是前面提到的3個可選模型： oops：(One Occurrence Per Sequence) 意思是每種motif在每條序列中只都出現(xiàn)一次。 zoops：(Zero or One Occurrence Per Sequence)，意思是每種motif在中至多出現(xiàn)一次。 anr：(Any Number of Repetitions)，意思是每種motif可以出現(xiàn)任意數(shù)目。這個其實是前面提到的tcm模型和zoops功能的一個組合。 搜索參數(shù)類 MEME預(yù)測的motif的過程其實是一個尋找給定目標(biāo)函數(shù)最優(yōu)的過程，通過目標(biāo)函數(shù)計算每個可能的motif的對數(shù)似然比，近似用E值表示，所以候選種子的E值越小，表示這個種子是真正的motif的可能性越大，所以結(jié)果文件會按E值從小到大的順序輸出。 所以下面的一些有關(guān)輸入輸出的參數(shù)是可以根據(jù)我們實際需要，進(jìn)行調(diào)整。 -nmotifs ：選擇輸出motif個數(shù)n，表示meme的運行結(jié)果文件輸出n個不同motif，這個參數(shù)的默認(rèn)值是1。 -evt ：E-value的閾值，當(dāng)這個參數(shù)為p，表示我們只輸出E值比p小的那些motif。所謂E值是一個描述motif模型可信度的一個標(biāo)志，一般而言，E值越小，結(jié)果可信度越高。這個參數(shù)的默認(rèn)值是無窮大，也就是在默認(rèn)狀態(tài)下，這個參數(shù)不起作用。 -nsite：表示每個motif出現(xiàn)的期望。 -minsites ：表示每個motif出現(xiàn)的期望的下限。 -maxsites ：表示每個motif出現(xiàn)的期望的上限。 這3個參數(shù)不是相互獨立的，只能取-nsite和其他兩個參數(shù)不能同時選擇。當(dāng)我們?nèi)?nsite時，只有那個給定數(shù)目的種子作為程序候選的motif進(jìn)行運算。而-minsites和-maxsites參數(shù)則給出這個期望數(shù)目的一個下限和上限約束的范圍。只有落在這個范圍內(nèi)的種子才給予考慮。注意，在OOPS參數(shù)下，這些參數(shù)的設(shè)置為無效設(shè)置。而在其他兩個模型參數(shù)的條件下，默認(rèn)值為-minsites：sqrt（序列數(shù)） -maxsites：當(dāng)ZOOPS時，默認(rèn)為序列數(shù)n 當(dāng)anr時，默認(rèn)為min（序列數(shù)5，50）。 -wnsites : 表示nsite初始的權(quán)重。它是控制種子滿足nsites或minsites和maxsites限制的一個貢獻(xiàn)值，這個權(quán)值是0,1)的一個值，這個值越大，表示motif傾向與滿足限制的可能越大。默認(rèn)值是用0.8。 -w ：表示motif種子的寬度。 -minw ：表示motif種子寬度的下限 -maxw ：表示motif種子寬度的上限 也就是說，當(dāng)-w給定，程序只嘗試寬度為-w的種子，否則，嘗試motif種子寬度在-minw和-maxw之間的值。默認(rèn)的-minw是8 -maxw是50。 -text：輸出文本格式（默認(rèn)的輸出格式是HTML文件） -maxsize ：最大的數(shù)據(jù)集大小 -sf ：打印作為序列文件名 -nostatus：不要在終端輸出打印報告 -time ：運行在時間前自動退出程序MEME高級參數(shù)技巧篇-nomatrim：這個參數(shù)的選用令MEME跳過用多序列聯(lián)配截短的步驟。-wg ：空位罰分，默認(rèn)值是11-ws ：空格罰分，默認(rèn)值是1-noendgaps：對最后的空位不罰分，默認(rèn)值是對最后的空位罰分。-bfile ：背景概率分布文件。-revcomp：翻轉(zhuǎn)序列，即motif可以在本身序列上，也可以在補鏈上-pal：考慮回文結(jié)構(gòu)-maxiter ：表示EM迭代次數(shù)，默認(rèn)是50次-distance ：表示收斂規(guī)則，表示兩個連續(xù)的頻率矩陣的歐幾里得距離小于這個值是，跳出迭代。默認(rèn)值是0.001。-prior ：先驗概率模型參數(shù)dirichlet：簡單Dirichlet先驗，作為-dna和-alph的默認(rèn)先驗。dmix：混合Dirichlets先驗，-protein的默認(rèn)分布。mega：方差非常小的混合Dirichlets分布，方差除以數(shù)據(jù)集變換。megap：在最后一次迭代用dmix參數(shù)，其他使用mega參數(shù)。addone：給每一個觀察值+1。（laplace 法則）。-b ：先驗?zāi)Ｐ偷膹姸取?= 0：表示使用固有的先驗強度。默認(rèn)值通常：0.01：對dirichlet先驗，=0.01。0 if：對于dmix而言，默認(rèn)=0。-plib ：存放Dirichlet先驗的文件，文件格式如prior30.plib所示。-spfuzz ：作模糊操作具體如spmap；-spmap ：使用映射函數(shù)，uni：增加一個先驗，并把一個子串轉(zhuǎn)為theta的估計，這里的theta是EM算法中的一個估計參數(shù)。默認(rèn)的參數(shù)是-spfuzz ：0.5pam：把PAM矩陣的子串轉(zhuǎn)化為theta的估計。默認(rèn)的參數(shù)-spfuzz ：120 (PAM 120)-cons ：廢除采樣的起始點，使用給定的字符串作種子。這個參數(shù)在一直特定motif的基本形式時會很有用。系統(tǒng)參數(shù)： 最后幾個參數(shù)是和大型機(jī)系統(tǒng)有關(guān)，MEME程序可以在服務(wù)器上支持并行運算，因此最后幾個參數(shù)只有在一定的服務(wù)器環(huán)境下可以使用。 -p ：用個CPU并行運算 以上是MEME各參數(shù)的說明，對于參數(shù)的選擇，需要大家根據(jù)具體問題靈活調(diào)整，而對于這里提到的高級參數(shù)部分，如果大家對MEME的概率模型不熟悉的話，建議使用默認(rèn)參數(shù)即可。 例1：假設(shè)我們有5條水稻的TSS上游調(diào)控序列，序列文件rice9311_sample.fa，我們需要尋找他們的motif，并把結(jié)果保存在rice_meme.html里面。 在命令行中輸入命令如下： $ meme rice9311_sample.fa -dna -nmotifs 4 -mod zoops -minw 5 -maxw 15 rice_meme_htmlFormat.html 或 $ meme rice9311_sample.fa -dna -nmotifs 4 -mod zoops -minw 5 -maxw 15 -textrice_meme_txtFormat.out 下表列出rice_meme_txtFormat.out的一部分文件如下，我們將分段解析結(jié)果文件包含的內(nèi)容，從整體上看，MEME的結(jié)果文件主要由文件頭、基本數(shù)據(jù)信息、輸入命令信息、motif結(jié)果條目以及motif綜合信息五大模塊組成。下面我們以例1運行的結(jié)果為例，解析各區(qū)段的作用。 （1）頭文件： 我們可以看到，MEME結(jié)果的文件頭可以分兩大部分，一部分包含MEME版本信息當(dāng)前版本是（v3.5.4）以及官方主頁。*MEME - Motif discovery tool*MEME version 3.5.4 (Release date: 3.5.4)For further information on how to interpret these results or to geta copy of the MEME software please access .This file may be used as input to the MAST algorithm for searchingsequence databases for matches to groups of motifs. MAST is availablefor interactive use and downloading at .* 頭文件的另一部分是關(guān)于文獻(xiàn)引用的說明，在論文中如果使用MEME的結(jié)果，可以應(yīng)用下面給出的T.Bailey1994的文獻(xiàn)。希望使用MEME的讀者能在論文中正確的引用MEME的文獻(xiàn)。*REFERENCE*If you use this program in your research, please cite:Timothy L. Bailey and Charles Elkan,Fitting a mixture model by expectation maximization to discovermotifs in biopolymers, Proceedings of the Second InternationalConference on Intelligent Systems for Molecular Biology, pp. 28-36,AAAI Press, Menlo Park, California, 1994.*（2）數(shù)據(jù)集基本信息。 這段數(shù)據(jù)保存了輸入數(shù)據(jù)集的一些基本的信息。包括輸入序列文件名和字符表以及每條輸入序列的序列名稱, 序列的權(quán)重，和序列長度。其中這里的序列權(quán)重是表示這條序列在預(yù)測motif中發(fā)揮作用大小的一個度量。*TRAINING SET*DATAFILE= rice9311_sample.faALPHABET= ACGTSequence name Weight Length Sequence name Weight Length- - - - - -BGK02457 1.0000 600 BGK03149 1.0000 600BGK04166 1.0000 600 BGK01655 1.0000 600BGK03537 1.0000 600 BGK02838 1.0000 600BGK04759 1.0000 600* 我們可以看到在這里我們的輸入序列文件名：rice9311_sample.fa。因為輸入是單純的DNA序列，序列中并沒有出現(xiàn)DNA的簡并形式，所以這里的字符表是A,C,G,T，我們還可以看到我們輸入的序列名稱，BGKXXXX，一般情況下，序列的權(quán)重是1，現(xiàn)在我們的輸入序列的長度、都是600bp。 （3）命令行參數(shù)。這一部分主要打印運行MEME時的人工參數(shù)以及程序默認(rèn)的參數(shù)。*COMMAND LINE SUMMARY*This information can also be useful in the event you wish to report a problem withthe MEME mand: meme rice9311_sample.fa -dna -nmotifs 4 -mod zoops -minw 5 -maxw 15 -dir/home/goodgoodmodel: mod= zoops nmotifs= 4 evt= inf objectfunction= E-value of product of p-valueswidth: minw= 5 maxw= 15 minic= 0.00width: wg= 11 ws= 1 endgaps= yesnsites: minsites= 2 maxsites= 7 wnsites= 0.8theta: prob= 1 spmap= uni spfuzz= 0.5em: prior= dirichlet b= 0.01 maxiter= 50 distance=1e-05data: n= 4200 N= 7strands: +sample: seed= 0 seqfrac= 1Letter frequencies in dataset:A 0.267 C 0.277 G 0.211 T 0.245Background letter frequencies (from dataset with add-one prior applied):A 0.267 C 0.277 G 0.211 T 0.245* 我們可以看到上表列出的參數(shù)中,在我們命令中，我們使用的是比較常用的參數(shù)，除了文件名和序列類型（-dna/-protein）以外，我們還使用了模型參數(shù)-mod zoops，在這里值得注意的是，這個模型參數(shù)的指定是MEME結(jié)果能否與生物學(xué)意義緊密聯(lián)系的重要參數(shù)之一。一般通常沒有特殊假設(shè)的情況下使用ZOOPS模型比較合適?？傊‘?dāng)模型的選擇將是我們解決問題捷徑。 另一個十分重要的參數(shù)是-revcomp，這個參數(shù)是用來指定程序是否讓序列的補鏈上也參與motif預(yù)測。因為這個參數(shù)的選擇是和我們的數(shù)據(jù)集以及我們期望解決的生物學(xué)問題緊密關(guān)聯(lián)，所以這個參數(shù)是不可以默認(rèn)的。在我們這個例子的問題中，由于我們的數(shù)據(jù)是水稻基因5端上游啟動子序列，所以我們并不需要補鏈信息，所以我們沒有使用-revcomp參數(shù)。關(guān)于其他的程序參數(shù)如何選擇，很多時候并沒太嚴(yán)格的規(guī)定，希望大家能在使用中仔細(xì)揣摩，因為參數(shù)的選擇問題有些時候更象是一門藝術(shù)。結(jié)果： 這一區(qū)段主要是MEME的預(yù)測得到的motif的形式和描述和評估。在這個部分，結(jié)果文件可以根據(jù)找到的motif的E值由小到大的順序，顯示我們在參數(shù)中指定輸出的最顯著的若干的motif。在這個區(qū)段中，MEME返回了多種形式的motif，并且給出了每個motif的統(tǒng)計顯著性的E值和P值。 對于這個區(qū)段里列出每個motif，都會包含motif的基本信息和描述、位置p 值的排序、塊的圖解和motif的塊形式、motif的位置特異的打分矩陣和位置特異概率矩陣形式以及正則表達(dá)式形式。1 motif 1 的基本信息。*MOTIF 1 width = 13 sites = 7 llr = 88 E-value = 7.9e-001* 這里給出這個motif的最基本信息，包括motif的寬度、位點的匹配數(shù)目、對數(shù)似然比(llr)和這個motif的E值。這里llr的基本公式可以用llr = log (Pr(sites | motif) / Pr(sites | back)計算，其中 Pr(sites|motif)條件概率，表示在motif條件下匹配的概率；Pr(sites | back)表示隨機(jī)背景條件下出現(xiàn)位點匹配的概率。這里的llr值和E值都是反映這個motif是否是真正motif的指標(biāo)，其中l(wèi)lr越大，E值越小，則說明這個motif越可靠。2motif 1 的描述- Motif 1 Description-Simplified A a143a:a:4:69:pos.-specific C :13:1probability G :637:a:966419matrix T :1:1:4:bits 2.2 * 2.0 * 1.8 * * 1.6 * * *Information 1.3 * * *content 1.1 * *(18.0 bits) 0.9 * * 0.7 * * 0.4 * 0.2 * 0.0 -Multilevel AGAGAGAGGGAAGconsensus CA ATGsequence G- 從上面的描述片斷，我們可以看到三種不同形式的motif的表示形式反映motif不同位置的保守性，一種是用簡單的頻數(shù)概率矩陣，一種是用motif各位置信息熵值，熵值越大，說明這個位點的保守性越高。我們看到motif從左往右的第5位、第6位和第7位的信息熵較大分別是是1.8、2.2和1.8，從頻度矩陣上也可以看到這幾個位置最保守。這里motif表示的第三種形式是多水平的一致序列形式，可以直觀的看到那些位置的保守堿基是哪些。 3位置P值的排序 這一部分列出了這個motif在序列上的匹配位置，并根據(jù)它的匹配的P值大小，按P值從小到大的順序排列序列，也就是根據(jù)匹配的可靠程度，排列序列，反映序列和motif的關(guān)系。這里的P值可以理解為假陽性率，P值越小，結(jié)果越可靠。-Motif 1 sites sorted by position p-value-Sequence name Start P-value Site- - - -BGK03149 547 1.21e-08 CTAGCGAGCG AGAGAGAGGGAAG GGTTGCGACTBGK04759 510 3.58e-08 GGATACAGGT AGAGAGAGGTGAG AAGGCAGTGGBGK04166 68 9.35e-07 CAATAGAGTT ATAGAGAGATAAG AGAAGAGGTABGK02457 558 1.40e-06 TTTCGTGAGC AGGAAGAGGGGGG GGGGGGGGGGBGK03537 102 1.62e-06 ACACCAAGTG ACCGAGAGATAAG GCTCATACAGBGK01655 278 2.51e-06 GACTAGACGG AGCAAGATGGAAG TACACAGTCABGK02838 315 4.48e-06 TTAAATAAAA AAGGAGAGAG