Incomplete parsing of entire genbank file using python/biopython
Question: The main goal of my script is to convert a genbank file to a gtf file. My
problem pertains to extracting CDS information (gene, position (e.g., CDS
2598105..2598404), codon_start, protein_id, db_xref) from **_all_** CDS
entries. My script should open/parse a genbank file, extract information from
each CDS entry, and write the information to another file. The script produces
no errors, but only writes information from the first 1/2 of the genbank file
before terminating. Here is my code...

    
    
    import Bio
    from Bio import GenBank
    from Bio import SeqIO
    
    fileList = ['data_files/e_coli_ref_BA000007.2.gb']
    qualies = ['gene', 'protein_id', 'db_xref']
    
    #######################################################DEFINITIONS################################################################
    def strip_it(string_name):
        stripers = ['[', ']', '\'', '"']
        for s in stripers:
            string_name = string_name.replace(s, '')
            string_name = string_name.lstrip()
        return string_name
    
    def strip_it_attributes(string_name):
        stripers = ['[', ']', '\'', '"', '{', '}',',']
        for s in stripers:
            string_name = string_name.replace(s, '') 
            string_name = string_name.lstrip() 
            string_name = string_name.replace(': ', '=')
            string_name = string_name.replace(' ', ';')
        return string_name
    #---------------------------------------------------------------------------------------------------------------------------------
    
    #######################################################################################################################
    for f in fileList:
        nameOut = f.replace('gb', 'gtf')
    
        with open(f, 'r') as inputFile:
            with open(nameOut, 'w') as outputFile:
                record = next(SeqIO.parse(f, 'genbank'))
                seqid = record.id
                typeName = 'Gene'
                source = 'convert_gbToGFT.py'
                start_codon = 'NA'
                attribute = 'NA'    
    
                featureCount = 0
                for f in record.features:
                    print(f.type)
                    string = ''
                    if f.type == 'CDS':
                        dic = {}
                        CDS = record.features[featureCount]
    
                        position = strip_it(str(CDS.location))
                        start = position.split(':')[0]
                        stop = position.split(':')[1].split('(')[0]
                        strand = position.split(':')[1].split('(')[1].replace(')', '')
                        score = '.'
    
                        for q in qualies:
                            if q in CDS.qualifiers:
                                if q not in dic:
                                    dic[q] = ''
                                dic[q] = strip_it(str(CDS.qualifiers[q]))
    
                        attribute = strip_it_attributes(str(dic))
    
                        if 'codon_start' in CDS.qualifiers:
                            start_codon = str(int(str(CDS.qualifiers['codon_start'][0]))-1) #need string when finished so it can be added to variable 'string'
    
                        string = '\t'.join([seqid, source, typeName, start, stop, score, strand, start_codon, attribute])
                        if attribute.count(';') == 2:
                            outputFile.write(string + '\n')
    
                        featureCount+=1
    
    #---------------------------------------------------------------------------------------------------------------------------------
    

The last line of the output file is:

    
    
    BA000007.2  convert_gbToGFT.py   Gene  2598104  2598404  .  +  0  protein_i     d=BAB36052.1;db_xref=GI:13362097;gene=ECs2629
    

The location of gene ECs2629 appears on line 36094 in the genbank file, but
the total number of lines in this file is 73498. I have re-downloaded the file
multiple times to see if there was a downloading issue and I have visually
inspected the file (I find no fault with it). I have also tried this script on
another equally large genbank file and was met with identical issues.

Can anyone offer some suggestions as to why the entire genbank file is not
parsed, how I could modify my code to remove this issue, or point me to
another possible solution?

(you can see the format of a genbank file from here:
<http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html>), however, I am
working with an E. _coli_ genbank file (Escherichia _coli_ O157:H7 str. Sakai
DNA, complete genome) which can be found here:
<http://www.ncbi.nlm.nih.gov/nuccore/BA000007.2>

I am using the following: Centos 6.7, Python 3.4.3 :: Anaconda 2.3.0 (64-bit),
Biopython 1.66

[EDIT] @Gerrat suggestions worked for the file in question, but not for other
files. Using <http://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3> with the
suggested edit yields ~28 lines of output where my original code output 2084
lines (however, there should be 4332 lines of output).


Answer: Change this line:

    
    
    CDS = record.features[featureCount]
    

to:

    
    
    CDS = f
    

You're skipping records by accessing them via the `featureCount' index (since
there are probably 1/2 as many feature Counts as records).

**EDIT: To elaborate on your comment:**

Your original script is just wrong (w.r.t. the way you're using
`featureCount`). My correction is necessary. If you have further issues, there
is something else wrong. In this case, there appear to be 28 CDS records with
an attribute count of 2. (I know nothing about gene sequencing, I'm just going
by the variable names in the script). When you switch back to using
`featureCount`, you're now looking at records where the "type" is not "CDS".
It is "gene", or "repeat_region". You're checking the type of the record, `f`
to see if it is `CDS`, but then using a completely different record,
`record.features[featureCount]`. These don't refer to the same record (check
the CDS.type of this record - it's no longer "CDS" in most cases).

