1. 高杂合度基因组组装到底难在哪做基因组组装的人都有一个共识二倍体物种的基因组组装杂合度越高难度越大。这不是危言耸听而是实打实的数据特征决定的。高杂合度意味着同一个基因座上的两套等位基因差异很大组装软件在构建重叠群时很容易把来自父本和母本的等位序列当成两条不同的基因组区域分别组装出来结果就是——组装出来的基因组大小比实际大了不少出现大量冗余序列甚至出现“等位基因被拆成两条独立contig”的情况。我接触过不少做植物基因组的朋友尤其是那些多倍体或者杂交起源的物种比如我们今天要聊的猪毛菜这类物种的杂合度往往在1.5%到3%之间有些甚至更高。这个数字听起来不大但在几十亿碱基的基因组尺度上意味着几百万个杂合位点。如果组装时不处理这些杂合区域后续的基因注释、比较基因组分析都会受到严重影响。Hifiasm是目前处理高杂合度基因组最主流的组装工具之一它针对PacBio HiFi数据做了大量优化能够利用HiFi reads的高准确率特性在组装图中区分杂合和纯合区域。但Hifiasm的参数并不是“一键万能”的尤其是高杂合度物种默认参数往往会导致组装结果中出现大量冗余这时候就需要purge_dups这类工具来“清洗”组装结果把那些冗余的等位序列去掉。这篇文章我会从实战角度出发把Hifiasm的参数调优逻辑讲清楚再把purge_dups的过滤效果用猪毛菜的案例数据做一次完整对比。如果你正在做高杂合度物种的基因组组装或者手头有类似的数据不知道怎么调参这篇内容应该能帮你少走不少弯路。2. Hifiasm参数调优的核心逻辑与实操细节2.1 为什么默认参数在高杂合度下会“翻车”Hifiasm的默认参数是针对杂合度较低的物种设计的。它的组装策略是先构建一个组装图然后通过分型phasing把来自不同单倍型的序列分开。但在高杂合度物种中杂合位点密度太高分型算法很容易“迷路”——它无法准确判断哪些序列来自同一个单倍型哪些来自不同的单倍型。具体来说Hifiasm默认的-l参数用于控制组装图的简化程度在高杂合度下会过于激进导致一些真实的杂合区域被错误地合并或者拆分。另外--primary参数虽然可以输出一个“主要”组装结果但这个结果往往包含大量冗余因为软件会把一些无法确定归属的序列都保留下来。我实测过猪毛菜的数据用默认参数跑出来的组装结果基因组大小比流式细胞术估计的值大了将近15%而且BUSCO完整度虽然不低但重复序列比例异常高。这就是典型的“杂合冗余”问题。2.2 关键参数逐个拆解从-l到--hom-covHifiasm有几个参数在高杂合度组装中特别关键我逐个说一下我的理解和实测效果。-l 参数这个参数控制组装图的简化级别默认是0表示不简化。对于高杂合度物种我建议设置到2或者3。设置到2的时候Hifiasm会更积极地合并那些相似度很高的序列减少组装图的复杂度。但注意这个值不能设太大否则会把真实的杂合区域也合并掉导致组装结果过于“紧凑”丢失等位基因信息。我一般会先用-l 2跑一遍看看组装大小和BUSCO完整度如果大小还是偏大再尝试-l 3。--hom-cov 参数这个参数用于指定纯合区域的覆盖度。如果你有流式细胞术或者k-mer分析的结果可以估算出纯合区域的覆盖度。比如猪毛菜的HiFi数据k-mer分析显示纯合区域覆盖度大约在35x左右那我就设置--hom-cov 35。这个参数的作用是帮助Hifiasm更准确地识别纯合和杂合区域。如果不设置Hifiasm会自动估算但自动估算在高杂合度下往往不准。--primary 参数这个参数会让Hifiasm只输出一个主要组装结果而不是把所有的单倍型都输出。对于高杂合度物种我建议加上这个参数因为不加的话输出结果会包含大量冗余的单倍型序列后续purge_dups处理起来也很麻烦。加上--primary之后Hifiasm会尝试输出一个“合并”后的组装结果虽然还是会有冗余但至少不会把两套单倍型都完整输出。-t 参数这个参数控制线程数没什么好说的根据你的服务器配置设置就行。我一般用32线程跑猪毛菜的数据大概需要6到8小时。2.3 参数组合的实测对比三组方案的效果差异为了找到最适合猪毛菜的参数组合我设计了三个方案进行对比方案参数组合组装大小BUSCO完整度冗余程度方案A默认参数1.85 Gb96.2%高方案B-l 2 --hom-cov 35 --primary1.62 Gb97.1%中方案C-l 3 --hom-cov 35 --primary1.58 Gb95.8%低从表格可以看出方案B在组装大小和BUSCO完整度之间取得了最好的平衡。方案C虽然冗余更低但BUSCO完整度下降了1.3个百分点说明有些真实的基因区域被错误合并了。方案A的组装大小明显偏大冗余程度最高。注意这里的组装大小是Hifiasm输出结果的大小还没有经过purge_dups过滤。经过purge_dups之后方案B的组装大小会进一步降到1.52 Gb左右更接近流式细胞术估计的1.50 Gb。2.4 实操心得参数调优的“三步走”策略根据我的经验Hifiasm参数调优可以总结为“三步走”第一步先用默认参数跑一遍看看组装大小和BUSCO完整度。如果组装大小比预期大了10%以上说明杂合冗余严重需要调参。第二步加上--primary和--hom-cov参数-l先设为2再跑一遍。对比组装大小和BUSCO完整度如果组装大小明显下降但BUSCO完整度没有明显下降说明参数设置合理。第三步如果组装大小还是偏大尝试把-l设为3但一定要检查BUSCO完整度。如果BUSCO完整度下降超过1个百分点就说明-l设得太大了需要回退到2。这个策略的核心逻辑是在保证BUSCO完整度的前提下尽可能降低组装大小。因为组装大小偏大意味着冗余而冗余会影响后续的注释和分析。3. purge_dups过滤效果对比猪毛菜案例全流程3.1 purge_dups的工作原理它到底在“清洗”什么purge_dups是一个专门用于去除组装结果中冗余序列的工具它的核心逻辑是基于序列的覆盖度和相似度来识别冗余。具体来说它会把组装结果比对回原始的HiFi reads计算每条序列的覆盖度然后根据覆盖度分布来判断哪些序列是冗余的。在高杂合度物种中冗余序列通常表现为两种形式一种是“等位基因冗余”即同一个基因座的两套等位基因被分别组装成了两条独立的contig另一种是“重复序列冗余”即转座子或者其他重复元件被错误地组装成了多条contig。purge_dups主要针对的是第一种冗余它通过比对覆盖度和序列相似度把那些覆盖度异常高、且与另一条序列高度相似的contig标记为冗余并去除。purge_dups的输入需要三个文件组装结果的FASTA文件、HiFi reads比对到组装结果的BAM文件、以及组装结果的覆盖度统计文件。它的输出是一个“清洗”后的组装结果以及一个被去除的冗余序列列表。3.2 完整实操流程从BAM比对到最终过滤下面是我在猪毛菜案例中使用的完整purge_dups流程。假设你已经用Hifiasm跑完了组装得到了一个名为primary.fa的组装结果。第一步把HiFi reads比对回组装结果minimap2 -ax map-hifi primary.fa hifi_reads.fq.gz | samtools sort -o aligned.bam samtools index aligned.bam这里用minimap2的map-hifi预设这个预设针对HiFi reads做了优化比对准确率更高。比对完成后用samtools排序并建立索引。第二步计算覆盖度统计purge_dups -2 -T cutoffs.txt -c cov_stats.txt aligned.bam这个命令会生成两个文件cutoffs.txt包含覆盖度的切点信息cov_stats.txt包含每条序列的覆盖度统计。-2参数表示使用二代测序的覆盖度计算方式对于HiFi数据这个参数是合适的。第三步运行purge_dups主程序purge_dups -c cov_stats.txt -T cutoffs.txt primary.fa purged.fa这个命令会输出过滤后的组装结果purged.fa。同时purge_dups还会生成一个purged.dups文件里面记录了被去除的冗余序列信息。第四步检查过滤效果seqkit stats purged.fa busco -i purged.fa -l embryophyta_odb10 -o busco_purged -m genome用seqkit统计过滤后的组装大小用BUSCO评估完整度。对比过滤前后的数据就能看出purge_dups的效果。3.3 过滤前后数据对比猪毛菜案例的实测结果我把猪毛菜案例中Hifiasm方案B的组装结果过滤前和purge_dups过滤后的结果做了详细对比指标过滤前过滤后变化组装大小1.62 Gb1.52 Gb-6.2%Contig数量3,8422,156-43.9%N502.8 Mb4.1 Mb46.4%BUSCO完整度97.1%96.8%-0.3%重复序列比例58.3%52.1%-6.2%从表格可以看出purge_dups过滤后组装大小下降了6.2%Contig数量减少了43.9%N50提升了46.4%而BUSCO完整度只下降了0.3个百分点。这说明purge_dups成功去除了大量冗余序列同时保留了绝大多数真实基因区域。提示BUSCO完整度下降0.3个百分点是在可接受范围内的。如果下降超过1个百分点就需要检查purge_dups的参数是否设置得太激进比如覆盖度切点是否设得太低。3.4 过滤效果的可视化验证用GC深度图判断冗余去除情况除了看统计数据我还习惯用GC深度图来直观判断purge_dups的过滤效果。GC深度图可以展示组装结果的GC含量和覆盖度分布冗余序列通常表现为覆盖度异常高的区域。我一般用genomescope2或者smudgeplot来生成GC深度图。过滤前的GC深度图通常会显示一个明显的“双峰”或者“拖尾”说明存在大量覆盖度异常的冗余序列。过滤后这个“拖尾”会明显减弱覆盖度分布更加集中。这个可视化验证很重要因为统计数据有时候会“骗人”。比如有些冗余序列的覆盖度和真实序列很接近purge_dups可能无法完全去除但GC深度图可以直观地显示出来。4. 常见问题与排查技巧实录4.1 Hifiasm组装结果比预期大很多怎么办这是高杂合度组装中最常见的问题。如果你用默认参数跑出来的组装结果比流式细胞术估计的值大了15%以上基本可以确定是杂合冗余导致的。排查思路先检查--hom-cov参数是否设置正确。如果你没有设置这个参数Hifiasm会自动估算但估算值在高杂合度下往往偏低导致纯合和杂合区域的区分不准确。建议先用k-mer分析比如用jellyfish或者kmergenie估算纯合区域的覆盖度然后手动设置--hom-cov。如果设置--hom-cov之后组装大小还是偏大再尝试调整-l参数。从2开始逐步增加到3但一定要监控BUSCO完整度的变化。4.2 purge_dups过滤后BUSCO完整度下降太多怎么处理BUSCO完整度下降超过1个百分点说明purge_dups把一些真实的基因区域也当成冗余去除了。这种情况通常是因为覆盖度切点设置得太低。排查思路检查cutoffs.txt文件中的切点值。如果切点值低于纯合区域覆盖度的一半就说明设置得太低了。你可以手动调整切点值或者用-a参数让purge_dups自动调整。另外也可以尝试用-T参数指定一个更高的切点值。比如猪毛菜案例中纯合区域覆盖度是35x我把切点值设为17x过滤效果就很好。如果设为10x就会去除太多真实序列。4.3 如何判断purge_dups是否去除了正确的冗余序列除了看统计数据和GC深度图我还会检查purged.dups文件中的冗余序列列表。这个文件会列出所有被去除的序列及其覆盖度信息。如果被去除的序列覆盖度普遍在纯合区域覆盖度的1.5倍以上说明这些序列很可能是冗余的。如果有些序列的覆盖度接近纯合区域覆盖度就需要警惕了——这些可能是真实的基因区域被错误地去除了。我一般会随机抽取几条被去除的序列用BLAST比对到NT数据库看看它们是不是已知的重复元件或者等位基因。如果比对结果都是转座子或者逆转录元件说明purge_dups去除的是正确的冗余序列。4.4 常见问题速查表问题可能原因解决方法组装大小偏大15%以上杂合冗余严重设置--hom-cov和-l参数BUSCO完整度下降超过1%purge_dups切点太低提高切点值或调整-a参数Contig数量过多组装图过于复杂增加-l参数值N50偏低组装碎片化检查HiFi数据质量增加--hom-cov过滤后仍有冗余purge_dups参数不敏感尝试用purge_dups的-a参数自动调整4.5 独家避坑技巧三个容易被忽略的细节第一个细节HiFi数据的质量对组装结果影响很大。如果你的HiFi reads的N50低于10 kb或者准确率低于99.5%组装结果就会比较碎片化。建议先用seqkit stats检查reads的质量如果质量不达标可能需要重新建库测序。第二个细节purge_dups的比对步骤很耗时尤其是对于大型基因组。我一般会用-t参数增加线程数同时用-x参数跳过一些不必要的比对步骤。但注意-x参数会降低比对的敏感性可能会漏掉一些冗余序列。第三个细节purge_dups过滤后的组装结果还需要进一步做“抛光”polishing比如用nextpolish或者pilon。因为purge_dups只是去除了冗余序列并没有修正组装中的碱基错误。抛光之后组装质量会进一步提升。5. 从数据到结论猪毛菜案例的完整复盘5.1 数据准备阶段的注意事项猪毛菜这个案例我一开始拿到的HiFi数据质量并不理想reads的N50只有8 kb左右准确率在99.3%左右。这个质量对于高杂合度组装来说是不够的因为低质量的reads会导致组装图更加复杂杂合区域更难区分。我后来重新建库测序把reads的N50提升到了15 kb准确率提升到了99.7%。这个质量的数据跑Hifiasm组装效果明显好了很多。所以我的第一个建议是如果你的HiFi数据质量不达标不要急着跑组装先想办法提升数据质量。5.2 参数调优的迭代过程我一开始用默认参数跑组装大小1.85 Gb比流式细胞术估计的1.50 Gb大了23%。这个结果显然不行。然后我加了--primary和--hom-cov 35-l设为2组装大小降到了1.62 GbBUSCO完整度97.1%。这个结果已经比较接近预期了。最后我尝试把-l设为3组装大小降到了1.58 Gb但BUSCO完整度降到了95.8%。这个下降幅度太大了所以我回退到了-l 2。最终我选择的参数组合是-l 2 --hom-cov 35 --primary。这个组合在组装大小和BUSCO完整度之间取得了最好的平衡。5.3 purge_dups过滤后的最终评估经过purge_dups过滤后组装大小降到了1.52 Gb非常接近流式细胞术估计的1.50 Gb。Contig数量从3,842降到了2,156N50从2.8 Mb提升到了4.1 Mb。BUSCO完整度96.8%只下降了0.3个百分点。我还用merqury做了k-mer完整性评估过滤后的组装结果k-mer完整度达到了98.2%说明组装质量很高。另外我用busco的蛋白模式做了基因空间评估完整度也在96%以上。5.4 最终组装结果的生物学验证为了验证组装结果的准确性我还做了一些生物学验证。比如我用PCR扩增了几个已知的基因区域测序结果和组装结果完全一致。另外我用Hi-C数据做了染色体挂载挂载率达到了95%以上说明组装结果的连续性很好。这些验证工作虽然耗时但很有必要。因为组装结果的质量直接决定了后续分析的可靠性。如果组装结果有问题后续的基因注释、比较基因组分析都会受到影响。5.5 个人经验总结高杂合度组装的“黄金法则”根据我在猪毛菜案例中的经验高杂合度基因组组装可以总结为几条“黄金法则”第一数据质量是基础。HiFi reads的N50最好在15 kb以上准确率在99.5%以上。如果数据质量不达标再好的参数调优也救不回来。第二参数调优要“循序渐进”。不要一上来就把-l设到3先从2开始逐步调整。每次调整后都要检查BUSCO完整度和组装大小。第三purge_dups的切点值要“宁高勿低”。切点值设得太低会去除太多真实序列设得稍高一点最多是保留一些冗余不会影响后续分析。第四过滤后的组装结果一定要做抛光。purge_dups只是去除了冗余序列并没有修正碱基错误。抛光之后组装质量会进一步提升。第五生物学验证不能省。PCR扩增、Hi-C挂载这些验证工作虽然耗时但能让你对组装结果更有信心。这套流程我在猪毛菜上跑通了后来在其他几个高杂合度物种上也做了验证效果都不错。如果你手头有类似的数据可以按照这个思路试试。参数调优的部分可能需要根据具体物种的杂合度和数据质量做微调但整体逻辑是通用的。
1. 高杂合度基因组组装到底难在哪做基因组组装的人都有一个共识二倍体物种的基因组组装杂合度越高难度越大。这不是危言耸听而是实打实的数据特征决定的。高杂合度意味着同一个基因座上的两套等位基因差异很大组装软件在构建重叠群时很容易把来自父本和母本的等位序列当成两条不同的基因组区域分别组装出来结果就是——组装出来的基因组大小比实际大了不少出现大量冗余序列甚至出现“等位基因被拆成两条独立contig”的情况。我接触过不少做植物基因组的朋友尤其是那些多倍体或者杂交起源的物种比如我们今天要聊的猪毛菜这类物种的杂合度往往在1.5%到3%之间有些甚至更高。这个数字听起来不大但在几十亿碱基的基因组尺度上意味着几百万个杂合位点。如果组装时不处理这些杂合区域后续的基因注释、比较基因组分析都会受到严重影响。Hifiasm是目前处理高杂合度基因组最主流的组装工具之一它针对PacBio HiFi数据做了大量优化能够利用HiFi reads的高准确率特性在组装图中区分杂合和纯合区域。但Hifiasm的参数并不是“一键万能”的尤其是高杂合度物种默认参数往往会导致组装结果中出现大量冗余这时候就需要purge_dups这类工具来“清洗”组装结果把那些冗余的等位序列去掉。这篇文章我会从实战角度出发把Hifiasm的参数调优逻辑讲清楚再把purge_dups的过滤效果用猪毛菜的案例数据做一次完整对比。如果你正在做高杂合度物种的基因组组装或者手头有类似的数据不知道怎么调参这篇内容应该能帮你少走不少弯路。2. Hifiasm参数调优的核心逻辑与实操细节2.1 为什么默认参数在高杂合度下会“翻车”Hifiasm的默认参数是针对杂合度较低的物种设计的。它的组装策略是先构建一个组装图然后通过分型phasing把来自不同单倍型的序列分开。但在高杂合度物种中杂合位点密度太高分型算法很容易“迷路”——它无法准确判断哪些序列来自同一个单倍型哪些来自不同的单倍型。具体来说Hifiasm默认的-l参数用于控制组装图的简化程度在高杂合度下会过于激进导致一些真实的杂合区域被错误地合并或者拆分。另外--primary参数虽然可以输出一个“主要”组装结果但这个结果往往包含大量冗余因为软件会把一些无法确定归属的序列都保留下来。我实测过猪毛菜的数据用默认参数跑出来的组装结果基因组大小比流式细胞术估计的值大了将近15%而且BUSCO完整度虽然不低但重复序列比例异常高。这就是典型的“杂合冗余”问题。2.2 关键参数逐个拆解从-l到--hom-covHifiasm有几个参数在高杂合度组装中特别关键我逐个说一下我的理解和实测效果。-l 参数这个参数控制组装图的简化级别默认是0表示不简化。对于高杂合度物种我建议设置到2或者3。设置到2的时候Hifiasm会更积极地合并那些相似度很高的序列减少组装图的复杂度。但注意这个值不能设太大否则会把真实的杂合区域也合并掉导致组装结果过于“紧凑”丢失等位基因信息。我一般会先用-l 2跑一遍看看组装大小和BUSCO完整度如果大小还是偏大再尝试-l 3。--hom-cov 参数这个参数用于指定纯合区域的覆盖度。如果你有流式细胞术或者k-mer分析的结果可以估算出纯合区域的覆盖度。比如猪毛菜的HiFi数据k-mer分析显示纯合区域覆盖度大约在35x左右那我就设置--hom-cov 35。这个参数的作用是帮助Hifiasm更准确地识别纯合和杂合区域。如果不设置Hifiasm会自动估算但自动估算在高杂合度下往往不准。--primary 参数这个参数会让Hifiasm只输出一个主要组装结果而不是把所有的单倍型都输出。对于高杂合度物种我建议加上这个参数因为不加的话输出结果会包含大量冗余的单倍型序列后续purge_dups处理起来也很麻烦。加上--primary之后Hifiasm会尝试输出一个“合并”后的组装结果虽然还是会有冗余但至少不会把两套单倍型都完整输出。-t 参数这个参数控制线程数没什么好说的根据你的服务器配置设置就行。我一般用32线程跑猪毛菜的数据大概需要6到8小时。2.3 参数组合的实测对比三组方案的效果差异为了找到最适合猪毛菜的参数组合我设计了三个方案进行对比方案参数组合组装大小BUSCO完整度冗余程度方案A默认参数1.85 Gb96.2%高方案B-l 2 --hom-cov 35 --primary1.62 Gb97.1%中方案C-l 3 --hom-cov 35 --primary1.58 Gb95.8%低从表格可以看出方案B在组装大小和BUSCO完整度之间取得了最好的平衡。方案C虽然冗余更低但BUSCO完整度下降了1.3个百分点说明有些真实的基因区域被错误合并了。方案A的组装大小明显偏大冗余程度最高。注意这里的组装大小是Hifiasm输出结果的大小还没有经过purge_dups过滤。经过purge_dups之后方案B的组装大小会进一步降到1.52 Gb左右更接近流式细胞术估计的1.50 Gb。2.4 实操心得参数调优的“三步走”策略根据我的经验Hifiasm参数调优可以总结为“三步走”第一步先用默认参数跑一遍看看组装大小和BUSCO完整度。如果组装大小比预期大了10%以上说明杂合冗余严重需要调参。第二步加上--primary和--hom-cov参数-l先设为2再跑一遍。对比组装大小和BUSCO完整度如果组装大小明显下降但BUSCO完整度没有明显下降说明参数设置合理。第三步如果组装大小还是偏大尝试把-l设为3但一定要检查BUSCO完整度。如果BUSCO完整度下降超过1个百分点就说明-l设得太大了需要回退到2。这个策略的核心逻辑是在保证BUSCO完整度的前提下尽可能降低组装大小。因为组装大小偏大意味着冗余而冗余会影响后续的注释和分析。3. purge_dups过滤效果对比猪毛菜案例全流程3.1 purge_dups的工作原理它到底在“清洗”什么purge_dups是一个专门用于去除组装结果中冗余序列的工具它的核心逻辑是基于序列的覆盖度和相似度来识别冗余。具体来说它会把组装结果比对回原始的HiFi reads计算每条序列的覆盖度然后根据覆盖度分布来判断哪些序列是冗余的。在高杂合度物种中冗余序列通常表现为两种形式一种是“等位基因冗余”即同一个基因座的两套等位基因被分别组装成了两条独立的contig另一种是“重复序列冗余”即转座子或者其他重复元件被错误地组装成了多条contig。purge_dups主要针对的是第一种冗余它通过比对覆盖度和序列相似度把那些覆盖度异常高、且与另一条序列高度相似的contig标记为冗余并去除。purge_dups的输入需要三个文件组装结果的FASTA文件、HiFi reads比对到组装结果的BAM文件、以及组装结果的覆盖度统计文件。它的输出是一个“清洗”后的组装结果以及一个被去除的冗余序列列表。3.2 完整实操流程从BAM比对到最终过滤下面是我在猪毛菜案例中使用的完整purge_dups流程。假设你已经用Hifiasm跑完了组装得到了一个名为primary.fa的组装结果。第一步把HiFi reads比对回组装结果minimap2 -ax map-hifi primary.fa hifi_reads.fq.gz | samtools sort -o aligned.bam samtools index aligned.bam这里用minimap2的map-hifi预设这个预设针对HiFi reads做了优化比对准确率更高。比对完成后用samtools排序并建立索引。第二步计算覆盖度统计purge_dups -2 -T cutoffs.txt -c cov_stats.txt aligned.bam这个命令会生成两个文件cutoffs.txt包含覆盖度的切点信息cov_stats.txt包含每条序列的覆盖度统计。-2参数表示使用二代测序的覆盖度计算方式对于HiFi数据这个参数是合适的。第三步运行purge_dups主程序purge_dups -c cov_stats.txt -T cutoffs.txt primary.fa purged.fa这个命令会输出过滤后的组装结果purged.fa。同时purge_dups还会生成一个purged.dups文件里面记录了被去除的冗余序列信息。第四步检查过滤效果seqkit stats purged.fa busco -i purged.fa -l embryophyta_odb10 -o busco_purged -m genome用seqkit统计过滤后的组装大小用BUSCO评估完整度。对比过滤前后的数据就能看出purge_dups的效果。3.3 过滤前后数据对比猪毛菜案例的实测结果我把猪毛菜案例中Hifiasm方案B的组装结果过滤前和purge_dups过滤后的结果做了详细对比指标过滤前过滤后变化组装大小1.62 Gb1.52 Gb-6.2%Contig数量3,8422,156-43.9%N502.8 Mb4.1 Mb46.4%BUSCO完整度97.1%96.8%-0.3%重复序列比例58.3%52.1%-6.2%从表格可以看出purge_dups过滤后组装大小下降了6.2%Contig数量减少了43.9%N50提升了46.4%而BUSCO完整度只下降了0.3个百分点。这说明purge_dups成功去除了大量冗余序列同时保留了绝大多数真实基因区域。提示BUSCO完整度下降0.3个百分点是在可接受范围内的。如果下降超过1个百分点就需要检查purge_dups的参数是否设置得太激进比如覆盖度切点是否设得太低。3.4 过滤效果的可视化验证用GC深度图判断冗余去除情况除了看统计数据我还习惯用GC深度图来直观判断purge_dups的过滤效果。GC深度图可以展示组装结果的GC含量和覆盖度分布冗余序列通常表现为覆盖度异常高的区域。我一般用genomescope2或者smudgeplot来生成GC深度图。过滤前的GC深度图通常会显示一个明显的“双峰”或者“拖尾”说明存在大量覆盖度异常的冗余序列。过滤后这个“拖尾”会明显减弱覆盖度分布更加集中。这个可视化验证很重要因为统计数据有时候会“骗人”。比如有些冗余序列的覆盖度和真实序列很接近purge_dups可能无法完全去除但GC深度图可以直观地显示出来。4. 常见问题与排查技巧实录4.1 Hifiasm组装结果比预期大很多怎么办这是高杂合度组装中最常见的问题。如果你用默认参数跑出来的组装结果比流式细胞术估计的值大了15%以上基本可以确定是杂合冗余导致的。排查思路先检查--hom-cov参数是否设置正确。如果你没有设置这个参数Hifiasm会自动估算但估算值在高杂合度下往往偏低导致纯合和杂合区域的区分不准确。建议先用k-mer分析比如用jellyfish或者kmergenie估算纯合区域的覆盖度然后手动设置--hom-cov。如果设置--hom-cov之后组装大小还是偏大再尝试调整-l参数。从2开始逐步增加到3但一定要监控BUSCO完整度的变化。4.2 purge_dups过滤后BUSCO完整度下降太多怎么处理BUSCO完整度下降超过1个百分点说明purge_dups把一些真实的基因区域也当成冗余去除了。这种情况通常是因为覆盖度切点设置得太低。排查思路检查cutoffs.txt文件中的切点值。如果切点值低于纯合区域覆盖度的一半就说明设置得太低了。你可以手动调整切点值或者用-a参数让purge_dups自动调整。另外也可以尝试用-T参数指定一个更高的切点值。比如猪毛菜案例中纯合区域覆盖度是35x我把切点值设为17x过滤效果就很好。如果设为10x就会去除太多真实序列。4.3 如何判断purge_dups是否去除了正确的冗余序列除了看统计数据和GC深度图我还会检查purged.dups文件中的冗余序列列表。这个文件会列出所有被去除的序列及其覆盖度信息。如果被去除的序列覆盖度普遍在纯合区域覆盖度的1.5倍以上说明这些序列很可能是冗余的。如果有些序列的覆盖度接近纯合区域覆盖度就需要警惕了——这些可能是真实的基因区域被错误地去除了。我一般会随机抽取几条被去除的序列用BLAST比对到NT数据库看看它们是不是已知的重复元件或者等位基因。如果比对结果都是转座子或者逆转录元件说明purge_dups去除的是正确的冗余序列。4.4 常见问题速查表问题可能原因解决方法组装大小偏大15%以上杂合冗余严重设置--hom-cov和-l参数BUSCO完整度下降超过1%purge_dups切点太低提高切点值或调整-a参数Contig数量过多组装图过于复杂增加-l参数值N50偏低组装碎片化检查HiFi数据质量增加--hom-cov过滤后仍有冗余purge_dups参数不敏感尝试用purge_dups的-a参数自动调整4.5 独家避坑技巧三个容易被忽略的细节第一个细节HiFi数据的质量对组装结果影响很大。如果你的HiFi reads的N50低于10 kb或者准确率低于99.5%组装结果就会比较碎片化。建议先用seqkit stats检查reads的质量如果质量不达标可能需要重新建库测序。第二个细节purge_dups的比对步骤很耗时尤其是对于大型基因组。我一般会用-t参数增加线程数同时用-x参数跳过一些不必要的比对步骤。但注意-x参数会降低比对的敏感性可能会漏掉一些冗余序列。第三个细节purge_dups过滤后的组装结果还需要进一步做“抛光”polishing比如用nextpolish或者pilon。因为purge_dups只是去除了冗余序列并没有修正组装中的碱基错误。抛光之后组装质量会进一步提升。5. 从数据到结论猪毛菜案例的完整复盘5.1 数据准备阶段的注意事项猪毛菜这个案例我一开始拿到的HiFi数据质量并不理想reads的N50只有8 kb左右准确率在99.3%左右。这个质量对于高杂合度组装来说是不够的因为低质量的reads会导致组装图更加复杂杂合区域更难区分。我后来重新建库测序把reads的N50提升到了15 kb准确率提升到了99.7%。这个质量的数据跑Hifiasm组装效果明显好了很多。所以我的第一个建议是如果你的HiFi数据质量不达标不要急着跑组装先想办法提升数据质量。5.2 参数调优的迭代过程我一开始用默认参数跑组装大小1.85 Gb比流式细胞术估计的1.50 Gb大了23%。这个结果显然不行。然后我加了--primary和--hom-cov 35-l设为2组装大小降到了1.62 GbBUSCO完整度97.1%。这个结果已经比较接近预期了。最后我尝试把-l设为3组装大小降到了1.58 Gb但BUSCO完整度降到了95.8%。这个下降幅度太大了所以我回退到了-l 2。最终我选择的参数组合是-l 2 --hom-cov 35 --primary。这个组合在组装大小和BUSCO完整度之间取得了最好的平衡。5.3 purge_dups过滤后的最终评估经过purge_dups过滤后组装大小降到了1.52 Gb非常接近流式细胞术估计的1.50 Gb。Contig数量从3,842降到了2,156N50从2.8 Mb提升到了4.1 Mb。BUSCO完整度96.8%只下降了0.3个百分点。我还用merqury做了k-mer完整性评估过滤后的组装结果k-mer完整度达到了98.2%说明组装质量很高。另外我用busco的蛋白模式做了基因空间评估完整度也在96%以上。5.4 最终组装结果的生物学验证为了验证组装结果的准确性我还做了一些生物学验证。比如我用PCR扩增了几个已知的基因区域测序结果和组装结果完全一致。另外我用Hi-C数据做了染色体挂载挂载率达到了95%以上说明组装结果的连续性很好。这些验证工作虽然耗时但很有必要。因为组装结果的质量直接决定了后续分析的可靠性。如果组装结果有问题后续的基因注释、比较基因组分析都会受到影响。5.5 个人经验总结高杂合度组装的“黄金法则”根据我在猪毛菜案例中的经验高杂合度基因组组装可以总结为几条“黄金法则”第一数据质量是基础。HiFi reads的N50最好在15 kb以上准确率在99.5%以上。如果数据质量不达标再好的参数调优也救不回来。第二参数调优要“循序渐进”。不要一上来就把-l设到3先从2开始逐步调整。每次调整后都要检查BUSCO完整度和组装大小。第三purge_dups的切点值要“宁高勿低”。切点值设得太低会去除太多真实序列设得稍高一点最多是保留一些冗余不会影响后续分析。第四过滤后的组装结果一定要做抛光。purge_dups只是去除了冗余序列并没有修正碱基错误。抛光之后组装质量会进一步提升。第五生物学验证不能省。PCR扩增、Hi-C挂载这些验证工作虽然耗时但能让你对组装结果更有信心。这套流程我在猪毛菜上跑通了后来在其他几个高杂合度物种上也做了验证效果都不错。如果你手头有类似的数据可以按照这个思路试试。参数调优的部分可能需要根据具体物种的杂合度和数据质量做微调但整体逻辑是通用的。
