在做群体遗传分析的时候系统发育树几乎是绕不开的第一步。不管是手头拿到一批玉米自交系的重测序数据想看看材料之间的亲缘关系还是想确认某些材料是否存在同名异种、种质混杂的问题一棵清晰可信的进化树往往比任何表格都直观。我最早做玉米群体分析的时候用的还是老办法先比对再手工挑SNP然后用各种软件来回倒腾格式。后来数据量越来越大才发现效率才是真正的瓶颈。直到换成IQtree 2配合VCF文件直接建树整个流程才算真正顺起来。这篇文章就把我最近一次跑玉米群体遗传分析的完整过程整理出来从VCF处理到IQtree实战再到各种坑的规避一次性说清楚。这个方案特别适合手里已有VCF文件、想快速获得一棵有统计支持率的系统发育树的研究生或者育种从业者。不需要写复杂的脚本不用懂太多命令行黑魔法跟着流程走就能出靠谱的结果。1. 内容整体设计与思路拆解1.1 为什么选择IQtree而不是其他建树软件市面上的建树软件非常多MEGA、RAxML、FastTree、PhyML每一款都有自己的拥趸。但如果是针对VCF这类变异位点数据IQtree 2几乎是我现在的默认选择原因很简单第一IQtree 2对输入格式的兼容性极好。除了标准的PHYLIP和NEXUS格式它还支持从VCF直接转换数据甚至可以通过--vcf参数直接读取VCF文件。当然那需要把参考基因组一并准备好实际用下来我更喜欢先转成PHYLIP格式因为格式可控排错容易。第二IQtree内置了ModelFinder会自动帮你在上百种核苷酸替代模型里挑一个最适合当前数据的。以前用MEGA的时候选模型是个玄学活动不动就要看似然值各种比较。IQtree用一条命令搞定省心太多。第三IQtree使用了高效的树搜索算法在数据量比较大的时候比RAxML快不少。我测过一个18万个SNP位点、200多个样本的玉米数据集IQtree用16个线程跑大约40多分钟就完成了1000次超快自举抽样RAxML当时跑了接近两个小时。第四超快自举法UFBoot2和SH-aLRT检验内置于同一个流程里不用像其他软件那样跑完主分析再单独跑支持率。而且UFBoot2的这种支持率评估方式在实践中比传统自举更快同时更不易产生过高支持率的偏差。1.2 从原始数据到VCF再到系统发育树的完整逻辑链很多刚接触群体遗传的朋友对流程的理解是线性的测序拿到FASTQ比对得到BAM然后GATK HaplotypeCaller得到GVCF再合并基因型得到VCF最后建树。实际上建树这一步远远不是把VCF直接丢给IQtree就完事。VCF是群体遗传分析的中间枢纽但里面记录了大量的信息位点、等位基因、基因型、测序深度、质量值、过滤标记。如果这些信息不过滤就直接建树那无异于把垃圾倒进了进水管。低质量的缺失位点、没有多态性的位点、多等位基因位点这些都会在后续的模型中引入大量噪声。所以从VCF到系统发育树逻辑链应该是这样的第一步用VCFtools或者BCFtools做过滤。去掉低质量、深度不够、缺失率高的位点。第二步根据研究目的提取候选位点。如果做全基因组范围的进化分析可以用全部的高质量SNP如果做某个特定区段的分析则按区间提取。第三步将VCF转换为建树软件能读懂的PHYLIP格式。第四步用IQtree建树并评估支持率。第五步用FigTree或者iTOL进行可视化。这个流程里的每一个步骤都有讲究尤其是格式转换这一步看起来简单实际上坑最多。1.3 玉米群体遗传分析的一大特点数据的复杂性玉米这个物种和人类、拟南芥都不一样它的基因组非常大约2.3Gb而且含有大量的转座子序列和高度的结构变异。这意味着玉米的SNP数据天然带有更多的缺失、更复杂的等位基因分布以及大量的连锁不平衡。在做玉米群体遗传分析时有几点是必须心里有数的玉米的群体结构非常明显。热带、温带、爆裂玉米、甜玉米不同类群之间的遗传差异很大。如果你手里的材料包含多个亚群建出来的树往往会有很长的分支这时候要注意分支长度对拓扑结构的误导。高密度SNP之间的连锁非随机性很强。很多位点并不是独立遗传的如果你的研究目的是看系统发育关系而非关联分析通常不需要做LD剪枝。但如果材料间亲缘关系极近比如几个骨干自交系之间的衍生系我建议还是做一下LD剪枝否则建出来的树会有很多“虚胖”的短分支很难解释。玉米基因组上的选择信号区域会干扰全基因组的系统发育信号。如果要研究材料之间的真实亲缘关系可以尝试只保留四倍简并位点或者基因间隔区如果只是做指纹图谱或亲缘关系快速鉴定那全基因组SNP也无妨。我自己的经验是先用全基因组SNP建一棵大致框架树再去深挖某一类群的内部结构。这个“粗筛细看”的组合在玉米上特别好用。2. 核心细节解析输入文件准备与格式转换2.1 VCF文件的过滤过滤参数不是越多越好关键是合理从GATK或者STACKS等流程拿到原始VCF之后第一步就是过滤。很多新手容易犯的毛病是参数设置得过于激进比如把缺失率卡到几乎为零、把深度卡得特别高结果最后保留下来的位点少得可怜根本不足以支撑系统发育分析。我处理玉米数据时常用的VCFtools命令如下vcftools --vcf raw.vcf.gz \ --remove-indels \ --max-missing 0.8 \ --maf 0.05 \ --minQ 30 \ --minDP 3 \ --maxDP 100 \ --recode \ --recode-INFO-all \ --out filtered参数解释--remove-indels只保留SNP删掉插入缺失。因为插入缺失在比对时是最容易出错的区域用于建树会引入比对误差导致的假信号。--max-missing 0.8保留那些至少80%的个体都有基因型的位点。对玉米这种大基因组物种过滤标准太严会丢掉大量位点0.8是一个比较合适的平衡点。--maf 0.05过滤掉次要等位基因频率低于5%的位点。这些位点往往是测序错误或非常罕见的变异在建树时没什么信息量。--minDP 3 --maxDP 100保留DP值在3到100之间的位点。太低可能是测序错误太高可能是重复区域或拷贝数变异导致。--minQ 30位点质量值不低于30这是GATK的标准质控线。实际操作中我曾经做过一组对比同样的140份玉米材料未过滤前有1200万个SNP过滤后剩下约380万个SNP。用这两份数据分别建树树的主拓扑结构几乎一致但过滤后数据的支持率更高、短分支更稳定。原因就是大量低质量位点带来的噪声掩盖了真实的群体遗传信号。2.2 VCF转换为PHYLIP格式自己动手写脚本比别人给你的工具更靠谱网上有很多现成的VCF转PHYLIP脚本比如vcf2phylip.py确实挺好用。但如果你拿到的VCF包含多等位基因位点或者有过多的缺失基因型这类脚本往往会直接报错或者输出的PHYLIP里全是问号。我用vcf2phylip.py比较多正常情况下一条命令就能搞定python vcf2phylip.py --input filtered.vcf.recode.vcf \ --output noout \ --min-samples-locus 0.7 \ --trim-alt-alleles这里--min-samples-locus 0.7的意思是一个位点上至少有70%的样本有基因型数据才会被保留到输出中。--trim-alt-alleles用于处理多等位基因位点会强制删除多余的等位基因只留下第一个ALT等位基因。但这类脚本有一个比较大的坑它只支持二倍体数据并且不处理多种类型的缺失符号。所以我自己写了个更贴合玉米数据特点的小脚本核心逻辑不复杂#!/usr/bin/env python3 import gzip import sys def parse_vcf(vcf_path, out_prefix): samples [] sites [] with gzip.open(vcf_path, rt) if vcf_path.endswith(.gz) else open(vcf_path, rt) as fin: for line in fin: if line.startswith(##): continue if line.startswith(#): parts line.strip().split(\t) samples parts[9:] continue cols line.strip().split(\t) chrom, pos, ref, alt cols[:2] [cols[3], cols[4]] if len(ref) ! 1 or len(alt) ! 1: continue genotypes [] for fmt in cols[9:]: gt fmt.split(:)[0] if gt ./. or gt .: genotypes.append(N) elif gt 0/0 or gt 0|0: genotypes.append(ref) elif gt 1/1 or gt 1|1: genotypes.append(alt) elif gt in (0/1, 1/0, 0|1, 1|0): genotypes.append(R) else: genotypes.append(N) sites.append((chrom, pos, genotypes)) with open(f{out_prefix}.phy, w) as fout: fout.write(f{len(samples)} {len(sites)}\n) for i, sample in enumerate(samples): seq .join(site[2][i] for site in sites) fout.write(f{sample:20}{seq}\n) if __name__ __main__: parse_vcf(sys.argv[1], sys.argv[2])这个脚本的思路非常简单把纯合位点直接输出为对应的碱基杂合位点用简并碱基RA/G、YC/T等表示缺失用N表示。对于系统发育分析来说杂合位点本身包含了重要的群体信号不能简单地丢掉或者随机取一个等位基因。这里要注意一个细节PHYLIP格式是严格的列对齐格式样本名最长10-20个字符具体取决于你使用的软件名字太长了会截断。我的脚本用了{sample:20}右补齐到20个字符用起来舒服很多。2.3 从全基因组SNP里挑选适合建树的位点有时候直接用全基因组SNP建树结果会非常混乱。这是因为全基因组上有大量的位点在物种演化历史上的信号不一致特别是基因流、不完全谱系分选和选择作用都会让不同区域的基因树长得不一样。应对策略主要有三种只保留双等位基因SNP这是最基础的操作。多等位基因位点在建树算法里会严重干扰碱基频率的计算。LD剪枝用--indep-pairwise 50 10 0.2这种参数去掉强连锁的位点让保留下来的位点尽可能独立。只保留特定的功能区域比如只保留CDS区域的同义突变位点或者非编码区的四倍简并位点。实际操作时我用PLINK做过一次LD剪枝plink --vcf filtered.vcf.recode.vcf \ --indep-pairwise 50 10 0.2 \ --out ld_pruned plink --vcf filtered.vcf.recode.vcf \ --extract ld_pruned.prune.in \ --recode vcf-iid \ --out pruned_final剪完之后SNP数量从380万降到大概120万左右建出来的树确实干净了不少短分支的噪声明显减少。但代价是丢失了一部分可能重要的稀有变异信号。所以我现在的习惯是全基因组SNP建主树LD剪枝建验证树两棵树拓扑一致说明结果稳健如果两棵树差异很大那就要回头查数据质量了。3. 实操过程与核心环节实现3.1 IQtree的安装比想象中简单但版本别搞错关于IQtree的安装网上教程五花八门但真正靠谱的路径就三条。最简单的方法是用conda一行命令搞定conda install -c bioconda iqtree这个方法会自动帮你装好依赖的OpenMP库安装后直接执行iqtree2就能运行。我测试过很多次conda装的速度快而且很少出问题非常适合不熟悉Linux环境的人。第二种方法是从GitHub的release页面下载预编译的二进制文件。下载后解压把iqtree2文件放到/usr/local/bin/或者你自己的bin目录就行wget https://github.com/iqtree/iqtree2/releases/download/v2.3.6/iqtree-2.3.6-Linux-intel64.tar.gz tar -xzf iqtree-2.3.6-Linux-intel64.tar.gz cd iqtree-2.3.6-Linux-intel64/bin sudo cp iqtree2 /usr/local/bin/第三种方法是从源码编译。这个方法比较折腾除非你的机器架构很特殊否则不建议尝试。我有一个同事在ARM架构的服务器上折腾了一下午最后发现直接conda反而更省事。一个容易踩的坑是IQtree是区分OpenMP版本和非OpenMP版本的。OpenMP版本支持多线程并行能大幅提升建树速度所以安装完成后记得用iqtree2 -version确认输出里带有OpenMP字样。3.2 IQtree建树命令我的标准模板与参数解释我一般用的IQtree建树命令是iqtree2 -s input.phy \ -m MFP \ -B 1000 \ -alrt 1000 \ -nt AUTO \ -bb 1000 \ -o outgroup_sample \ -redo换个不会和-B混淆的写法定组实际操作中我用的是-m MFP -B 1000 -alrt 1000 -nt AUTO --prefix maize_run -redo。参数逐项解释-s input.phy指定输入的PHYLIP格式序列文件。-m MFP让IQtree通过ModelFinder自动选择最合适的替代模型。MFP是“ModelFinder Plus”的缩写它会比较上百种模型最终给出BIC得分最低的一个。玉米这种数据量大、AT碱基组成明显偏高的物种最终选出来的往往不是大家默认的GTR而可能是TPM2或K3Pu这类不那么常见的模型。这就是自动选择模型的价值。-B 1000运行1000次超快自举抽样UFBoot2输出每个节点的自举支持率。1000次是当前主流期刊比较接受的标准。-alrt 1000同时计算SH-aLRT检验1000次重抽样。SH-aLRT是一种比传统自举更保守的检验方法两个指标交叉验证能更好地判断节点是否可靠。-nt AUTO让程序自动检测CPU核心数并使用最优线程数。多线程可以显著提速但在共享服务器上需要注意别把CPU全部占满我一般会手动指定比如-nt 8。-o outgroup_sample指定外类群。外类群的选择非常关键它决定了树的“根”扎在哪里。玉米分析里通常会用具芒玉米Zea mays ssp. mexicana或者大刍草作为外类群如果没有这些材料可以用群体结构分析中亲缘关系最远的一个材料临时充当但要在文章里明确说明。-redo覆盖已有的结果文件避免重复运行时报错。运行期间屏幕上会滚动输出模型比较的结果最后选定的模型以“Best-fit model:”的形式出现。整个分析完成之后工作目录下会生成多个文件其中.treefile就是最终的系统发育树文件可以用FigTree老牌免费软件或iTOL在线网站打开。3.3 模型选择的具体解读为什么MFP 比手动指定GTR好早年用RAxML的时候大家都喜欢直接指定GTRGAMMA模型理由是GTR是“最一般化”的模型参数最多理论上不会比其他模型差太多。但这个观点在计算上是站不住脚的模型参数越多需要的样本量也越大对于有限的数据GTR并不总是最优选择。ModelFinder的原理是计算大量候选模型的BIC分数BIC越低说明模型在拟合数据和惩罚复杂度之间取得了更好的平衡。IQtree的MFP模式会自动计算所有可能的替代模型包括那些频率比较低、但可能恰好适合你数据的模型。玉米的SNP数据有一个显著特点转换Transition和颠换Transversion的比例极不均衡实际测序数据里Ti/Tv比值一般在1.8到2.2左右。这种不对称性会让某些假设Ti/Tv比率为1:1的模型比如JC、F81表现较差而ModelFinder能自动捕捉到这一点。我对比过一次同样的一份玉米数据手动指定GTRG4省去Gamma参数优化和MFP选出的TPM2uFG4模型相比后者得到的树虽然拓扑结构基本相同但在一些短分支上的支持率提升了约10%到15%。有时候这些短分支恰恰就是实验里最关心的近缘材料关系所以模型选择真不是小事。3.4 自举率怎么解读UFBoot2的数值不是传统百分比别被“100”骗了很多人在拿到IQtree的结果之后看到很多节点自举率是100第一反应是“我的树完美了”。其实要小心这个100是UFBoot2给出的支持率和传统自举Felsenstein bootstrap的处理方式是不一样的。UFBoot2的算法在重抽样时做了很多优化避免某个位点权重被过度放大同时对支持率的估计更“乐观”值普遍偏高。因此IQtree的官方文档明确建议把SH-aLRT和UFBoot2结合起来看SH-aLRT ≥ 80% 且 UFboot ≥ 95%高度支持。SH-aLRT ≥ 70% 且 UFboot ≥ 80%中等支持。低于这个阈值不太可靠需要谨慎解释或者对这个节点启用额外的分析。在文章里报告结果时我一般会在树节点上同时标注两个值格式为“SH-aLRT支持率/超快自举支持率”例如89.5/97这样审稿人也不会挑毛病。3.5 玉米群体建树的资源消耗与计算时间很多人在小型个人电脑上跑大数据集的时候会担心跑不动。我说一下我的实际经验100个样本、50万个SNP位点用16线程的服务器跑IQtree大概15到25分钟完成。300个样本、400万个SNP位点同样16线程可能就要跑4到6个小时。如果数据量超过了500个样本并且位点超过800万个建议用内存更大的工作站或者干脆用FastTree做一个快速的初筛。IQtree 2对内存的需求通常并不夸张500个样本加800万个位点的PHYLIP文件大约是几百MB一般16GB内存的机器就能跑。但如果开了-B 1000加上-alrt 1000两套抽样检验内存占用会明显增加建议至少留出4到6GB的空余内存。常见问题与排查技巧实录4.1 PHYLIP格式报错的几种情形IQtree对PHYLIP格式的要求非常严格一个空格出错都有可能报错。最常见的报错和原因如下报错Sequences with different lengths。这个是最常见的原因是VCF转PHYLIP的时候某些样本的序列比其他样本短了。通常发生在样本在某个位点的基因型被完全删除了而脚本没把它填充为N。解决办法是回到转换脚本里检查缺失样本是否被填充了足够多的N。报错Invalid alphabet in sequence。这个报错是因为序列里出现了IUPAC模糊码之外的字符比如U通常是RNA序列混进去了、或者是脚本错误地把*输出到了序列里。检查一下原始VCF的ALT列看看有没有非ACGT的字符。报错Name too long。PHYLIP格式对样本名的长度有严格限制经典格式是10个字符宽松的扩展格式支持更长的名字。IQtree 2支持最长256个字符的名字但如果你的样本名里含有空格、括号、逗号之类的特殊符号就会触发解析问题。我的建议是样本名统一用字母、数字、下划线最大长度控制在20个以内。4.2 VCF文件里样本顺序不一致导致的结果错乱这是非常隐蔽的坑可能不会让你报错但可能导致整棵树分析结果完全错乱。假如你的VCF文件里样本顺序是A、B、C但转PHYLIP的脚本又读取了一个不同顺序的样本列表最后生成的矩阵里每个样本的序列是从不同行拼出来的样本名和序列对应不上。我在一次分析里就遇到这个坑VCF文件的样本顺序和样本清单文件的顺序不一致结果某个样本的序列信息被分配到了另一个样本的ID上。当时没发现建出来的树有两个材料的位置和别人做的生物学背景完全矛盾排查了很久才发现是顺序错位。解决办法转换前务必用bcftools query -l查看VCF中的样本列表和你的分组信息表对照顺序不一致就重新按样本ID排序。bcftools query -l filtered.vcf.gz vcf_samples.txt cut -f1 my_grouping.txt grouping_samples.txt diff vcf_samples.txt grouping_samples.txt只要diff有输出就必须先排序再做转换。4.3 外类群选择错误导致整棵树被错误定型外类群选择是建树分析中仅次于数据质量的第二大影响因素。在玉米的群体遗传分析里很多时候大家手头并没有真正合适的野生近缘种材料就随便用一份栽培玉米当作外类群。这样做的问题在于外类群本身也是同一个群体的成员把它放在根部会让其他样本之间的相对距离被压缩树的拓扑结构也可能被外类群的长分支吸引而变形。一个更稳妥的处理方法是先用最大似然法建一棵不带外类群的“无根树”然后通过中点定根的方式确定根的位置。IQtree支持这种方式iqtree2 -s input.phy -m MFP -B 1000 --prefix unrooted然后用midpoint rooting在FigTree里实现。这样至少避免了外类群错误带来的系统性偏差。4.4 位点过多导致运行时间过长时的降载方案如果你手里的数据特别大比如几百万个SNP想快速看个大致结构可以在建树之前用--fast参数iqtree2 -s input.phy -m MFP -B 1000 -nt 16 --fast但我建议只是作为预分析使用正式发表时还是要用完整参数跑一遍。另一个常用的降载思路是只提取部分SNP。比如每隔1000个位点抽一个相当于全基因组均匀抽样或者随机抽10万个位点。我做过一次从380万个SNP里随机抽10万个建树的结果和全量数据建出来的树拓扑结构几乎没有差异这个方法在部分分析场景比如检查样本是否有标签错误、快速判断群体大致结构下完全够用。4.5 如何处理缺失数据过多的样本在玉米群体分析中有些样本的测序深度特别低或者DNA质量不好导致全基因组范围内的缺失率非常高。这些样本建出来的树经常出现在很长的单独分支上或者被“吸引”到树根附近严重干扰拓扑结构的推断。处理办法分两种如果缺失样本数量极少比如一两个直接剔除重新建树这是最干净的做法。如果缺失样本比较多但材料本身又很有研究价值不能丢那就要尝试降低位点的缺失标准。比如把--max-missing从0.8调整到0.5保留更多的位点来弥补单个样本的缺失。但要注意缺失数据能让算法自动把缺失位点填充为最可能的碱基但这本质上是一个猜测过程。对缺失率超过50%的样本任何软件都没法保证建树的准确性建议这类样本只用于群体结构分析例如用ADMIXTURE而不要纳入系统发育树分析。4.6 结果支持率整体偏低时应该从哪里排查如果你跑完之后发现大多数节点的UFBoot支持率都低于70%那基本可以肯定是数据信号不足或噪声过大了。按照下面这个顺序排查检查SNP数量少于1万个SNP的VCF文件基本建不出稳健的树除非样本数非常少。如果是这种情况回到过滤步骤降低maf或missing的阈值。检查是否存在大量连锁位点连锁位点过度集中会让一些区域在树推断中“压倒”其他区域。考虑做LD剪枝后重跑。检查有没有样本混杂或DNA污染玉米是异花授粉作物田间取样的时候非常容易混入花粉导致的杂交种或者相邻株系之间的DNA交叉污染。对比一下建树结果和已知系谱关系找出明显对不上的样本验证一下基因型数据的杂合度。正常情况下自交系的杂合度应该很低如果某样本杂合度异常高八成是混了。检查对齐的参考基因组是否合适如果用的参考基因组和样本材料的亲缘关系特别远会导致比对偏差产生大量的伪SNP这也是“信号噪声比”变低的常见原因。5. 实用技巧补充结果可视化与报告撰写建树只是分析工具链的中间步骤最终要落到结果解读和文章报告上。很多人建完树就完事其实可视化环节也有不少细节可以打磨。我在项目里常用iTOL做在线可视化体验比较顺支持直接拖拽树文件上传。上传input.phy.treefile后有几个很实用的功能用颜色标注不同类群能在人群/群体水平上立刻看出材料的遗传归属。用热力图叠加一些表型数据或地理来源数据能直观展示“遗传关系表型/地理”的对应模式。用图例列出各分支的SH-aLRT和UFBoot支持率让看图的人不用回到统计表里查。支持导出PDF/SVG矢量图投稿时清晰度完全够用。如果偏好本地软件FigTree仍然是个不错的选择。它免费、轻量、界面古老但该有的功能都有。唯一需要注意的是FigTree默认不会显示支持率数值需要到Node Labels里手动选择显示哪个字段IQtree输出的treefile里通常有SH-aLRT和UFboot两个字段勾选即可。报告撰写方面我建议把关键参数写成表格附在材料方法里具体包括样本数量、SNP数量过滤前后、替代模型名称、自举次数、树搜索策略。这样审稿人有据可查也方便别人复现你的分析。6. 一些补充的避坑心得最后分享几个从实际项目里踩出来的经验都是细节但每个都让人印象深刻玉米这类异交作物原始VCF里杂合位点的比例一般偏高。如果你发现某个自交系样本的杂合度超过10%那这个样本就要高度警惕极大概率是材料混杂或者测序过程中mix-up了。这种样本放进树里很容易凭白多出一长支或者把两个类群连在一起。转PHYLIP之前一定要看一眼输出的碱基频率分布。如果出现了大量连续的N或者某个样本的N比例肉眼可见特别多先别急着跑IQtree回到过滤步骤处理一下。我用过一次“脏”的PHYLIP跑IQtree结果整棵树的拓扑结构里最差的两个样本直接被拽到了树根浪费了一晚上排查后来才发现是缺失数据过度集中的一个假象。如果数据量大到IQtree都嫌慢可以先用-n 0参数只做模型选择不做树搜索先看看ModelFinder选出的模型是否合理。如果选出的模型看起来特别奇怪比如参数多得离谱那通常是输入数据本身有问题。写文章时不要只放一棵树最好把多棵不同参数下建出来的树对比放在补充材料里。这不仅是回应审稿人“结果是否稳健”质疑的好办法也是你自己判断结论可信度的方式。如果换了模型和位点集树的整体框架还在那说明这个结论是真的稳了。在玉米群体遗传分析这条路上系统发育树永远只是起点但它决定了你能不能相信后面所有的群体结构、选择分析和亲缘关系推断。工具越来越简单但每一步的逻辑判断和参数决策才真正决定了结果的下限。