Phylogenomic_Tutorial || ML_Tree inference

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*
IQTREE 最大似然法ML法建树；推荐2.06版本以上； https://github.com/Cibiv/IQ-TREE/releases
FigTree系统树的可视化软件；FigTree
BEAST2分子进化；分子时间计算经典软件；https://www.beast2.org
TracerBayesian贝叶斯树的后验计算；Tracer
Blast+ 经典序列比对blast软件；Blast
PAML 系统发育推断软件；主要利用里面的codeml子程序计算dN/dS；PAML (Phylogenetic Analysis by Maximum Likelihood)
ASTRAL 根据基因树来推断物种树的软件；ASTRAL
aTRAM2 组装相关软件；aTRAM github repository
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

结果确实会比单基因建树的支持率效果更好。多基因串联建树也是后期由基因树推断物种树的一个重要的方法。

参考资料

Substitution Model Selection
Maximum-Likelihood Phylogenetic Inference