搜索
查看: 3125|回复: 4

【直播】我的基因组(十五):提取未比对的测序数据

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-2 10:08:51 | 显示全部楼层 |阅读模式
【直播】我的基因组(十五):提取未比对的测序数据

之前我们说了比对上的数据,那么会有人想到有没有没有比对上的数据呢?

既然是从我的血液里面提取到的DNA进行测序的,那么理论上测序仪出来的所有测序reads都应该是我的,也应该都可以比对到人类的参考基因组,但是实际过程中的确存在着未必对上的数据。

现有人类基因组毕竟只是个参考,也许我有某些独特的DNA序列呢?而且也不一定测序数据就都是人类的,也许我血液里面会有那么些微不纯粹呢?也有可能是某些片段变异的太多了,超过了常见的比对软件的承受能力,可以试用SHRIMP这个软件来提取一下。当然,这不是本讲的重点。



前面我们已经详细讲解了sam文件的格式,就是为了给这个做铺垫,如果还不清楚的,可以回过头再仔细阅读([url=]http://genome.sph.umich.edu/wiki/SAM[/url])。sam格式文件的第3和第7列,可以用来判断某条reads是否比对成功到了基因组的染色体,左右两条reads是否比对到同一条染色体。有两个方法可以提取未比对成功的测序数据,sam文件的第3列是*的(如果是PE数据,需要考虑第6,7列),或者sam文件的flag标签包含0x4的,代码如下:


[mw_shl_code=perl,true]samtools view -f4 sample.bam > sample.unmapped.sam

samtools view sample.bam |perl -alne '{print if $F[2] eq "*" or $F[5] eq "*" }' > sample.unmapped.sam[/mw_shl_code]


虽然上面两个方法得到的结果是一模一样的,但是这个perl脚本运行速度远远比不上上面的samtools自带的参数。

sam文件的说明书里面有这样一句话;https://samtools.github.io/hts-specs/SAMv1.pdf

An unmapped segment without coordinate has a ‘*’ at this field.  However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting. If RNAME is ‘*’, no assumptions can be made about POS and CIGAR.

(其实也不一定要自己写脚本,我们前面讲到的用来把巨大的bam文件按照染色体分割的小软件bamtools也可以完成这个需求,用 bamtools -split -in my.bam -mapped 即可!)


小写的f是提取,大写的F是过滤。因为我们测序数据的双端的,那么sam文件的第3列是reads1的比对情况,第6列是reads2的比对情况。所以未比对成功的测序数据可以分成3类,仅reads1,仅reads2,和两端reads都没有比对成功。

也可以用下面的代码分步提取这3类未比对成功的reads:

[mw_shl_code=perl,true]samtools view -u  -f 4 -F264 alignments.bam  > tmps1.bam samtools view -u -f 8 -F 260 alignments.bam  > tmps2.bam samtools view -u -f 12 -F 256 alignments.bam > tmps3.bam
samtools merge -u - tmps[123].bam | samtools sort -n - unmapped
bamToFastq -bam unmapped.bam -fq1 unmapped_reads1.fastq -fq2 unmapped_reads2.fastq[/mw_shl_code]

可以简单的统计一下未比对成功的reads有多少:
[mw_shl_code=perl,true]
cut -f 3,6 P_jmzeng.unmapped.sam |sort |uniq -c >unmapped.counts  [/mw_shl_code]

结果如下

ab17e6017863114387e4ef8199a6ea7b.png

如果对bamtools软件的结果来统计:

[mw_shl_code=perl,true]samtools view P_jmzeng.final.REF_unmapped.bam |cut -f 3,7 |sort |uniq -c >unmapped.counts[/mw_shl_code]
得到的结果只有7492014 *        * 说明它只考虑了PE reads均为比对成功的情况。

很奇怪,看起来我的未比对成功的测序数据里面竟然没有右端成功,而左端失败的情况,这个我没办法解释。 我也还需要学习才能搞明白这件事。
(其实之前我也搞错了,如果PE reads的左右两端均没有比对成功,那么第3,6,7列都是*,4,5,8,9都是0,第2列flag只有77,141这两种情况。(77代表PE,而且PE的两条reads都是unmanned的,141跟77一样,只是它们分别指代unmanned的的PE的reads的两端,结合[url=]https://broadinstitute.github.io/picard/explain-flags.html[/url]来理解)

0c5d8a53b92af4277dfaf3a7c17940d9.jpg

如果是左右两端reads只有一个比对成功,另一个reads没有比对上,如果是read1比对了,read2失败了,那么第3列应该是read1的染色体,第7列应该是*号表明read2没有比对成功。同理,如果read2比对成功,read1失败,按照道理,我们应该看得第7列有染色体,第3列是*号,但是我们在提取的unmapped文件里面,没有发现这种情况。

但其实不管是左端还是右端,第3列都是有染色体的,第7列是=号,但这并不能说明左端跟右端有着同样的比对结果。而第6列CIGAR是*,这个才是判断左右端是否匹配失败的标准。

sam文件的说明书里面有这样一句话;https://samtools.github.io/hts-specs/SAMv1.pdf
对于第6列CIGAR来说,An unmapped segment without coordinate has a ‘*’ at this field. However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting. If RNAME is ‘*’, no assumptions can be made about POS and CIGAR.
这也就是我为什么没有发现第7列有染色体,第3列是*号的reads。即使PE reads的右端匹配,左端未匹配,它只会把这个read比对的染色体写在第3列,而不是第7列!所以说要想探究它是左端还是右端未比对成功,得看flag。

da66ec1196a861291befeec0de35556a.png

这样就提取出来了未比对的测序数据,但是还需要做进一步分析,看看这些reads究竟是何方大神!具体要在第25讲之后了,敬请期待!


文:Jimmy、吃瓜群众

图文编辑:吃瓜群众






上一篇:必须要了解的美国国立综合癌症网络NCCN
下一篇:做一个癌症相关基因检测产品需要什么?
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-2-2 10:17:17 | 显示全部楼层
我看你更新的挺快的,不知道你看懂了木有呢?可以把你的疑问留在下面,也记录一下自己的学习经历
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
 楼主| 发表于 2017-2-2 11:27:04 | 显示全部楼层
Jimmy 发表于 2017-2-2 10:17
我看你更新的挺快的,不知道你看懂了木有呢?可以把你的疑问留在下面,也记录一下自己的学习经历 ...

主要看思路和逻辑,具体代码还是不懂,毕竟用python
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

633

主题

1189

帖子

4054

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4054
发表于 2017-2-2 11:32:43 | 显示全部楼层
zckoo007 发表于 2017-2-2 11:27
主要看思路和逻辑,具体代码还是不懂,毕竟用python

(⊙o⊙)…代码就不用看了,我的代码太丑了,思路和逻辑是重点
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
 楼主| 发表于 2017-2-2 11:44:46 | 显示全部楼层
Jimmy 发表于 2017-2-2 11:32
(⊙o⊙)…代码就不用看了,我的代码太丑了,思路和逻辑是重点

的确,你代码对我来说就是天书
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2020-3-31 10:34 , Processed in 0.030963 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.