Github/mmatschiner的phylogenetic & phylogenomic學習教程記錄【一】多序列比對；核算替換模型的選擇；最大似然法建樹的學習

[TOC]

軟件準備Preparation

Basics

Bash

Ruby 2; 版本> 2

Python 3

R

安裝包Libraries

Python3的安裝包：包括有numpy, scipy msprime

R安裝包：包括有ape, coda，用于introgression分析。

軟件Programs

MAFFT: 多序列比對軟件；MAFFT

AliView： 多序列比對的可視化軟件；AliView

BMGE： 多序列比對后保守序列區段鑒別軟件，polish軟件；類似熟悉的Gblocks。BMGE

PAUP: 系統發育分析軟件MP法；PAUP*

FigTree系統樹的可視化軟件；FigTree

BEAST2分子進化；分子時間計算經典軟件；https://www.beast2.org

TracerBayesian貝葉斯樹的后驗計算；Tracer

Blast+ 經典序列比對blast軟件；Blast

ASTRAL 根據基因樹來推斷物種樹的軟件；ASTRAL

Abyss 同組裝相關軟件；Abyss

BWA 短序列比對軟件；BWA

bcftools vcf文件操作軟件；bcftools

Dsuite 計算Introgression的D-statistics數值軟件；Dsuite

Relate 大規模樣本數據計算譜系歷史的軟件；較新的軟件；Relate

多序列比對

Outline

多序列比對是系統發育分析的基礎；其目的是為了確定序列之間同源的區段(homologous nucleotides可比較的序列)；本節為了解決以下問題：

利用MAFFT軟件做多序列比對；

鑒定并排除可能錯誤比對的序列區段；

利用公共數據庫NCBI GeneBank補全已有的數據集合；

Datasets

數據來源于Matschiner et al. 2017。此為研究cichlid(慈鯛科魚種)在各大洋的歷史演化及分布。物種包括41species，分為不同地理的groups(Non-cichlid, Neotropical cichlids, African; Indian; Malagasy)。此處利用到兩個基因：線粒體mitocondrial 16S, 核基因RAG1

Softwares

利用到的軟件包括：

MAFFT: 多序列比對軟件；MAFFT

AliView： 多序列比對的可視化軟件；AliView

BMGE： 多序列比對后保守序列區段鑒別軟件，polish軟件；類似熟悉的Gblocks。BMGE

利用MAFFT比對，AliView可視化

利用MAFFT網絡版或本地版做序列比對；其中比對策略包括有不同全局比對和局部比對(Global and Local alignment)的算法策略，對于大多數人來說采用默認的--auto 軟件即會自動匹配最適算法。

## 1

mafft --auto 16s.fasta >16s_aln.fasta

## 2 Gap penalty

mafft --auto --op 2 16s.fasta >16s_op2_aln.fasta

利用AliView直接做序列比對的可視化

利用BMGE去除低質量比對區域

多序列比對軟件對于一些區域比對質量較差，可能會出現比對錯誤的情況，對后續系統發育分析會產生影響，因此需要將其剔除。我們利用BMGE軟件去除這些低質量比對的區域。BMGE(block mapping and gathering with entropy)由JAVA所寫。

## 查看幫助

java -jar BMGE.jar -?

## 運行

java -jar BMGE.jar -i 16s_aln.fasta -t DNA -of 16s_filtered.fasta -oh 16s_filtered.html

程序before和after會告訴排除了多少問題比對的序列。后續同樣的對核基因rag1也可以做相同的操作；

image

核酸替換模型的選擇Substitution Model Selection

當完成多序列比對之后，進行似然法likelihood的系統發育分析之前，隨后要進行的就是核酸替換模型的選擇，包括有Jukes-Cantor(JC)模型，HKY模型，GTR模型等等。通常是根據核酸替換模型軟件計算得出最佳的AIC值(Akaike Information Criterion)的模型，再根據此模型進行后續的系統發育分析。

Preparation

此篇教程是利用PAUP軟件進行的模型篩選。http://phylosolutions.com/paup-test/。開發于1980s，是系統發育分析最古老軟件之一。雖然目前很多功能已被其它更快更新的軟件所替代，但里面目前仍包括較為全面的likelihood-based系統發育分析的功能；目前在用的為版本4.0；

PAUP接受的格式為nex格式文件；

利用PAUP模型篩選

下載GUI圖形化版本的PAUP軟件，打開文件16s_filtered.nex

PAUP計算核酸替換模型需要首先利用NJ法快速建樹；

根據"Automated Model Selection"選項進行模型篩選；其中模型中包含有較多的術語：通常篩選標準包括有AIC(Akaike information criterion); AICc(Akaike information criterion corrected for small sample sizes); BIC(Bayesian information criterion); DT(decision-theoretic criterion)。通常是根據AIC值最小的Model選為最佳Model。

Model for among-site rate variation中G(GAMMA)分布; I(Invariable sites)

