IQtree 是一款用于构建系统发育树的软件,支持 Phylip、Fasta、Nexus、Clustalw 等多序列比对后的数据格式。但随着 SNP 数据费用的大幅下降,研究人员利用 SNP 数据构建系统发育树的情况越大越多。本文为用 vcf 格式的 SNP 数据构建系统发育树的教程,主要包括:Vcf 文件转 Phylip 文件、Phylip 文件输入 IQtree 建树 两部分。
数据集
本教程使用的数据集矩阵大小为 577 × 10000。样本取自 577 个玉米品系,包含 70 个大刍草(par)、208 个热带品系(tst)、172 个温带品系(temp)、127 个混合品系(mixed)。SNP 取自 3 号染色体上 743w SNP,经过过滤(maf > 0,miss = 0)与随机抽样后,最终得到 1w SNP。数据已上传(链接如下),免费下载(0 积分)。
IQtree:使用SNP数据(vcffile)构建系统发育树(数据)
1. Vcf 文件转 Phylip 文件
python vcf2phylip.py --input SNP.vcf--input FILENAME Name of the input VCF file, can be gzipped (but the extension must be .vcf.gz)
输出文件:SNP.min4.phy
使用 vcf2phylip.py 脚本将 vcf 文件转化为 phylip 格式。vcf2phylip.py 脚本格式转换速度较快,例如,处理20GBVCF(约 300w SNP x 650 个体)约 27 分钟。除了默认的输出格式 PHYLIP 外,还支持将 VCF 转换为 FASTA、NEXUS 或二进制 NEXUS 格式。
vcf2phylip.py 脚本的下载、详细介绍参见:/edgardomortiz/vcf2phylip
2. Phylip 文件输入 IQtree 建树
iqtree -s SNP.min4.phy -keep-ident -st DNA -m GTR+F+R6+ASC -fast -nt 50-s输入文件-keep-ident保留序列相同的样本-st指定输入序列的类别(当 iqtree 识别错误时使用)-m选择替换模型(substitution model)-fast加速建树-nt最大线程数
程序大概率会报错(Error)并输出文件:SNP.varsites.phy
ASC(ascertainment bias correction)是 IQtree 处理SNP数据时需要添加的参数,可以矫正因 SNP 数据所导致的系统发育树分枝长度过长与结构错误等问题(参见附录 ASC)。ASC模式要求输入的 SNP 必须是variable sites,但 IQtree 对 constant、singleton、variable 概念定义较为特殊(参见附录 Ambiguous site),可能 phy 文件中有大量的位点被识别为 constant sites。如本教程数据集 1w SNP 中有 6466 SNP 被识别为 constant sites。若 phy 文件中包含 constant sites,IQtree 会报错并输出varsites.phy文件,包含输入数据中的 singleton 与 variable sites。IQtree 默认在建树过程中会删除相同序列的样本,待建树完成后再添加,以节省分析时间。但这也会使 varsites.phy 输出文件中的样本数不足 577。虽然本数据集中不存在样本序列相同,但如果 SNP 较少且样本间亲缘关系较近,则群体中可能包含相同序列。-keep-ident参数可以保留序列相同的样本,使输出样本数与输入样本数保持一致。-fast参数通过仅使用最大简约(maximum parsimony)和BIONJ2 种方法构建系统发育树,仅进行 2 轮树的优化(-nstop=2),且无 OPTIMIZING CANDIDATE TREE SET 步骤。否则(不使用 -fast 参数)使用 100 种构建方法(default -ninit=100),并至多迭代 100 次(default -nstop=100)进行树的优化。经过 ModelFinder 检索(-m MF),本数据集的最佳替换模型为 GTR+F+R6+ASC,所以这里直接指定-mGTR+F+R6+ASC 模型。
iqtree -s SNP.varsites.phy -m GTR+F+G4+ASC -fast -nt 50
将输出文件 SNP.varsites.phy.treefile 放入 iTOL 中即可观察所构建出的系统发育树结构。
附录
Ambiguous site
IQtree 将位点(site)分为 3 类:constant(invariant)、singleton、informative(variable)。constant site 指位点在所有样本中只包含 1 种碱基;singleton site 指位点包含 2 种类型碱基,但突变仅出现 1 次;informative site 指位点包含 2 种类型碱基,且每种碱基出现至少 2 次。从算法角度来说,IQtree 无法利用 singleton 对样本分群,即无法提供分群信息,所以 constant 与 singleton 位点属于 uninformative site 。
IQtree 对歧义位点(Ambiguous site)位点采用交集(intersection)形式分析:如果交集非空则忽略歧义碱基,若为空则视为 variant,同列中歧义碱基不堆叠。
site 1234567species_1 CCCCCCCspecies_2 RTCYRCCspecies_3 CYTYYRTspecies_4 CCCYCRTsite1 is singleton, because R={A or G}, which exclude C, uninformative.site2 is singleton, because Y={C or T}, which include C, ignore Y, uninformative.site3 is singleton, uninformative.site4 is constant, because Y={C or T}, which include C, ignore Y, uninformative.site5 is singleton, because R={A or G}, which exclude C, Y={C or T}, which include C, ignore Y, uninformative.site6 is singleton, because R={A or G}, which exclude C, uninformative.site7 is informative.
由于歧义位点的交集策略,当同列中仅包含非空歧义碱基时,位点会被判定为 constant site(如上例中 site 4)。因此 SNP 数据从 vcf 转 phylip 后存在部分 site 被判定为 constant 的情况。如本教程 1w SNP 中包括 1807 parsimony-informative sites,1867 singleton sites,6466 constant sites。
ASC
IQtree 构建系统发育树时,informative sites 用于推断树的拓扑结构与分枝长度,uninformative site 则用于矫正。由于个体间固定长度区间内 SNP 越多,说明个体间积累的突变越多,分离的时间越长。所以SNP构成的 Fasta 数据大幅增加了突变的密度,使个体间分离时间大幅提前,导致分枝过长。
同时,分枝过长会压缩祖先之间的分枝长度,降低祖先之间拓扑结构推断的准确性。因为祖先与个体之间时间跨度越长,祖先的可能性就越多,如果祖先对应一个概率分布,则祖先与祖先之间的分布交集就越大。
例如 A / B 的祖先是 D,C 的祖先是 E,D / E 的祖先是 F,如果 A / B / C 与 D / E 的距离越远,D / E 的可能序列就越多,D / E 分布重合的区间就越大;因为 D / E 间相似度高,所以 D / E 与 F 的距离会缩短,即压缩 D / E 与祖先 F 之间的分枝长度,降低树结构预测的准确性。
为了矫正输入 SNP 数据导致发育树分枝过长的情况,IQtree 提供 +ASC 参数(ascertainment bias correction)。使用本教程数据,分别评测了 GTR+F+G4 和 GTR+F+G4+ASC 两种替换模型下的树结构,结果参见下节(Test 2&3)。+ASC 虽然对 Log-likelihood 影响较小(1.2%=1-79848/80786),但对树结构影响较大(23.3%=1-5.536/7.216),且能够大幅提高稳定性(148%->0%)。
不同参数的性能测试
Test 1&2:使用 IQtree 过滤 SNP 数据集中的 constant site 后再输入 IQtree,不仅可以使用 +ASC 提高模型精度,还能减少建树时间(min4 299s,varsites 50s)。Test 2&3:ASC使树的分枝长度大幅降低(23.3%=1-5.536/7.216),且能够大幅提高稳定性(148%->0%),但对似然值改变较小(1.2%=1-79848/80786)。Test 3&4:更换模型使树的分枝长度大幅降低(75%=1-1.328/5.536),但似然值改变较小(3%=1-77426/79848)。Test 4&5:多轮迭代使树的分枝长度大幅降低(20%=1-1.079/1.328),但似然值改变微小(3‰=1-77239/77426)。Test all:所有模型中中间分枝(internal branch)长度占总枝长比例基本固定(约 20%),near-zero 与 tree length 相关性较高,当发育树整体分枝长度增加时,内部分枝长度也相应增加,其中近 0 分枝数量减少。中间分枝长度占比 20 % 说明80%的枝长分布在末端分枝,个体间距离较大,分离时间较早。部分祖先间的进化关系太紧密(near-zero),用多叉树描述可能更符合这些节点。
重复
iqtree -s SNP.varsites.phy -m GTR+F+R6+ASC -p SNP.varsites.nex -fast -nt 50
Log-likelihood 变化较小,但 tree length 浮动较大。near-zero 与 Sum of internal 与 tree length 相关性较高。树的拓扑结构存在差异。
参考
/doc/Substitution-Models
Binary and morphological modelsAscertainment bias correction
/doc/Command-Reference
/doc/Frequently-Asked-Questions
Why does IQ-TREE complain about the use of +ASC model?What are the differences between alignment columns/sites and patterns?