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比对，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
网友评论