pdbqt에 대한 설명은 아래 글을 참고 하세요.
https://novelism.tistory.com/32
pdbqt를 pdb 로 문제 없이 복원하기 위하선 reference 를 사용하는 것이 제일 좋습니다.
reference를 사용하려면 3d conformer pdb 파일을 생성할 때부터 atom name 을 명확시해야 합니다.
다음 코드는 3d conformer를 생성한 후, pdb 파일에서 atom_name 을 수정하는 코드입니다.
pdb 파일포멧에서 하나의 ligand 나 residue에 속하는 모든 원자들은 독립된 이름을 배정 받습니다.
예를들면 " CA ", " CB ", " CG " 등이 있습니다.
그래서 이런 atom name을 key로 사용할 수 있습니다. 하지만 obabel로 생성한 구조는 이런 이름이 붙지 않고, 원자기호만 붙습니다. 그래서 아래 코 드에서 구조를 생성한 후 원자마다 독림된 이름을 배정합니다.
이 코드를 실행하기 위해선 openbabel 3 이상의 버전이 필요합니다.
import subprocess
def fix_ligand_atom_idx(line_list):
atom_dict = dict()
total_line_out = str()
for line in line_list:
line = line.rstrip('\n')
if line[0:6] == 'ATOM ' or line[0:6] == 'HETATM':
atom = line[12:16]
at = atom[0:2]
chain = line[21]
line = 'HETATM%s%s UNK %s 1 %s' % (line[6:12], at,
chain, line[27:])
if line[0:6] == 'HETATM':
atom = line[12:16]
at = atom[0:2]
if at not in atom_dict:
atom_dict[at] = 0
# chain = line[21]
atom_dict[at] += 1
idx = atom_dict[at]
line_out = line[0:12] + '%s%-2d' % (at, idx) + line[16:]
else:
line_out = line
total_line_out += line_out + '\n'
return total_line_out
def gen_3d_pdb(smi, ligand_file, timeout=10):
result = subprocess.run(run_line.split(), capture_output=True,
check=True, universal_newlines=True,
timeout=timeout)
result_lines = result.stdout.strip('\n').split('\n')
total_line_out = fix_ligand_atom_idx(result_lines)
fp = open(ligand_file, 'w')
fp.write(total_line_out)
fp.close()
smi = 'Cc1cccc(n1)c1nc(c2ccccc2n1)Nc1ccncc1'
ligand_file = 'M855.pdb'
gen_3d_pdb(smi, ligand_file, timeout=10)
아래는 생성된 pdb 파일을 pdbqt 로 변환하는 코드입니다.
AutoDock tools의 prepare_ligand4.py 를 사용합니다.
import subprocess
def pdb_to_pdbqt(pdb_file, pdbqt_file):
run_line = 'prepare_ligand4.py -l %s -o %s' % (pdb_file, pdbqt_file)
run_line += ' -U nphs_lps'
e = None
try:
subprocess.check_output(run_line.split(), stderr=subprocess.STDOUT,
universal_newlines=True)
except Exception as e:
return e
return e
pdb_file = 'M855.pdb'
pdbqt_file = 'M855.pdbqt'
pdb_to_pdbqt(pdb_file, pdbqt_file)
위에서 생성된 M855.pdbqt 파일을 M855.pdb 파일을 reference로이용해서 pdb 포멧으로 변환하는 코드입니다.
pdbqt는 bond 정보 (CONECT) 가 제거되었고, 원자의 순서도 pdb와 달라집니다. 그리고 극성 수소만을 남기고 나머지 수소는 제외됩니다.
bond length를 이용해서 bond 정보를 만들어내긴 하지만, aromatic ring에 대한 결합 정보를 잘못 만드는 경우도 있습니다.
그래서 원본 pdb 를 이용해서 bond 정보를 복원해내는 것이 안전합니다.
atom name 을 key로 사용하여 reference pdb 파일과 pdbqt 파일에서 동일한 원자를 찾습니다.
import subprocess
def read_ref_pdb_ligand(pdb_file):
fp = open(pdb_file)
lines = fp.readlines()
fp.close()
atom_dict = dict()
conect_dict = dict()
for line in lines:
if line[0:6] == 'HETATM':
atom_num = int(line[6:11])
# atom_name = line[12:16]
atom_dict[atom_num] = line
if line[0:6] == 'CONECT':
conect_list = []
for i in range(0, 8):
ini = i * 5 + 6
fin = (i + 1) * 5 + 6
atom_num = line[ini:fin].strip()
if len(atom_num) > 0:
conect_list += [int(atom_num)]
conect_idx = conect_list[0]
if conect_idx not in conect_dict:
conect_dict[conect_idx] = conect_list[1:]
else:
conect_dict[conect_idx] = conect_dict[conect_idx] + \
conect_list[1:]
return atom_dict, conect_dict
def read_pdbqt_file(pdbqt_file):
model_dict = dict()
model_num = 0
fp = open(pdbqt_file)
lines = fp.readlines()
fp.close()
for line in lines:
if line[0:6] == 'MODEL ':
model_num = int(line[6:].strip())
if model_num not in model_dict:
model_dict[model_num] = dict()
model_dict[model_num]['REMARK'] = list()
model_dict[model_num]['HETATM'] = dict()
if line[0:6] == 'REMARK':
model_dict[model_num]['REMARK'] += [line]
if line[0:6] == 'HETATM':
atom_name = line[12:16]
pos = line[30:54]
model_dict[model_num]['HETATM'][atom_name] = pos
return model_dict
def write_pdb_one_ref(model, ref_atom_dict, ref_conect_dict):
total_line_out = ''
remark_list = model['REMARK']
for line in remark_list:
total_line_out += line
coor_dict = model['HETATM']
total_atom_list = list()
keys = ref_atom_dict.keys()
for atom_num in keys:
atom_line = ref_atom_dict[atom_num]
atom_name = atom_line[12:16]
if atom_name in coor_dict:
total_atom_list += [atom_num]
line_out = '%s%s%s' % (
atom_line[:30], coor_dict[atom_name], atom_line[54:])
total_line_out += line_out
keys = ref_conect_dict.keys()
for atom_num in keys:
if atom_num not in total_atom_list:
continue
ans = ref_conect_dict[atom_num]
ans2 = list()
for an in ans:
if an in total_atom_list:
ans2 += [an]
num_conect = len(ans2)
line_out = ''
for i_con in range(num_conect):
if i_con % 4 == 0:
line_out += 'CONECT%5d' % (atom_num)
line_out += '%5d' % (ans2[i_con])
if i_con % 4 == 3:
line_out += '\n'
if len(line_out.strip()) < 1:
continue
if line_out[-1] != '\n':
line_out += '\n'
total_line_out += line_out
return total_line_out
def pdbqt_to_pdb_ref(input_pdbqt_file, output_pdb_file, ref_pdb_file):
ref_atom_dict, ref_conect_dict = read_ref_pdb_ligand(ref_pdb_file)
model_dict = read_pdbqt_file(input_pdbqt_file)
model_list = model_dict.keys()
num_model = len(model_list)
fp_out = open(output_pdb_file, 'w')
for model_id in model_list:
total_line_out = write_pdb_one_ref(
model_dict[model_id], ref_atom_dict, ref_conect_dict)
if num_model > 1:
line_out = 'MODEL %8d\n' % model_id
fp_out.write(line_out)
fp_out.write(total_line_out)
if num_model > 1:
line_out = 'ENDMDL\n'
fp_out.write(line_out)
line_out = 'END\n'
fp_out.write(line_out)
fp_out.close()
run_line = 'obabel %s -h -O %s' % (output_pdb_file, output_pdb_file)
subprocess.check_output(run_line.split(), stderr=subprocess.STDOUT,
universal_newlines=True)
return
input_pdbqt_file = 'M855.pdbqt'
output_pdb_file = 'M855_new.pdb'
ref_pdb_file = 'M855.pdb'
pdbqt_to_pdb_ref(input_pdbqt_file, output_pdb_file, ref_pdb_file)
'Drug > Computer-Aided Drug Discovery' 카테고리의 다른 글
| 인공지능 신약개발과 데이터의 문제 (6) | 2022.09.06 |
|---|---|
| 인공지능 신약개발 온라인 교육 플랫폼 LAIDD 2.0 open (0) | 2022.09.01 |
| 알파폴드가 예측한 단백질 구조는 완벽한가? (3) | 2022.08.05 |
| chimera GUI에서 docking (vina) 사용법: 3. chimera로 도킹하기 (0) | 2022.07.03 |
| chimera GUI에서 docking (vina) 사용법: 2. pdb 파일 구하기 (0) | 2022.07.03 |