使用 python/biopython 对整个 genbank 文件的解析不完整

Question

我的脚本的主要目标是将 genbank 文件转换为 gtf 文件。我的问题涉及从 中提取 CDS 信息（基因、位置（例如，CDS 2598105..2598404）、codon_start、protein_id、db_xref）所有 个 CDS 条目。我的脚本应该 open/parse 一个 genbank 文件，从每个 CDS 条目中提取信息，并将信息写入另一个文件。该脚本不会产生任何错误，但只会在终止前从 genbank 文件的前 1/2 部分写入信息。这是我的代码...

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

#---------------------------------------------------------------------------------------------------------------------------------

输出文件的最后一行是：

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

基因ECs2629的位置出现在genbank文件的第36094行，但是这个文件的总行数是73498，我重新下载了好多次文件，看是不是下载问题，我已经目视检查了文件（我发现它没有问题）。我还在另一个同样大的 genbank 文件上尝试过这个脚本，但遇到了同样的问题。

谁能就整个 genbank 文件未被解析的原因、我如何修改我的代码以解决此问题或向我指出其他可能的解决方案提出一些建议？

（您可以从此处查看 genbank 文件的格式：http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html），但是，我正在使用大肠杆菌coli genbank 文件（大肠杆菌coli O157:H7 str. Sakai DNA, complete genome) which can be found here: http://www.ncbi.nlm.nih.gov/nuccore/BA000007.2

我正在使用以下内容： Centos 6.7，Python 3.4.3 :: Anaconda 2.3.0（64 位），Biopython 1.66

[编辑] @Gerrat 建议对有问题的文件有效，但对其他文件无效。使用 http://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3 和建议的编辑会产生约 28 行输出，其中我的原始代码输出 2084 行（但是，应该有 4332 行输出）。

Answer 1

出于好奇，如果您通过更改遍历每一行会发生什么：

with open(f, 'r') as inputFile:

至

with open("file") as infile:
    for line in infile:
        do_something_with(line)

在循环遍历文件中的行并每次执行 variable += 1 以查看行号是否符合您的预期之前，将一些变量设置为零也很有趣

Answer 2

更改此行：

CDS = record.features[featureCount]

至：

CDS = f

您正在通过“featureCount”索引访问记录，从而跳过这些记录（因为特征计数可能是记录的 1/2）。

编辑：详细说明您的评论：

您的原始脚本是错误的（w.r.t。您使用的方式 featureCount）。我的纠正是必要的。如果您还有其他问题，则还有其他问题。在这种情况下，似乎有 28 个 CDS 记录，属性计数为 2。（我对基因测序一无所知，我只是按照脚本中的变量名称进行操作）。当您切换回使用 featureCount 时，您现在正在查看 "type" 不是 "CDS" 的记录。它是 "gene"，或 "repeat_region"。您正在检查记录的类型 f 以查看它是否为 CDS，但随后使用了完全不同的记录 record.features[featureCount]。这些不引用相同的记录（检查此记录的 CDS.type - 在大多数情况下它不再是 "CDS"）。

Answer 3

感谢@Gerrat 的评论。我重新编写了脚本，它运行得很流畅。

import Bio 
from Bio import GenBank
from Bio import SeqIO

fileList = ['F1.gb', 'F2.gb']

for f in fileList:
      with open(f, 'rU') as handle:
         for record in SeqIO.parse(handle, 'genbank'):
            for feature in record.features:
               if feature.type=='CDS':
                  #[extract feature values here]
                  count+=1

   print('You parsed', count, 'CDS features')

使用 python/biopython 对整个 genbank 文件的解析不完整

Incomplete parsing of entire genbank file using python/biopython

python

biopython

genbank