# Splitting a GenBank file into smaller files

I have a genbank file that I have been trying to run in HMMER using the following command:

hmmscan --domtblout NC_000962.3.final.gbk.out -E 0.000001 prot_profile.hmm NC_000962.3.final.gbk


However, since the file is too big for HMMER, I get the following error:

Fatal exception (source file p7_pipeline.c, line 697):
Target sequence length > 100K, over comparison pipeline limit.
(Did you mean to use nhmmer/nhmmscan?)
Aborted (core dumped)


I need to split this file into smaller files for the HMMER run. The problem is, I would need the ORIGIN field (containing the nucleotide sequence) at the bottom of the file split into relevant chunks into the subfile as well. This is where I am getting stuck.

I have already looked at GenBank Parser and seqretsplit to no avail.

I would greatly appreciate some help with this. Thank you.

• Can do it with biopython – Liam McIntyre Jan 28 at 20:42

You can use NCBI EntrezDirect to download data in GenBank format for a specific region of the sequence as follows:

efetch -db nuccore -id NC_000962.3 -format gb -seq_start 1 -seq_stop 99000 > seq.gb

• The genome fragments from NCBI do not fit my purpose. The input file is obtained from another program and I need to split that file as mentioned in the question. – BhushanDhamale Jan 28 at 12:36

I don't have time to give you a polished fully functional version but this will get you started. The next step will be a bit harder as you will need to capture the start of the first CDS in the chunk and end of the last CDS instead of just hardcoding as I have here.

fout = open('/path/first3.gbk','w')
for record in SeqIO.parse('/path/NC_000962.3.final.gbk','genbank'):
for i, feature in enumerate(record.features):
new_record = record
new_record.seq = new_record.seq[0:4437]
new_record.features = new_record.features[:3]
if i > 2:
SeqIO.write(new_record, fout, 'genbank')
fout.close()
break