Повторный доступ к БОЛЬШИМ файлам fasta. Самый запоминающийся метод?

Я использую Biopython, чтобы открыть большой файл fasta с одной записью (514 мегабаз), чтобы я мог извлечь последовательность ДНК из определенных координат. Возвращать последовательность достаточно медленно, и мне просто интересно, есть ли более быстрый способ выполнить эту задачу, которую я не понял. Скорость не будет проблемой для одного или двух попаданий, но я перебираю список из 145 000 координат, и это займет несколько дней :/

import sys
from Bio import SeqIO
from Bio.Seq import Seq

def get_seq(fasta, cont_start, cont_end, strand):
  f = fasta
  start_pos = cont_start
  end_pos = cont_end
  for seq_record in SeqIO.parse(f, "fasta"):
   if strand == '-' :
    return seq_record.seq[int(start_pos):int(end_pos)].reverse_complement()
   elif strand == '+':
    return seq_record.seq[int(start_pos):int(end_pos)]
   else :
    print ' Invalid syntax! 
    sys.exit(1)

python performance biopython fasta dna-sequence
person user1995839    schedule 10.06.2013    source источник
comment
Вы хотите получить последовательность с этими координатами для всех записей в файле FASTA? Если это так, вы, вероятно, захотите использовать yield- это просто вернет часть последовательности для первого seq_record.   -  person David Cain    schedule 10.06.2013
comment
@JoachimIsaksson: нет, это правильный синтаксис. Второй параметр — это строка, определяющая формат ввода. (Кроме того, f=fasta здесь, так что это приведет к тому, что одно и то же значение будет передано дважды).   -  person David Cain    schedule 10.06.2013
comment
@davidcain Ах, упс, вот что, по-видимому, получается при чтении кода перед кофе по утрам :)   -  person Joachim Isaksson    schedule 10.06.2013

Ответы (1)


arrow_upward
2
arrow_downward

Ваша функция анализирует весь файл каждый раз, когда вы хотите найти одну цепочку. В этом нет необходимости — «лучшим способом» было бы разобрать все последовательности один раз и сохранить их в памяти для последующего доступа. Самый простой способ сделать это — просто преобразовать генератор, возвращаемый SeqIO.parse, в list или аналогичную структуру данных.

В качестве альтернативы вы можете сохранить проанализированные SeqRecord объекты в базе данных: ZODB, shelve или простое использование pickle может достичь этого.

Однако, судя по внешнему виду вашей функции, вы всегда просто возвращаете результаты из первого SeqRecord, найденного в файле (первая итерация через SeqIO.parse(f, "fasta") либо вернет, либо вызовет sys.exit(1)). Вы имеете в виду yield вместо этого (я предполагаю, что вы делаете?).

Вот как бы я справился с этим:

# Parse once, store in a list (alternatively, place DB "load" command here
all_seq_records = list(SeqIO.parse(f, "fasta"))

def get_seq(cont_start, cont_end, strand):
    assert strand in ["+", "-"], "Invalid strand parameter '%s'" % strand
    for seq_record in all_seq_records:
        segment = seq_record.seq[int(start_pos):int(end_pos)]
        yield segment if strand == "+" else segment.reverse_complement()
person David Cain    schedule 10.06.2013
comment
Спасибо за ваши комментарии и предложенный код Дэвид! Я анализирую только один файл fasta и извлекаю последовательность, указанную координатами, из анализа файла аннотаций gff. Итак, я вызываю get_seq, чтобы получить последовательность для начальной и конечной позиции. Я поиграю с вашим предложением. Там определенно есть что-то новое для меня, чтобы подумать :) - person user1995839; 10.06.2013