手里的玉米VCF文件终于整理完了,几百份材料、几十万个变异位点躺在服务器上,接下来最想干的一件事,就是先建一棵系统发育树,快速看群体之间的亲疏远近。这类分析在玉米群体遗传里太常见了,无论是自交系谱系梳理、地方品种群体归属,还是育种材料亲缘关系排查,都会用同一套思路:从SNP数据出发,用最大似然法把样本聚成有层次的树,而IQtree是最常用的建树工具之一。
在这篇文章里,我会把从VCF到系统发育树的整个流程讲透,包括IQtree安装、VCF过滤、模型选择、跑树命令、结果解读,以及一堆文档里不会写的坑。核心不是机械地跑完一条命令,而是理解每个参数为什么这么设,遇到问题能从原理层面去排查,而不是到处搜答案。
适合谁看呢?适合手里已经有一份VCF文件、想用SNP数据做玉米或其他作物群体遗传分析,但还没正式跑过IQtree的读者。不管你是刚接触群体遗传的学生,还是想换工具重新处理旧数据的从业者,按这套流程走基本不会跑冤枉路。如果你只是想找现成命令,可以直接跳到第四部分,但我还是建议把预处理看完,因为这种项目一半以上的坑都埋在VCF处理环节。
1. 项目概述与整体思路
1.1 为什么群体遗传分析要先建树
群体遗传分析里,系统发育树是最直观的一张图。它的本质是基于样本间的遗传差异程度,把所有个体重新排列成一颗有层次的树,亲缘关系近的样本聚在一起,关系远的自然分离。在玉米这种驯化作物里,数据通常包含不同热带、温带亚群甚至野生近缘种,一棵树扫一眼就能看出样本是否按预期分组、有没有标注错误、有没有混入异源材料。我自己的做法是拿到处理干净的VCF之后,不给任何预设信息,先让树自己说话,再结合群体结构分析去验证。
但这里有一个前提必须说清楚:系统发育树适合描述分化关系,不适合描述存在大量基因流和重组的群体结构。玉米是异交为主的作物,群体内部重组信号很强,树的分枝不一定完全代表真实亲缘关系,尤其在亚群内部。我的习惯是把树当作“第一眼筛选工具”,而不是最终结论。它可以快速暴露异常样本,也能为后续PCA、admixture分析提供方向;但反过来,这些分析之间是互补关系,不能互相替代。建树和PCA的视角不同,树更强调历史的谱系关系,PCA更擅长展示连续的遗传梯度,两者结合才能给出完整判断。
1.2 为什么选IQtree而不是其他建树工具
建树工具很多,老牌的RAxML、FastTree、MEGA、MrBayes各有拥趸。我选IQtree 2,核心原因是它在速度、准确性和易用性之间平衡得最好。RAxML-NG很快,但不支持直接读VCF;FastTree在大数据量下很能打,但模型相对粗糙,也没有VCF输入接口;MEGA适合教学和小数据量,跑几十万位点会非常吃力。IQtree 2有几个非常实用的特点:内置ModelFinder自动选模型、支持VCF直接输入、自带UFBoot2和SH-aLRT两个分支支持度检验、多线程和内存控制都很成熟。
| 工具 | 最大亮点 | 直接支持VCF | 适合场景 |
|---|---|---|---|
| IQtree 2 | 自动模型选择、VCF直读 | 支持 | 群体规模SNP建树,速度与准确度平衡 |
| RAxML-NG | 高并行、大规模ML树 | 不支持 | 基因组级数据,传统PHYLIP输入 |
| FastTree | 极快,近似ML | 不支持 | 超大数据量快速看树形 |
| MEGA | 图形界面友好 | 不支持 | 教学、小数据集 |
这种对比不是说要否定其他工具,而是提醒你在选工具时先看输入格式。如果你的数据管线已经习惯生成PHYLIP或FASTA,RAxML也完全能用。但对我这种日常工作要处理大量VCF文件的人来说,少一步格式转换、少踩一步坑,就是实打实省时间。IQtree还有一个好处是模型选择集成在同一个程序里,不用像RAxML那样先跑ModelTest再跑建树。
1.3 完整流程预览
我自己的标准流程固定五步:拿到VCF后先做位点和样本层面的过滤;按需做LD剪枝控制位点间冗余;决定是直接输入VCF还是转成PHYLIP格式;跑IQtree,包含模型选择、建树、支持率评估;最后解读输出文件并可视化。后面每一部分都会展开。如果只是想先跑通,可以直接跳到第四部分复制命令;但我还是建议把预处理看完,因为大多数“树很奇怪”的问题都出在VCF没处理干净,而不是建树软件本身。
整个项目从拿到VCF到出图,通常一周内可以完成,其中大头时间往往耗在VCF过滤和反复调参上。真正跑IQtree的时间反而不长。把握好这个时间分配,你的效率会明显提升。
2. 环境准备:把IQtree装稳再谈建树
2.1 三种安装方式怎么选
IQtree安装方式常见有三种:conda/mamba环境、官方预编译二进制、源码编译。我推荐前两种,源码编译除非你对自己编译依赖很有把握,否则没必要碰。
用mamba创建独立环境最省事:
mamba create -n phylo -c bioconda -c conda-forge iqtree conda activate phylo iqtree2 -vconda环境的好处是所有依赖自动打包,不会污染系统里的其他软件。这在集群上尤其重要,因为很多服务器管理员不会给你装系统级软件,你自己装一个conda环境就能解决权限问题。如果你机器上还没有mamba,可以用conda install -n base -c conda-forge mamba先装上。
官方预编译二进制适合没有conda环境或者想快速部署的场景:
wget https://github.com/iqtree/iqtree2/releases/download/v2.3.6/iqtree-2.3.6-Linux.tar.gz tar -xzf iqtree-2.3.6-Linux.tar.gz iqtree-2.3.6-Linux/bin/iqtree2 -v解压后把二进制路径加进PATH,或者建立软链接。没有sudo权限就把下面这行写进~/.bashrc,然后source ~/.bashrc:
export PATH=/path/to/iqtree-2.3.6-Linux/bin:$PATH需要注意版本。IQtree从2.0开始才支持VCF直接输入,1.x版本完全没有这个功能。如果你在集群里只找到了老版本,用VCF输入会直接报错或不识别参数,所以起手先确认iqtree2 -v输出的至少是2.x版本。
2.2 安装自检与常见报错
安装完成后第一件事就是跑一条最小命令验证:
iqtree2 -v如果提示command not found,多半是PATH没设置好;如果提示没有执行权限,用chmod +x iqtree2修复即可。
我遇到比较多的是这个报错:
iqtree2: error while loading shared libraries: libgomp.so.1: cannot open shared object file这是缺少OpenMP运行库。在conda环境里基本不会出现,系统裸装时可以通过安装gcc或libgomp解决。Debian/Ubuntu用apt install libgomp1,CentOS/Rocky用yum install libgomp。如果服务器权限不够,就别折腾系统环境,直接回到conda方案,把依赖一起装进去,省事得多。
还有一类问题是二进制平台不一致。下载Linux版本放到macOS上跑,或者反过来,都会报各种奇怪的错误。下载前先确认uname -a的系统类型。这一步虽然简单,但在集群环境里经常搞错。
3. VCF预处理:这步决定建树成败
3.1 过滤标准:质量、MAF、缺失率一个都不能少
VCF文件里的位点不是都能直接喂给IQtree。原始VCF里混有低质量位点、稀有变异、高缺失位点,甚至多等位位点。其中多等位位点尤其麻烦,IQtree的VCF模块在很多版本里会明确报错,要求输入双等位位点,所以第一步就是过滤。
我自己常用的bcftools过滤命令:
bcftools view --threads 4 \ -m2 -M2 -v snps \ -i 'QUAL>30 && MAF>0.05 && F_MISSING<0.2' \ input.vcf.gz -Oz -o clean.vcf.gz bcftools index -t clean.vcf.gz逐个解释这些参数:-m2 -M2表示只保留双等位位点,多等位位点会直接排除;-v snps只保留SNP,把INDEL丢掉,因为INDEL在进化模型里的处理方式和SNP不一样;QUAL>30过滤掉低质量的变异;MAF>0.05要求次等位基因频率大于5%,用来剔除极稀有变异;F_MISSING<0.2要求位点在样本中的缺失率低于20%。
注意:阈值不是铁律。如果你的样本量只有几十份,MAF>0.05会丢掉大量有效信息,可以放宽到0.02或者干脆不做MAF过滤。玉米这种大基因组、多亚群的材料,原始SNP动辄上百万,过滤后剩十万到几十万都很正常。
样本层面的缺失率也得看。我会用vcftools看每个样本的缺失情况:
vcftools --gzvcf clean.vcf.gz --missing-indv输出的out.imiss文件里F_MISSING过高的样本要及时标记,建树时可以考虑剔除。一个缺失率70%的样本挂在树上,不仅自己形成一条不真实的长枝,还会拉低周围分支的bootstrap支持率,这种问题在预处理阶段解决成本最低。
3.2 LD剪枝:避免个别基因组区域绑架整棵树
连锁不平衡是群体遗传分析里绕不开的问题。如果基因组上某个区域因为选择、参考基因组组装不完整或测序偏好贡献了大量连锁SNP,那么这些位点会在建树时被重复加权,导致树被这个区域“绑架”。玉米是大基因组,不同染色体区域的LD衰减速度差异很大,不做剪枝很容易让树反映的是某个染色体片段的历史,而不是全基因组水平的群体关系。
LD剪枝常用plink:
plink --vcf clean.vcf.gz --double-id --allow-extra-chr \ --indep-pairwise 50 10 0.2 --out clean plink --vcf clean.vcf.gz --double-id --allow-extra-chr \ --extract clean.prune.in --make-bed --out clean_ld解释一下关键参数:50 10 0.2的意思是窗口大小为50个SNP,每次滑动10个SNP,窗口内两两位点的r²超过0.2时剔除其中一个;--allow-extra-chr是因为玉米数据里除了1到10号染色体,还可能有大量scaffold,不加这个参数plink会直接报错;--double-id解决样本ID里有空格或FID/IID拆分的问题。
剪枝后的数据可以转回VCF,也可以直接从plink生成的clean.prune.in文件里提取保留位点。如果说清楚适用场景:如果你的目标是看全基因组水平的群体分化,LD剪枝后的数据更稳妥;如果你想尽量保留所有信息做精细的样本亲缘关系梳理,不剪枝也可以,但最好两个版本都跑一遍,对比结果。
3.3 输入格式选择:直接读VCF还是转PHYLIP
IQtree 2直接支持VCF输入,命令里用-vcf指定文件,这是最方便的路径。不过要注意几点:第一,IQtree会把VCF里的杂合基因型转成IUPAC简并碱基,相当于把二倍体的不确定性编码进序列里,但它并不真正解析单倍型,所以树实际反映的是样本在SNP位点上的综合相似度;第二,缺失基因型会被当作未知状态处理,缺失率过高会让树的支撑率崩掉;第三,VCF文件最好用bgzip压缩并用bcftools索引,防止读取大文件时出错。
如果你更习惯先转成PHYLIP或FASTA,我常用的工具是vcf2phylip.py:
python vcf2phylip.py --input clean_ld.vcf.gz --output-prefix maize转格式最大的好处是文件结构透明,出了问题更容易定位;坏处是多一道工序、多占一点磁盘。在几百份样本、几十万SNP的规模下,直接输入VCF和转成PHYLIP在运行时间上没有本质差别,按个人习惯选择即可。有一点我建议保持固定:无论走哪条路,都要用过滤剪枝后的VCF,不要把原始VCF直接灌进去。
4. IQtree建树核心参数与实战命令
4.1 最稳妥的启动命令
我的标准建树命令长这样,假设直接用过滤剪枝后的VCF:
iqtree2 -vcf clean_ld.vcf.gz \ -m MFP+ASC \ -B 1000 -alrt 1000 \ -T AUTO --prefix maize转成PHYLIP之后,把-vcf换成-s maize.phy即可,其余参数一样。逐项解释:-m MFP+ASC让ModelFinder自动选择含ASC校正的最佳模型;-B 1000跑1000次UFBoot2超快自助法;-alrt 1000同时做1000次SH-aLRT似然比检验;-T AUTO自动检测可用CPU;--prefix maize让所有输出文件都以maize开头,不会污染目录。
这个命令会把模型选择、最大似然建树、两种分支支持率一次跑完。对于几百份样本、几十万SNP的数据,通常几十分钟到几小时内可以完成。如果位点数量特别大,比如超过百万,跑之前先看一眼模型选择阶段的时间,因为模型选择往往比建树本身更耗时。
4.2 模型选择:MFP还是直接指定GTR+ASC
-m MFP是IQtree 2的招牌功能,全称ModelFinder,它会在输入数据上评估上百种替代模型,挑出最适合的一个。但从实际经验看,位点越多,模型选择越慢。对于50万甚至上百万SNP,直接让MFP完整跑完,可能在模型选择阶段就卡很久。
两种策略:数据量在十万位点以内,直接用-m MFP+ASC省心;数据量很大,先用-m GTR+ASC固定模型跑一版,或者抽样一部分位点跑MFP看最佳模型——通常结果都是GTR+ASC这类复杂模型——然后全量数据用这个固定模型重新跑。
关于ASC,我必须单独提醒一下。绝大多数从VCF构建系统发育树的场景,位点是只含有变异位点的SNP,不变位点根本不在数据里。这种情况下如果不用ASC校正,模型会低估不变位点的相对比例,结果就是枝长被系统性拉长、模型参数出现偏差,严重时影响拓扑结构。+ASC让模型在已知所有位点都可变的前提下重新计算似然,近似校正这种由于SNP筛选造成的偏差。我早期没有加ASC时,树长得离奇,尤其是根部枝长爆炸,加上+ASC之后整体才回归正常。你可以自己对比两个版本,效果非常直观。
4.3 分支支持率:UFBoot2和SH-aLRT怎么配合
IQtree跑完后,-B 1000的UFBoot2和-alrt 1000的SH-aLRT会同时出现在树文件和报告文件里。前者是基于自举重采样的频率支持度,后者是每个分支的似然比检验值。两者互补:SH-aLRT对单分支的替代信号更敏感,UFBoot2更看重整体重采样稳定性。
我自己的经验判断标准:SH-aLRT大于等于80且UFBoot大于等于95,可以算强支持;SH-aLRT大于等于70且UFBoot大于等于90,算中等支持;低于这个范围的分支,写论文时我会谨慎描述,不能硬撑。
但这套标准是在正常分化尺度下用的。玉米不同亚群之间差异大,支持度通常很高;亚群内部的个体之间SNP差异很小,很多短枝的支持率天然就低,就算把所有位点都用上也很难改变。这时候别硬调参数把所有分支都撑到99,那往往是过拟合。合理做法是承认分辨率边界,用PCA或亲缘关系矩阵辅助解读近缘个体之间的位置。
4.4 并行与内存控制的实际调优
IQtree默认会用满机器上所有核心,-T AUTO很方便,但在共享服务器上我推荐手动指定核心数,避免把节点资源全抢光:
iqtree2 -s maize.phy -m GTR+ASC -B 1000 -alrt 1000 -T 16 -mem 60G-mem限制最大内存使用量。几十万SNP、几百样本的情况下,60G通常够用。如果内存不够,先检查VCF是否压缩、是否包含大量无效位点,实在不行就减少样本或位点量。我试过在内存只有32G的机器上跑30万SNP,指定-mem 30G后勉强能跑,但速度会明显变慢,所以有条件还是优先保证CPU和内存的平衡。
还有一个小建议:建树之前先df -h看一眼磁盘剩余空间。IQtree在大数据量下会生成临时文件,如果写满磁盘,运行到一半报错或段错误,排查起来非常痛苦。建议留出至少数据文件体积十几倍的空间。
5. 结果解读与可视化:从输出文件到发表级图片
5.1 IQtree输出文件到底该怎么看
一次完整的IQtree运行会生成一堆文件,第一次跑的人容易看懵。我按重要性逐个说明:
maize.iqtree:主报告,纯文本。包含运行参数、最佳模型、碱基频率、有效位点数、树的分支支持率统计,分析前必看;maize.treefile:最大似然树,Newick格式,每个分支都带着UFBoot和SH-aLRT支持值;maize.contree:基于UFBoot2得到的共识树,如果想展示一个相对可信的拓扑,一般用这个文件;maize.model.gz:压缩后的模型参数,里面记录替换速率、均衡频率等;maize.log:运行日志,排查问题用。
打开.iqtree报告时,我最先看三处:开头的“Model of substitution”确认最终模型,中间的“Number of informative sites”确认有效位点数,树的部分看支持率分布。如果有效位点只有一两千个,那后面描述得再漂亮都是虚的,得回到过滤那一步放宽条件。
5.2 三种常用的可视化路线
拿到maize.contree之后,可视化方式很多。
最简单的是FigTree,图形界面操作,直接把Newick文件拖进去,在Node Labels里选择显示支持率,再导出PDF。适合快速看树形、检查样本分组。
如果想做文章级别的图,我推荐用R语言的ggtree:
library(ape) library(ggtree) tr <- read.tree("maize.contree") p <- ggtree(tr, layout="circular") + geom_tiplab(size=3) + geom_nodelab(aes(label=label), size=2) ggsave("maize_tree.pdf", p, width=10, height=10)如果你有样本对应的群体信息,可以映射到颜色上,把不同亚群标出来,比单纯黑白的树直观很多。还有在线工具iTOL,上传树文件后可以通过网页调整颜色、形状、标签,很多期刊里精致的进化树就是这么做的。
必须提醒一句:可视化解决的是展示问题,不是验证问题。无论图上多好看,结论都要回到.iqtree报告里的支持率和模型参数。如果支持率很低,图上颜色再丰富也不能改变结论的可靠性。
6. 避坑指南与问题排查
6.1 我踩过的五个真实坑
第一个坑是直接拿原始VCF跑IQtree,没过滤多等位位点。当时跑了三个多小时快结束时报错VCF has multi-allelic sites,整个人都很崩溃。后来把-m2 -M2放进预处理命令,再也没出过这个问题。
第二个坑是ASC校正。有一版树整体结构很奇怪,所有个体之间的枝长都明显偏长,根部距离大到不真实。群里的朋友提醒可能是SNP数据没有ASC校正,我加上+ASC重跑后树才正常。从那以后,我默认只要是VCF里的多态SNP建树,就至少跑一个含ASC的版本作对照。
第三个坑是高缺失样本。我有一份玉米数据,某个样本测序质量差,缺失率接近70%。跑出来的树里这个样本被挂在很长的孤立枝条上,接着这个长枝又拉低了周围几个分支的支持率。后来用vcftools查看了样本缺失率,删掉这个样本后整棵树都稳定了。现在我在预处理阶段一定先看样本缺失率,不让明显质量差的样本进树。
第四个坑是模型选择耗时失控。早期拿到40万SNP,直接-m MFP+ASC上,模型选择阶段跑了快三个小时还没结束。后来改成先随机抽样3万位点、固定GTR+ASC模型做一个快速预跑,确认拓扑合理之后再放全量数据跑,效率提升非常明显。
第五个坑是不做LD剪枝。当时用的SNP里有一段染色体区域因为参考基因组重复序列导致成百上千个连锁位点富集,树的亚群结构严重偏向这个区域。剪枝之后重新跑,树形才和已知系谱吻合。所以现在正式分析我都会同时跑剪枝和不剪枝两个版本,防止单版本误导结论。
6.2 常见报错速查表
| 现象 | 主要原因 | 处理建议 |
|---|---|---|
command not found: iqtree2 | PATH没有配置好 | 检查~/.bashrc或使用绝对路径 |
libgomp.so.1: cannot open shared object file | 缺少OpenMP运行库 | 装libgomp,或直接用conda环境 |
VCF has multi-allelic sites | 多等位位点未过滤 | bcftools view -m2 -M2 |
segmentation fault | 内存不足、文件损坏 | 检查内存和磁盘,用-mem限制,重新生成VCF |
| 模型选择速度过慢 | 位点太多 | 抽样跑MFP,或直接固定GTR+ASC |
| 几乎所有bootstrap都很低 | 近缘个体多或信息位点不足 | 增加位点、剔除高缺失样本,调整对短枝支持率的预期 |
| 树中某个样本长枝不自然 | 混合样本、样本鉴定错误或缺失过高 | 核查样本信息,必要时剔除问题样本 |
| 树拓扑与已知群体关系不符 | 未加ASC、未做LD剪枝或样本标签错乱 | 逐项排查,先跑剪枝加ASC的版本对照 |
最后再分享一个小技巧:我在这类项目里从不只跑一棵树。同一份VCF,我会同时生成一个“全位点+GTR+ASC”版本和一个“剪枝后+MFP+ASC”版本,对比哪些分支稳定、哪些分支对位点选择敏感。如果两者差异很大,说明数据里可能存在地区特异的信号或样本结构干扰,需要进一步排查。这种多重对照的习惯,比把宝押在一次运气好的运行上靠谱得多。