PDBをFASTAに変換

PDB 形式のファイルの ATOM レコードまたは SEQRES レコードにはアミノ酸配列の情報が記載されている。このページでは、PDB ファイルから FASTA ファイルを作成するスクリプトの例を示す。

ディクショナリを利用した方法

ファイルを開いて ATOM 行を探し、Cα を見つけたら、そのアミノ酸の 3 文字コードを取得し、1 文字に変換する。最後に、それらを出力する。

import sys
argvs = sys.argv

pdb_file = '1ALK.pdb'  # PDBファイルの位置
fa = {}                # 読み込んだデータを保存

# 3文字アミノ酸コードを1文字アミノ酸コードに変換するディクショナリー
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'
}


# PDBファイルを1行ずつ処理する
with open(pdb_file) as fh:
	for buff in fh:
		# ATOM レコードでなければ次の行に進む 
		if (buff[0:4] != 'ATOM'):
			continue
		
		# ATOM レコードならばチェーン名、残基番号、3文字アミノ酸記号を取得
		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] = []

		# アミノ酸1文字表記をとりあえずXとする
		aa = 'X'
		# ディクショナリーから検索し、もしあれば3文字表記を1文字表記に変換
		if (amino_acid in amino_code):
			aa = amino_code[amino_acid]
		
		# 残基番号が飛び飛びになっていれば、その間をXで補完する
		if (len(fa[chain_name]) != res_number):
			fa[chain_name] += ['X'] * (res_number - len(fa[chain_name]))
		
		# 1文字表記のアミノ酸を変数に保存する
		fa[chain_name][res_number - 1] = aa

# 取得したアミノ酸配列をFASTA形式で画面に出力
for k, v in sorted(fa.items()):
	print ('>' + k)
	print (''.join(v))

PDB から 1ALK.pdb をダウンロードし、上のスクリプトを実行させると、FASTA 形式が画面に出力される。

ディクショナリを利用した方法 2

もう少し簡単なスクリプト。

import sys
argvs = sys.argv

pdb_file = argvs[1]
pdb_id = ''
chain_name = ''
chain_names = []
aa_index = 0
seq = ''

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:6] == 'HEADER'):
            pdb_id = buff[62:66]

        if (buff[0:4] != 'ATOM'):
            continue

        chain_name = buff[21:22]
        amino_acid = buff[17:20]
        res_number = int(buff[22:26])

        if not (chain_name in chain_names):
            print('\n>' + pdb_id + '|' + chain_name)
            chain_names.append(chain_name)

        aa = ''
        if (amino_acid in amino_code):
            aa = amino_code[amino_acid]

        if (aa_index != res_number):
            print(aa, end = '')
            aa_index = res_number

このスクリプトはターミナルで次のように実行する。PDB ファイルのパスはコマンド引数として与える。上のスクリプトを pdb2fa.py の名前で保存して以下のように実行する。

python3 pdb2fa.py 1ALK.pdb
## 
## >1ALK|A
## TPEMPVLENRAAQGNITAPGGARRLTGDQTAALRNSLSDKPAKNIILLIGDGMGDSEITAARNYAEGAGGFFKGIDALPLTGQYTHYALNKKTGKPDYVTDSAASATAWSTGVKTYNGALGVDIHEKDHPTILEMAKAAGLATGNVSTAELQDATPAALVAHVTSRKCYGPSATSQKCPGNALEKGGKGSITEQLLNARADVTLGGGAKTFAETATAGEWQGKTLREEAEARGYQLVSDAASLNSVTEANQQKPLLGLFADGNMPVRWLGPKATYHGNIDKPAVTCTPNPQRNDSVPTLAQMTDKAIELLSKNEKGFFLQVEGASIDKQDHAANPCGQIGETVDLDEAVQRALEFAKKEGNTLVIVTADHAHASQIVAPDTKAPGLTQALNTKDGAVMVMSYGNSEEDSQEHTGSQLRIAAYGPHAANVVGLTDQTDLFYTMKAALGLK
## >1ALK|B
## TPEMPVLENRAAQGNITAPGGARRLTGDQTAALRNSLSDKPAKNIILLIGDGMGDSEITAARNYAEGAGGFFKGIDALPLTGQYTHYALNKKTGKPDYVTDSAASATAWSTGVKTYNGALGVDIHEKDHPTILEMAKAAGLATGNVSTAELQDATPAALVAHVTSRKCYGPSATSQKCPGNALEKGGKGSITEQLLNARADVTLGGGAKTFAETATAGEWQGKTLREEAEARGYQLVSDAASLNSVTEANQQKPLLGLFADGNMPVRWLGPKATYHGNIDKPAVTCTPNPQRNDSVPTLAQMTDKAIELLSKNEKGFFLQVEGASIDKQDHAANPCGQIGETVDLDEAVQRALEFAKKEGNTLVIVTADHAHASQIVAPDTKAPGLTQALNTKDGAVMVMSYGNSEEDSQEHTGSQLRIAAYGPHAANVVGLTDQTDLFYTMKAALGLK

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)