【直播】我的基因组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咯。

http://bedtools.readthedocs.io/en/latest/content/tools/bamtofastq.html

https://gsl.hudsonalpha.org/information/software/bam2fastq

命令如下：

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、阿尔的太阳

图文编辑：吃瓜群众