samtools的coverage和bedcov结果不一致

Log

分别用samtools coverage和bedtools bedcov计算区域内的Reads数，得到的结果不一致。

软件版本

  • samtools 1.18

问题

  • samtools coverage的-r参数只能输入1个区域。
  • samtools bedcov可以统计整个bed文件的所有区域。
  • 但是coverage和bedcov -c，统计同一个bam文件的同一个区域，得到的Reads数不一样。

原因

  • bedcov和mpileup一样，受-d depth限制，只统计到设定的深度上限，再多的Reads就不统计了。
  • coverage默认的-d是1000000。bedcov的帮助文档没写是多少。

具体参数

bedcov帮助文档

  • -d参数：设置此参数后，结果会多一列，是参考碱基超过-d设定的深度的数量
  • -c参数：设置此参数后，结果会多一列，是这个区域内的Reads数量
    bedcov

coverage帮助文档

  • -d参数：最大深度阈值，默认1000000（一百万），如果设置为0则阈值是整形变量的最大值，相当于没有任何深度限制。
    coverage

详细说明

  • bedcov的-d参数统计有多少个参考碱基超过-d设定的深度，没有注明默认值是多少，当-d有设置值时，会在第4列的总碱基数量后面多1列结果，是参考碱基超过-d设定的深度的数量。但实际上结果中第4列的总碱基数量以及-c统计出的Reads数，都会受到-d限制。
  • 在-d不设置的状态下，-c统计出的Reads数最高只会有64020，但是当-d设置为1000000时，Reads数会统计出更高的结果，而且数量与coverage模块的默认参数（-d 1000000）一样。
  • coverage模块的-d参数如果设置为0，会去除深度限制。但是如果bedcov模块的-d设置为0，不会去除深度限制，只是总碱基数量后面会多1列，统计有多少个参考碱基深度超过0。
  • 查看bedcov的源码，samtools/bedcov.c at 1.18 · samtools/samtools (github.com)
    • 默认最大深度阈值应该是DEFAULT_DEPTH 64000（第43行）
    • -d输入的参数应该是min_depth（第120行）
    • 如果min_depth > 64000则用min_depth作为最大深度阈值，否则64000就是最大深度阈值（第223~226行）

测试代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
echo "# samtools coverage #"
/xxx/Software/samtools-1.18/bin/samtools coverage -r chr2:29446329-29446350 Test.bam
/xxx/Software/samtools-1.18/bin/samtools coverage -r chr6:117650534-117650556 Test.bam
echo ""

echo "# samtools bedcov -c #"
/xxx/Software/samtools-1.18/bin/samtools bedcov -c Test.bed Test.bam
echo ""

echo "# samtools bedcov -c -d 1000000 #"
/xxx/Software/samtools-1.18/bin/samtools bedcov -c -d 10000000 Test.bed Test.bam
echo ""

echo "# samtools bedcov -c -d 0  #"
/xxx/Software/samtools-1.18/bin/samtools bedcov -c -d 0 Test.bed Test.bam
echo ""

测试结果

result