Извлечение последовательностей CDS в Biopython

Всем добрый день,

Я начинаю программировать на Biopython и мне интересно, как извлечь последовательности генов и идентификаторы белков из файла GenBank генома (*.gb), имеющего координаты всех признаков.

Моя идея состоит в том, чтобы создать текстовый файл, содержащий идентификаторы белков, координаты генов и последовательности генов.

Если у вас есть какие-либо идеи, я был бы очень признателен за них.

Я пробовал это до сих пор:

for seq_record in seq_record.features: 
    if seq_record.type == 'CDS':
       x=seq_record.qualifiers['protein_id']
       i=seq_record.location._start.position
       f=seq_record.location._end.position
       sq = seq_record.seq
       FEAT_LIST.append('START END STRAND ID')
       FEAT_LIST.append(str(((i, f), s, x, sq)))
       print(FEAT_LIST)

Однако я получаю эту ошибку: sq = seq_record.seq AttributeError: объект «SeqFeature» не имеет атрибута «seq».

Спасибо за помощь.


biopython sequences
person user3579670    schedule 28.04.2014    source источник
comment
Добро пожаловать в StackOverflow! Что ты пробовал? Вы смотрели учебник/вики [Biopython](www.biopython.org)? Короткий ответ: сделайте синтаксический анализатор, если нет такого, который делает то, что вы хотите.   -  person llrs    schedule 28.04.2014

Ответы (3)


arrow_upward
5
arrow_downward

FeatureLocation имеет хороший метод extract, который берет родительскую последовательность и дает вам новый объект SeqRecord. Над этим объектом вы можете использовать обычный .seq, чтобы получить последовательность:

from Bio import SeqIO

for rec in SeqIO.parse("sequence.gb", "genbank"):
    if rec.features:
        for feature in rec.features:
            if feature.type == "CDS":
                print feature.location
                print feature.qualifiers["protein_id"]
                print feature.location.extract(rec).seq
person xbello    schedule 28.05.2014

arrow_upward
1
arrow_downward

Я бы порекомендовал вам изучить документацию Biopython для объектов SeqIO и SeqRecord, таких как parse и read. Формат genbank реализован в парсере, так что у вас не должно возникнуть проблем с чтением файла. Действительно, вам нужно только указать genbank в качестве параметра.

Здесь у вас даже есть пример чтения файлов генбанка.

Редактировать. Насколько я понимаю, у вас возникли проблемы при просмотре записей. Проблема, которую я вижу, заключается в том, что существует путаница между объектами SeqRecord и объектами SeqFeature. Вы не можете сделать:

for seq_record in seq_record.features:

потому что тогда seq_record является объектом SeqFeature, НЕ SeqRecord. Когда вы впервые анализируете файл GenBank, вы можете перебрать SeqRecord объектов:

for record in SeqIO.parse('my_file.gbk','genbank'):
    print "Record %s has %i features and sequence: %s" % (record.id, len(record.features), record.seq)

Каждый объект SeqRecord имеет атрибут seq и список объектов SeqFeature в атрибуте features. Если вы хотите что-то сделать с функциями, вам нужно перебирать их для каждой записи.

person cnluzon    schedule 28.04.2014
comment
Спасибо за информацию. Однако моя проблема возникает, когда я пытаюсь запустить этот код: for seq_record in seq_record.features: `if seq_record.type == 'CDS':` ` x=seq_record.qualifiers['protein_id']` ` i=seq_record.location._start.position` ` f=seq_record.location._end.position` ` sq = seq_record.seq` ` FEAT_LIST.append('START END STRAND ID')` ` FEAT_LIST.append(str(((i, f), s, x, sq) ))` print(FEAT_LIST) Я получаю эту ошибку: sq = seq_record.seq AttributeError: объект 'SeqFeature' не имеет атрибута 'seq' - person user3579670; 29.04.2014
comment
Я отредактировал свой ответ, чтобы попытаться быть более пояснительным. Я также думаю, что код, который вы разместили, должен быть виден в вашем вопросе. Опубликовал правку по этому поводу. - person cnluzon; 29.04.2014

arrow_upward
0
arrow_downward

Объект seqfeature имеет метод извлечения, который избавляет вас от необходимости копаться в FeatureLocation. Он возвращает новую последовательность.

from Bio.Seq import Seq
from Bio.Alphabet import generic_protein
from Bio.SeqFeature import SeqFeature, FeatureLocation
seq = Seq("MKQHKAMIVALIVICITAVVAAL", generic_protein)
f = SeqFeature(FeatureLocation(8, 15), type="domain")
f.extract(seq)
# returns: Seq('VALIVIC', ProteinAlphabet())
person cts    schedule 29.03.2018