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

欢迎关注”生信修炼手册”!

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

1. samtools depth

命令如下

samtools depth input.bam > sample.depth

2. samtools mpileup

命令如下

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

3. bdetools genomecov

命令如下

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

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

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

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

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

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

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

第二列的flag含义如下

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

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

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

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

·end·

—如果喜欢&＃xff0c;快分享给你的朋友们吧—

原创不易&＃xff0c;欢迎收藏&＃xff0c;点赞&＃xff0c;转发&＃xff01;生信知识浩瀚如海&＃xff0c;在生信学习的道路上&＃xff0c;让我们一起并肩作战&＃xff01;

本公众号深耕耘生信领域多年&＃xff0c;具有丰富的数据分析经验&＃xff0c;致力于提供真正有价值的数据分析服务&＃xff0c;擅长个性化分析&＃xff0c;欢迎有需要的老师和同学前来咨询。

  更多精彩

• KEGG数据库&＃xff0c;除了pathway你还知道哪些

• 全网最完整的circos中文教程

• DNA甲基化数据分析专题

• 突变检测数据分析专题

• mRNA数据分析专题

• lncRNA数据分析专题

• circRNA数据分析专题

• miRNA数据分析专题

• 单细胞转录组数据分析专题

• chip_seq数据分析专题

• Hi-C数据分析专题

• HLA数据分析专题

• TCGA肿瘤数据分析专题

• 基因组组装数据分析专题

• CNV数据分析专题

• GWAS数据分析专题

• 2018年推文合集

• 2019年推文合集
  

  写在最后

转发本文至朋友圈&＃xff0c;后台私信截图即可加入生信交流群&＃xff0c;和小伙伴一起学习交流。

扫描下方二维码&＃xff0c;关注我们&＃xff0c;解锁更多精彩内容&＃xff01;

一个只分享干货的

生信公众号