Extract 16S rRNA gene from GenBank files?

Question

[MSRS]

Hi, I have found a number of posts about 16S gene extract from genank files. But those are not working till now. Looks error (IndentationError: expected an indented block )

Any updated code or programs available?

TIA

Answer 1

[Mensur Dlakic]

The code is correct, but is not properly indented. The problem is with these two lines:

for genome in gbank:
for gene in genome.features:

You can't have a for loop and not indent after that. So all the lines after for genome in gbank: need to be right-indented by 4 spaces:

from Bio import SeqIO, SeqFeature
import sys

gbank=SeqIO.parse(open(sys.argv[1],"rU"),"genbank")

for genome in gbank:
    for gene in genome.features:
        if gene.type=="rRNA":
            if 'product' in gene.qualifiers:
                if '16S' in gene.qualifiers['product'][0]:
                    start = gene.location.nofuzzy_start
                    end = gene.location.nofuzzy_end
                    if 'db_xref' in gene.qualifiers:
                        gi=[]
                        gi=str(gene.qualifiers['db_xref'])
                        gi=gi.split(":")[1]
                        gi=gi.split("'")[0]
                        print(">GeneId|%s|16S rRNA|%s\n%s" % (gi,genome.description,genome.seq[start:end]))
                    else:
                        print(">GeneId|NoGenID|16S rRNA|%s\n%s" % (genome.description,genome.seq[start:end]))

[MSRS]

Thank you so much.

Almost done.

Little error at line 20

print ">GeneId|%s|16S rRNA|%s\n%s" % (gi,genome.description,genome.seq[start:end])
          ^
SyntaxError: invalid syntax

[Mensur Dlakic]

There is no error at my end. Line 20 in the script is the very last line, not the one you listed above. I would make sure that everything is properly pasted.

If you are using python3, there should be parentheses in print lines, like this:

                print(">GeneId|%s|16S rRNA|%s\n%s" % (gi,genome.description,genome.seq[start:end]))
            else:
                print(">GeneId|NoGenID|16S rRNA|%s\n%s" % (genome.description,genome.seq[start:end]))

[MSRS]

Okay. Thank you so much for your valuable time. Clear to me.