ちょっと調べたらいけそうな気がしてきた。Biopythonを使うといろいろ簡単にできる模様。
まずはmultifastaを開いて配列を順番に読み込む
import sys from Bio import SeqIO fasta_file = sys.argv[1] for record in SeqIO.parse(fasta_file, 'fasta'): ids = record.id desc = record.description seq = record.seq print('id:', ids) print('desc:', desc) print('seq:', seq)
とりあえずこんなスクリプトをfasta_in.pyという名前で作れば
$ python3 fasta_in.py <fasta_file.fasta>
という感じに投げるとidと説明書きとシークエンスを取り出すことができる。
bioinformatics - Python find longest ORF in DNA sequence - Stack Overflow
そんでもって最長ORFはここに出ている例を使って取り出せるので
import sys from Bio import SeqIO import re fasta_file = sys.argv[1] for record in SeqIO.parse(fasta_file, 'fasta'): for strand, seq in (1, record.seq), (-1, record.seq.reverse_complement()): for frame in range(3): index = frame while index < len(record) - 6: match = re.match('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(seq[index:])) if match: orf = match.group() index += len(orf) if len(orf) > 1300: pos = str(record.seq).find(orf) + 1 print(">{}, pos {}, length {}, strand {}, frame {}".format\ (record.id, pos, len(orf), strand, frame )) print(orf) else: index += 3
これでいいかな?と思ったが、なんかおかしい。
そもそもcDNAのfastaを見ているのでstrandは1方向固定でいいから
for strand, seq in (1, record.seq), (-1, record.seq.reverse_complement()):
このforループはいらないな。
matchじゃなくてfindallのほうがいいかな?
もうちょっと調べる。
import sys from Bio import SeqIO import re fasta_file = sys.argv[1] for record in SeqIO.parse(fasta_file, 'fasta'): match = max(re.findall('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(record.seq)), key = len, default=0) if match: print(">" + record.id) print(match)
これで最長ORFを取り出せた。
あとはこれをアミノ酸に置き換える。
#!/usr/bin/python3 import sys, re from Bio import SeqIO from Bio.Alphabet import IUPAC from Bio.Seq import Seq fasta_file = sys.argv[1] for record in SeqIO.parse(fasta_file, 'fasta'): match = max(re.findall('(ATG(?:\S{3})*?T(?:AG|AA|GA))', str(record.seq)), key = len, default=0) if match: seq = Seq(match, IUPAC.ambiguous_dna) print(">" + record.id) print(seq.translate())
できた。これで一気に処理できる。
問題があるとすると、stopコドンがなくて尻切れトンボのORFを持つcDNAからはアミノ酸を読み出せないことかな。