BiopythonでFASTAファイル操作 その２：reverse complementをとる

BiopythonのSeqIOを用いて、fastaファイル中の配列の相補配列を一気に書き出してみます。

まずは、標準出力に打ち出してみます。区別のため、def lineの頭に'RC_'が付くようにします。

# revComp_stdout.py

import sys
from Bio import SeqIO

mfa_in = sys.argv[1]

for record in SeqIO.parse(mfa_in, 'fasta'):
    rc_desc = 'RC_' + record.description # add 'RC_' to the def line
    rc_seq = record.seq.reverse_complement()
    print('>' + rc_desc)
    print(rc_seq)

結果をファイルに保存するには、リダイレクトを使います。

python revComp_stdout.py input_fasta > RC_result.mfa


これでも便利なのですが、毎回結果ファイルに名前を与えなければいけないので、'RC_'を入力ファイル名につけたものをカレントディレクトリに出力するようにしてみます。

# revComp_file.py

import sys
import os
from Bio import SeqIO

records = []
for record in SeqIO.parse(mfa_in, 'fasta'):
    record.id = 'RC_' + record.id
    record.description = record.id + ' ' + record.description.split(' ', 1)[1]
    record.seq = record.seq.reverse_complement()
    records.append(record)


rc_file = 'RC_' + os.path.basename(mfa_in) # remove upper path

SeqIO.write(records, rc_file, 'fasta')

ここでは、SeqIOの'O'のほうも使いました。SeqRecord objectをため込んで、一気に書き込んでいます。同名のファイルがある場合には上書きされてしまいますので、注意が必要です。
IUPACのmixed base一文字表記（R=A/G, M=A/C, B=C/G/Tなど）やgap ('-')もうまく扱ってくれるので、便利です。