PAUP結果中-lnL(-log likelihood); K(the number of free parameters in the model), 代表branch length

結果顯示最佳AIC模型是GTR模型；GTR模型是系統發育分析中最常用的model。

pic

已有的核酸替換模型的篩選都是認為一個模型適用于所有位點，所有的位點進化速率都是相同的。但實際上一些位點的進化速率是不同的，例如cds序列中的third-codon位點是相較于第一第二位點速率是更快的。根據"Automated Partioning"結果顯示

最大似然法建樹ML phylogenetic Inference

Preparation

利用熱門的軟件IQTREE做系統發育分析，Figtree做基本的可視化展示。Datasets利用16s_filtered.nex and rag1_filtered.nex

利用IQTREE做ML系統發育分析

IQTREE軟件可以直接快捷進行系統發育分析：程序會自動計算使用的threads線程數；每條序列內的gap數目；選擇的最佳模型(根據BIC模型所選擇)；以及最大似然樹ML treefile(Newick格式文件)

iqtree --help

## iqtree 默認最簡參數；

iqtree -s 16s_filtered.nex

結果會生成Newick格式的treefile；((Ambassispcxxxx:0.04,Synbranmarmora:0.24):0.02,Mugilxxcephalu:0.13)；其中

根據figtree顯示：

pic

數字代表枝長branch(表示一定的遺傳距離)，括號內的兩個species為姊妹支。

選擇我們預設已知的某一支系作為外類群，reroot進行設置。隨后即可展示各個類群之間的系統發育關系。

利用bootstrap法獲取ML的各個node支持率

根據默認參數的-s所得到的系統發育樹，我們尚無法確定各個支系的支持率程度，可靠度(reliability)，因此我們需要借用bootstrap法獲取ml樹的支持率；

根據--help 我們利用"ULTRAFAST BOOTSTRAP/JACKKNIFE"的參數-B最少以1000次bootstrap做計算

## -B 為ultrafast的bootstrap值，--prefix為預設輸出名字

iqtree -s 16s_filtered.nex -B 1000 --prefix 16s_filtered.bs.nex

結果利用figtree做展示，根據Note label下的diasplay Label即會展示每個節點的bootstrap值；

pic1

對于小數據而言，不同次的bootstrap運算，以及設置不同的model都會對支持率有一定的影響(slight difference)。而一般對于iqtree的ultra bootstrap值來說，>90才被認為是可靠的節點。一些節點的支持率較低，多是因為序列數不多，信息位點數不多，不能很好的區分開序列之間的差別。

對cds核基因分段(Partitioned)的ML系統推斷

由于前面做模型選擇時我們發現，對于cds的codon(1,2,3)序列來說，不同位置會用不同的替換模型，在IQTREE中也可以同樣根據此進行Partitioned inference；

#

#NEXUS

BEGIN SETS;

CHARSET codon1 = 1-1368\3;

CHARSET codon2 = 2-1368\3;

CHARSET codon3 = 3-1368\3;

END;

### 將以上輸出到partitions.txt文件中

## run iqtree

iqtree -s rag1_filtered.nex -p partitions.txt -B 1000 --prefix rag1_filtered.bs.nex

比較不同ML樹之間的差異

通過不同的基因我們會獲得不同的系統發育樹。而為檢查兩個樹之間的差異以及評估整體的差別，我們可以利用Robinson-Foulds distance來鑒定樹之間拓撲結構差異。在IQTREE中的-t和--tree-dist2

iqtree -t 16s_filtered.bs.nex.treefile --tree-dist2 rag1_filtered.bs.nex.treefile

同樣也可以計算不同樹文件中的平均bootstrap值判斷支持率。

串聯比對(concatenated alignments)下的系統推斷

一個基本前提是，序列位點越多，最后獲得各個支系的支持率也就越大，所以我們可以將多個基因串聯起來統一進行系統發育推斷。當然另一個前提是不同的基因會隨著物種進行同樣的進化歷程，這個前提在對一些近緣物種推斷時是會存在問題的。

#NEXUS

BEGIN SETS;

CHARSET 16S = 16s_filtered.nex: *;

CHARSET rag1_codon1 = rag1_filtered.nex: 1-1368\3;

CHARSET rag1_codon2 = rag1_filtered.nex: 2-1368\3;

CHARSET rag1_codon3 = rag1_filtered.nex: 3-1368\3;

END;

## 根據組合方式進行統一iqtree計算 -p

iqtree -p partitions.txt -B 1000 --prefix concatenated.bs.nex

結果確實會比單基因建樹的支持率效果更好。多基因串聯建樹也是后期由基因樹推斷物種樹的一個重要的方法。

參考資料

總結

以上是生活随笔為你收集整理的figtree需要在JAVA下运行吗_Phylogenomic_Tutorial || ML_Tree inference的全部內容，希望文章能夠幫你解決所遇到的問題。

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