反复访问 LARGE fasta 文件。最节省内存的方法？

Question

我正在使用 Biopython 打开一个大的单项 fasta 文件（514 兆碱基），这样我就可以从特定坐标中提取 DNA 序列。返回序列的速度相当慢，我只是想知道是否有更快的方法来执行我还没有想到的这项任务。速度不会只是一两次点击的问题，但我正在遍历 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)

Answer 1

您的函数在每次要查找单链时解析整个文件。没有必要这样做——“更好的方法”是一次性解析所有序列，并将它们存储到内存中以供以后访问。最简单的方法是将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()