拿到 Hi-C 辅助组装的结果,有人会对照参考从左到右逐个看片段名。预期是甲、乙、丙,程序却给出丙、乙、甲,是否意味着整条染色体拼反了?如果整条双链 DNA 连同每个片段的方向一起翻转,两种写法可以表示同一个连接。我们先在一个已知顺序的小基因组上用 YaHS 做真实挂载,再从输入中拿掉一组特定的接触,看看哪个连接真的失去支持。两次的片段序列分毫未动,变化发生在哪个位置可以直接从 AGP 重建。YaHS 作者的算法说明
让一个片段在输入时就反向
模拟三条各二十万碱基的独特序列,原始染色体从左到右依次为甲、乙、丙。写入待挂载的 FASTA 时,甲和丙保持原方向,乙写成反向互补。因此按输入文件给部件标方向,期望连接是正甲、反乙、正丙。在长分子上的两个接触点按随距离衰减的简化分布抽样,随后才映射回三个输入片段自己的坐标。第一个点若落在反向写出的乙上,坐标也要用片段长度减一减去原坐标。只换名字、不换坐标的模拟,会给挂载程序制造出一个错误的方向。
共生成十二万条人工接触对,给两端相同的高比对质量。输入采用每行一对的路径格式,保存读对名、两个片段名、各自位置和两端比对质量。我们实际安装了 YaHS 1.2.2,并用索引好的相同序列运行两次。第一轮保留全部接触,第二轮只去掉连接甲与乙的六千二百一十六条,剩十一万三千七百八十四条。乙与丙之间的六千三百一十二条全部保留,每条片段内部的接触也原样保留。为把问题集中在挂载,两轮都关闭片段纠错和后续挂载纠错,使用相同分辨率。YaHS 输入与参数
从 AGP 重建连接,不靠 FASTA 标题猜方向
第一轮结果是一条长六十万零二百碱基的 scaffold。依次读取三个序列部件,是反丙、正乙、反甲。把整个结果倒过来读,同时对每个部件取反,便得到正甲、反乙、正丙,与设计完全一致。若只反转片段顺序而不切换各自符号,就会错判乙相对两侧的方向。对于没有外部染色体方向标记的挂载,先看哪两个端点相邻,再看整条染色体该以哪一端作为显示起点。后一项可以随展示约定改变,前一项直接影响跨断点基因与结构解释。
三个序列部件的长度总和始终是六十万。第一轮在两个接头各添加一百个未知碱基作占位,AGP 用 U 表示长度未确定的间隔,连接证据写明来自邻近连接实验。因此最后比真实序列多二百,并不意味着程序测到了两段各一百碱基的新序列。要统一多个版本的组装总长,可以同时报告序列部件碱基数、间隔占位数以及包含间隔的挂载长度。遇到连续的未知碱基时,再查它们是原始片段内部本来存在的,还是此次挂载新添的。NCBI AGP 格式规范
| 接触对输入 | 片段间的有效连接 | scaffold 数 | AGP 间隔 |
|---|---|---|---|
| 全部 120000 条 | 甲乙 6216、乙丙 6312 | 1 | 两个,各 100 bp |
| 去掉全部甲乙接触 | 甲乙 0、乙丙仍 6312 | 2 | 一个,100 bp |

还有一个不依赖参考序列的检查。两轮把所有标为序列部件的行按原始片段名汇总,确认甲、乙、丙各出现一次,每段起点为一、终点为二十万。再分开汇总不同挂载的未知间隔。这样若有一段突然重复计入两个位置,或一条短片段在处理时被悄悄过滤,即使最终挂载的条数完全符合预期,也能从原始部件清单中查出来。换用真实基因组时,先把被长度阈值排除的片段单列出来,用同一批输入定义分母。
实跑经验|判断一个分支要看指定两端的证据
samtools faidx contigs.fa
yahs --no-contig-ec --no-scaffold-ec -r 5000,10000,20000 -o all contigs.fa all.pa5
yahs --no-contig-ec --no-scaffold-ec -r 5000,10000,20000 -o reduced contigs.fa no_AB.pa5
awk -F '\t' '$5=="W"{print $1,$6,$9}' reduced_scaffolds_final.agp
去掉甲乙接触后,第二轮输出一条由反丙、正乙组成的挂载,另有一条独立的正甲。乙丙之间的相对顺序和方向与第一轮一致,原来甲连到乙的位置消失。两条挂载加起来依旧包含全部六十万已知碱基,只有占位减少了一百。若只看总长几乎没变,很容易错过一次结构变化。反过来,若把两份 AGP 中同一个片段的正负号逐行求差,又会把单独存在的甲换方向也统计成一次生物学倒位。可以先将每个接头写成“哪个片段的哪一端连接到哪一端”,再比较两套接头集合。
真实数据的困难在于,观测并不都像本例这样单一。重复片段会出现多个落点,过度扩增可能带来多个几乎重合的接触,染色体之间也有背景联系。评估可疑连接时,可以把它放回不同分辨率的接触图,在候选断点两侧分别考察邻近距离的支持是否连续,再核对配对质量和唯一定位。一次局部删减导致断开,只说明这个模拟中的那组证据对该连接重要,并不构成“达到多少接触数就一定可连”的通用阈值。若同一物种还有长读段、遗传连锁图或同源共线性,可以把这些证据按接头编号留在同一份审核表中。
最后保存 AGP 和最初的输入 FASTA。只保留挂载后的序列会失去部件方向与原始接头,后续改动一个区段时难以判断坐标属于哪个版本。利用 AGP 的一基闭区间可逐段校验长度,再把候选基因映射到新坐标,对靠近未知间隔的多外显子模型单独复核。这样即使以后取消甲乙连接,也能迅速找到哪些注释需要重新检查。
当多个样本要反复调整接头并保存接触图及坐标映射时,可以在云生信的生物信息分析服务器共同复核各版挂载。
