0

特定の座標から DNA シーケンスを引き出すことができるように、Biopython を使用して大きな単一エントリの fasta ファイル (514 メガ ベース) を開きます。シーケンスを返すのはかなり遅いので、私が理解していないこのタスクを実行するより高速な方法があるかどうか疑問に思っています。1 回か 2 回のヒットであれば速度は問題になりませんが、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 に答える