第一次用 hifiasm 组装二倍体基因组,结果目录里会出现 hap1 和 hap2。如果一个候选基因在两份文件里各有一条序列,能把它们分别叫作父本和母本拷贝吗?如果两条序列一模一样,又是组装失败了吗?
这两个问题都要追溯到同一件事:手里的读段,究竟把哪些差异连在了同一条 DNA 分子上? 下面以 hifiasm 0.25.0、同一个二倍体个体的 PacBio HiFi 数据为例,沿着这条证据读懂组装结果。hifiasm 仓库与版本
一条长读段,怎样帮助区分两套序列?
假设某段 DNA 有两个杂合位点。第一个位置有 A 和 G,第二个有 C 和 T。只知道这四个碱基,还无法确定 A 与 C 在同一条染色体上,还是 A 与 T 在一起。如果多条跨越两个位点的 HiFi 读段反复读到 A…C,另一批反复读到 G…T,两套局部组合就有了直接的分子证据。这个过程叫分相。

HiFi 的长度让一条读段有机会跨过多个差异,高准确度则有助于区分真实差异与偶发测序错误。hifiasm 寻找读段之间的重叠,并在校错时参考重叠读段在信息位点是否一致。这样,另一套单倍型上反复出现的真实碱基就有机会保留下来。接着,程序用相容的重叠关系建立组装图,把能连续延伸的路径拼成更长的序列,即 contig。这个过程可以从读段出发,无需先有参考基因组。hifiasm 原始方法论文
把组装图想成沿序列延伸的路径:两套拷贝相同的部分可以共用一段路径,遇到杂合差异时,路径可能分开,随后再汇合。这种局部结构常被称为“气泡”。读图时要留意,重复序列和残余错误也可能造成分支。看到气泡,要回到读段支持和邻近序列,才能判断它是否对应一对等位序列。
哪些数据能把局部分相延伸到更远?
假设左边一段已经分清 A/G,右边一段也分清 C/T,但中间是一段很长的纯合序列。若没有单条读段或连续重叠的读段链把两端的信息位点连起来,单靠局部序列,就可能无法确定哪一条左侧路径应该接上哪一条右侧路径。“两边各自分得开”和“从左到右始终属于同一套”需要不同的证据。
| 手里有什么数据 | 额外信息从哪里来 | 拿到两套输出后怎样解释 |
|---|---|---|
| 只有子代 HiFi | 同一读段及相互重叠的读段把杂合位点串成局部区块。覆盖足够、两套拷贝都有读段支持,且信息位点能连续衔接时,分相可以延伸得更远。 | 默认输出为部分分相的两套 contig。跨过缺乏信息的区域,亲本归属可能切换;hap1 的编号本身没有父源含义。 |
| 子代 HiFi,加父母双方的短读段 | 用 yak 分别统计双亲的短序列标记(k-mer),借助能区分双亲的标记为子代组装图上的路径定向。 | 若把父本标记交给 -1、母本标记交给 -2,两套输出就有相应的亲本解释。亲本样本身份、标记质量及目标区域有没有有效标记,都影响这一解释。 |
| 同一个体的 HiFi 和 Hi-C | 成对的 Hi-C 读段提供远距离接触信息,帮助判断相隔较远的区块该归到哪一套。 | 可以获得长距离分相。Hi-C 本身不给两套序列冠以父源或母源,分相跨度也应和 contig 连续性、染色体级挂载分别检查。 |
同一段在两套结果里完全一致,可以符合真实的纯合状态。若长重复超过了可桥接的范围,也可能出现两份拷贝被折叠的情况;这时要查覆盖、拷贝数和跨越重复边界的读段。单凭“相同”或“不同”都很难下结论。以上讨论面向二倍体,两套输出不能直接代表多倍体的全部拷贝。三种模式的官方用法 · Hi-C 长距离分相方法
先认清输出文件,再导出 FASTA
下面是只使用 HiFi 的命令示例,输入文件名和线程数按实际环境替换。第二条命令只导出 hap1,方便说明 GFA 与 FASTA 的区别。
hifiasm -o sample.asm -t 16 reads.hifi.fastq.gz 2> sample.asm.log
awk -F '\t' '$1=="S" {print ">"$2; print $3}' sample.asm.bp.hap1.p_ctg.gfa > sample.hap1.fa
GFA 中的 S 行存放序列片段,L 行记录片段间的连接。导出的 FASTA 可以用于比对和注释;原始 GFA 要留下,复查断点或分支时还能看到图结构。
常见文件(省略 sample.asm 前缀) |
适合回答的问题 |
|---|---|
.bp.p_ctg.gfa |
查看主组装的 contig。主组装从不同区段选择代表路径,不能直接当作某一个亲本的完整基因组。 |
.bp.hap1.p_ctg.gfa、.bp.hap2.p_ctg.gfa |
对照仅用 HiFi 得到的两套部分分相结果。想找一个基因的两条等位序列,应从这对文件的对应区域入手。 |
.dip.hap1.p_ctg.gfa、.dip.hap2.p_ctg.gfa;Hi-C 模式对应 .hic.hap1.p_ctg.gfa、.hic.hap2.p_ctg.gfa |
核对所选分相模式的输出。比较两套组装时,保留模式、输入样本和版本记录,避免把不同运行生成的文件配成一对。 |
*.bin、运行日志 |
二进制文件保存校正读段、重叠等中间结果,兼容的后续运行可以复用。日志可用于确认这一次实际读入了什么、是否复用了已有结果。 |
拿到结果后,按这三个问题验收
一、该有的序列找回来了吗? 分别统计 hap1 和 hap2 的总长、contig 数与 N50。对常规二倍体常染色体,每套的总长可先与独立估计的单倍体基因组大小对照,性染色体、真实结构变异等要另行解释。再用可信的高质量读段 k-mer 检查完整度和碱基准确度(QV)。两套合起来覆盖了多数可信 k-mer,并不能告诉你每一套是否各自完整。N50 只描述连续性,错误连接也可能把它抬高。Merqury 的 k-mer 评估方法
二、相位沿序列保持一致吗? 有双亲数据时,可以查看亲本特异 k-mer 在每条 contig 上的分布。若一段持续富集父本标记,后面转为持续富集母本标记,交界处值得检查是否出现了相位切换。标记稀少的区段证据较弱,可结合跨区读段,或在有 Hi-C 时检查接触图。即使一条 contig 很长,它内部有充分证据、能保持同一相位的区块也可能短得多。Merqury 的亲本标记与分相评估
三、候选基因的“两份”真是等位序列吗? 回到开头的问题。先把候选基因和两侧邻近序列一起比对到 hap1、hap2。若两边都有对应的邻近基因、顺序一致,关键差异又有多条 HiFi 读段跨过并支持,解释为等位差异才更有把握。若只能匹配基因内部一个保守结构域,或同一套组装在多个位置都能匹配,先排查旁系同源拷贝和重复;若候选恰在 contig 末端,则先排查组装缺口。这样才能区分“这套组装里没有找到”和“这个个体确实缺失一份”。
正式运行时可从作业记录中查看峰值内存,并在组装过程中记录工作盘最低剩余空间,把需要保留的 *.bin 与两套 GFA 一起算进去。若本地资源不足,带上预估基因组大小、HiFi 数据量和已有运行记录,再对照云生信的内存与工作盘配置选择租期即可。云生信高内存服务器租赁
