0
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Can do it with biopython $\endgroup$ Jan 28 at 20:42
  • $\begingroup$ Could you please expand on your answer, Liam? $\endgroup$ Jan 28 at 20:45
1
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Jan 28 at 12:36
1
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.