如何从bam文件提取unique mapping reads呢？

  • Q1: 如何从比对结果bam文件里面提前unique mapping的reads呢？

  • Q2: 什么样的比对软件或者参数可以控制只输出unique mapping情况呢？

这类问题我已经被问烦了，不得不写文来说明这件事，而且我实在是不明白为什么一个简单的提取多比对情况的答案都搜索不到呢！！！

如果你的比对工具输出了多比对情况，也就意味着双端序列会输出大于2的记录，在bam文件里面，很简单写一个脚本查看几个这样的片段即可。

  1. cat  bwa.sam |perl -alne '{$h{$F[0]}++;}END{foreach(keys %h){print if $h{$_}>2}}'  | head

然后再从bam文件里面搜索其中一个，看看情况，如下：

  1. ERR1698194.3036751    83  3   5742854 60  21S80M  =   5742854 -80

  2. ERR1698194.3036751    163 3   5742854 60  28S73M  =   5742854 80

  3. ERR1698194.3036751    419 3   5742613 60  32M69H  =   5742854 321

很明显，一条片段应该是测了2条reads，所以比对情况应该是只有2个，但是这个片段有3个，说明发生了多比对情况。那么就需要仔细看，如何去区分它们！！！

发现第二列flags有一点奇怪，去 https://broadinstitute.github.io/picard/explain-flags.html 解释一下：

  1. Summary:

  2.    read paired (0x1)

  3.    read mapped in proper pair (0x2)

  4.    mate reverse strand (0x20)

  5.    second in pair (0x80)

  6.    not primary alignment (0x100)

又是很明显， (0x100) 代表着多比对情况，所以直接用 samtools view -f 0x100可以提取 multiple比对的 情况

不得不提的是 samtools flagstat 会输出bam文件的行数，所以多比对情况会造成比reads数量还要多的比对情况出来。

而不少比对工具，都是有参数可以控制不输出多比对情况的，这个时候就需要详细阅读比对工具的readme了。

当然，真实情况是，我也试图用关键词去谷歌搜索了，发现居然还真的就搜索不到，interesting！！！

赞 (0)

相关推荐