玉米群体进化树构建实战从577个品系1万SNP到有根系统发育树玉米作为全球最重要的粮食作物之一其遗传多样性研究对育种和进化生物学具有重要意义。本文将详细介绍如何利用IQtree软件基于577个玉米品系的1万个SNP标记数据构建有根系统发育树的全流程。这个实战案例不仅适用于玉米研究也可为其他作物的群体遗传分析提供参考。1. 数据准备与预处理构建高质量系统发育树的第一步是确保输入数据的可靠性。我们的数据集包含577个玉米品系涵盖70个大刍草(parviglumis)、208个热带品系(tropical)、172个温带品系(temperate)和127个混合品系(mixed)。这些品系代表了玉米的广泛遗传多样性。原始数据来自3号染色体上的743万个SNP位点经过以下过滤步骤质量过滤最小等位基因频率(MAF) 0缺失率(missing rate) 0SNP抽样# 使用vcftools进行随机抽样 vcftools --vcf input.vcf --snps 10000 --recode --out sampled_10k提示SNP抽样数量需要平衡计算资源和信息量。1万个高质量SNP通常足以反映群体结构同时保持计算可行性。经过处理后我们获得了一个577×10000的SNP矩阵保存为VCF格式。为了适应IQtree的输入要求需要将其转换为PHYLIP格式# 使用vcftools转换为phylip格式 vcftools --vcf sampled_10k.vcf --phylip --out SNP.varsites2. IQtree模型选择与优化IQtree提供了丰富的进化模型选择功能。对于SNP数据特别是当我们需要构建有根树时模型选择尤为关键。2.1 Lie Markov模型遍历IQtree从v2.1.3开始支持使用Lie Markov模型进行无外群生根。我们首先遍历所有Lie Markov模型以寻找最适合我们数据的模型iqtree -s SNP.varsites.phy -mset liemarkov -nt 100 -fast --prefix SNP.varsites.liemarkov.fast这个命令会测试所有Lie Markov模型并记录每个模型的log-likelihood值。在我们的案例中发现了几个关键现象当遍历到6.7a模型时会出现警告WARNING: Numerical underflow for non-rev lh-branch Noname WARNING: evec not invertible遍历到9.2b模型时程序会报错终止2.2 最优模型确定经过比较6.7aASCR5模型表现出最高的log-likelihood值(-78103)因此被选为最终建树模型模型Log-likelihood稳定性6.7a-78103警告6.7aASC-78102警告6.7aR5-78101警告6.7aASCR5-78100稳定注意当输入数据不包含constant sites时IQtree在遍历某些Lie Markov模型时容易报错。这是当前版本的一个已知问题。3. 有根系统发育树构建确定了最优模型后我们可以正式构建有根系统发育树iqtree -s SNP.varsites.phy -m 6.7aR5ASC -nt 100 -fast --prefix SNP.varsites.6.7_R5_ASC.fast这个命令使用了100个线程(-nt 100)来加速计算。构建完成后我们会得到几个重要文件.treefileNewick格式的树文件.log运行日志文件.iqtree详细的分析报告3.1 结果可视化使用iTOL等工具对结果树进行可视化时我们可以按品系类型进行颜色标注红色大刍草品系蓝色热带品系黑色温带品系绿色混合品系这种可视化有助于直观理解不同玉米群体的进化关系。4. 实战中的挑战与解决方案在实际操作中我们遇到了几个关键挑战并总结了相应的解决方案4.1 IQtree官方推荐流程的问题IQtree官方文档建议采用两步法构建有根树先用常规模型构建无根树iqtree2 -s SNP.varsites.phy -p SNP.varsites.nex -T AUTO --prefix rev_snp再用Lie Markov模型合并并生根iqtree2 -s SNP.varsites.phy -p rev_snp.best_scheme.nex --model-joint 12.12 -T AUTO --prefix nonrev_snp然而我们发现这种方法在处理SNP数据时经常报错特别是第二步的合并过程。4.2 替代方案评估我们测试了另一种策略先构建无根树再以其为模板进行生根# 第一步构建无根树 iqtree -s SNP.varsites.phy -m GTRFR6ASC -fast -nt 100 --prefix rev_snp # 第二步以无根树为模板生根 iqtree -s SNP.varsites.phy -m 6.7aASCR5 -te rev_snp.treefile -fast -nt 100 --prefix nonrev_snp实验结果发现IQtree不会修改-te参数提供的树的拓扑结构只会调整分枝长度。三次重复实验的结果比较重复实验Log-likelihood总树长内部分枝长度比例1-774411.13619.94%2-774421.07719.94%3-773971.24819.91%这表明该方法虽然稳定但并不能真正优化树的拓扑结构。4.3 功能兼容性问题我们发现Lie Markov模型与某些IQtree功能存在兼容性问题--model-joint参数iqtree -s SNP.varsites.phy -p rev_snp.best_scheme.nex --model-joint 6.7aASCR5 -nt 100 -fast --prefix nonrev_snp经常导致程序崩溃似然映射(--lmap)iqtree -s SNP.varsites.phy -m 6.7aASCR5 -lmap 10000 -fast -nt 100 --prefix nonrev_snp.lmap1w无法正常工作根检验(--root-test)iqtree -s SNP.varsites.phy -m 6.7aASCR5 --root-test -zb 1000 -te SNP.varsites.6.7_R5_ASC.fast.treefile -fast -nt 100 --prefix nonrev_snp.root_test会产生意外结果5. 生物学意义解读与验证构建系统发育树后我们需要从生物学角度解读结果。在我们的案例中树形结构显示大刍草品系(红色)位于树的基部支持玉米驯化于大刍草的理论热带品系(蓝色)和温带品系(黑色)形成明显分支反映地理隔离导致的遗传分化混合品系(绿色)分布在多个位置表明存在广泛的基因交流为了验证结果的可靠性我们建议重复实验使用不同的SNP抽样集重建树检查拓扑结构的稳定性参数敏感性分析测试不同模型参数对结果的影响与已知标记比较将系统发育树与已知的功能标记或表型数据进行对比在577个品系×1万SNP的数据规模下整个分析流程在100线程服务器上约需6-8小时完成。主要的计算瓶颈在于模型测试阶段特别是遍历Lie Markov模型时。
给玉米群体进化树‘安家’：577个品系1万SNP的IQtree有根树实战
玉米群体进化树构建实战从577个品系1万SNP到有根系统发育树玉米作为全球最重要的粮食作物之一其遗传多样性研究对育种和进化生物学具有重要意义。本文将详细介绍如何利用IQtree软件基于577个玉米品系的1万个SNP标记数据构建有根系统发育树的全流程。这个实战案例不仅适用于玉米研究也可为其他作物的群体遗传分析提供参考。1. 数据准备与预处理构建高质量系统发育树的第一步是确保输入数据的可靠性。我们的数据集包含577个玉米品系涵盖70个大刍草(parviglumis)、208个热带品系(tropical)、172个温带品系(temperate)和127个混合品系(mixed)。这些品系代表了玉米的广泛遗传多样性。原始数据来自3号染色体上的743万个SNP位点经过以下过滤步骤质量过滤最小等位基因频率(MAF) 0缺失率(missing rate) 0SNP抽样# 使用vcftools进行随机抽样 vcftools --vcf input.vcf --snps 10000 --recode --out sampled_10k提示SNP抽样数量需要平衡计算资源和信息量。1万个高质量SNP通常足以反映群体结构同时保持计算可行性。经过处理后我们获得了一个577×10000的SNP矩阵保存为VCF格式。为了适应IQtree的输入要求需要将其转换为PHYLIP格式# 使用vcftools转换为phylip格式 vcftools --vcf sampled_10k.vcf --phylip --out SNP.varsites2. IQtree模型选择与优化IQtree提供了丰富的进化模型选择功能。对于SNP数据特别是当我们需要构建有根树时模型选择尤为关键。2.1 Lie Markov模型遍历IQtree从v2.1.3开始支持使用Lie Markov模型进行无外群生根。我们首先遍历所有Lie Markov模型以寻找最适合我们数据的模型iqtree -s SNP.varsites.phy -mset liemarkov -nt 100 -fast --prefix SNP.varsites.liemarkov.fast这个命令会测试所有Lie Markov模型并记录每个模型的log-likelihood值。在我们的案例中发现了几个关键现象当遍历到6.7a模型时会出现警告WARNING: Numerical underflow for non-rev lh-branch Noname WARNING: evec not invertible遍历到9.2b模型时程序会报错终止2.2 最优模型确定经过比较6.7aASCR5模型表现出最高的log-likelihood值(-78103)因此被选为最终建树模型模型Log-likelihood稳定性6.7a-78103警告6.7aASC-78102警告6.7aR5-78101警告6.7aASCR5-78100稳定注意当输入数据不包含constant sites时IQtree在遍历某些Lie Markov模型时容易报错。这是当前版本的一个已知问题。3. 有根系统发育树构建确定了最优模型后我们可以正式构建有根系统发育树iqtree -s SNP.varsites.phy -m 6.7aR5ASC -nt 100 -fast --prefix SNP.varsites.6.7_R5_ASC.fast这个命令使用了100个线程(-nt 100)来加速计算。构建完成后我们会得到几个重要文件.treefileNewick格式的树文件.log运行日志文件.iqtree详细的分析报告3.1 结果可视化使用iTOL等工具对结果树进行可视化时我们可以按品系类型进行颜色标注红色大刍草品系蓝色热带品系黑色温带品系绿色混合品系这种可视化有助于直观理解不同玉米群体的进化关系。4. 实战中的挑战与解决方案在实际操作中我们遇到了几个关键挑战并总结了相应的解决方案4.1 IQtree官方推荐流程的问题IQtree官方文档建议采用两步法构建有根树先用常规模型构建无根树iqtree2 -s SNP.varsites.phy -p SNP.varsites.nex -T AUTO --prefix rev_snp再用Lie Markov模型合并并生根iqtree2 -s SNP.varsites.phy -p rev_snp.best_scheme.nex --model-joint 12.12 -T AUTO --prefix nonrev_snp然而我们发现这种方法在处理SNP数据时经常报错特别是第二步的合并过程。4.2 替代方案评估我们测试了另一种策略先构建无根树再以其为模板进行生根# 第一步构建无根树 iqtree -s SNP.varsites.phy -m GTRFR6ASC -fast -nt 100 --prefix rev_snp # 第二步以无根树为模板生根 iqtree -s SNP.varsites.phy -m 6.7aASCR5 -te rev_snp.treefile -fast -nt 100 --prefix nonrev_snp实验结果发现IQtree不会修改-te参数提供的树的拓扑结构只会调整分枝长度。三次重复实验的结果比较重复实验Log-likelihood总树长内部分枝长度比例1-774411.13619.94%2-774421.07719.94%3-773971.24819.91%这表明该方法虽然稳定但并不能真正优化树的拓扑结构。4.3 功能兼容性问题我们发现Lie Markov模型与某些IQtree功能存在兼容性问题--model-joint参数iqtree -s SNP.varsites.phy -p rev_snp.best_scheme.nex --model-joint 6.7aASCR5 -nt 100 -fast --prefix nonrev_snp经常导致程序崩溃似然映射(--lmap)iqtree -s SNP.varsites.phy -m 6.7aASCR5 -lmap 10000 -fast -nt 100 --prefix nonrev_snp.lmap1w无法正常工作根检验(--root-test)iqtree -s SNP.varsites.phy -m 6.7aASCR5 --root-test -zb 1000 -te SNP.varsites.6.7_R5_ASC.fast.treefile -fast -nt 100 --prefix nonrev_snp.root_test会产生意外结果5. 生物学意义解读与验证构建系统发育树后我们需要从生物学角度解读结果。在我们的案例中树形结构显示大刍草品系(红色)位于树的基部支持玉米驯化于大刍草的理论热带品系(蓝色)和温带品系(黑色)形成明显分支反映地理隔离导致的遗传分化混合品系(绿色)分布在多个位置表明存在广泛的基因交流为了验证结果的可靠性我们建议重复实验使用不同的SNP抽样集重建树检查拓扑结构的稳定性参数敏感性分析测试不同模型参数对结果的影响与已知标记比较将系统发育树与已知的功能标记或表型数据进行对比在577个品系×1万SNP的数据规模下整个分析流程在100线程服务器上约需6-8小时完成。主要的计算瓶颈在于模型测试阶段特别是遍历Lie Markov模型时。
