ちょっと調べたらいけそうな気がしてきた。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からはアミノ酸を読み出せないことかな。
