构建系统发育树

生信喵

背景
      构建系统发育树属于群体遗传学分析范畴，随着时间和地理位置的变化，新冠病毒经过多次迭代，在基因组上会累积不同的突变，已经与祖先产生明显的不同。通过对多个序列进行系统发育分析，不仅可以厘清不同物种之间的亲缘关系，而且可以重塑新冠病毒的演化过程，具有重要的现实意义。例如某地新发疫情，可以对样本快速测序，构建全基因组序列，然后对其进行系统发育分析，快速定位到系统发育树中，可以快速鉴定新发菌株的亲缘关系，对于疫情防控溯源具有重要的指导作用。
      构建系统发育树本质上是一种聚类分析，通过不同基因组之间两两比对，构建距离矩阵，然后进行聚类。
      首先，将多个样品基因组合并为一个文件，然后进行多序列比对。比对之后就可以根据两两样品之间序列的差别构建距离矩阵，然后进行聚类，构建系统发育树。本节中我们将比较新冠病毒各个突变株以及 SARS 等已有序列，构建系统发育树，比较各个基因组之间的亲缘关系。

一、下载序列
      下载新冠病毒多个突变类型的序列，并加上参考系列以及 SARS 病毒参考菌株，蝙蝠分离株参考序列一起构建系统发育树。利用 efetch 批量下载各个序列，直接输入序列的Accession Number 即可。
      Accession Number 与 Name 对应表
      

1. 下载测序数据
   
2. mkdir 32.tree;cd 32.tree;
   
3. /share/Software/miniconda3/envs/test/bin/efetch -db nuccore -format fasta -id "MN908947 MG772933 KT444582 MZ310552 MZ202314 MZ169911 MZ318159 MZ373479 MZ169912 MW852494 MZ257684 MZ310903 MZ310580">ncov.fasta
   

复制代码

     为了后面便于识别，可以将 Accession Number 修改为常见的名字，例如将 KT444582 修改为 SARS，MG772933 修改为 batSARS，当让也可以后面直接在系统发育树中修改 ID。

二、多序列比对
      构建系统发育树的基础是多序列比对。所谓多序列比对（Multiple sequence alignment；MSA），是对三个以上的生物学序列，如蛋白质序列、DNA 序列或 RNA 序列所作的序列比对。一般来说，是输入一组假定拥有演化关系的序列。从多序列比对的结果可推导出序列的同源性，而种系发生关系也可引导出这些序列共同的演化始祖。如图可视化展示的多序列比对，例如单碱基的错配，或是如删除突变与插入突变。都可以显示出来。多序列比对的结果常常用于系统发育分析。
      多序列比对的工具有很多，可以使用 muscle，clustalW，mafft 等，比对之后直接就可以用于构建系统发育树
      mega 多序列比对截图
      

1. #muscle 多序列比对
   
2. muscle -in ncov.fasta -out ncov.clw -clw
   
3. #treebest 画树
   
4. treebest nj ncov.clw >ncov.nwk
   

复制代码

     
      比对时间26min

三、不同算法之间的差别
      构建系统发育树，本质上是一种聚类分析。聚类分析属于一种数据规约技术。对于层次聚类，最常见的算法包括单联动，全联动，平均联动，也就是 UPGMA，质心和 Ward 法。对于划分聚类，最常用的算法是 K 均值，也就是常见的 kmean，和围绕中心的划分（PAM），有不少进化树的软件名字中都包含 PAM，用的就是围绕中的划分算法。
      那么这些算法是如何使用的呢？主要用来计算距离上。也是画树之中最重要的步骤。假设10 个样品，两两之间都要进行比较，比较完了需要给定一个量化指标，这就是距离矩阵。

距离如何计算呢，其实原理并不难。假设有三个样品 ABC。
      计算距离矩阵，最简单的就是欧几里得距离，也称为欧氏距离。差值的平方和在开根号。除了欧氏距离，还有马氏距离，闵可夫斯基距离，切比雪夫距离等等。欧氏距离适合连续型变量，比如上面都是数字，那么系统发育树中是序列，应该属于名义型变量，不用欧氏距离了。

