关于map当中的unique mapped reads问题 18

当我们使用bwa或者bowtie map reads之后,当我们察看日志文件,我们会看到诸如:

20164115 reads; of these:
  20164115 (100.00%) were unpaired; of these:
    255188 (1.27%) aligned 0 times
    14142345 (70.14%) aligned exactly 1 time
    5766582 (28.60%) aligned >1 times
98.73% overall alignment rate

这样的报告。很多人会告诉我们,对于下游分析,需要去除map至多个位点的reads,只保留所谓的unique mapped reads。那么这个问题是怎么产生的,为什么要只保留unique mapped reads,现在人们又是如何来处理这个问题的呢?

unique mapping 这个概念是由可能最早的reads对比工具Eland提出的。当时,在Eland的比对结果中,它会报告一个如'[UR][0-2]NM’的标签(tag),来指示这个map结果是uniquely, repetitively还是unmapped。从那以后,我们就开始谈论’mapping uniqueness’了。对于Eland,mapping uniqness是一个清晰的概念。对Eland来说,其比对结果是不允许gap的,如果一个reads在相同数量的mismatch时会有两个或者两个以上的位点都吻合,那么,它就是比对结果不唯一。

那人们怎么只保留unique mapped的结果呢?对于bowtie1,我们可以简单的设置m 1这样的参数来对结果进行过滤。但是对于bowtie2以及bwa,并不存在这样的参数。人们就使用bwa输出中的XT:A:U(用于标记unique hit)以及tophat中的NH:i:1(用于标记unique hit)来进行过滤。

对于Single End的测序来说,这是很有效的,但对于pair end的测序,就不一定了。有一种情况是一端是unique hit,而另一端的mate不是(标记为XT:A:M),怎么处理呢?有人说,应该考虑,因为它们之间的间距可以用于筛选。也有人反对,认为直接一刀切比较稳妥。

随着历史的发展,人们对于测序结果中每一个位点的质量也加入的比对,于是,即使相同的mismatch时,一个质量高的碱基mismatch和一个质量低的碱基mismatch还是有区别的,在此情况下,我们完全可以相信测序质量低的那个错配位点是因为测序质量太低而引起的mismatch,而不是真实的mismatch.

而实际的情况是,无论是bwa还是bowtie它们的比对结果都是由特定的算法来给出比对评分,然后从中选择一个最佳的。但事实上很有可能这个评分高的并不一定是最佳匹配。再加上引入的插入,空隔,翻转等比对概念,对于mapping uniqueness的概念,就更加模糊了。再加上测序的长度增加了,得到相同评分这种情况更加不可能,大多数是一个相近的评分,这个时候如何判断它是否是mapping uniqueness呢?

所以现在人们已经开始避免使用unique mapped reads这个概念了,而转向使用mapq值来保留高质量的比对结果。因为mapq值反应了一组比对结果发生的可能性,MapQ = -10 log10(P), 比如结果为10,那就是1/10的概率会出现这个比对结果,如果我们认为0.05%是一个小概率的话,那个mapq值为15就可以用于筛选了, 如果认为0.01%是个小概率的话,mapq值为20就可以用于筛选了。但是人们往往从30这个值开始试起(1/1000的概率),如果它的筛选结果符合你的测序要求,就可以使用它了。如果不行,可以适当的调整这个筛选值。

但是科学是不允许比较合适这样的字眼的。有什么原则来确定这个筛选值呢?因为每个软件给出的mapq值不一定相等,所以得依照不同的软件进行调整。因为现在还没有明确的筛选标准,我个人理解为筛选后所得应该全部为unique mapped reads。但有些时候,在很高的筛选值也出现不unique mapped reads怎么办?这时候要分情况,如果是single-end的ChIP-seq或者RNA-seq,我们可以不理彩,如果是paired-end,对于SNP,INDEL的研究,最好还是去除,可以查看比例,如果比例低于1%,可以直接使用reads name来对其筛选,对于重复的reads name直接使用Picard MarkDuplicates去除就可以了。

下面是samtools对mapq值筛选的示例命令,其中-q参数就是依据mapq值筛选:

 samtools view -bhS -q 30 input.sam > output.bam

