如下图,是tophat2比对完之后的结果
其中,有aligned >1 times的比例。unique mapping 这个概念是由可能最早的reads对比工具Eland提出的。当时,在Eland的比对结果中,它会报告一个如'[UR][0-2]NM’的标签(tag),来指示这个map结果是uniquely, repetitively还是unmapped。从那以后,我们就开始谈论’mapping uniqueness’了。对于Eland,mapping uniqness是一个清晰的概念。对Eland来说,其比对结果是不允许gap的,如果一个reads在相同数量的mismatch时会有两个或者两个以上的位点都吻合,那么,它就是比对结果不唯一。
所以我们在比对完之后都需要抽出unique mqpped reads,至今(2020.3.6)学习过的RNA-seq、ChIP-seq,MNase-seq都需要抽取unique reads。比对软件会根据reads的比对情况赋予reads不同的tag。
1、bowtie2
XS:i: Alignment score for the best-scoring alignment found other than the alignment
reported. Can be negative. Can be greater than 0 in --local mode (but not in --end-to-endmode). Only present if the SAM record is for an aligned read and more than one alignment was found for the read. Note that, when the read is part of a concordantly-aligned pair, this score could be greater than AS:i.
方法:grep “AS:” aligned.sam | grep –v“XS:” >unique_alignments.sam #AS:Align≥1 time XS:出现第二个最佳比对的地方
2、tophat2/hasit2
"One possible solution is to use the next best alignment score (ZS:i:, previously XS:i:). If the next alignment score is the same as the primary alignment score (AS:i:), that means there are at least two equivalent alignments found by HISAT. This additional SAM field (ZS:i:) is only available at HISAT github, which I’ll include in the next release of HISAT."
NH:i:: N=1 表示为unique。常用于tophat/hisat2产生的sam文件unique read筛选。
ZS:i: : 如果有多余2个的比对,则会出现该标签。
方法:grep 'NH:i:1' out.sam >unique.sam
samtools view xxxx.bam | grep 'NH:i:1' >xxxx.bam
随着历史的发展,人们对于测序结果中每一个位点的质量也加入的比对,于是,即使相同的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
转载自:关于map当中的unique mapped reads问题