【直播】我的基因组(17):初步分析一下multiple mapping 的情况

分析multiple mapping的情况,还是先用公司已经提供的bam文件来吧,后面我会用自己比对得到的bam文件再重新分析一次。

公司给我的bam数据是把有效测序数据通过 BWA(Li H et al.)比对到参考基因组 (b37),这是目前使用率很高的一个软件。在比对过程中,如果一个或一对 read(s) 在基因上可以有多个比对位置,BWA 的处理策略是从中选择一个最好的,如果有两个或以上最好的比对位置,则从中随机选择一个。这种多重比对(multiple hits)的处理对 SNP、INDEL 以及 CNV 等的检测有重要影响。

通常检测 SNP 或 INDEL 的时候要使用高质量的比对结果(alignment), 即比对质量值[MAPQ]大于0或更高。

前面我们详细讲了sam文件格式,在sam文件的第五列是MAPQ值,对某些软件,比如BWA,直接通过它就可以区分unique mapping和mutiple mapping的情况,就是判断这条reads是否只比对到参考基因组的唯一的位点。所以对于这种情况下的sam文件,提取unique mapping非常简单,用samtools view -q 1 (过滤掉MAPQ值低于1的情况)。

正好公司提供的bam文件就是这样的情况,因为他们在使用bwa这个软件来把测序数据比对到参考基因组的时候并没有加上-a这个参数,那么输出的sam文件里面,bwa会对每一个有multiple mapping情况的reads的MAPQ值设置为0,所以提前multiple mapping的reads是非常容易的。

但其它软件,其它参数产生的sam比对文件,如果我们要提取multiple mapping的reads就不那么容易了,因为MAPQ为0,并不一定就是multiple mapping。如下:

可以看到上面很多序列都是150M,但是mapping的quality都是0,也有很多其它比对序列也是quality为0,它们在IGV显示如下:

每一个矩形箭头都是一条比对好的reads,如果是白色的说明比对质量值为0,如果是灰色的,说明是正常的reads。

那么这个地方为什么会集中一大堆的multiple mapping的reads呢?

一些软件在根据bam文件来选择变异位点的时候会忽略掉这些mapping quality为0的reads。

也就是说对于这部分软件来说,这些mapping quality为0的reads是没有用的,相当于损失掉了,假设整体基因组的覆盖深度是很平均的,那这些MAPQ为0的位置的覆盖深度相当于降低了。这很有可能影响SNV位点的可信度。

那有multiple mapping情况的reads都集中在基因组的哪些区域呢?在哪些基因附近呢?

我们可以先用公司提供的bam文件提取出MAPQ为0的reads[里面包含大多数 multiple reads],看一下大致的分布,后面再用我自己比对得到的bam作进更加准确的分析。

 

请扫描以下二维码关注我们,获取直播系列的所有帖子!

菜鸟团公众号二维码

菜鸟团公众号二维码

Comments are closed.