搜索
查看: 5636|回复: 3

【直播】我的基因组39:从bam中提取我们的原始测序数据

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-4 22:28:38 | 显示全部楼层 |阅读模式
【直播】我的基因组39:从bam中提取我们的原始测序数据
公司给了raw data, clean data,还有alignment的bam文件.在这之前我的博客提到,虽然公司给了比对好的bam文件,但我还是想要自己再比对一下,这就需要把fastq文件上传到服务器上。

我把55G的bam文件上传到服务器,因为网速太慢,不想再重新上传fastq文件了,所以打算从bam文件里面提取fastq文件,这可以节省很多时间。

老规矩:

一个是自己写脚本,就是重复造轮子;

一个是利用现成的工具。


自己写脚本,本质上,就是考验对sam/bam格式文件的理解能力。同样也可以锻炼我们对生物信息数据的处理能力。bam文件是sam文件的二进制,所以bam文件和sam文件内容本质上没有区别的。

我们前面之前也详细的讲解了sam文件的格式(【直播】我的基因组(十三):了解sam格式比对结果),sam格式的文件是要比原fastq文件要大的,因为sam不仅包含了fastq文件的信息更包含了比对的很丰富的信息。

它的第1,10,11列可以提取出来还原成我们的测序数据fastq格式。因此这就是我们能够从bam文件中提取fastq文件的基本原理。

由于我们的数据是配对reads[双端测序的PE reads],首先要把bam文件根据reads对的名字排序,现在如下图,同一个reads的第一列应该是一致的,但是下面的bam是按照染色体坐标排序,而双端的两条reads往往是比对到不同的位置的,这就把reads pair给分开了,所以我们要根据reads的名称重新排序。



samtools sort -@ 20 -n  -o  P_jmzeng.Nsort.bam P_jmzeng.final.bam
这一步非常耗时 , 而且文件会有所增加,55G变成了99G。


先别着急提取,细心的朋友可能会发现一个问题,就是在上面的图片中,序号尾号是39704的一对reads本应出现两次,但是图中为什么出现了三次呢??

由于数据处理太慢了  我在写文章的时候使用了  hisat2的比对结果作为范例,hisat2比对对于同一条reads是允许输出多个比对位置的。而bwa、bowtie等,对于reads比对到的多个位置也只允许输出一个最佳的位置。因此,不同的比对软件输出的结果是有差异的,大家需要注意。

然后我们就要动手处理数据了,稍微明白fastq格式和sam格式,都知道该怎么写了吧!脚本如下:


samtools view P_jmzeng.Nsort.bam | head | perl -alne 'BEGIN{open FH1,">1.fq";open FH2,">2.fq"}{if($.%2==0){print FH1 "$F[0]\n$F[9]\n+\n$F[10]" }else{ print FH2  "$F[0]\n$F[9]\n+\n$F[10]"}}'

我用了head命令测试脚本正确与否,你运行的时候去掉就可以啦!


需要注意的是,双端测序数据的sam,有些reads可能是缺失了配对序列,需要注意。然后就是有些sam格式,一条reads出现多条记录,在sam文件。



随意Google一下,就能拿到现成的工具。这里挑选大名鼎鼎的bedtools咯。
[url=]http://bedtools.readthedocs.io/en/latest/content/tools/bamtofastq.html[/url]
[url=]https://gsl.hudsonalpha.org/information/software/bam2fastq[/url]


命令如下:

bedtools=~/biosoft/bedtools/bedtools2/bin/bamToFastq

nohup time $bedtools -i control_1.Nsort.bam -fq tmp1.fq -fq2 tmp2.fq &

这一步也会非常耗时:
8.9亿的150bp长的reads的fastq文件,这个文件大小很容易就算出了,参加前面的帖子哈。【直播】我的基因组(四):计算资源的准备

到这里,我们就成功从bam中提取到了fastq文件,省去了很多上传时间。



文:Jimmy、阿尔的太阳

图文编辑:吃瓜群众



本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:【直播】我的基因组 38:我得了艾滋病?我是暴躁狂?
下一篇:【直播】我的基因组 40:不同lane的bam文件的比较
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

0

主题

15

帖子

88

积分

注册会员

Rank: 2

积分
88
发表于 2017-2-19 01:00:33 | 显示全部楼层
“samtools sort -@ 20 -n  -o  P_jmzeng.Nsort.bam P_jmzeng.final.bam
这一步非常耗时 , 而且文件会有所增加,55G变成了99G。”
给点小建议,samtools似乎没有多线程版本,可以用sambamba这个软件,操作和samtools差不错,有线程参数,可以省的bam split等操作了。。
回复 支持 反对

使用道具 举报

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
 楼主| 发表于 2017-2-19 19:41:18 | 显示全部楼层
lefroyqiu 发表于 2017-2-19 01:00
“samtools sort -@ 20 -n  -o  P_jmzeng.Nsort.bam P_jmzeng.final.bam
这一步非常耗时 , 而且文件会有所 ...

高手,你看的好快啊
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

0

主题

15

帖子

88

积分

注册会员

Rank: 2

积分
88
发表于 2017-2-19 20:57:01 | 显示全部楼层
zckoo007 发表于 2017-2-19 19:41
高手,你看的好快啊

谢谢版主回复!
老菜鸟一个。。出于好奇,昨天夜里看完了目前论坛和微信里的58条直播,确实学习到很多啊。
虽然是注重生物信息学注重方法的论坛,我其实更好奇的是生物学意义,即单个个体基因组能挖掘到什么程度的生物学意义或故事。例如直播中涉及到位点变异和潜在疾病表型的关联(直播36-38),利用PCA挖掘自己和已有人群数据的亲缘关系(直播55-58)等都很有意思~ 以上是给Jimmy大神生信技能树直播一点小建议哈。
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-9-16 15:39 , Processed in 0.047110 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.