一、引言

??? GENBANK是目前最大而權威的分子序列數據庫，調用其中數據可以進行分子系統發育樹的構建。

1、序列數據獲取（以皿蛛系統發育樹為例）

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

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

????? 以上說得元數據在更新時一定要小心，以免造成破化。建議采用自動手動結合的方式進行更新。

?實戰：Linyphiidae分析系統發育樹構建。

???????????? 元數據：taxaid——DB中獲得

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

????????????? 序列數據：得到GI列表后，通過R語言腳本read_genbank.r獲得序列數據，并轉換為fasta格式文件存儲待用。

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

??? 例子：為了運算速度快，只選擇4種皿蛛進行分子系統發育樹構建。分別是187178Bolyphantes alticeps，187180Neriene variabilis，187181Neriene radiate，233285Pimoa sp. X131，Bolyphantes屬于皿蛛亞科，兩種蓋蛛屬于M-EClade，外群為Pimoa。

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

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

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

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

????????? （1）首先是定義校準點。在taxa標簽下，將具有最近共同祖先MRCAs的taxa組成一組，表示可以設定校正點。在本例中，除Pimoa外，所有皿蛛為MRCA，強制單系，以Pimoa為外群。在實際應用中，可設置多個tMRCA，但taxa之間不能有重復，否則Beast時會報錯。

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

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

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

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

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

???????????? （4）Operators。全部默認設置

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

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

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

?

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

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

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

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

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

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

??? C. 運行Tracer分析結果。用Tracer打開log文件分析。這個過程可以看出跑樹結果的好壞。

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

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

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

? b. 和議完成后，可以通過PSRF參數評價樹的好壞。?

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

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

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

代碼中的2048為jvm可調用內存的限值，運算量越大，這個值應調的越高。

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

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



總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。