一、引言

??? GENBANK是目前最大而權(quán)威的分子序列數(shù)據(jù)庫，調(diào)用其中數(shù)據(jù)可以進(jìn)行分子系統(tǒng)發(fā)育樹的構(gòu)建。

1、序列數(shù)據(jù)獲取（以皿蛛系統(tǒng)發(fā)育樹為例）

???? 在GenBank中，每一個(gè)物種或階元都有一個(gè)taxid，他是taxa的ID。而且taxa之間存在父子關(guān)系。我們的研究對象是蜘蛛目（Aranaea），其taxaid為6893，其父級階元是蛛形綱（Arachnida），taxaid為6854。按照遞歸查詢的原則，只要有自6893以下所有taxa的父子關(guān)系對照表，就能查詢到目前在GenBank中記錄的所有蜘蛛的名錄。事實(shí)上，這個(gè)想法已經(jīng)可以實(shí)現(xiàn)了！NCBI有一個(gè)public的ftp（ftp.ncbi.nlm.nih.gov），從中的/pub/taxonomy/taxdmp.zip壓縮包中可以下載到相應(yīng)的信息。taxadmp.zip含有9個(gè)文件，其中比較重要的是其中的names.dmp和nodes.dmp。通過查看readme.txt文檔，得知其中names.dmp為genbank中所有taxaid的基本信息，nodes.dmp為taxaid有關(guān)父子關(guān)系的信息。將兩個(gè)文件通過NCBI_TAXDMP_NAMES.nopi和NCBI_TAXDMP_NODES.nopi導(dǎo)入文件導(dǎo)入到DB中，再遞歸查詢并篩選出屬于6893以下階元的文件，便可得到GenBank中有記錄的所有spider的名錄。在我的DB中，導(dǎo)入的這兩個(gè)表稱為NCBI_TAXDMP_NAMES和NCBI_TAXDMP_NODES。為了和DB中其他表格協(xié)同，特新建了MV_SIG_TAXDMP_NAMES和MV_SIG_TAXDMP_NAMES_ID兩個(gè)物化視圖。

???????? 在Gen其每一條序列數(shù)據(jù)都有唯一的GI號，它就是序列的身份證。在得知TaxaID后，建立GI、Gene和TaxaID的對應(yīng)關(guān)系就十分重要。在DB中，這個(gè)表稱為SIG_NUCL，它有4個(gè)字段，分別是tax_id,GI,Gene和Flag。其中，Gene的外鍵是SIG_Gene表中值。由于同一個(gè)taxa的某一個(gè)基因可能有多個(gè)序列，且長短不一，最終要選哪個(gè)必須由研究者手動(dòng)來定，所以表中增設(shè)了Flag開關(guān)，1代表選取某序列，0代表暫存但在計(jì)算中不使用。

????? 以上說得元數(shù)據(jù)在更新時(shí)一定要小心，以免造成破化。建議采用自動(dòng)手動(dòng)結(jié)合的方式進(jìn)行更新。

?實(shí)戰(zhàn)：Linyphiidae分析系統(tǒng)發(fā)育樹構(gòu)建。

???????????? 元數(shù)據(jù)：taxaid——DB中獲得

??????????????????????????? GI ——NCBI提供的Entrez在線檢索器用腳本代碼生成結(jié)果xml文件，再從xml文件中獲得某一個(gè)taxaid的GI列表。見entrez.txt

????????????? 序列數(shù)據(jù)：得到GI列表后，通過R語言腳本read_genbank.r獲得序列數(shù)據(jù)，并轉(zhuǎn)換為fasta格式文件存儲(chǔ)待用。

????????????? 其他工具：序列對齊工具ClustalX2，跑樹工具BEAST工具包。

??? 例子：為了運(yùn)算速度快，只選擇4種皿蛛進(jìn)行分子系統(tǒng)發(fā)育樹構(gòu)建。分別是187178Bolyphantes alticeps，187180Neriene variabilis，187181Neriene radiate，233285Pimoa sp. X131，Bolyphantes屬于皿蛛亞科，兩種蓋蛛屬于M-EClade，外群為Pimoa。

???? 1. 首先，用R腳本生成待所有taxa序列的fasta文件，再導(dǎo)入ClustalX2進(jìn)行自動(dòng)對齊，再用NOTEPAD++進(jìn)行手動(dòng)對齊，之后導(dǎo)出為aln文件。

