kuroの覚え書き

96の個人的覚え書き

multi FASTA (DNA)からmulti FASTA (Amino Acid)を機械的に作成する(その2)

ちょっと調べたらいけそうな気がしてきた。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)
    if match:
        print(">" + record.id)
        print(match)

これで最長ORFを取り出せた。
あとはこれをアミノ酸に置き換える。

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)
    if match:
        seq = Seq(match, IUPAC.ambiguous_dna)
        print(">" + record.id)
        print(seq.translate())

できた。これで一気に処理できる。
問題があるとすると、stopコドンがなくて尻切れトンボのORFを持つcDNAからはアミノ酸を読み出せないことかな。