I am trying to extract the protein sequence of specific genes from gene bank like format file obtained from antismash part of which looks like this

LOCUS       scaffold_10            47160 bp    DNA     linear   UNK 01-JAN-1980
DEFINITION  scaffold_10.
ACCESSION   scaffold_10
VERSION     scaffold_10
COMMENT     ##antiSMASH-Data-START##
            Version      :: 6.0.1-25e43a7(changed)
            Run date     :: 2021-11-05 17:39:45
            NOTE: This is a single cluster extracted from a larger record!
            Orig. start  :: 2049841
            Orig. end    :: 2097001
FEATURES             Location/Qualifiers
     protocluster    1..47160
                     [2070175:2070449](+), [2070507:2071973](+),
                     [2072023:2074178](+), [2074232:2074784](+),
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
     proto_core      join(20001..20287,20335..20608,20667..22132,22183..24337,
                     /detection_rule="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS
                     or hyb_KS or itr_KS or tra_KS))"
     cand_cluster    1..47160
                     /detection_rules="cds(PKS_AT and (PKS_KS or ene_KS or
                     mod_KS or hyb_KS or itr_KS or tra_KS))"
     region          1..47160
                     /rules="cds(PKS_AT and (PKS_KS or ene_KS or mod_KS or
                     hyb_KS or itr_KS or tra_KS))"
     CDS             1899..2720
                     /NRPS_PKS="Domain: PKS_KR (26-149). E-value: 7.6e-08.
                     Score: 24.2. Matches aSDomain:
                     /NRPS_PKS="type: other"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1001:short-chain dehydrogenase/reductase SDR (Score:
                     161.6; E-value: 3.9e-49)"
     aSDomain        1977..2345
                     /specificity="KR activity: inactive"
                     /specificity="KR stereochemistry: C1"
     aSModule        1977..2345
     CDS_motif       1980..2042
     CDS             join(3252..3263,3321..3443,3493..4428)
                     /NRPS_PKS="Domain: PKS_ER (28-277). E-value: 6.4e-13.
                     Score: 40.2. Matches aSDomain:
                     /NRPS_PKS="type: other"
                     /gene_functions="biosynthetic-additional (smcogs)
                     SMCOG1040:alcohol dehydrogenase (Score: 275.7; E-value:

Now from this I want to extract the protein sequence of gens having PKS domain (PKS_ER/PKS_AT/PKS_KS ...). So far I have tried the following biopython code but it is giving me error

from Bio import SeqIO
from Bio.SeqFeature import SeqFeature

gbk = "myfile_predicted_SMclusters.gbk"
fa = "myfile_PKS_protein.fa"
input_handle = open(gbk, "r")
output_handle = open(fa, "w")

for record in SeqIO.parse(input_handle, "genbank"):
    features = [feature for feature in record.features if feature.type == "CDS"]
    for feature in features:
        if feature.qualifiers["NRPS_PKS"][1] == "PKS":
        assert len(seq_feature.qualifiers['translation'])==1
        output_handle.write(">%s %s\n%s\n" % (


The error is as follows:

if feature.qualifiers["NRPS_PKS"][1] == "PKS":
KeyError: 'NRPS_PKS'

The link to the sample inputfile is below: sample my_file.gbk

I'm using python 3.8.8

Please help!!!


1 Answer 1


I see 2 problems.

First, you need to check if the qualifier your are searching for exists in the given CDS feature.

Second, the qualifier holds more then just the string PKS for example Domain: PKS_AT (537-846). E-value: 1.5e-98.

You could modify your code like this:

feature_type = "CDS"
feature_qualifier = "NRPS_PKS"
target = "PKS"

def has_pks_value(feature, qualifier, part):
    return any(part in target for target in feature.qualifiers[qualifier])

for record in SeqIO.parse(input_handle, "genbank"):
    for feature in record.features:
        if feature.type == feature_type and feature_qualifier in feature.qualifiers.keys() and has_pks_value(feature, feature_qualifier, target) and len(feature.qualifiers['translation'])==1:
            output_handle.write(f">{record.description} {feature.location.start} {feature.location.end} {feature.location.strand} \n {SeqFeature(feature.location).extract(record.seq)}")

I've created some variables which for your string. This might prevent you from making some typos.

Also instead of a list comprehension, I just run over all features and checked with a single if statement if the asked qualifier exists and contains the given value.

I also changed the string formatting to f-string. But I think you should create the Fasta Files using the SeqIO module of Biopython.

As result I got 37 sequences

Update - I added a missing method to check for the correct value of the qualifier

  • $\begingroup$ Thank you for your suggestions, But I god the following error now. output_handle.write(f">{record.description} {feature.location.start} {feature.location.end} {feature.location.strand} \n {SeqFeature(feature.location).extract(record.seq)}") ^ SyntaxError: invalid syntax could you please Help, I recently started using python (last 3 days) just for this purpose only trying to make sense of resources available online. $\endgroup$
    – ambika
    Commented Nov 15, 2021 at 21:08
  • $\begingroup$ Could you try run it without the f-string part? Just use your original version for the file writing $\endgroup$
    – Mr_Z
    Commented Nov 15, 2021 at 21:36
  • $\begingroup$ It has another error ` output_handle.write(">%s %s\n%s\n" % (record.description, feature.location, seq_feature.qualifiers['translation'][0])) NameError: name 'seq_feature' is not defined ` $\endgroup$
    – ambika
    Commented Nov 15, 2021 at 23:15
  • $\begingroup$ The variable seq_feature doesn't exist anymore. Use feature instead $\endgroup$
    – Mr_Z
    Commented Nov 16, 2021 at 18:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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