测序深度的计算，你真的掌握了吗

计算测序深度是数据分析中的一个常规操作，常用的工具有以下几种

命令如下

代码语言：javascript复制

samtools depth input.bam > sample.depth

命令如下

代码语言：javascript复制

samtools mpileup -A  -Q 0 chrY.bam | cut -f1,2,4 > sample.depth

命令如下

代码语言：javascript复制

bedtools  genomecov -ibam input.bam -d > sample.depth

尽管方法不同，但是输出结果的格式是相同的，示意如下

第一列是染色体名称，第二列是染色体上的坐标，第三列是对应的测序深度。原本以为计算测序深度就是这么轻松的一件事，但是在比较不同方法的输出结果时，却发现部分区域samtools计算的结果和bedtools的结果对应不上，结果如下

第三列为samtools软件计算出来的结果，第四列为bedtools软件计算出来的结果，可以看到，samtools的结果比bedtools少了很多。

为了确定这段区域真实的测序深度，将bam文件导入IGV软件之中进行查看，对应区域的结果如下

从reads来看，确实应该是10和16条，那么samtools计算出来的结果又是什么，总不可能是samtools的bug吧，作为一个应用这个广泛的软件，不可能有如此低级的错误。

从结果来看，samtools在计算depth的过程中对reads进行了过滤，那么它过滤的原则是什么呢？通过查看bam文件中第二列的flag我找到了规律，143-150bp的reads有以下10条

第二列的flag含义如下

代码语言：javascript复制

0x91    145 PAIRED,REVERSE,READ2
0x163   355 PAIRED,PROPER_PAIR,MREVERSE,READ1,SECONDARY
0x63    99  PAIRED,PROPER_PAIR,MREVERSE,READ1

可以看到，其中有一个SECONDARY比对，对应数字355，这样的reads有4条，去除这4条之后，就是samtools计算的最终结果了。

作为SAM文件的官方操作工具，samtools内置了对flag的过滤，而bedtools默认没有进行这样的过滤，只是简单统计了该区域比对上的reads数目。

在使用IGV展示测序深度时，可以用igvtools先计算出tdf文件再导入IGV中，这样比直接导入bam文件快很多。在tdf文件中，其测序深度和bedtools软件的计算结果是一致的，也是只统计了read数目，没有做过滤。

0 人点赞