【直播】我的基因组(十四):bam文件给按照染色体给分割成小文件
昨天,我们了解了一下SAM格式的比对结果,不知道大家理解的怎么样。但是全基因组测序数据实在是太大了,即使比对后把sam文件压缩成二进制的bam文件也还有55G(如何压缩转换可查看直播十二),如果完整的导入IGV查看会略微考验计算机配置。
如果按照染色体(chr1-chr22,chrX,chrY,chrMT)来分割写一个脚本其实很容易,无非是效率的高低而已。但是我Google了一下,发现有现成的工具,也顺便试用一下这个软件bamtools。
如果需要手动切割,用下面的脚本,其中$BAM是需要传进去的参数。
for chrom in `seq 1 22` X Y MT do samtools view -bh $BAM chr${chrom} | samtools sort - chr${chrom} samtools index chr${chrom}.bam done
如果需要使用现成的工具bamtools的话,该软件的github地址是: 。安装也是非常容易,因为没有二进制可执行版本,所以需要下载源码自己编译。
## Download and install variationtoolkit
## https://github.com/pezmaster31/bamtools/wiki/Building-and-installing
cd ~/biosoft
mkdir bamtools && cd bamtools
git clone git://github.com/pezmaster31/bamtools.git
cd bamtools
cmake --version ## BamTools requires CMake (version >= 2.6.4).
mkdir build && cd build
cmake ../
make
~/biosoft/bamtools/bamtools/bin/bamtools
与我以前安装的软件不太一样的是要先cmake然后再make,而且保证cmake的版本不低于2.6.4
用法非常简单:
bamtools split -in file.bam -reference
我的代码如下:
~/biosoft/bamtools/bamtools/bin/bamtools split \
-in /data/project/myGenome/bamFiles/P_jmzeng.final.bam \ -reference
## 这里指定按照reference来分离bam文件
还可以指定 -tag RG 来把这个bam文件按照原来的测序上样品的lane给分离开(因为本身测序文件就是多个,比对后merge的bam)
也可以指定-mapped来分离比对成功与否的bam文件!
默认split后的小bam文件,就在原来的大的bam文件目录下,这个55G的文件,运行了近8个小时。
上面的脚本也好,这个bamtools工具也好,都是一个个染色体依次运行,所以速度很慢,其实可以同时开25个文件句柄,一次读入,全部写出!!!
最后呢,留个问题给大家,对于PE reads,如果左端的reads比对到1号染色体, 但是右端比对到2号染色体,这个应该归于哪个染色体的比对情况呢?欢迎大家评论区留言!
文:Jimmy、吃瓜群众
图文编辑:吃瓜群众