搜索
查看: 6595|回复: 0

怎么获取unique mapping read?

[复制链接]

365

主题

512

帖子

1713

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1713
发表于 2017-3-10 09:47:47 | 显示全部楼层 |阅读模式
这个非常简单,就是考验你对sam文件的理解,每个做alignment工具给出的sam有略微的区别:

在NGS数据中,Unique reads指的是只能map到一个位置的reads。要判断一条read是否unique,需要看SAM文件的optional fields。optional fields的格式是TAG:TYPE:VALUE,在SAM文件的第12列。多个optionalfields由空格分开。与unique read相关的主要TAG是“NH”,代表一条read贴在了基因组的多少个位置,数据类型为i,即整数,所以"NH:i:1"就代表此条read只贴到了基因组上的一个位置。(但align的目标未必总是基因组,也有可能是转录组,小RNA组等,取决于实验设计。)以"X"(或"Y","Z")开头的TAG是保留给用户的。如tophat结果的“XS:A:+”就代表产生这条reads的来源是正链转录本,“XS:A:-”则是负链转录本来源。BWA的结果比Tophat(hisat2)的结果在mapping信息上更为详细。在BWA结果中,"X0"代表打开gap的个数(BWA不支持splicejunction,但支持gap),数据类型为i。"XT"代表mapping类型,数据类型为"A",即可打印字符,其值可以是Unique/Repeat/N/Mate-sw,有"XT:A:U"就代表这一条是unique read。由于BWA还提供了"AS"(alignment score)和"XS"(suboptimalalignment score)这两项,所以即使有些reads没有被Uniquely mapped,只要AS的值和XS的值相差足够大,也有一定价值。
附:
1、Sam文件各标签含义(tophat/hisat2)
•NH:i:<N>: N=1时 为unique。常用于tophat/hisat2产生的sam文件unique read筛选。
•CC:Z: 当为‘=’为map到同一条基因上,一般在map基因组时由于内含子存在而容易出现,他只代表两种不同的方式,计数时应记为1。此处一般为其他基因的名字。CP:i 和HI:i标签为map到第i条基因及起始位置。
•YT:Z:<S>:代表的含义与bowtie产生的sam也不同。具体还未知!其他标签AS,XN,XM,XO,XG,NM,MD等如下图可以看出都相同。
•对于tophat/hisat2比对产生的sam文件我们可以直接筛选NH标签。
grep‘NH:i:1’ out.sam >unique.sam
2、Sam文件各标签含义(bowtie2)
•AS:i:<N>Alignmentscore.可以为负的,在local下可以为正的。 只有当Align≥1 time才出现
•XS:i:<N>Alignmentscorefor second-best alignment. 当Align>1 time出现
•YS:i:<N>Alignmentscorefor opposite mate in the paired-end alignment. 当该read是双末端测序中的另一条时出现
•XN:i:<N>Thenumber of ambiguous bases in the reference covering this alignment.(推测是指不知道错配发生在哪个位置,推测是针对于插入和缺失,待查证)
•XM:i:<N>错配碱基的数目
•XO:i:<N>Thenumberof gap opens(针对于比对中的插入和缺失)
•XG:i:<N>Thenumberof gap extensions(针对于比对中的插入和缺失延伸数目)
•NM:i:<N>Theeditdistance。(edits:插入/缺失/替换数目)
•YF:Z:<S>该reads被过滤掉的原因。可能为LN(错配数太多,待查证)、NS(read中包含N或者.)、SC(match bonus低于设定的阈值)、QC(failing quality control,待证)
•YT:Z:<S>值为UU表示不是pair中一部分、CP是pair且可以完美匹配、DP是pair但不能很好的匹配、UP是pair但是无法比对到参考序列上。
•MD:Z:<S>比对上的错配碱基的字符串表示。
由于bowtie2产生的sam文件并没有NH标签,所以提取uniqueread可能比较麻烦。首先提取“AS”标签表示能比对上的read(>=1 time),然后利用grep反正则表达式过滤掉XS标签得到我们需要的unique read。
grep “AS:” aligned.sam | grep –v“XS:” >unique_alignments.sam
对于双端测序用bowtie2比对筛选unique concordant pair时则需要在上一步的基础上增加如下命令:
grep ‘YT:Z:CP’ unique.sam>pair-end_unique.sam
3、Bwa获取unique
samtools view bwa.bam | grep "XT:A:U"



上一篇:血液肿瘤诊断和治疗那些事儿
下一篇:第六章 方差分析(生物统计学基础)
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 09:08 , Processed in 0.027229 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.