0

我正在使用 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)
4

1 回答 1

2

每次您想查找单链时，您的函数都会解析整个文件。没有必要这样做——“更好的方法”是一次性解析所有序列，然后将它们存储到内存中以供以后访问。最简单的方法是将返回的生成器转换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()
于 2013-06-10T04:52:30.240 回答