关于unique mapped reads

如下图,是tophat2比对完之后的结果


mapping rate.png
   其中,有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问题

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,504评论 6 496
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,434评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,089评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,378评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,472评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,506评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,519评论 3 413
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,292评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,738评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,022评论 2 329
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,194评论 1 342
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,873评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,536评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,162评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,413评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,075评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,080评论 2 352

推荐阅读更多精彩内容