???? 2. 對齊之后，為了滿足BEAST對于格式的要求，再M(fèi)esquite中將aln文件轉(zhuǎn)換為nex格式。方法是：（1）在Mesquite中打開aln文件，在彈出的對話框中選擇Clustal（DNA/RNA）,將nex保存在aln同一目錄。但是，由于BEAST軟件并不能識別這種復(fù)雜格式的nex，還需要在Mesquite中將其導(dǎo)出為簡單格式的nex文件。做法是打開gi_seq_aln.nex文件，（忽略提示錯(cuò)誤），選擇export->Simplified_Nexus導(dǎo)入。導(dǎo)出后的nex只包含文件頭，taxa標(biāo)簽和對齊后的序列信息，之前nex的參數(shù)文件都不見了。

???? 3. 使用BEAST跑樹。?????????????????????

??????? A. 首先，打開BEAUti，生成跑BEAST所需的xml文件。具體做法參考BEAST使用說明。

????????? （1）首先是定義校準(zhǔn)點(diǎn)。在taxa標(biāo)簽下，將具有最近共同祖先MRCAs的taxa組成一組，表示可以設(shè)定校正點(diǎn)。在本例中，除Pimoa外，所有皿蛛為MRCA，強(qiáng)制單系，以Pimoa為外群。在實(shí)際應(yīng)用中，可設(shè)置多個(gè)tMRCA，但taxa之間不能有重復(fù)，否則Beast時(shí)會(huì)報(bào)錯(cuò)。

?????????? （2）設(shè)置分子進(jìn)化模型Site Models。Substitution Model：GTR。Base frequeencies: Empirical. Site Heterogeneity Model選為Gamma

??????????? (3) 設(shè)置分子鐘模型Clock Models。 將Molecular Clock Model選為“Relaxed Clock:Uncorrelated Lognormal”, estimate不選擇，fix mean rate of molecular clock model不選擇。

??????????? （4）設(shè)置Trees。選擇Tree Prior為Speciation:Birth-Death Process,starting tree處選擇User defined.此處可以import一顆ML樹作為BEAST的起始樹。具體做法是在file菜單下點(diǎn)import data，將nex格式的ML樹導(dǎo)入。需要注意的是，這一nex樹不能用mesquite直接生成的simplified nex，格式要更加簡化。translate標(biāo)記要取消，而且腳本以begin tree開始，end結(jié)束，以newick格式記錄樹形及支長。在引入自定義的Starting tree后，tMRCA的設(shè)定便有了限制，tMRCA的taxa分組不能和Starting tree拓?fù)浣Y(jié)構(gòu)沖突，否則在BEAST的時(shí)候會(huì)報(bào)錯(cuò)。

???????????????????????

?????????????(5)?? 設(shè)置Prior。在tmrca（Linyphiidae）中選擇Lognormal，參數(shù)為0，2，125（參考Dimitrov，2012）。

???????????? （4）Operators。全部默認(rèn)設(shè)置

???????????? （5）設(shè)置MCMC。鏈長設(shè)為1000000，Echo設(shè)為10000，Log設(shè)為200（這些參數(shù)需要根據(jù)研究需要自定義）。

???????????? （6）生成XML文檔（忽略）。

????? B.運(yùn)行BEAST。在BEAST中打開xml文件，運(yùn)行后生成后綴為trees和log的文件。

?

????????????????? 由于BEAST跑樹計(jì)算量非常大，而且要持續(xù)數(shù)天，因此不建議在PC上運(yùn)行，最好在服務(wù)器上運(yùn)算。從運(yùn)算效率上考慮，建議使用Linux平臺的操作系統(tǒng)，我選擇的是CentOS。BEAST在CentOS系統(tǒng)中的安裝步驟：

????????????????? 1.從BEAST官網(wǎng)下載運(yùn)行在Linux/Unix系統(tǒng)下的tgz包，之后在PC上用7zip解壓，將解壓好的文件夾拷到root目錄下。

?????????????????? 2.變更best文件夾的權(quán)限。最簡單的辦法就是輸入如下命令，使所有文件夾和文件對所有用戶都有最大權(quán)限。

