拿到一张组装图，看到一个节点标着三倍覆盖，最自然的追问是有几条分子经过它。两个数字之间并没有简单的一一对应。一个比对可以只走过节点的一半，一条读段也可以有多条候选比对。Verkko 的图处理工具里有一个很短的覆盖计算脚本，正适合把这两个问题拆开核算。我们取官方发布的 2.3.2 版，在本地实际运行其中的节点覆盖计算组件，用一个可以手算的图与五行路径比对得到下面的结果。此次实验从已有图和比对开始，验证的范围是这一项统计。Verkko 版本与输入输出说明
先把图上路径画成有长度的坐标轴
设四个节点分别长一百、二百、五十和四十个碱基，前三个依次相连，第四个没有比对经过。节点之间的重叠特意设为零，图里保留每段真实序列。第一条比对沿正向走过前三个节点，整条路径长三百五十，实际使用从四十到三百三十的区域，按零基左闭右开解释。第一个节点从四十到一百用了六十个碱基，第二个完整用了二百，最后一个用了三十。一共二百九十个碱基。除以各自节点长度，贡献分别为零点六、一和零点六，你可以先把这三个数写在纸上。
第二条读段从相反方向走同一组节点。路径的起点变成五十碱基节点的另一头，沿这个方向从二十走到三百一十，仍然覆盖两个端点各自的六成。若看见反向标记就直接把第一条路径的起止数字原样抄过去，会错置剪裁的方向。在图上找到路径首节点和尾节点，按照行进方向算剩下多少碱基，再以不带方向的节点名汇总，第二条的贡献恰好与第一条相同。图的节点深度因此可以出现小数，即使测序分子的数量总是整数。GAF 路径坐标与方向字段规范
接着放入一条只比到中间节点的读段。它在二百个碱基里覆盖了五十到一百五十这一百个位置，贡献零点五。为演示记录筛选，我们人为放入第二行同名、同区间的比对，给它的比对质量设为零，原来的那行是六十。最后再放入一条单独覆盖第三节点全长的读段。整个文件中只有四个不同读名，却有五条数据行。第二行同名比对是刻意制造的对照，实际项目中的多重比对还可能出现在不同路径或不同方向上，应按自己的分析目的定义取舍。
| 节点 | 长度 | 五行全部参与 | 去掉低质量行 | 哪些比对贡献了覆盖 |
|---|---|---|---|---|
| A | 100 | 1.2 | 1.2 | 两个端部覆盖各占 0.6 |
| B | 200 | 3.0 | 2.5 | 两条路径全长经过，加上一行或两行半节点命中 |
| C | 50 | 2.2 | 2.2 | 两个 0.6 与一个全长 |
| D | 40 | 0.0 | 0.0 | 节点存在，输入中没有比对经过 |
核对结果时，把未观察到的节点也放在表里。第四个节点有明确的长度且出现在图中，所以覆盖为零有可追溯的分母。如果某节点根本没有出现在这一版图里，左连接后应先标记为缺失，不能直接补零并解释为这个样本在该区域没有读段。两版组装对齐节点名时尤其容易混淆这两种情况。
实跑经验｜去掉一行后，用碱基总数验算
python3 calculate_coverage.py nodes.gfa < reads.gaf > all.tsv
awk -F '\t' '$12 >= 1' reads.gaf | python3 calculate_coverage.py nodes.gfa > filtered.tsv
cat all.tsv
cat filtered.tsv
直接运行官方组件，五行输入得到表格左侧的覆盖值。删除比对质量为零的那行再运行，中间节点从三变为二点五，两端节点与没有读段的节点都保持原值。两轮不同读名数完全一样。这项程序从标准输入逐行处理，没有在内部按读名去重，也没有在该组件中使用比对质量门槛。项目的其他阶段可以提供已筛选的输入，因此复核别人交付的覆盖表时，最值得先找的是当时实际传给这一项程序的路径文件，以及它在前一步经历了哪种筛选。Verkko 2.3.2 覆盖组件源码
还有一个不依赖脚本输出外观的验算办法。将每个节点的平均覆盖乘以自身长度，再把四项相加。第一轮是八百三十个累计经过的碱基，第二轮是七百三十，差值恰好等于移去的单条比对覆盖的一百个位置。也可以直接把每行路径已用的区间长相加，两种算法在这个零重叠图上应当完全一致。如果屏幕里看起来数值很合理，按长度加权却对不上，就要检查是否误用查询坐标代替路径坐标，或者在反向路径的末端少减了一个未覆盖区间。
这里的覆盖还和“节点上至少有一条读段覆盖的位置比例”不同。中间节点原本已有两条通穿路径，那一行额外的半节点比对会增加平均覆盖，却不会扩大其已覆盖的位置并集。反过来，独立读段数又是第三种指标。想用长读段证明一条图边是否真的可走，应另外统计从上一个节点跨到下一个节点、在两侧都有足够独特锚定的分子。单个节点深度较高并不能告诉你它的每一个出口都得到了同样多的支持，也不能单靠此值给重复片段定拷贝数。
小图中各段没有重叠，是为了直接检验上述碱基守恒。对真实组装图，先查节点序列是否完整、路径是否使用同一版图的节点名、连接是否有重叠、端点剪去多少，以及是否含特殊间隙节点。若 GFA 已将序列替换成星号只留下长度标签，本版本组件会从序列字段直接计长，这样的文件不适合作为本实验的输入。留下一条能够手工复核的往返路径，通常比凭全图平均值追查一个异常节点来得清楚。
当多个大型基因组需要共同保存组装图、长读段比对与覆盖计算版本时，可以了解云生信的生物信息分析服务器。
