bam文件急速depth计算工具

WGS, exome, or targeted sequencing的数据分析的时候,经常需要计算全基因组整体的平均测序深度或者 target 区域的测序深度的分布情况以考察试剂盒的捕获效率。

通常samtools depth调用的是pileup的API来计算每个base的 depth,这个工具其实运算速度挺慢的,因为对于每条read都要遍历每个碱基。

我几年前开发了这个工具,在合理增加内存使用量的情况下,极大地加速了这项分析的性能。原理其实很简单,先来看图。

我们将每条read的start坐标作为key,其出现的次数作为value,存入哈希hashStart;每条read的end坐标作为key,其出现的次数作为value,存入哈希hashEnd,如果在哈希中出现过这个key,则在其value值累加1,如果没有出现过这个key,则value为1。

然后将hashStart 和 hashEnd的所有keys 合并到一个集合UnionSet,并按坐标顺序排序。然后逐一遍历UnionSet的每个key,如果hashStart中出现了这个key 则depth 加上对应的value,如果hashEnd中出现了这个key,则depth 减去对应的value。

我们对sorted bam文件按每条染色体逐条处理。就可以计算出每条染色体的depth分布情况了。

对于WES或者target Sequencing的数据,则只需要 将targetRegion.bed与以上计算出来的depth.bed取overlap,就可得到target区域的 depth分布情况了,这样还能统计 off-target的depth或者percentage分布。

受别人的项目启发( https://github.com/brentp/mosdepth ),再不开源,别人还以为是他自己的原创了。这个工具(或相似工具)本早在三四年前就已实现,现本人已将此工具开源,欢迎star

https://github.com/xiongxu/bam2bedGraph2

  • 发表于:
  • 原文链接http://kuaibao.qq.com/s/20171222G0YS9100?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 yunjia_community@tencent.com 删除。

扫码关注云+社区

领取腾讯云代金券