PDB 形式のファイルの ATOM レコードまたは SEQRES レコードにはアミノ酸配列の情報が記載されている。このページでは、PDB ファイルから FASTA ファイルを作成する Python スクリプトの例を示す。
ディクショナリを利用した方法
ファイルを開いて ATOM 行を探し、Cα を見つけたら、そのアミノ酸の 3 文字コードを取得し、1 文字に変換する。最後に、それらを出力する。
import sys
argvs = sys.argv
pdb_file = '1ALK.pdb'
fa = {}
amino_code = {
'ALA':'A', 'ARG':'R', 'ASN':'N', 'ASP':'D',
'CYS':'C', 'GLN':'Q', 'GLU':'E', 'GLY':'G',
'ILE':'I', 'LEU':'L', 'LYS':'K', 'MET':'M',
'PHE':'F', 'PRO':'P', 'SER':'S', 'THR':'T',
'TRP':'W', 'TYR':'Y', 'VAL':'V', 'HIS':'H',
'ASX':'B', 'GLX':'Z', 'UNK':'K'
}
with open(pdb_file) as fh:
for buff in fh:
if (buff[0:4] != 'ATOM'):
continue
chain_name = buff[21:22]
res_number = int(buff[22:26])
amino_acid = buff[17:20]
if not (chain_name in fa):
fa[chain_name] = []
aa = 'X'
if (amino_acid in amino_code):
aa = amino_code[amino_acid]
if (len(fa[chain_name]) != res_number):
fa[chain_name] += ['X'] * (res_number - len(fa[chain_name]))
fa[chain_name][res_number - 1] = aa
for k, v in sorted(fa.items()):
print ('>' + k)
print (''.join(v))
PDB から 1ALK.pdb をダウンロードし、上のスクリプトを実行させると、FASTA 形式が画面に出力される。
BioPython を利用した方法
BioPython にある PDBParser などのモジュールを利用する場合、以下のように FASTA 形式のファイルを作成できる。
from Bio.PDB import *
amino_code = {
'ALA':'A', 'ARG':'R', 'ASN':'N', 'ASP':'D',
'CYS':'C', 'GLN':'Q', 'GLU':'E', 'GLY':'G',
'ILE':'I', 'LEU':'L', 'LYS':'K', 'MET':'M',
'PHE':'F', 'PRO':'P', 'SER':'S', 'THR':'T',
'TRP':'W', 'TYR':'Y', 'VAL':'V', 'HIS':'H',
'ASX':'B', 'GLX':'Z', 'UNK':'K'
}
p = PDBParser()
structure = p.get_structure('X', '1ALK.pdb')
for model in structure:
for chain in model:
chain_id = chain.get_id()
chain_seq = ''
for residue in chain:
aaa = residue.get_resname()
a = ''
if (aaa in amino_code):
a = amino_code[aaa]
chain_seq += a
print('>' + chain_id)
print(chain_seq)