生物信息学生物信息数据科学

64.《Bioinformatics Data Skills》之

2021-09-03  本文已影响0人  DataScience

一般来说除了查看文件,我们很少用到sam文件,大多数程序设计出来都是直接读取二进制的bam文件而不是sam文件字符串,例如samtools的子命令sortindexdepthmpilup等为了效率都要求bam文件作为输入。这一章节的内容主要涉及使用samtools工具对bam文件进行操作。

本节内容数据下载

sam与bam文件转换

使用samtools的子命令view结合-b参数可以将sam文件转换为bam输出:

$ samtools view -b celegans.sam > celegans.bam

bam文件也可以转回sam文件:

$ samtools view celegans.bam > celegans_copy.sam

然而如此转换后sam文件便无法再转回bam文件,因为头部信息已经丢失:

$ samtools view -b celegans_copy.sam > celegans_copy.bam
[E::sam_parse1] missing sam header
[W::sam_read1] parse error at line 1
[main_samview] truncated file.

因此将bam文件转为sam文件需要加参数-h代表包括头部内容:

$ samtools view -h celegans.bam > celegans_copy.sam
$ samtools view -h celegans.sam > celegans_copy.bam

bam文件排序

很多程序为了效率起见还会要求bam文件是排完序的(例如突变calling)。可以使用子命令sort 对 bam文件进行排序:

$ samtools sort celegans_unsorted.bam -o celegans_sorted.bam

面对大型的文件,排序操作将会耗费大量时间与算力。sort命令方便地提供了归并排序算法(即将文件拆分为多个临时文件最后合并)与并行线程,例如:

$ samtools sort -m 50G -@ 2 celegans_unsorted.bam -o celegans_sorted.bam

其中-m参数可以设置每个拆分文件的大小(单位为K,M或G),越大效率越高,-@接线程数目。这里由于文件太小,不会有明显的作用。

bam文件建立索引

整个bam文件可能非常大,如果我们只关注很小的一段区域而将整个序列都读进内存是非常低效的,为文件建立索引可以方便针对性的提取特定区域。bam文件建立索引与fasta文件类似,使用index子命令,要求输入为排序好的bam文件:

$ samtools index celegans_sort.bam
# 当前目录下会生成下面文件
$ ls celegans_sort.bam.bai
celegans_sort.bam.bai
上一篇下一篇

猜你喜欢

热点阅读