同事说一份文件里有十条,有人过滤后剩七条,还有人按名字去重只看到四个。核查覆盖之前,值得先问一句,每个人究竟在数比对记录、已比对的读段端,还是测序模板?我们准备了只有十行的比对文件,让程序逐行筛选,还用两个重叠的区间查了两次。这样,输出里的每一个差值都可以追到某个标志位或某条重复返回的记录,读者可以照这个办法给自己的结果选择合适的分母。
一条测序模板可以在 BAM 里占几行
双端测序通常给同一个模板的两个端相同的查询名称,用标志位区分第一端和第二端。一个端可能落在多个候选位置,除主记录外还有次要记录。同一个端若分段比对,还可能有补充记录,描述同一查询在另一个位置或方向的部分。一行对应一个比对记录。完整遍历文件数行很容易,但这一数字既包括已定位的记录,也可能包含未定位的占位行。同名两行有时是同一模板的两端,有时是同一端的两个位置,这两种情形需要读标志位才能区分。SAM 规范对模板、主次与补充比对的定义
本例有两个双端模板,合计四行。另有一个单端查询含主、次两行,一个拆分查询含主、补充两行,一条标为重复的单端主记录,以及一条未比对记录。所有名字在这份构造数据中各有明确归属,所有记录都通过平台质量标记。用 SAMtools 1.19.2 对这份表排序、建索引,再分别查询。实验的输入标记由我们设定,因此结果只能说明筛选逻辑,不能用来评价某个比对器识别重复或断点的准确性。
实跑经验|按顺序加上三个排除条件
先看完整文件,十条记录对应六个不同的查询名。只去掉未比对记录后,剩下九条,对应五个不同名字。再排除次要和补充记录,剩七条,仍有五个名字。刚才减掉的两条各有一条同名主记录,所以行数减少,名字数不变。最后去掉标为重复的那条单端主记录,还剩六条,对应四个名字。其中两个名字是成对的,各占两行,另外两个名字各占一行。
| SAMtools 筛选 | 保留的比对记录 | 不同查询名 |
|---|---|---|
| 不筛选 | 10 | 6 |
| 排除未比对,标志掩码 0x4 | 9 | 5 |
| 只计已比对的主记录,0x904 | 7 | 5 |
| 再排除已标记重复,0xD04 | 6 | 4 |

标志位可以并存,适合用位掩码检查。排除掩码表示只要记录命中其中任一位便去掉。已比对主记录的掩码合并未比对、次要、补充三位。再加上重复位则得到最后一行。没有标为重复只表示输入里相应位未被设置。若上游根本没做重复标记,这个筛选照样可以正常执行,却不会去掉任何记录。若实验使用分子标签,按查询名统计也无法直接得到分子数。先读文件头中上游程序和实验设计,再决定是否有依据使用某个标记。SAMtools view 的标志与过滤选项
samtools sort -o sorted.bam fixture.sam
samtools index sorted.bam
samtools flagstat -O json sorted.bam
samtools view -c -F 0x904 sorted.bam
samtools view -F 0xD04 sorted.bam | cut -f1 | sort -u | wc -l
同一文件的汇总表给出已比对百分之九十,已比对主记录百分之八十七点五。第一项用九除以全部十行,第二项用七除以八条主记录。未比对那条也是主记录,因为它既没有次要位,也没有补充位。两项都按各自分母算对了。如果把七除以十后拿去与第二项对账,永远差十七点五个百分点。双端部分另报四端均正确配对,比例百分之百,它又有自己的分母,不应再拿全文件十行去除。SAMtools flagstat 对各项百分比的定义
这份小表里,汇总显示质量未通过的一栏全为零,因此已比对占比没有混入另一组记录。遇到真实数据,该工具把通过与未通过平台质量检查的记录分栏,读出两个数后,应先明确是否计划合并。若要计入后者,也要让相应分母一并改变。把两栏的分子相加,却沿用其中一栏的分母,得到的比对率甚至可能超过百分之百。
两个探针区间相交,记录会被返回几次
我们查询一号序列上一百到一百三十,以及一百一十五到一百四十五两个区间。把两个位置参数直接传入时,程序返回十行。改用多区间迭代器访问它们的并集,实际只有六条不同的比对记录,对应四个查询名。普通查询中另四行来自同时与两个区间相交的记录。坐标排序和索引能加快查找,没有消除不同查询区间的重叠。这里的重复属于检索行为,和哪条记录带有重复标志是两个独立问题。
samtools view -c sorted.bam chr1:100-130 chr1:115-145
samtools view -c -M sorted.bam chr1:100-130 chr1:115-145
若想分别评价每个外显子能捕获多少记录,同时与两个区间相交的读段可以各记一次,输出应保留区间身份。若想给整套捕获面板报一次总记录数,就应取区间并集,并保留所需的主次过滤条件。若还要统计两端都落在面板内的完整片段,仅对查询名去重不能区分另一端落在面板外、未比对或被过滤,要回到成对记录逐端检查。命令行区间使用一基、两端包含的坐标,使用 BED 文件时还要按其零基、右开格式准备边界,尤其注意只覆盖一个碱基的探针。
跨批次比较前,把每个分母写进列名。例如分别保存总记录、已比对主记录、排除重复后的主记录,以及对应的不同查询名数。再留一份带有全部标志的手核小样,不要只保存聚合后的百分比。若合并多条泳道或多个样本,先核查名字在不同读组之间是否可能重用,必要时使用读组和名字的组合键。这样在覆盖度突然变化时,可以先定位是哪一道条件改变了进入分析的记录。
固定统计口径后,需要在多个队列上批量核对这些指标时,可使用云生信的生物信息分析服务器集中处理 BAM 与索引文件。