前面介绍过层次聚类有五种算法，主要差别就是在计算距离上的不同。
      单联动是一个类与另一个类中点的最小距离。
      全联动是一个类与另一个类中点的最大距离。
      平均联动，顾名思义，只一个类与另一个类中的点的平均距离。
      质心，两类中质心之间的距离。
      ward 法，两个类之间所有变量的方差分析的平方和。
这是层次聚类，划分聚类主要就是 K 均值法和基于中心点的划分 PAM 方法。
那么在实际过程中该使用哪种方法呢，这个要根据具体的数据特点。
下面我们总结一下几种画树软件的算法的使用。
      1、NJ 法计算速度较快，适合序列相似度较高的序列。对于相似度很低的序列，会出现所谓长枝吸引现象。
      2、数据是连续型时，比如一个数字矩阵，差异表达数据属于这种，可以使用 UPGMA 法计算距离。


四、利用 MEGA 构建系统发育树
      mega 的全称是（Molecular Evolutionar Genetics Analysis），是一款非常好用的软件。也是一款非常经典的进行系统发育分析的软件。mega 不像 GATK 那种软件，比如找snp，后者操作特别繁琐。
      Mega 已经有很多年的历史了，应该是 93 年诞生的。mega 的引用率非常高。我们在文献中看到很多系统发育树的图，都是用 mega 来做的。最近几年基因数据越来越多，软件更新的也比较频繁。之前 mega 至于 windows 版本和 Mac 版本，都是图形化界面。也是因为数据越来越大，后来又有了命令行版的 mega-cc 工具。现在也有 Linux 的图形化工具了。
      以前软件版本是 Mega7，最近又出了最新版的 mega x，目前最新版本的 mega 软件为mega11。
      mega 的网站：https://www.megasoftware.net/

      mega 系统发育树结果
      
4.1 MegaGUI
      mega 提供图形化的 GUI 版本以及命令行的 megacc。图形化的版本使用起来更方便，里面集成了多序列比对，计算距离矩阵以及构建系统发育树等功能。使用 mega 比对之后直接够可以用于构建系统发育树了。但是图形化版本的 mega 安装在个人电脑上，需要较大的计算资源，这种情况下就需要命令行版本的 megacc。
      
      这就是图形化的在比对了，我已经听到电脑的风扇在疯狂旋转了。
      
      比对了接近50min（内存设置1g多）

4.2 megacc
      在线帮助文档：https://www.megasoftware.net/web_help_10/index.htm#t=MEGACC_Overview.htm
4.2.1 生成.mao 文件
      1、切换到 prototype 模式

      
      2、选择分析数据类型

      
      3、点击菜单上对应的分析模式，比对或者建树等；
      
      4、设置好参数，保存文件

      
4.2.2 使用案例

1. #megacc 比对
   
2. megacc -a muscle_align_nucleotide.mao -d ncov.fasta -o mega/ncov -f MEGA

复制代码

     选项参数：
       -a 配置文件，mao 格式
       -d 输入文件，meg 格式或者 fasta 格式
       -o 输出结果文件

1. #mega 画树
   
2. megacc -a infer_UPGMA_nucleotide.mao -d mega/ncov.meg -o mega/ncov

复制代码

4.3 bootstrap
      boostrap 值，也叫作自展值，就是图中枝节点是行的数字，100,99,94 等。它是用来检验进化树分支可信度的。简单地讲就是把序列的位点都重排，重排后的序列再用相同的办法画树，如果原来树的分枝在重排后的树中又出现了，就给这个分枝加 1。
类似于选举的时候画正字。一般需要设置重排 1000 次，也就是 bootstrap 值设置为 1000，最终重现的次数占重排次数的百分比就是自展值。99 就表示百分之 99%的频率。通常认为自展值大于 75%是可靠的，再小点 70%也凑合，如果自展值过低，说明这个树的分枝关系并不可靠。


4.4 nwk 文件
      Newick 格式，windows 系统下扩展名为 nwk，这种格式是一种比较通用的树文件格式，有个网站专门介绍这种格式，其实并不难，都是一些纯文本，表示拓扑结构，看起来稍微有点乱。这个就是系统发育树了。在 treeview 软件中，就会变成图片来显示了。

http://evolution.genetics.washington.edu/phylip/newicktree.html