Parabricks FQ2BAM 验收｜用一条官方比对记录算清读长、覆盖范围和片段长度

拿到一份双端比对结果，同事问这条记录有多长，三个回答可能都对。测序仪交出了多少个碱基，这一端横跨参考上的多少个位置，这一对读段来自多长的片段，分别要读不同的字段。NVIDIA 的 FQ2BAM 教程恰好列出一条同时带有插入、配对坐标和另一端比对摘要的记录。我们用 SAMtools 把它做成可查询的文件，逐项验证三个长度，最后再考察一个很容易误解的现象。明明标记成了正常配对，在另一端的坐标上却可以一条记录也取不到。NVIDIA FQ2BAM 教程

收到结果先看一条能算清楚的记录

FQ2BAM 接收参考序列和成对的原始读段，生成按参考坐标排序的比对文件。核对时先跟对一个片段从哪里来，把原始读段名、样本名、测序批次、使用的参考版本留在记录中。不同来源的样本都可能覆盖同一个位置，因此取到一条看似合理的记录，并不足以说明交接的是预期样本。生产数据可以按读组追到样本和文库，浏览器显示的染色体前缀也要与文件头一致。接下来所有坐标都采用 SAM 的一基、两端包含约定。

我们取的是官方教程展示的首条数据，在本地补一个只含所需参考条目的最小文件头，用实际安装的 SAMtools 1.19.2 转成 BAM 并建索引。保留官方记录原有的方向、碱基、质量、位置和标签，不为配偶编造第二条记录。读者可以对自己的验收样本做同样的事，事先写下预期的四个边界再查询。教程里的 GPU 程序和阶段时间出自 NVIDIA 的示例，本次本地运行的是下游格式检查与区间查询。SAM 格式规范

读长为 115，参考跨度为 110

该记录从一号染色体的 10027 开始。把比对串按操作拆开，依次是二十四个参与匹配的碱基、五个相对参考的插入、八十六个参与匹配的碱基。要回答测了多长，就把消耗查询序列的三段相加，得到一百一十五。要回答在参考上走了多远，只加两段匹配操作，得到一百一十。起点本身占第一个位置，所以末坐标还要减一。直接数这条记录的序列字符串和质量字符串，两者长度也都是一百一十五。

匹配操作允许碱基相同或不同。若想知道参考与读段实际差多少，需要再读编辑距离与错配描述。例子里的错配描述显示六次参考碱基替换，其中两处之间相隔零个匹配碱基。这个零是两个连续错配之间的间隔。再从比对串加入插入的五个碱基，正好与编辑距离十一相合。把编辑距离直接记成十一处单碱基替换，或者把两段匹配操作全部算作精确相同，都会让错误类型统计失真。遇到参考缺失时还要识别错配描述里以尖号开头的删除串，不能把其中字母逐一当作替换。

配偶的软剪切去了哪里

十进制标志九十九经工具解码，表明这是成对测序中第一端，软件将其标为正常配对，另一端落在反向链。这个标志只说记录自带的状态。要追踪另一端，先读配偶参考名称和起点，再看配偶比对摘要。例子中的另一端有三十四个软剪切碱基，它们保留在那一端的查询序列中，在参考上不向前推进。还有一个插入也只占查询。由此可以独立算出配偶本身有一百一十五个碱基，横跨的参考却只有八十个位置。

本例两端位于同一参考，外侧端点相距二百三十一个位置，恰好等于第一端记录里给出的正片段长度。如果把两条读长直接相加，得到二百三十，数值只差一，理由却完全不同。两个读段之间留有未被这两端覆盖的区间，配偶前部又发生软剪切。本例能从配偶标签交叉验证外侧端点，是因为标签齐全。实际文件可能没有这个可选标签，此时要取回配偶记录，从它自己的比对串重算。遇到跨染色体和补充比对时也要重新确认片段长度字段是否有可解释的含义。

实跑经验｜用两个刚好落空的坐标检查算法

samtools flags 99
samtools view -b -o one.bam one.sam
samtools index one.bam
samtools view -c one.bam chr1:10027-10027
samtools view -c one.bam chr1:10136-10136
samtools view -c one.bam chr1:10137-10137
samtools view -c one.bam chr1:10178-10178

分别查询起点前一位、起点、末位、末位后一位，返回数依次是零、一、一、零。配偶起点处也返回零。原因可以直接从实验设计找到，重建的文件只包含第一端，索引按文件里真实存在的比对记录建立。一个标签里提到另一个位置，不会在该处凭空生成第二条记录。以后有人用一个窄窗口截出双端数据，再拿它评估成对恢复率，先问截取方法是否同时找回了落在窗口以外的配偶。整库与小窗口的分母往往就在这里发生变化。SAMtools 区间读取说明

区域查询返回的是比对记录的重叠关系。若下一步要交付覆盖深度，需要另外定义碱基质量、比对质量、删除位置和重叠双端各如何计数。尤其在有较大缺失或跨内含子比对时，一条记录从区间左边延伸到右边，区间中央却可能没有任何被测到的碱基。所以这组边界测试适合检查坐标换算，深度验收还要取同一坐标查看堆叠结果。把这两个问题分开记录，后续才能解释浏览器里一个长条和深度表里一个零为什么可以同时成立。

官方日志另有三个阶段分别用时五十、五和十七秒，相加为七十二秒，从首段开始到末段结束的时钟跨度却是七十四秒。做流程验收时，将阶段统计和经过时间分别记在两栏，并注明数据规模、是否包括排序与重复标记。这样同事既能复核结果文件，又能看懂时间是从哪一步开始计的。

把这套小窗口验收用到整批全基因组数据时，如果还要长期保留原始读段和比对结果供多人核查，可以了解云生信的生物信息分析服务器。