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?
