Неполный анализ всего файла genbank с использованием python/biopython
Основная цель моего сценария - преобразовать файл genbank в файл gtf. Моя проблема связана с извлечением информации CDS (ген, позиция (например, CDS 2598105..2598404), codon_start, protein_id, db_xref) из всех записей CDS. Мой скрипт должен открывать / анализировать файл genbank, извлекать информацию из каждой записи CDS и записывать информацию в другой файл. Скрипт не выдает ошибок, а только записывает информацию из первой половины файла genbank перед завершением. Вот мой код...
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 отображается в строке 36094 в файле genbank, но общее количество строк в этом файле составляет 73498. Я повторно загрузил файл несколько раз, чтобы увидеть, была ли проблема с загрузкой, и я визуально проверил файл (Я не вижу в этом вины). Я также пробовал этот скрипт на другом, не менее большом файле genbank, и у меня возникли идентичные проблемы.
Может ли кто-нибудь предложить некоторые предложения относительно того, почему весь файл genbank не анализируется, как я могу изменить свой код для устранения этой проблемы или указать мне другое возможное решение?
(Вы можете увидеть формат файла genbank здесь: http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html), однако я работаю с файлом E.bank colin (Escherichia coli O157: H7 str. Sakai DNA, полный геном), который можно найти здесь: http://www.ncbi.nlm.nih.gov/nuccore/BA000007.2
Я использую следующее: Centos 6.7, Python 3.4.3:: Anaconda 2.3.0 (64-bit), Biopython 1.66
[EDIT] Предложения @Gerrat работали для рассматриваемого файла, но не для других файлов. Использование http://www.ncbi.nlm.nih.gov/nuccore/NC_000913.3 с предлагаемым редактированием дает ~28 строк вывода, где мой исходный код выводит 2084 строки (однако должно быть 4332 строки вывода).
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')
Измените эту строку:
CDS = record.features[featureCount]
чтобы:
CDS = f
Вы пропускаете записи, получая доступ к ним через индекс featureCount (поскольку количество объектов, вероятно, вдвое меньше, чем количество записей).
РЕДАКТИРОВАТЬ: Чтобы уточнить ваш комментарий:
Ваш оригинальный сценарий просто неверен (в зависимости от того, как вы используете featureCount). Мое исправление необходимо. Если у вас есть другие проблемы, значит что-то не так. В этом случае, кажется, есть 28 записей CDS с количеством атрибутов 2. (Я ничего не знаю о последовательности генов, я просто использую имена переменных в сценарии). Когда вы вернетесь к использованию featureCountТеперь вы просматриваете записи, в которых "тип" не является "CDS". Это "ген" или "repeat_region". Вы проверяете тип записи, f чтобы увидеть, если это CDS, но затем, используя совершенно другую запись, record.features[featureCount], Они не относятся к одной и той же записи (проверьте CDS.type этой записи - в большинстве случаев это больше не "CDS").
Из любопытства, что произойдет, если вы перебираете каждую строку, изменяя:
with open(f, 'r') as inputFile:
в
with open("file") as infile:
for line in infile:
do_something_with(line)
Также было бы интересно установить некоторую переменную в ноль, прежде чем перебирать строки в файле и делать variable += 1 каждый раз, чтобы увидеть, является ли номер строки, что вы ожидаете