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

Question

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

Я начинаю программировать на 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'

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

Answer 1

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

Answer 2

Я бы порекомендовал вам ознакомиться с документацией по Biopython для SeqIO а также SeqRecord объекты, такие как parse а также read, genbank Формат реализован в парсере, поэтому у вас не должно возникнуть проблем с чтением файла. На самом деле, вам нужно только указать 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 приписывать. Если вы хотите что-то сделать с функциями, вам нужно выполнить их итерацию для каждой записи.

Answer 3

У объекта seqfeature есть метод extract, который избавляет вас от необходимости копаться в 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())