前言
如果對較大的bam文件進行操作時梳凛，往往相當耗時险毁，所以將一個很大的bam文件分割為幾個較小的bam文件效床，在并行進行分析可以大大縮短所需要的時間身堡。通常來說分割的方法是用bamtools邓尤，bamtools可以很方便的將bam文件按照染色體分割為若干個小的bam，代碼如下：
bamtools split -in test.bam -reference
其中test.bam就是所需要分割的文件盾沫，-reference是按照染色體進行分割裁赠。
但是問題來了，bamtools的操作是逐個進行赴精，且相當緩慢佩捞，當一個很大的bam文件進行分割時，速度慢的令人發(fā)指蕾哟。有沒有比較快的方法呢一忱？答案當然是有的，我們可以利用samtools進行分割谭确，操作為：
samtools view -@ 12 -b test.bam chr1 > chr1.bam
方案
通過如上代碼我們就可以將1號染色體的Bam給分割出來帘营，且可以通過參數(shù)-@調(diào)整所需要的線程數(shù)。
但是按照上面的操作仍然有一個很大的問題逐哈，如果我不知道這個bam的文件具體有那些染色體芬迄，這些染色體又叫做什么名字，且需要一次性進行分割該如何操作呢昂秃？這里我給大家寫了一個簡單的python腳本禀梳，復(fù)制如下代碼保存為splitbam.py杜窄，運行命令：python splitbam.py -input test.bam -p 12 -outdir ./即可。
其中-input為輸入的bam文件算途，-p為指定使用的線程數(shù)塞耕，-outdir為輸出分割文件的路徑，如果不存在會自動創(chuàng)建嘴瓤，可選參數(shù)--index扫外，可以為每個分割的bam文件創(chuàng)建索引。
注意：使用該工具前務(wù)必正確安裝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
#創(chuàng)建輸出文件路徑
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!")
最后：如果想了解更多和生信或者精品咖啡有關(guān)的內(nèi)容歡迎關(guān)注我的微信公眾號：生信咖啡筛谚，更多精彩等你發(fā)現(xiàn)！
