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

[TOC]

軟件準(zhǔn)備Preparation

Basics

Bash

Ruby 2; 版本> 2

Python 3

R

安裝包Libraries

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

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

軟件Programs

MAFFT: 多序列比對軟件；MAFFT

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

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

PAUP: 系統(tǒng)發(fā)育分析軟件MP法；PAUP*

FigTree系統(tǒng)樹的可視化軟件；FigTree

BEAST2分子進(jìn)化；分子時(shí)間計(jì)算經(jīng)典軟件；https://www.beast2.org

TracerBayesian貝葉斯樹的后驗(yàn)計(jì)算；Tracer

Blast+ 經(jīng)典序列比對blast軟件；Blast

ASTRAL 根據(jù)基因樹來推斷物種樹的軟件；ASTRAL

Abyss 同組裝相關(guān)軟件；Abyss

BWA 短序列比對軟件；BWA

bcftools vcf文件操作軟件；bcftools

Dsuite 計(jì)算Introgression的D-statistics數(shù)值軟件；Dsuite

Relate 大規(guī)模樣本數(shù)據(jù)計(jì)算譜系歷史的軟件；較新的軟件；Relate

多序列比對

Outline

多序列比對是系統(tǒng)發(fā)育分析的基礎(chǔ)；其目的是為了確定序列之間同源的區(qū)段(homologous nucleotides可比較的序列)；本節(jié)為了解決以下問題：

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

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

利用公共數(shù)據(jù)庫NCBI GeneBank補(bǔ)全已有的數(shù)據(jù)集合；

Datasets

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

Softwares

利用到的軟件包括：

MAFFT: 多序列比對軟件；MAFFT

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

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

利用MAFFT比對，AliView可視化

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

## 1

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

## 2 Gap penalty

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

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

利用BMGE去除低質(zhì)量比對區(qū)域

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

## 查看幫助

java -jar BMGE.jar -?

## 運(yùn)行

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

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

image

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

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

Preparation

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

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

利用PAUP模型篩選

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

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

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

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

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

結(jié)果顯示最佳AIC模型是GTR模型；GTR模型是系統(tǒng)發(fā)育分析中最常用的model。

pic

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

最大似然法建樹ML phylogenetic Inference

Preparation

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

利用IQTREE做ML系統(tǒng)發(fā)育分析

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

iqtree --help

## iqtree 默認(rèn)最簡參數(shù)；

iqtree -s 16s_filtered.nex

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

根據(jù)figtree顯示：

pic

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

選擇我們預(yù)設(shè)已知的某一支系作為外類群，reroot進(jìn)行設(shè)置。隨后即可展示各個類群之間的系統(tǒng)發(fā)育關(guān)系。

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

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

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

## -B 為ultrafast的bootstrap值，--prefix為預(yù)設(shè)輸出名字

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

結(jié)果利用figtree做展示，根據(jù)Note label下的diasplay Label即會展示每個節(jié)點(diǎn)的bootstrap值；

pic1

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

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

由于前面做模型選擇時(shí)我們發(fā)現(xiàn)，對于cds的codon(1,2,3)序列來說，不同位置會用不同的替換模型，在IQTREE中也可以同樣根據(jù)此進(jìn)行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樹之間的差異

通過不同的基因我們會獲得不同的系統(tǒng)發(fā)育樹。而為檢查兩個樹之間的差異以及評估整體的差別，我們可以利用Robinson-Foulds distance來鑒定樹之間拓?fù)浣Y(jié)構(gòu)差異。在IQTREE中的-t和--tree-dist2

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

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

串聯(lián)比對(concatenated alignments)下的系統(tǒng)推斷

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

#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;

## 根據(jù)組合方式進(jìn)行統(tǒng)一iqtree計(jì)算 -p

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

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

參考資料

總結(jié)

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

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