比较两版染色体组装时,你可能看到组件从“甲在前、乙在后”变成“乙在前、甲在后”,每条组件的方向也都反了。单看坐标表,这很像剪开后交换了位置。可如果整条 scaffold 只是改为从另一端开始读,它的内部接头、序列成员和长度可以完全相同。我们运行 HapHiC 的参考排序命令,生成两套真实的坐标表及序列,逐碱基区别两种解释。本篇专门验证交付前的参考排序步骤,输入已有组件布局。HapHiC 官方项目
哪些文件决定了一份排序结果
原始 FASTA 保存 contig 的碱基。AGP 文件逐段描述每条 scaffold 由哪些组件组成、放置方向如何、间隔留多少未知碱基。第三份输入是原始组件与参考染色体的比对,对应哪个染色体、在哪个方向、覆盖多长,都从这里来。关键是使用和 AGP 组件栏一致的序列名。若把成品 scaffold 当作查询序列去比对,再把结果交给期待原始 contig 名的步骤,文件看上去都有碱基和坐标,却没法把证据正确关联到指定组件。作者的参考排序和输入说明
读这类数据先选一个有方向的接头。原布局若是正向的甲,接一百个未知碱基,再接正向的乙,那么从另一端读,顺序必然变成反向互补的乙、同样的一百个未知碱基、反向互补的甲。整条翻转会同时改变组件次序和方向。如果两版交换了组件的次序却只翻转其中一条,或把中间的接头换到别处,简单的整条翻转便解释不了。这条判断只用到序列本身,适合在看热图或共线图前先做。
实跑经验|五个原始组件,两条参考染色体
我们固定随机种子构造五条唯一序列。甲组由两万和一万七千碱基的两个组件组成,乙组由两万二千和一万八千碱基的两个组件组成,另有四千五百碱基的一条单独放置。两个已拼接组各留一百个未知碱基。制作模拟参考时,把甲组整体反向互补,乙组保持原方向,单独组件不放进参考。用 minimap2 为原始 contig 实际生成比对表,四个可匹配组件都端到端映射,映射质量均为六十。然后安装 HapHiC 1.0.7 的依赖,运行两次 refsort。
| 条件 | 输出的 scaffold 顺序 | 甲组的组件和方向 | 乙组和单独组件 |
|---|---|---|---|
| 原始布局 | scB、scA、unplaced | a 正向、b 正向 | c 与 d 正向,e 单独 |
| 按参考默认顺序输出 | scA、scB、unplaced | b 反向、a 反向 | 与原始相同 |
| 指定先列 chrB 再列 chrA | scB、scA、unplaced | b 反向、a 反向 | 与原始相同 |
把输出的甲组逐碱基与原始甲组的反向互补比较,三万七千一百个位置全部一致,中间仍是一百个未知碱基。乙组四万零一百个碱基逐字相同,未放置组件也完整保留。两次排序的文件顺序不同,但按 scaffold 名分别计算序列校验值,三个条目在两次输出间一一相同。这两次运行既没有引入新接头,也没有把未放置组件锚到参考上。可以从逐条碱基、组件身份和真实输出同时复核。

整条翻转后,旧注释坐标怎样换算
对这样一条仅整体翻转、且长度已核实为三万七千一百的序列,可以选一个小区间手算。原来一开始计数并包含两端的第一百零一到第二百位,翻转后落在第三万六千九百零一到第三万七千位,同时要交换链方向。计算过程是用总长度减去原来的右端再加一,得到新的左端,再用总长度减原来的左端加一作为右端。若采用从零开始的半开区间,公式的加一应相应改变。实际转换前先把原区间提取出来,对新位置拿到的序列做反向互补,两端和方向都能逐字核对。
minimap2 -x asm5 -c --secondary=no ref.fa contigs.fa > ref.paf
haphic refsort raw.agp ref.paf --fasta contigs.fa \
--fout default.fa > default.agp
haphic refsort raw.agp ref.paf --fasta contigs.fa --ref_order chrB,chrA \
--fout reordered.fa > reordered.agp
如果只有多条 scaffold 在 FASTA 里的排列顺序变了,很多序列索引工具仍根据名称找数据,但旧索引里的文件偏移量可能已失效。重新生成与新文件对应的索引,再核对染色体名、长度和方向。若标注依附于翻转的甲组,采用相同名字并不会自动修正原来的基因坐标。同理,变异位置和参考等位序列也要对新参考逐项核实。为了追溯,可以保留旧新两版布局、原始组件文件和位点转换表,避免仅交付同名序列让下游误以为坐标可直接沿用。
收到两个版本时,我会先按 scaffold 名比较长度、未知碱基数和原向序列校验值,只有原向核对不一致的才继续查反向互补。如果换向后仍然不一致,再沿组件表逐一比较成员、组件采用的起止坐标、间隔长度和接头。两个同名组件若使用了不同的截取范围,即使图上顺序和正负号都一致,也已经改变了交付序列。把核查一步步收窄到具体组件,能直接指出应该请哪一方回到上游文件核实。
对于未放置序列,最好再记录它在本次参考比对中有没有足够长的匹配、是否存在多个相近的命中,以及排序前后的序列是否一致。这里的第五个组件从未出现在模拟参考里,所以输出末尾仍有这条序列完全可以预期。换成真实材料,它可能来自参考缺少的插入片段、细胞器、污染或一段难以唯一定位的重复,需要靠样本自己的读段与序列特征继续判断。
什么情况下还需要看原始 Hi-C 证据
排序后如果希望宣称解决了染色体误拼,应具体指出哪个内部邻接发生了变化,以及哪些独立数据支持新接头。核对接头两侧的 Hi-C 联系、跨越分子、近缘材料中的共线关系,同时保留可能的真实重排。在多倍体材料中,近缘参考与多个同源拷贝可能都能很好匹配,不能仅凭统一朝向确定它们属于同一单倍型。此次实验的模拟序列能唯一比对,参考就是按已知答案制作的,读者可以在这一条件下练习核对顺序,却不能用它估计多倍体自动分组的准确率。参考排序实现与输出规则
下载本篇实验脚本与结果(12 篇合集 ZIP),解压后打开 080-haphic 目录。
需要让组装、注释和比较基因组人员反复核查同一批大文件时,可在云生信的生物信息分析服务器保留原始布局和交付坐标,供各环节按版本引用。
