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。


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 脚本格式转换速度较快，例如，处理 20GB VCF（约 300w SNP x 650 个体）约 27 分钟。除了默认的输出格式 PHYLIP 外，还支持将 VCF 转换为 FASTA、NEXUS 或二进制 NEXUS 格式。

vcf2phylip.py 脚本的下载、详细介绍参见：https://github.com/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）和 BIONJ 2 种方法构建系统发育树，仅进行 2 轮树的优化（-nstop=2），且无 OPTIMIZING CANDIDATE TREE SET 步骤。否则（不使用 -fast 参数）使用 100 种构建方法（default -ninit=100），并至多迭代 100 次（default -nstop=100）进行树的优化。
  • 经过 ModelFinder 检索（-m MF），本数据集的最佳替换模型为 GTR+F+R6+ASC，所以这里直接指定 -m GTR+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        1234567
  species_1   CCCCCCC
  species_2   RTCYRCC
  species_3   CYTYYRT
  species_4   CCCYCRT

  site1 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

input

model

optimize

Log-likelihood

Total tree length

Sum of internal branch lengths (ratio)

near-zero internal branches

Lmap not informative ratio

1

min4.phy

GTR+F+G4

fast

-95636

3.751

0.799 (21%)

93

8.48%

2

varsites.phy

GTR+F+G4

fast

-80786

7.216

1.549 (21%)

97

14.82%

3

varsites.phy

GTR+F+G4+ASC

fast

-79848

5.536

1.185 (21%)

101

0%

4

varsites.phy

GTR+F+R6+ASC

fast

-77426

1.327

0.269 (20%)

263

0%

5

varsites.phy

GTR+F+R6+ASC

complete

-77239

1.079

0.210 (19%)

316

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

Repeat

Log-likelihood

Total tree length

Sum of internal branch lengths (ratio)

near-zero internal branches

1

-72435

1.383

0.281 (20%)

287

2

-72046

1.892

0.385 (20%)

217

3

-72416

1.724

0.345 (20%)

234

4

-72531

1.607

0.332 (20%)

220

  • Log-likelihood 变化较小，但 tree length 浮动较大。near-zero 与 Sum of internal 与 tree length 相关性较高。
  • 树的拓扑结构存在差异。


参考

http://www.iqtree.org/doc/Substitution-Models

  1. Binary and morphological models
  2. Ascertainment bias correction

http://www.iqtree.org/doc/Command-Reference

http://www.iqtree.org/doc/Frequently-Asked-Questions

  1. Why does IQ-TREE complain about the use of +ASC model?
  2. What are the differences between alignment columns/sites and patterns?