【直播】我的基因组（十七）:初步分析一下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作进更加准确的分析。