??????????????????????????????? chmod -R 777 [beast目錄]

?????????????????? 3.由于BEAST是基于Java虛擬機(jī)運(yùn)行的程序，因此在CentOS中首先要安裝并運(yùn)行Java。由于BEAST在運(yùn)行中需要CentOS的圖形界面，因此要在gnome桌面模式下運(yùn)行。

????????????????? 4. 在/beast/bin目錄下找到beast程序，運(yùn)行BEAST。彈出的對話框打開編輯好的xml文件，就會(huì)自動(dòng)運(yùn)行beast，結(jié)果trees和log文件存儲(chǔ)在bin目錄下。需要說明一點(diǎn)，如果在windowsPC下編輯生成的xml文件，在CentOS平臺下不能直接使用，由于EOF在UNIX和DOS/WINDOWS下編碼不同，直接使用會(huì)報(bào)錯(cuò)。解決方法是在NOTEPAD++的EDIT菜單下將文件的EOF模式設(shè)置為UNIX格式，保存后將其上傳至服務(wù)器，再用BEAST打開就OK了。

??? C. 運(yùn)行Tracer分析結(jié)果。用Tracer打開log文件分析。這個(gè)過程可以看出跑樹結(jié)果的好壞。

?????? (1)Tracer顯示了BEAST后得到的各種參數(shù)，包括likelihood,meanRate(平均進(jìn)化速率)，coefficientOfVariation(每一分支分子進(jìn)化的變異系數(shù)，表示分子變化的速率在不同分支上是否等速),root和每一個(gè)tMRCA分化時(shí)間的估計(jì)值等等。同時(shí)選中root.height和tMRCA，可以通過箱線圖和邊界密度圖顯示預(yù)設(shè)的化石校正點(diǎn)估計(jì)值是否一致。

? （2）Tracer作為一個(gè)通用工具，不僅可以分析BEAST跑出的bayers樹，還可以分析MrBayers等跑出來的樹。對于MrBayers跑出來的樹，可以做如下分析：

???????? a. 確定Burnin的值，在tracer中加載兩輪.p文件，分析trace圖，改變Burnin的值，直到trace圖顯示lnl值能夠圍繞一個(gè)值小幅波動(dòng)為止。一般來說，bayers分析burnin值不超過15%。

? b. 和議完成后，可以通過PSRF參數(shù)評價(jià)樹的好壞。?

sump burnin=250（在此為1000個(gè)樣品，即任何相當(dāng)于你取樣的25％的值），參數(shù)總結(jié)summarize the parameter，程序會(huì)輸出一個(gè)關(guān)于樣品（sample）的替代模型參數(shù)的總結(jié)表，包括mean，mode和95 % credibility interval of each parameter，要保證所有參數(shù)PSRF（the potential scale reduction factor）的值接近1.0，如果不接近，分析時(shí)間要延長.RSRF值在*.vstat文件中。 c. 和議完成的Bayers樹，每個(gè)節(jié)點(diǎn)還會(huì)有probs值，該值表示某一結(jié)點(diǎn)的支持度。

??? D. 運(yùn)行TreeAnnotator得到合議樹。用TreeAnnotator打開trees文件，設(shè)定參數(shù)得到合議樹。特別要指出的是，TreeAnnotator中Output File文件默認(rèn)是讓選擇，其實(shí)直接取一個(gè)輸出的文件名點(diǎn)ok即可。需要注意的是，由于jvm有內(nèi)存限制，對于超過10000棵樹需要和議時(shí)，不能直接運(yùn)行treeannotator，而是進(jìn)入beast的lib目錄，運(yùn)行以下代碼啟動(dòng)treeannotator。

java -Xmx2048m -Xms2048m -classpath beast.jar dr.app.tools.TreeAnnotator

代碼中的2048為jvm可調(diào)用內(nèi)存的限值，運(yùn)算量越大，這個(gè)值應(yīng)調(diào)的越高。

? ? E.最后用FIGTREE打開tree文件，分析和議得到的Dated Phylogeny。 為了使結(jié)果更為直觀，可以在Figtree中顯示tMRCA中node的值及95%HPD

??????????????



總結(jié)

以上是生活随笔為你收集整理的使用GENBANK数据进行分子系统发育树的构建的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。