hifiasm 实跑｜两套组装都长了一倍，怎样判断序列拼对了、单倍型却分错了？

拿到两个单倍型文件，里面都找到了候选基因的两种等位。研究人员可能马上讨论基因加倍，但还存在另一种很容易验证的解释。两条等位序列都被正确组装出来，随后在分配到两个文件时各保留了一份。我们用两条已知模板实跑 hifiasm，把输出的每个片段和模板逐碱基核对，然后只改变纯合覆盖参数再组装一次。两轮都找回全部设计的标记，只有一轮把两套标记分到了不同文件。这个区别决定了下游基因拷贝数到底能不能这么数。hifiasm 版本与用法

分清两个层级的正确

先生成一条十万碱基的序列，再在相隔二百二十七个碱基的位置替换一个字母，构成另一条等长模板。总共三百八十八个差异，其他位置都相同。为每个差异在每条模板上各取一个以它为中心、长三十一个碱基的窗口，共两套标记。逐一验证每个标记及其反向互补在两条模板合在一起的序列中只对应一个位置。如此一来，结果里同时检出某条标记的正向和反向表示时，只记同一个等位，不因为组装方向不同而漏算或重算。

从每条模板各取五百条无错长读段，每条长度介于六千八百和九千六百之间，交替使用正向和反向。少量读段特意触及两端，其余起点随机。将两组混在一起交给程序，不提供亲本序列、染色体构象数据或标记列表。实验设计者知道两条模板是谁，组装器看到的是同一份混合长读段。验收时分别问某一条组装片段里有没有另一模板的标记，以及一整个单倍型文件是否容纳了本该分开的两种片段。一个片段内没有混杂，并不自动保证整个文件的分组正确。

同时看片段、文件和两套文件的并集

安装并运行 hifiasm 0.25.0，在为小基因组关闭布隆过滤器的同一设定下，先让程序自行估计覆盖。第一份单倍型文件得到两条各十万碱基的片段。第一条命中模板甲的三百八十八个标记，模板乙的标记为零。另一条的计数倒过来。第二份单倍型文件也得到同样的两条，连主组装文件也如此。我们再把方向统一比较，两条片段分别与甲、乙全长逐碱基完全相同。这里确实保住了两种等位的序列，问题在于每个单倍型文件都含有两套完整序列。

若把两份单倍型文件拼成一个 FASTA，总长度会是四十万，而构造的二倍体合计只有二十万。只看“两套标记都找齐了”，两轮都会报满分。只看每条片段的长度或相似度，同样看不出重复分配。把文件名当作一列，片段名当作另一列，再分别统计甲、乙标记，就能看出重复发生在哪个层面。如果是真实基因所在区域，还应查看两条候选片段与侧翼是否一起共线、是否对齐到同一个等位位置。只凭同一个基因名出现两次，无法判断真实拷贝数。hifiasm 各类输出的含义

计数时还要防止输入文件本身制造假象。图文件可能同时有带序列的版本和省略序列的版本，后者的片段行以星号占位。应只从含实际序列的片段行导出 FASTA，连接行另存以便追查图结构。每次导出后核对片段条数、总长度和唯一名称，再去统计标记，可以避免将同一批序列的两种图文件一起加入结果。

实跑经验｜先估算覆盖峰应该落在哪里

hifiasm -f0 -t 2 -o auto reads.fa 2> auto.log
hifiasm -f0 -t 2 --hom-cov 82 -o fixed reads.fa 2> fixed.log
awk -F '\t' '$1=="S"{print $2, length($3)}' fixed.bp.hap1.p_ctg.gfa
grep 'coverage threshold' auto.log fixed.log

本例一千条读段平均约八千二百碱基，总产量除以十万碱基的单倍体长度，预计共享区域约八十二倍，每一种模板特有的区域约四十一倍。程序自动运行时在去重复阶段记录的纯合读段覆盖阈值是四十。第二轮先用输入规模做出八十二这个独立估计，再以新前缀显式设置并重跑，同一项日志变成八十二。在这一轮，两个单倍型文件各留下一条完整序列，主组装只保留其中一种。不能把主组装和两套单倍型再次简单相加作为总基因组量，它们提供不同的输出视角。hifiasm 关于组装膨胀和覆盖峰的答疑

这里的分母特别重要。若一开始用二十万碱基的二倍体总长度做除法，会把两套共同贡献的覆盖估成约四十一倍，数值恰巧很像一个单倍型自己的深度。野外样本还可能混有线粒体、性染色体和不均匀覆盖的区域，不能简单把任何项目的估值乘二照抄。先查看读段长度分布、原始产量、独立得到的单倍体基因组大小及覆盖直方图，同时挑一批有把握的单拷贝和杂合位点核对。把改动前后的多个指标一起保留，才能发现为缩短总长度而误删真实拷贝的情况。

这个小试验提供了理想模板，实际项目通常缺少全长真值，更依赖亲本特异标记、长程联系或独立遗传图谱。可以沿组装片段分窗口画出两套标记密度，观察一条片段是否在中途从甲占优变成乙占优，再看两份文件有没有整段重复出现。这样能分别记录局部相位切换和跨文件重复。没有亲本证据时，文件序号本身也不能直接赋予父源或母源身份。对于特定候选基因，把两侧序列、标记和跨越读段放进同一份表，会比给序列贴上亲本名称更有解释力。

需要为大型杂合基因组保存多轮组装图并比较标记分布时，可了解云生信的生物信息分析服务器。