18 thoughts on “关于map当中的unique mapped reads问题

  1. Reply Puriney 10月 8,2013 12:25 下午

    STAR,那个是找最长可map区域,会把质量值低的碱基踢掉,然后把reads完整地打补丁粘起来

  2. Reply zhenhua 1月 28,2014 9:11 上午

    博主,我最近正处理Multimapping的问题,不知博主有没有看过哪个比对软件,可以给出每个reads比对到基因组后的位置数量(即下面的locations)?可以给出如下的输出结果?
    #ID locations chromosome strand start end count
    s3.25mer.txt-1 1 chr12 + 105579297 105579321 1
    s3.25mer.txt-4 1 chr8 + 95642182 95642206 1
    s3.25mer.txt-7 6 chr13 + 66975161 66975185 1
    s3.25mer.txt-7 6 chr13 – 72592620 72592644 1
    s3.25mer.txt-7 6 chr14 – 46332831 46332855 1
    s3.25mer.txt-7 6 chr19 – 32540873 32540897 1
    s3.25mer.txt-7 6 chr1 – 113777719 113777743 1
    s3.25mer.txt-7 6 chr2 + 70297183 70297207 1

  3. Reply xuhao 6月 18,2014 2:23 上午

    20164115 reads; of these:
    20164115 (100.00%) were unpaired; of these:
    255188 (1.27%) aligned 0 times
    14142345 (70.14%) aligned exactly 1 time
    5766582 (28.60%) aligned >1 times
    98.73% overall alignment rate

    请问哪里可以看到这样的信息

  4. Reply xuhao 6月 19,2014 12:40 上午

    哦,我是用的bwa aln,所以想询问一下如何才能知道多少reads是unique map的,谢谢

    • Reply admin 6月 23,2014 8:19 上午

      对于bwa,原作者的观点是不赞成uniquely mapped reads的概念。但是如果你一定要知道的话,你可以通过过滤reads的名称出现的次数来获得uniquely mapped reads数。

  5. Reply 彭勇 9月 25,2014 8:33 上午

    欧老师好,我想问您几个问题
    1. 只能用unique mapped reads去做peak calling 吗?
    2. 若只保留unique mapped reads, 除了你说的那种方法,samtools view -bq 1 file.bam > file_unique.bam 能达到这种效果么?
    3. peak calling 的软件,除了MAC2 和Homer, 还有哪些软件比较好用且在持续更新? 新开发的也可以。 peak calling的软件,很多两三年都没有人维护了。在持续维护的、且比较得到公认的peak calling toolkit 我就找到2个: MACS2 and Homer.

    • Reply admin 9月 27,2014 8:05 上午

      请不要叫我老师。大家在这里是探讨问题。闻道有先后,术业有专攻。

      1. 不一定。这需要看情况。我想我的博文中有一定这方面的探讨,但是并不深入。如果深入的话,可能需要大量的实验验证。

      2. 不能。只能说去除了大部分的非unique mapped reads。现在的流行做法是用q参数来过滤。因为大家的观点没有统一,有些人认为可以mapped到多个位置的reads,只要map quality足够好,还是可以保留的。

      3. 现在对于ChIP-seq或者说ChIP-chip来说,最标准的还是MACS。因为算法已经很完善了。并不是说软件越新越好。对于一些宽峰的侦测,倒是有些新的算法。还有clip-seq等,也有些新算法。

  6. Reply 彭勇 9月 25,2014 8:49 上午

    欧老师, 第二个问题我问错了。老师您说的MAPQ>=30的mapped reads, 是这样选出来么:samtools view -bq 30 file.bam > file_unique.bam
    是不是SAM format file 的第5行的值?

  7. Reply Snower 12月 4,2014 9:48 下午

    请问tophat2 中-g的参数跟unique map 有什么关系呢?看了一些说明没有看懂。。

    • Reply admin 12月 5,2014 10:46 上午

      这个参数就是指reads最多可以map到多少个位置。多于它就不报了。少于它报最好的一个,或者随机最好的一个。类似bowtie1里的-M加-k 1参数,但不同于-m参数。所以它并不是报unique mapped reads。理解起来有点绕。简单地讲,就是如果一个read可以map到多个位置,使用biowtie1中的-m 1参数,它是不汇报的。但是在bowtie2和tophat2当中,是报最好的一个。这种做法,同人们对unique mapped reads的理解的变化相关。

  8. Reply Rocky Chen 12月 23,2014 9:38 下午

    博主,您好。我有两个问题想请教您一下:
    1、关于MapQ,原文:MapQ = -10 log10(P), 比如结果为10,那就是1/10的概率会出现这个比对结果。那不是说MapQ的值越小,越好吗?根据samtools -q [int]参数说明,是去除掉小于[int]的值。 而且我发现我的结果里。MutliMap的read MapQ的值反而很小。我觉得MapQ = -10 log10(1-P) 不知道对不对。
    2、关于Mate 我一直不太理解,指的是双短测序中的另一端,还是该序列的互补配对的序列呢?如果是是互补配对的序列,为什么单端不用考虑呢?

    • Reply admin 12月 26,2014 10:10 上午

      你的这两个问题很好。
      1. MapQ是越大质量越好。如果1/10会出现这个比对结果,那mapQ值就是-10 * log10(.1) = 10, 如果是1/100的话,那mapQ的值就是-10 * log10(.01) = 20。你需要注意就是,这个概率,是指这一比较结果是随机行为的可能性,所以它的概率越低,说明它越不可能是随机的结果。
      2. Mate是指双端测序中的另一端。就是一对测序结果。

Leave a Reply to xuhao Cancel Reply

  

  

  

%d 博主赞过: