ixxmu
commented
1 year ago 手把手教你绘制进化树（理论与实践） by 水瓶之翼
********目录*******
目的
通过系统发育进化分析，确定目标序列与其它同源序列的演化关系、遗传距离远近，或者对细菌/病毒序列进行基因分型等。（本文分享的主要是基于编码区DNA的系统发育分析）
方法与过程
下载 MEGA 11.0 软件
使用 blast 确定序列同源性，下载目标序列
使用 MEGA 11.0 进行多序列比对
使用 MEGA 11.0 构建进化树
将进化树导出，使用R包 ggtree 进行可视化
1. 软件下载
MEGA(Molecular Evolutionary Genetics Analysis)是一个开源软件，我们可以直接下载最新版用于windows系统的图形化界面(GUI)版 MEGA 11.0。
MEGA 11 对大数据集的分析进行了优化，界面也更加美观，其它功能与 MEGA 7 没什么区别，当然你也可以两个版本都下载比较一下，windows上可以同时安装不同版本的MEGA。下载网址为：
https://www.megasoftware.net/
下载时需要匿名填一些信息，下载完直接按照提示安装即可。
2. 目标序列获取
一般来说，我们自己至少会有1条感兴趣的目标序列，其格式应该是 FASTA 格式。我这里举了个例子，是一个细菌序列的某个基因，第一行为>开头，后面为该序列的注释信息，第二行之后的所有内容是该基因的核酸序列。
如果你已经有自己准备好的序列就不需要进行目标序列获取这一步啦。
>CP000702.1:1005082-1006524 Thermotoga petrophila RKU-1, complete genome
ATGCCATCTGTGAAGATCGGTATCATCGGTGCGGGGAGCGCGGTGTTTTCTCTGAGGCTTGTGAGTGATC
TTTGCAAAACGCCGGGACTCTCTGGCAGCACGGTCACCCTCATGGATATCGACGAAGAAAGACTCGACGC
TGTTCTGACCATCGCGAAAAAATACGTTGAAGAAGTGGGAGCGGATCTGAAATTCGAAAAAACCATGAAT
TTAGATGACGTCATCATCGACGCGGATTTTGTGATAAACACAGCGATGGTGGGTGGCCATACCTACTTGG
AGAAAGTCAGACAGATCAGTGAGAAATACGGCTACTACAGAGGAATAGACGCTCAGGAGTTTAACATGGT
CTCCGACTACTACACCTTCTCCAACTACAACCAGCTCAAGTACTTCGTTGAAATAGCAAGGAAGATAGAG
AAGCTCTCCCCAAAAGCCTGGTACTTGCAGGCAGCGAATCCCGTTTTCGAAGGAACAACCCTTGTGACAA
GAACGGTTCCCATAAAGGCAGTGGGATTCTGCCATGGACACTACGGCGTGATGGAGATCGTAGAGAAACT
GGGGCTGGAAGAAGAAAAAGTAGATTGGCAGGTCGCAGGAGTGAACCACGGTATCTGGCTGAATAGGTTC
AGATACAACGGGGAGAACGCGTATCCCCTCCTTGACAGGTGGATCGAGGAAAAATCAAAAGATTGGAAAC
CAGAGAACCCCTTCAACGATCAGCTCTCTCCCGCTGCGATAGATATGTACAGATTCTACGGTGTGATGCC
CATCGGTGACACCGTGAGAAACTCTTCGTGGAGGTACCACAGGGATCTTGAGACCAAGAAGAAATGGTAC
GGTGAACCCTGGGGAGGAGCAGATTCTGAAATAGGCTGGAAATGGTACCAGGACACACTTGGAAAGGTAA
CGGAGATCACAAAGAAGGTGGCAAAGTTCATCAAAGAAAATCCGTCCGCGAGGCTCTCCGACCTTGGAAG
TGTTCTGGGAAAAGACCTCTCAGAAAAGCAGTTTGTGCTCGAAGTAGAGAAAATCCTCGATCCAGAAAGA
AAGAGTGGAGAGCAGCACATCCCATTCATCGATGCGCTGTTGAACGATAACAAGGCAAGATTCGTGGTGA
ACATACCAAATAGGGGTATCATCCACGGAATAGACGATGACGTGGTCGTTGAAATCCCAGCCCTTGTGGA
CAAGAACGGAATCCATCCCGAGAAGATCGAACCACCGCTTCCAGATCGCGTGGTCAAGTACTACCTGAGA
CCAAGGATCATGAGAATGGAAATGGCTCTGGAGGCGTTCCTCACGGGGGACATAAGGATCATAAAAGAAC
TTCTCTACAGAGATCCAAGAACGAAGAACGACGAACAGGTAGAAAAGGTGATCGAGGAAATCCTCGCACT
CTCGGAAAACGAAGAGATGCGGAAACATTATCTGAAGAGATGA
这里需要注意一个概念，我们用于进化分析的所有序列理论上应该都是同源的，所以并不是说随便下载几条序列都可以进行进化分析的。有些序列比对后仅仅能证明它们相似而不一定能证明它们同源，而 NCBI的blast工具就可以进行显著性检验，帮助我们判断序列同源的可靠性。
需要区分同源与相似的概念。
如果两个序列享有一个共同的进化上的祖先,则这两个序列是同源(homologous)的。对这个定义需要注意的一点是，同源是个定性的概念，没有“度”的差异。对两个序列，它们或者同源或者不同源，不能说它们70%同源或80%同源。
与同源相关但不同的两个概念是相似(similarity)和相同(identity)，它们都是定量的概念，基于对序列中字符的精确比较，既可以说两个核酸序列高度相似，也可以说它们 70% 的碱基相同。
2.1 blast获取同源序列
进入ncbi在线网站 -> 点击右侧的blast -> 选择 Nucleotide BLAST -> 将前面的序列复制进去，其它参数默认即可 -> 然后点击下方的 BLAST，等待比对结果。
比对结果的解释 ：
Description 列，是对当前 subject 序列的简短介绍，它还有一个链接，鼠标左键点击它即可通往该序列比对的结果页面；
Scientific Name 列，不是我们关注的重点；
Max Score 列，代表比对到 subject 序列的片段的最高评分（因为 blast 比对时是将待比对序列拆分成一系列片段去与subject序列进行比对的）；Total Score 列，同理，即所有比对片段的评分之和。这两个评分越高，代表两个序列越相近；
Query Cover 列，代表与 query 序列的全长相比，有多少比例的 query 序列与 subject 序列能比对上（即序列比对时能覆盖到的比例）；
E value 列，即期望值。是 blast 软件进行统计显著性计算得到的值，E值越接近0，代表该 subject 序列越有可能与 query 序列同源；
Per.Ident 列，即两个序列能比对上的部分完全相同的碱基的比例；
Acc.Len 列，即 subject 序列的原始长度；Accession 列表明了 subject 序列的在 genbank 上的登录号。
序列同源性判断与下载：
我们点击 Description 列的第一条序列，可以通往比对结果的链接。
先看下图中红色方块中的内容，Score/Expect/Identities 与前面的一致（identities可以检查比对的序列是否有失配，包括插入缺失以及碱基的替换）；
Gaps代表空格数（序列比对时如果有插入/缺失，blast会通过插入空格来使序列对齐）；
Strand代表 query 序列与 subject 序列的方向，其中
Plus/Plus 代表 subject 序列与 query序列是同向的，Plus/Minus则代表是反向的。
绿色方块中显示了具体的比对情况，有gap的位置会通过
-表示。接着我们要判读该 subject 序列是不是可以下载下来作为我们进化分析的同源序列。事实上并没有一个严格的标准来界定序列是否同源，有文献把 Query Cover>=60%，E value <1e-3 的序列作为区分是否同源的依据，也有文献选择 E value <1e-5，所以需要我们根据自己的研究领域进行判断，如果符合我们的要求，就可以认为比对得到的这条 subject 序列与我们的 query 序列同源，就可以下载下来用于进化分析了。
在本例中(Query Cover=100%，Strand=Plus/Plus)，直接点击粉色方块中的 Download -> 选择 FASTA(aligned sequences)可以直接下载到完整的subject序列。
当 Query Cover < 100%，或者 Strand=Plus/Minus。比如我们比对到的下图中这条 subject序列，它的序列其实和我们输入的序列是反向互补的，如果直接点击 Download 其实下载的是反向互补的序列，如果我们没注意的话后面比对结果会不准确。
因此我们点击下面蓝色方块中的 Genbank 链接，会进入一个新的窗口（下图2），然后在右侧的 Customize view 中勾选 Show reverse complement，接着点击 Update View即可。
接着可以再手动确认一下更新的序列前几位是否与我们的 query 序列一致，然后即可下载这条 subject 序列，点击上方的 Send to 下载 FASTA 文件（下图3）
总结：
这样一条条序列检查是否同源，是否同向，然后下载是比较稳妥的方法，不会出错，但是序列多了就会比较浪费时间，我们也可以直接在比对页面直接选择多个序列批量下载 FASTA 文件(如下图所示)，但是下载后一定要注意有些序列与我们的 query 序列是反向互补的，需要将其反向互补后重新整合才能用于后面的序列比对（不过之后用MEGA也可以对序列进行反向互补操作）。
2.2 自己判断序列同源性
还有一种情况是我们自己知道下载的序列一定是同源的，比如下载某个病毒不同毒株的某个基因序列，就没有必要通过 blast 等软件检查序列是否同源。
比如我在 NCBI 的 Nucleotide 数据库中检索 PRRSV 病毒的 ORF5基因，如下图所示，然后下载其中一些序列用于进化分析，勾选想要的序列后还是和前面一样点击 Send to 下载 FASTA 文件即可。
2.3 序列检查
不管我们是如何获取到的序列，在将其用于比对之前一般应确保以下原则（1）尽量保证序列完整，既没有缺失也没有多余（2）用于进化分析的序列应同源（3）序列的方向应一致。
如果我们选择的是某个基因全长或者CDs序列(即编码序列)这种会比较好检查。获取到所有序列后一般将所有序列整合到一个FASTA文件中，接着推荐使用 Sublime 软件直接打开 FASTA 文件，可以方便地对序列进行批量检查。检查的方式可以根据我们的序列特征进行，比如下面几种方法：
检查
>的数量。因为1条序列由1个
>开头，>的数量应与我们实际需要比对的序列数量一致。比如我这次需要比对18条序列，在 Sublime 中按Ctrl+F进入搜索模式，然后搜索>，如下图所示，此时左下角会出现文档中的>的数量(绿色方框所示)，显示18个即说明序列数量正确。
检查序列方向。针对我感兴趣的序列，因为我要比对的这18条序列都是CDs序列，因此大部分的序列应该都是以 ATG 起始密码子起始（有些物种可能出现GTG），我们不能直接搜索以ATG开头的行，这样会与序列中间的以ATG开头的行混淆。再仔细看一下这些序列，发现都是
>所在的行以1个数字结尾，然后另起一行以ATG开头，我们通过精确匹配就能准确搜索到所有以ATG开始的序列。在 Sublime 中按
Ctrl+F进入搜索模式，如下图所示，点击下方最左侧的正则表达式搜索，搜索式为[1-9]\nATG，此时左下角显示18个即说明序列方向应该正确。（该搜索式中，[1-9]代表第1个字符匹配1-9之间的所有数字，\n代表换行符）
同样的，我用表达式 TAG\n\n搜索时，发现所有18条序列都是以TAG结尾（终止密码子一般是 TAA、TAG、TGA）,因此可以确认我的所有序列都是完整的，且方向也是对的。
批量修改注释行，一般从genbank上下载下来的序列
>所在的行不会像我的例子中这么干净，一般会有一串信息。因为之后构建进化树的时候软件一般会默认使用>后面的信息给树中的序列命名，所以我们可以提前用 Sublime 进行批量整理。在下面的例子中，如果我们想只保留
>后面的 genbank登录号，只需要搜索:，然后点击右侧的 Find All，此时光标会出现在所有匹配的:上，此时想批量删除:右侧的信息，只需要按Fn+Shift+右方向键选中，然后按退格键删除即可。
总结：
前面的序列检查方法并不唯一，针对我们自己的序列特征可以进行不同的检查（比如还可以检查文档的行数，如果1个基因序列一般都占10行的话，20条序列应该就是200行，Sublime中左侧显示了行数，可以防止下载成了别的基因）。总之，其主要目的都是为了确保所有序列的方向一致，所有序列同源且完整。
3. 使用 MEGA 进行多序列比对
3.1 导入待比对的 fasta 序列
依次点击：ALIGN -> Edit/Build Alignment -> Retrieve a sequence from a file -> OK -> 选择我们前面整理好的fasta文件 -> 出现 Alignment Explorer 窗口：
导入序列后如下图所示：（这时如果想更改某条序列的名称可以在该原名字处右键选中 Edit Sequence Name编辑）
3.2 多序列比对
MEGA提供了两种序列比对方法，Clustal W 与 MUSCLE ，有文献报道 MUSCLE 的准确性更高，且速度比Clustal W 更快，因此这里我们直接使用 MUSCLE 进行多序列比对演示。
按Ctrl + A选中所有序列 -> 点击 Alignment -> 选择 Align by MUSCLE -> 参数设置界面直接点击 OK 即可，一般不需要更改参数。（如果用Clustal W 进行核酸序列比对的话一般也不需要修改参数）
比对后结果：
3.3 保存比对结果
依次点击 Data -> Export Alignment -> MEGA Format，然后保存为.meg格式的文件即可（保存时会提示输入标题以及是否是编码序列），我这里命名为example.meg。另外我们这里也可以另外保存一份 FASTA格式的文件用于我们最后一步的可视化（我这里命名为example.fas）。
3.4 比对结果可靠性检验
首先，我们可以大致看一下比对之后有没有异常序列，如果发现比对后存在明显与其它序列不同的序列，可以用MEGA进行相应编辑或者考虑删除后重新比对。
其次可以检验比对结果的可靠性。
如果比对结果不可靠，那么构建的进化树也会不可靠。Thompson等人（1999）的一项研究比较了一些比对程序，表明当配对比对中氨基酸一致性(identity)的平均百分比太低时，多重比对的准确性低于它们能够产生可靠的系统发育树的水平。Thompson的研究表明，当氨基酸平均identity<20%时，则<50%的残基是被正确比对的。在20%-30% identity 之间的“过渡带”中，约80%的残基正确对齐，>30%的identity中，90%的残基正确比对。
因此我们可以通过计算序列之间的 p 距离(p-distance)来检验核酸或氨基酸比对的可靠性。比如在氨基酸比对中， p-distance 值为 0.6，代表氨基酸平均 identity为 1-0.6=40%，如前所述，氨基酸的 identity 一般需要 >30% 才比较可信(即 p-distance<0.7)，因此当 p-distance=0.6时，我们可以认为氨基酸的比对结果是比较可靠的。
对于核酸来说，检验方法类似，假如 p-distance 为 0.33，则 67% 的 identity 就可能不足以用于准确估计进化关系。p-distance<0.33 一般被认为结果比较可靠。
一般来说，假如核酸序列或者氨基酸序列的 p-distance 较高时，不应将多序列比对的结果用于构建进化树，而是应该删除其中不可靠的序列后重新比对，使得核酸序列或者氨基酸序列的 p-distance 满足 核酸：<0.33/氨基酸：<0.7。
在 MEGA 中依次点击：DISTANCE -> Compute Overall Mean Distance -> 选择前面导出的比对完的 .meg文件：
然后按照下图进行参数设置：
点击OK后，可以得到 1 个平均 p-distance， - 本例中是0.26<0.33，因此我们可以认为
- 本次比对结果比较可靠。
总结：
比对方法以及比对可靠性的检验方法有很多种，上面只是介绍了其中一种。比如 MAFFT、PRANK等比对方法就被认为准确性比较高，有文献报道 MAFFT准确性比Clustal W高、速度也更快，PRANK准确性被认为比MAFFT与Clustal W都高，但速度却相对较慢。大家在实际比对时可以尝试多种方法，并在准确性与速度方面进行一定的权衡。
4. 构建进化树
构建进化树的方法主要可以分为两类
基于距离的方法，基于所计算出的序列间遗传距离来构建进化树，快速，会产生单一的树。
邻接法：Neihbor Jioning(NJ) 非加权分组平均法(UPGMA)
基于字符的方法，基于描述遗传字符演化的数学模型构建进化树，速度较慢，会产生多棵树，然后根据其最优选择标准来找出最佳进化树。
最大简约法(Maximum Parsimony,MP) 最大似然法(Maximum Likelihood,ML) 贝叶斯马尔科夫链蒙特卡洛法(BMCMC)
根据《phylogenetic trees made easy》作者所述，MEGA是构建NJ树最好的程序，且其强烈反对使用UPGMA法构建系统发育树。
4.1 构建 NJ 树
邻接法：Neihbor Jioning(NJ)是使用得最广泛的其中一种方法。我们使用MEGA根据前面的比对结果构建NJ树。
邻接法基本介绍：
邻接法是一种基于最小进化原理的算法，基于距离矩阵构建树，它不检验所有可能的拓扑结构，能同时给出拓扑结构和分支长度。它的优点是重建的树相对准确，假设少，计算速度快，只得到一棵树。
NJ树适用于进化距离较小、信息位点较少的短序列，可用于大型数据集（速度很快）。其缺点主要表现在将序列上的所有位点等同对待，所分析序列的进化距离不能太大（对于距离过远的序列会出现长枝吸引现象）。
4.1.1 检验是否适合构建NJ树
假如 Jukes-Cantor(JC)距离 >1.0，则序列不适合构建NJ树，而应采用其它方法。
我们采用与 3.4 相同的步骤，在 MEGA 中依次点击：DISTANCE -> Compute Overall Mean Distance -> 选择步骤3.3多序列比对得到的.meg文件，然后在参数设置中将 Model/Method 改为 Jukes-Cantor 即可计算出 JC 距离。
在本例中计算出的 JC 距离为 0.34，因此我们认为适合构建 NJ 树。
4.1.2 NJ 树的估计-理论
为什么说是估计，类似于模型构建，因为事实上我们永远无法得到真实的物种演化的系统发育树，我们所谓的构建进化树只是在用各种算法尝试得到最接近真实进化树的一颗树。
我们估计进化树时一般会估计多次，首先是在没有 Bootstrap 的情况下先估计得到一颗树，初步查看是否有明显异常的序列被异常分类到某个子树中（关于什么是 Bootstrap 后面会介绍）。如果有这种情况应对序列进行检查，必要时还需要删除某些序列，然后重新进行多序列比对，再重新估计树。当我们认为树中没有异常序列时，再设置 Bootstrap 次数（一般为1000以上，至少100），估计最终的树。
对于最大似然法，这样估计树可以节约很多时间，因为最大似然法设置 Bootstrap 后速度会非常慢，对于 NJ 法，在序列不多的情况下（只有数百条），可以直接设置 Bootstrap，因为 NJ 法速度很快。
在 MEGA中依次点击：PHYLOGENY -> Construct/Test Neighbor Joining Tree -> 设置参数，我们可以先看一下参数设计界面的各个参数。
估计 NJ 树时各参数的意义与设置：
Test of Phylogeny：进化树的检验设置，我们可以根据前面说的设置为None，表明不检验（速度快），也可以设置为 Bootstrap method，然后在下面的 No.of Bootstrap Replications 设置检验次数，关于 Bootstrap 的介绍可见后面备注3。 Substitutions Type：本文只介绍比对核酸序列构建进化树，因此直接选择默认的 Nucleotide。 Model/Method：进化模型。我们这里选择 Maximum Composite Likelihood，算法原理在后面备注1中介绍。 Substitution to include：选择哪些替代类型用于计算，一般不需要改。transition 代表转换，嘌呤之间或者嘧啶之间，如 A<->G,C<->T；transversion 颠换，嘌呤与嘧啶之间，如A<->C/T,G<->C/T Rates among Sites ：位点之间的替换率分布。我们这里选择Gamma Distributed (G)，在 Gamma Parmamter中填的参数值可以通过MEGA计算得到。关于 Rates 的介绍可见后面备注2 Patterns among Lineages 默认即可，我们一般比对的都是同源序列。 Gaps/Missing Data Treatment：该参数设置如何处理空格(gaps)。几乎完全由空格组成的位点对估计树中的分支顺序或分支长度的贡献很小，但它们确实增加了计算时间。Complete deletion，代表直接删除所有包含gap的位点。Pairwise deletion代表配对删除，从而使缺少数据的位点仅在需要时被删除而不是总被删除。Partial deletion 代表部分删除，当 Site Coverage Cutoff设置为95%时，代表如果有95%的序列在此位点上是空格，则将该位点删除。 Select Condon Positions: 代表选择密码子的哪些位点来构建树，通常全选。除非你有特殊的需求，如只想根据密码子的第三位碱基构建进化树。 Number of threads：代表分析时电脑的线程数设置，根据自己的电脑情况可随意设置。
备注1 进化模型
在祖先序列进化为子序列的过程中，如果仅仅计算两个核苷酸之间的距离是不够的，比如某个位点由 A 替换成了 T，我们假设计算出的它们之间的距离是0.5，然而我们无法知道这中间是否发生了其它变化，比如可能先由 A 替换成 G，接着再从 G 替换成 T，这与 A 直接替换成 T 的替换率显然是不同的。又或者某个位点未发生替换，我们计算出的它们的距离为0，然而它中间也有可能经过了其它替换，比如A->G->A。因此我们需要选择进化模型来试图解释这些替换，以更准确地估计序列进化的过程。
单参数模型：
在单参数模型中，碱基替换可以通过以下矩阵表示，称之为Q矩阵，只有一种核酸替换率 α。
Jukes-Cantor model(JC)：该模型认为任何核苷酸改变为其它核苷酸的几率是相等的，只需要知道核苷酸替换率即可，因此称为 one-parameter model 或者 Jukes-Cantor model Felsenstein 81 (F81)模型拓展了 Jukes-Cantor model，允许碱基替换的基础替换率不同。
其它模型：
然而每种核苷酸的替换几率往往是不一样的
Kimura 2-parameter model(K2P) 拓展了Jukes-Cantor model，考虑核苷酸转换和颠换的替换率不一样。K2P模型假设 Q矩阵中有两种替换率，转换一种，颠换一种。
因此，K2P的Q矩阵如下，转换的替换率为α，颠换的替换率为b
Tamura 3-parameter model 添加了对成分偏差的修正，如果基础比率与相同频率有很大差异，可能是由于突变偏差，那么就需要考虑这种差异
Tamura-Nei model 延伸了前一个模型，把嘌呤之间的的转换和嘧啶之间的颠换的替换率进行了区分。
Felsenstein 84 (F84) model 和 HKY model也基于类似 Tamura-Nei model的假设
Maximum Composite Likelihood model 是一种基于似然性的 Tamura-Nei model的应用，提高了计算配对距离(pairse distance)的精确度
General Time-Reversible(GTR) model。该模型假设 Q矩阵有6种不同的率，认为碱基不同方向的替换概率相同(即A->C与C->A的替换率相同)。
总的来说，Kimura 2-parameter 与 Maximum Composite Likelihood 模型是目前构建NJ树时文献中用得比较多的模型。
备注2 位点之间的替换率
前面的进化模型中所说的碱基替换率指的是某个位点的碱基替换为其它碱基的频率（考察的是单个位点的碱基替换），这些模型假设所有位点的碱基替换率都是一样的，但实际上，不同位点的替换率可能会不同，比如编码序列的起始密码子处，通常为ATG（有时也为GTG），因此我们可以推测起始密码子处2，3位点的碱基替换率是0。
因此引入了序列中估计不同位点之间替换率的方法。Uniform Rates 假设不同位点之间的碱基替换率是相同的。还有一种假设是不同位点的碱基替换率的分布存在一种普遍的模式，常用的分布是伽马分布(Gamma distribution)，我们可以使用MEGA预测Gamma分布的参数。
备注3 进化树的可靠性检验
前面已经提过，事实上我们估计的所有树都是错误的，只是通过模型算法尽量找到一颗最接近真实的树。因此我们需要一种方法来评估通过这些模型构建出来的树的可靠性。其中 Bootstrap 方法是使用得最广泛的一种，NJ、ML等方法都使用 Bootstrap 检验进化树的可靠性，不过 BMCMC 通过进化枝后验概率计算树的可靠性。需要注意的是 Bootstrap 评估的是我们估计树的可重现性，而非树的准确性。
下面介绍 Bootstrap 的原理与过程。
首先，我们把未经过Bootstrap第一次估计出的树称为原始树，比如以下面的比对结果构建出一颗原始树：
上面的比对结果中包括10个碱基位点，5条序列。在第1次 bootstrap 时，会重新构建1个包含10个位点的比对结果，首先从10个位点中随机选取1个位点的序列作为伪比对第1个位点的碱基，接着重新从上面10个位点的序列中选取1个位点的序列作为伪比对第2个位点的碱基，如此直到形成一个新的包含10个位点的伪比对结果。这样得到的伪比对结果与原比对结果的不同在于有些位点可能会重复（有放回抽样），比如说第1次 bootstrap得到下面的伪比对结果，然后基于这个伪比对结果用相同的模型与参数估计一颗新的树：
接着按照相同的方式进行第2次、第3次…直到第1000次 bootstrap。最后会得到一颗最终的树，该树的每个节点旁都会记录一个自展值(Bootstrap value)。假如物种A的自展值为90，代表在1000次bootstrap中，有900次物种A在进化树中都被归为某一类(clade)，剩下的100次物种A被归类到别的clade。一般认为自展值>70的估计出来的分支是比较可靠的。
4.1.3 NJ 树的估计-实践
估计 gamma 参数：
按照前面的理论知识，如果我们认为自己的序列在不同位点具有不同的碱基替换率，可以设置gamma分布，所以先用MEGA得到gamma分布的参数：
依次点击，RATES -> Gamma Parameter for Site Rates (ML) -> 选择步骤3.3多序列比对得到的
.meg文件，不用更改参数，点击 OK即可
在结果页面，我们可以得到一个 gamma分布的参数为 0.499 ：
估计 NJ 进化树：
在 MEGA中依次点击：PHYLOGENY -> Construct/Test Neighbor Joining Tree -> 选择步骤3.3多序列比对得到的
.meg文件 -> 设置参数。根据前面的理论知识，我们按照下图设置参数，接着点击OK即开始估计进化树。
4.1.4 树的认识与MEGA中的可视化
进化树的基本结构
外部节点(绿色圆点)又被称为“叶节点”，表示采样及测序的实际生物体(例如，传染病研究中的病毒)，在进化生物学术语中又被称为“分类单元(taxa)”。内部节点(蓝色圆点)表示外部节点的假设祖先。根(红色圆点) 是进化树中所有物种的共同祖先。水平线条表示树的分支，又表示生物所发生的以时间或遗传分歧衡量的演变(灰色数字)。底部的线条表示这些分支长度的标尺。
MEGA中的进化树查看
本例中 MEGA 得到的进化树如下图1，先看红色框，这块区域显示了树的拓扑结构，与前一个图中稍有不同的是，在内部节点旁边还出现了1个数值（红色箭头所示），这是因为我们设置了 Bootstrap，该数值即表示自展值。可以看出现在的进化树是比较拥挤的，因为枝长按照实际计算出的值被画了出来，这样可以显示生物体演变的实际程度，此时树下方会有枝长的标尺。
蓝色框中是 MEGA 中绘制进化树的图形化界面选项，我们可以根据自己的需求定制不同风格的进化树，比如我们只想看一下各个分类单元的分类情况，可以点击蓝色框中 Layout -> Toggle Scaling of the Tree（粉红色箭头），就可以隐藏枝长的实际值（如下方第2张图所示），此时分支的长度不代表实际意义。
又比如我们前面提到自展值过小的估计出的分支一般不是很可靠，我们也可以通过蓝色框中的 Compute -> 设置自展值 Cutoff 值为50，则会合并自展值小于50的分支（如下方第3张图所示）。
另外我们还可以在蓝色框中设置是否显示自展值、枝长、显示为环形进化树、设置clade、提取子树等，只需要点击即可实时显示效果。
另外在绿色方框中会有两个树，Bootstrap Consensus Tree可以让我们更好地看清自展值与对应的节点，不过我们后续进化树可视化一般使用Original Tree 即可。
4.1.5 绘制有根树
关于无根树与有根树
我们构建的树（包括其它方法如ML法、MP法等）基本上得到的都是无根树。无根树没有方向性可言，能从一个节点追溯其父节点也可以追溯其子节点，如果移除无根树的某个分支，可以拆分到不同的分组，这些不同的分组称为 splits。有根树有方向性，不能从一个节点追溯其父节点，只能追溯其子节点，因为物种进化是有方向的，如果移除有根树的某个分支，得到的分组称为 clades。
虽然得到的是无根树，但是 MEGA 绘制出来的进化树看起来容易让人认为是有根树。我们在进化分析时往往希望得到的是有根树，因为只有有根树才能分析物种的演化关系。因此有时候通常会通过无根树人为绘制有根树。在一堆同源性较高的序列中去寻找它们的根往往不容易，因此我们往往需要借助额外的信息，常用的方法是借助外群(outgroup)。
外群(outgroup)可以被定义为1条或多条序列，这些序列与内群(ingroup)序列的进化距离比内群序列之间的距离都要远。
在前面的例子中，我选了18条 PRRSV 病毒的 ORF5 基因的序列，要绘制有根树，我需要另外加入1条外群序列，我选择的是与 PRRSV 病毒在同一个科下的另一种病毒的 ORF5 基因序列。在重新进行多序列比对、进化树估计后，得到下方的进化树，红框中的序列为新加入的外群序列，可以发现其与前面的18条序列相比形成了单独的一个分支：
接下来绘制有根树，首先在图中鼠标左键点击外群所在的分支，然后鼠标右键选择 Root tree 即可（如下图1所示），接着我们得到了1颗有根树（如下图2所示）。
虽然我们绘制的有根树看起来和无根树没什么区别，但实际上已经设置了1个根节点（下方图2中绿色箭头所示），而且仔细观察可以发现两颗树的分支长度发生了变化。上图无根树的分支看起来有两个分支分成了外群与内群，但实际上只有1个分支，其实际分支长度为这两个长度之和，MEGA中将其绘制成了看起来有两个分支的图形。下图2的有根树的分支确实存在两个分支，将其分成了内群与外群，因此分支长度也发生了改变。
接着我们就可以根据这个有根树添加不同的 clade 了。首先点击某个类别所有分类单元所在的共同分支（下图红色箭头），然后依次选择 Subtree -> Draw Options -> Format Subtree，即可进入设置 clade 的界面。（如果要删除设置的某个clade，还是这样选择，会出现一个Clear clade或者Clear All subtree options的选项，代表清除某个分类或者清除所有分类的设置）
另外有一系列选项可以进行子树的设置（下图1），大家可以自行探索。
我们可以发现 MEGA 自带的绘图功能其实也很强大，可以达到很美观的效果，设置 clade 后如下图2所示。
4.1.6 树的保存与导出
如果不需要对树进行进一步编辑，可以直接导出 MEGA 制作的进化树图片，选择 Image，可以保存为各种格式的图片文件。
如果想要将进化树数据导出，在其它软件中编辑，我们可以选择 File -> Export Current Tree(Newick)（下图1） -> 勾选两个选项，即同时导出枝长与自展值（下图2），接着会出现
.nwk格式文件的编辑窗口（下图3），直接点击保存然后命名为想要的名字即可。我这里保存为example.nwk等待步骤5使用。Newick格式是计算机可读形式的树文件标准格式，与 NEXUS、Phylip 等均是较为常见的储存系统发育树及其节点和分支相关数据的文件格式。我们通过MEGA得到树文件后，将其导出为 Newick 格式的树文件，然后可以使用 ggtree 等进行可视化。
我们也可以保存一份树文件，等待下一次使用MEGA时继续编辑，只需要点击 File -> Save Current Session -> 然后保存为 .mtsx或者.mts格式的文件即可。
4.2 构建 ML 树
最大似然法(Maximum Likelihood,ML)也是使用得比较普遍的一种进化树构建方法。因此我这里也介绍一下如何用MEGA构建MJ树。
最大似然法基本介绍：
最大似然法是一种基于进化模型和位点信息的树构建方法，其理论基础是基于两条基本假设：不同的性状进化是独立的、物种发生分歧后进化独立。
在ML法中，以一个特定的替代模型分析既定的一组序列数据，使所获得的每一个拓扑结构的似然率均为最大，挑出其最大似然率的最大拓扑结构为最终树。其主要优点是如果进化模型选择合理，ML树与真实进化的吻合度最高。然而，ML树的缺点在于计算强度大，非常耗时。
4.2.1 选择最佳模型
我们先看一下MEGA中估计ML树的参数，选择 PHYLOGENY -> Construct/Test Maximum Parsimony Tree(s) -> 可以进入参数设置页面，可以看到在 Model/Method 中提供了6个模型选项，Rates among sites提供了4个选项。
接着我可以们直接使用MEGA软件预测一个适合我们数据的上述两个参数的最佳组合。我们依次点击：MODELS -> Find Best DNA/Protein Models (ML)… -> 选择步骤3.3多序列比对得到的.meg文件 -> 参数设置页面不需要更改（如下图1），点击OK后一段时间会出现一个模型预测结果的窗口（如下图2）。
上图中红色框区域为MEGA为我们找到的模型结果，绿色框区域有一个对该结果的介绍。我们主要关注前5列（蓝色框所示）。
Model：替换模型与Rates参数简称的组合，从上到下排序为最好到最差。如 K2+G+I，代表 Model/Method 参数选择 Kimura 2-parameter 模型；Rates among Sites参数选择 Gamma Distributed With Invariant Sites(G+I)。 Paramaters：相应模型下必须要估计的参数的数量，在所有条件都相同的情况下，通常认为估计更少的参数会更好。这里的条件指的是用合适的模型来衡量（当然条件通常不相同） BIC、AICc、InL是三种方法计算出的分数，BIC分数计算采用的是贝叶斯信息准则(Bayesian Information Criterion)，分数越低代表模型越好，MEGA对Model列的排序就是基于BIC进行排序的。AICc分数计算采用的是校正后AIC信息准则(Akaike Information Criterion,correction)，分数越低代表模型越好。InL代表极大似然值(Maximum Likelihood value)，绝对值越低代表模型越好。
在本例中，虽然参数数量第一个模型组合比第二个模型组合大1，但其BIC、AICc、InL分数都是第一种模型组合更好一点，因此我们选择第一个组合（K2+G+I）。有时候第一个模型组合的BIC分数可能小于第二个模型组合，但是如果相差不大，且AIC、InL、参数数量都是第二个模型组合比较好的话，我们也可以考虑选择使用第二个模型组合。
4.2.2 估计 ML 树
如下图1，依次点击：PHYLOGENY -> Construct/Test Maximum Likelihood Tree -> 选择步骤3.3多序列比对得到的.meg文件，设置参数如下图2所示。点击OK即可开始比对。
关于参数的设置
Bootstrap：注意估计 ML 树设置 Bootstrap 检验会很耗时（电脑不太好的话有时候MEGA程序还会卡住或者闪退）。
Gaps/Missing Data Treatment与NJ树一致，其它参数我们直接选择默认即可，一般不需要修改。
最后我们可以得到最大似然法估计出的进化树，如下图所示，左下角红色箭头处会显示该树的log似然值（该值比真实树的似然值肯定会差一些）。其它关于树在MEGA中的可视化与NJ树一致，可见步骤4.1.4-4.1.6
5.使用 ggtree 进行可视化
前面步骤 4.1.6 我们得到了导出的 Newick 格式的树文件，接着就可以使用R包 treeio 处理以及 ggtree 进行美化了。可视化的部分并没有统一的规范，每个人的需求不同。建议有R语言与ggplot一定基础的同学去看一看Y叔的书《R实战：系统发育树的数据集成操作及可视化》，里面讲得很详细（了解ggplot2语法学起来会很快）。不想用R画图的同学可以借助一些软件如iTOL、Figtree等，或者直接使用前面介绍的MEGA中绘制出来的进化树图片。
下面我仅简单介绍一下如何使用treeio以及ggtree进行简单绘图，而不深入探索高级绘图的细节。我们一步一步探索ggtree的绘图过程：
如果没有安装ggtree，treeio，tidytree，需要先安装：
if(!require("ggtree",quietly = T)) BiocManager::install("ggtree")
if(!require("tidytree",quietly = T)) BiocManager::install("tidytree")
if(!require("treeio",quietly = T)) BiocManager::install("treeio")
5.1 treeio与ggtree使用介绍与练习
加载需要的包
library(treeio)
library(ggtree)
library(tidytree)
导入Newick格式的进化树
tree <- read.tree("example.nwk")
我们可以先看一下 read.tree() 函数读入树文件后储存的对象格式，phylo对象是储存系统发育数据常用的一种格式，系统发育分析中用到的大多数R包都依赖于phylo对象。
class(tree)
## [1] "phylo"
然后就可以使用ggtree画出最简单的树的拓扑结构了
ggtree(tree)
接着我们使用geom_tiplab()给每个分类单元添加名称，其实ggtree与ggplot函数的语法基本是一致的，如果要掌握ggtree最好先去了解一下ggplot语法。
ggtree(tree)+
geom_tiplab(size=1.5,color="purple") # size与color参数分别设置字体大小与颜色
我们前面画出的树中，其分支长度代表了物种实际的进化程度，因此可以使用geom_treescale()添加标尺，其中标尺x,y参数的设置只需要将整个图想象成第一象限的坐标轴，然后调整标尺的位置即可：
ggtree(tree)+
geom_tiplab(size=1.5,color="purple")+
geom_treescale(x=2,y=20)
前面的树不是很整齐，如果我们不需要体现分枝长度，只需要知道各分类单元属于哪个clade，那么可以设置 branch.length="none"参数去除分支长度，使树更加美观，layout参数可以指定树的布局方式。
ggtree(tree,branch.length = "none",layout="circular")+
geom_tiplab(size=1.5,color="purple")
这个phylo树对象在R语言中实际上是一个列表。前面画的矩形树图没有1个明显的根分支，我们可以给列表添加1个元素，即根分支的长度，我这里把长度设置为0.5，然后使用geom_rootedge()添加根分支图层
tree$root.edge <- 0.5
ggtree(tree)+
geom_tiplab(size=1.5,color="purple")+
geom_rootedge()
接下来我们探索如何对各个分类单元进行分组，即设置clade。在这之前我们先了解一下我们的树数据。因为phylo对象不容易直观理解，因此可以通过as_tibble()函数将phylo对象转换为tibble表格，与数据框差不多：
tree_tibble <- as_tibble(tree)
我们打开 tree_tibble 看一下就可以发现数据其实很简单，只有4列数据，node列保存了分类单元所在的节点编号，parent列保存了其父节点的编号，branch.length列保存了分支长度，label列保存了该节点对应的标签（叶节点的标签即分类单元的名称，内部节点的标签即自展值），关于节点、分支、分类单元等概念忘记的可以去步骤 4.1.4 进行回顾。
因为每个节点都有一个编号或者标签，因此我们可以通过该编号定位到某个节点。这样就可以实现对不同 clade 进行分类、标注颜色等操作。比如说下图中红色和蓝色圆圈所在的节点下所有分类单元分别可以被归为clade1与 clade2，我们只要往tree数据中添加基于此的分组信息即可方便地使用ggtree进行分组标记颜色与注释clade。
我们先看一下可以提取节点信息的几个函数，对于tree这个phylo对象的数据，parent()函数返回了第1个节点的父节点的节点编号
parent(tree,1)
## [1] 26
相应的还有child()函数可以提取某个节点子节点的编号，
child(tree,26)
## [1] 1 2
ancestor()函数提取某个节点所有父辈节点编号，直到追溯到根节点
ancestor(tree,1)
## [1] 26 25 24 23 22 21 20
offspring()函数提取某个节点所有子代节点编号
offspring(tree,20)
## [1] 21 31 19 22 11 23 28 24 27 25 4 26 3 1 2 5 6 7 29 8 30 9 10 12 32
## [26] 33 35 34 15 13 14 16 36 17 18
sibling()函数返回某个节点的兄弟节点
sibling(tree,1)
## [1] 2
上述函数也可以针对前面转换的tibble数据进行提取，可以得到更全的信息，比如：
child(tree_tibble,26)
## # A tibble: 2 × 4
## parent node branch.length label
## <int> <int> <dbl> <chr>
## 1 26 1 0.0102 MT268280.1
## 2 26 2 0.00610 FJ548853.1
回到前面讲的，我们的目的是获取红蓝色圈对应的节点编号，以实现分类，因此我们只需要找到根节点编号，然后根据根节点编号获取其两个子节点编号就是我们相要的两个节点编号。
有一种比较简单的方法，先将所有节点编号显示在图中，然后我们从图中找到直接用于后面的输入，使用geom_nodelab()函数映射节点编号即可，可以发现我们想要的编号是21与31：
ggtree(tree,branch.length = "none")+geom_nodelab(aes(label=node))
还有一种方法是直接提取，我们首先通过 ancestor(tree,1)得到任意节点（这里选择了节点1）的所有祖先节点，它返回的编号顺序是一级一级往前推的，因此最后一个编号就是根节点编号，然后通过rev()函数反向，再取其第一个编号就是根节点编号了。
( node_root <- rev(ancestor(tree,1))[1] )
## [1] 20
然后我们只需要使用child()函数即可得到想要的节点编号，示例的进化树中根节点分出了3个节点，前两个就是我们想要用于分组的节点。
( node <- child(tree,node_root) )
## [1] 21 31 19
接着通过groupClade函数，我们可以给树添加分组信息，只需要传入我们想要分组的节点编号的向量即可。下面的例子中，groupClade函数把节点21的所有子节点分为一组，把节点31的所有子节点分为另一组，转换为tibble数据查看一下，我们可以明显发现数据中添加了1列group信息，1代表第1组，2代表第2组，如果有别的节点不包含在我们输入的节点的子节点中，该函数会默认将这些节点标记为另外一组0。（需要注意groupClade只能对内部节点进行分组，想要将叶节点进行分组可以使用groupOTU函数）
groupClade(tree,c(21,31)) %>% as_tibble()
## # A tibble: 36 × 5
## parent node branch.length label group
## <int> <int> <dbl> <chr> <fct>
## 1 26 1 0.0102 MT268280.1 1
## 2 26 2 0.00610 FJ548853.1 1
## 3 25 3 0.00733 OQ883907.1 1
## 4 24 4 0.00288 OQ357724.1 1
## 5 27 5 0 NC_038291.1 1
## 6 27 6 0.00921 AF331831.1 1
## 7 28 7 0.0821 JN660150.1 1
## 8 29 8 0.0547 MW803134.1 1
## 9 30 9 0.0601 MZ416787.1 1
## 10 30 10 0.0816 JN654459.1 1
## # ℹ 26 more rows
分组以后我们需要将其保存为一个新的tree变量才可以进一步画图，这里保存为tree2。
tree2 <- groupClade(tree,c(21,31))
接着我们只需要使用类似于ggplot2的语法，使用aes函数将group分组信息映射到图中即可
p <- ggtree(tree2,aes(color=group))+
geom_tiplab()
p
接下来我们要进一步在分类单元右边添加clade标注。先构建一个数据框，node中还是填入刚才用于分类的编号形成的向量，name中填入对应的clade名称组成的向量，比如我想把节点21所在的分类命名为L1，把节点31所在的分类命名为L2,可以这样设置：
dat <- data.frame(node=c(21,31),
name=c("L1","L2"))
然后使用geom_cladelab函数添加一个图层，将前面制造的dat数据映射上去即可，offset参数代表设置文本到条带之间的距离。
p+geom_cladelab(data = dat,
mapping = aes(node=node,label=name),
offset=0.5)
如果我们想要移除某个分类单元怎么办呢？比如我现在不想NC_001639.1出现在图中，只需要使用drop.tip函数传入需要移除的分类单元名称组成的向量即可，然后对新的树tree3绘图。
tree3 <- drop.tip(tree2,"NC_001639.1")
ggtree(tree3)+
geom_tiplab()+
geom_cladelab(data = dat,
mapping = aes(node=node,label=name),
offset=0.5)
接下来，如果我们只想给其中几个分类单元标上颜色，只需要映射想要标记颜色的分类单元到颜色参数即可。
可以直接映射分类单元名称
mark <- c("MT268280.1","JN660150.1","OQ803248.1")
ggtree(tree3,aes(color=label %in% mark))+
geom_tiplab()+
geom_cladelab(data = dat,
mapping = aes(node=node,label=name),
offset=0.5)
也可以映射分类单元的节点编号：
ggtree(tree)+
geom_tiplab(size=2,aes(color=node %in% c(1,7,14)))
另外treeio也可以提取某个子树进行画图，比如我现在只想要L1所对应的子树，只需要使用subtree函数对子树进行提取即可。node参数指定要提取的子树的某个叶节点的编号或者标签均可以，levels_back参数指定子集树应该包括多少从选定节点返回的节点（即叶节点往前追溯多少数量的节点）
tree4 <- tree_subset(tree,node=1,levels_back = 6)
ggtree(tree4)+
geom_tiplab(size=2,color="purple")
5.2 进化树的可视化实践
下面我展示了我自己画出来的一个进化树：
# 载入包
library(treeio)
library(ggtree)
library(tidytree)
#导入树文件
tree <- read.tree("example.nwk")
tree <- drop.tip(tree,"NC_001639.1") # 移除不需要展示的分类单元
#设置根分支长度
tree$root.edge <- 0.3
#获取分类节点编号
tree_tibble <- as_tibble(tree)
node_root <- rev(ancestor(tree,1))[1] #获取根节点编号
node <- child(tree,node_root)[-3] #获取分类节点编号
tree2 <- groupClade(tree,node)
#定义clade名称
dat <- data.frame(node=node,
name=c("L1","L2"))
#ggtree可视化
p <- ggtree(tree2,
aes(color=group), #映射分组信息
branch.length = "none",
size=1)+ #分支粗细
geom_tiplab(size=3,color="purple",offset = 0.2)+ #叶节点字体大小与颜色
geom_tippoint(size=2)+ #给叶节点添加圆点
geom_cladelab(data = dat, #添加clade分类数据
mapping = aes(node=node,label=name),
offset=1, #线条与分类单元文本之间的距离
offset.text=0.02, #线条与clade文本之间的距离
geom="label", #为clade文本添加背景
fill="lightblue", #clade文本背景颜色
fontsize=3.5, #clade文本字体大小
barsize=0.7, #线条粗细，align参数即是否对齐
align=T,
barcolor="grey")+#线条颜色
geom_rootedge()+ #添加根分支图层
theme_tree("#FEE4E9")+ #设置主题
theme(legend.position = "none") #去除图例
p
#保存进化树图片
ggsave(p,file="plot1.png",width = 8,height = 12)
在得到上面的树后，如果我们想要调换某个分支的展示顺序，只需要简单地通过rotate函数即可完成,比如我对上图的节点20对应的分支调换展示顺序：
rotate(p,20)
我们还可以将进化树结果与多序列比对结果放在同一张图中进行展示，多序列比对结果（example.fas）我们已经从步骤3.3得到
#载入需要的包
library(treeio)
library(ggtree)
library(tidytree)
#导入Newick树文件与多序列比对FASTA文件
tree <- read.tree("example.nwk")
msa <- read.fasta("example.fas")
#因为read.fasta函数导入多序列比对结果后_会识别为空格，因此需要先将空格转换为_，否则会导致无法与进化树的分类单元名匹配而绘图失败。
library(stringr)
names(msa) <- str_replace(names(msa)," ","_")
#先绘制进化树
p1 <- ggtree(tree,branch.length = "none")+geom_tiplab(size=2)
#然后用msaplot函数将两个结果绘制在一起
msaplot(p1,msa,
offset = 3,width = 2) # offset为多序列比对结果与进化树之间的距离
最后，进化树的可视化就完成了，上面只展示了一部分比较基础的绘图方法，想要继续绘制更多高级的好看的图形还是需要自己去看书学习。
这次的分享就到这里啦，如果喜欢就点个赞吧。
参考资料
《Phylogenetic Trees Made Easy》第五版，Barry G. Hall，2018
《R实战 系统发育树的数据集成操作及可视化》第一版，余光创，2023
《生物信息学》第一版，李霞，2010
https://mp.weixin.qq.com/s/8LWMm7wwfKehGxTNCqEwkw