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
领取专属 10元无门槛券
私享最新 技术干货