前言
如果对较大的bam文件进行操作时，往往相当耗时，所以将一个很大的bam文件分割为几个较小的bam文件，在并行进行分析可以大大缩短所需要的时间。通常来说分割的方法是用bamtools，bamtools可以很方便的将bam文件按照染色体分割为若干个小的bam，代码如下：
bamtools split -in test.bam -reference
其中test.bam就是所需要分割的文件，-reference是按照染色体进行分割。
但是问题来了，bamtools的操作是逐个进行，且相当缓慢，当一个很大的bam文件进行分割时，速度慢的令人发指。有没有比较快的方法呢？答案当然是有的，我们可以利用samtools进行分割，操作为：
samtools view -@ 12 -b test.bam chr1 > chr1.bam
方案
通过如上代码我们就可以将1号染色体的Bam给分割出来，且可以通过参数-@调整所需要的线程数。
但是按照上面的操作仍然有一个很大的问题，如果我不知道这个bam的文件具体有那些染色体，这些染色体又叫做什么名字，且需要一次性进行分割该如何操作呢？这里我给大家写了一个简单的python脚本，复制如下代码保存为splitbam.py，运行命令：python splitbam.py -input test.bam -p 12 -outdir ./即可。
其中-input为输入的bam文件，-p为指定使用的线程数，-outdir为输出分割文件的路径，如果不存在会自动创建，可选参数--index，可以为每个分割的bam文件创建索引。
注意：使用该工具前务必正确安装samtools和python
代码如下：
# -*- coding: utf-8 -*-
import os
import argparse
#############################################################################
#Author：Terry Xizzy txizzy#gmail.com
#Describtion: splitbam.py is a mulit processes tool for spilt the big bam file.
#Version: V1.0
#Data: 2021/10/19
#Example: python splitbam.py -input test.bam -p 12 -outdir ./
#############################################################################
parser=argparse.ArgumentParser(description='splitbam.py is a mulit processes tool for spilt the big bam file.')
parser.add_argument('-input',type=str,help='Input your bamfile',required=True)
parser.add_argument('-outdir',type=str,help='Input your output bamfile directory',required=True)
parser.add_argument('-p',type=int,help='Input the thread counts',default=12)
parser.add_argument('--index',action="store_true",help='Make index for splited bam files')
args=parser.parse_args()
bamfile = args.input
thread_count = args.p
outdir = args.outdir
#创建输出文件路径
if not os.path.exists(outdir):
os.makedirs(outdir)
print("Split bam file, it may takes some time……")
def split_bam():
split_cmd = "samtools view -H %s | cut -f 2 | grep SN | cut -f 2 -d \":\" > %s/chr.txt"%(bamfile,outdir)
os.system(split_cmd) #得到该bam文件的所有染色体号
with open(outdir+'/chr.txt','r') as chr_file:
chr_n = chr_file.readlines()
for item in chr_n:
item = item.rstrip('\n')
#分割bam，并对bam文件排序
os.system("samtools view -@ {tn} -b {bam} {chr} | samtools sort -@ {tn} -o {out}/{chr}.bam".format(tn=thread_count,bam=bamfile,chr=item,out=outdir))
if args.index:
os.system("samtools index {out}/{chr}.bam {out}/{chr}.bam.bai".format(out=outdir,chr=item))
os.remove(outdir+'/chr.txt')
split_bam()
print("All done!")
最后：如果想了解更多和生信或者精品咖啡有关的内容欢迎关注我的微信公众号：生信咖啡，更多精彩等你发现！
