64.《Bioinformatics Data Skills》之
2021-09-03 本文已影响0人
DataScience
一般来说除了查看文件,我们很少用到sam文件,大多数程序设计出来都是直接读取二进制的bam文件而不是sam文件字符串,例如samtools
的子命令sort
,index
,depth
,mpilup
等为了效率都要求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