pdbqt 는 Autodock 계열 docking 프로그램 (Autodock, Autodock Vina, Quick vina Smina 등)
에서 사용하는 파일 포멧입니다. 기본적으로 pdb와 유사합니다.
아스피린에 대해서 pdb와 pdbqt 파일을 비교해보겠습니다.
아스피린의 SMILES O=C(C)Oc1ccccc1C(=O)O 입니다.
obabel로 3d 구조를 생성해봅시다.
$ obabel -:'O=C(C)Oc1ccccc1C(=O)O' -O a.pdb --gen3d
$ cat a.py
COMPND UNNAMED
AUTHOR GENERATED BY OPEN BABEL 3.1.0
HETATM 1 O UNL 1 2.800 -1.310 -0.132 1.00 0.00 O
HETATM 2 C UNL 1 2.323 -0.194 -0.276 1.00 0.00 C
HETATM 3 C UNL 1 0.854 0.108 -0.238 1.00 0.00 C
HETATM 4 O UNL 1 3.032 1.002 -0.458 1.00 0.00 O
HETATM 5 C UNL 1 4.446 1.040 -0.563 1.00 0.00 C
HETATM 6 C UNL 1 5.155 -0.001 -1.172 1.00 0.00 C
HETATM 7 C UNL 1 6.544 0.075 -1.291 1.00 0.00 C
HETATM 8 C UNL 1 7.227 1.197 -0.829 1.00 0.00 C
HETATM 9 C UNL 1 6.520 2.248 -0.241 1.00 0.00 C
HETATM 10 C UNL 1 5.123 2.184 -0.109 1.00 0.00 C
HETATM 11 C UNL 1 4.462 3.349 0.538 1.00 0.00 C
HETATM 12 O UNL 1 5.037 4.071 1.335 1.00 0.00 O
HETATM 13 O UNL 1 3.203 3.592 0.117 1.00 0.00 O
HETATM 14 H UNL 1 0.294 -0.809 -0.447 1.00 0.00 H
HETATM 15 H UNL 1 0.597 0.837 -1.011 1.00 0.00 H
HETATM 16 H UNL 1 0.580 0.484 0.750 1.00 0.00 H
HETATM 17 H UNL 1 4.654 -0.886 -1.556 1.00 0.00 H
HETATM 18 H UNL 1 7.094 -0.746 -1.747 1.00 0.00 H
HETATM 19 H UNL 1 8.309 1.255 -0.924 1.00 0.00 H
HETATM 20 H UNL 1 7.065 3.122 0.114 1.00 0.00 H
HETATM 21 H UNL 1 2.889 2.966 -0.562 1.00 0.00 H
CONECT 1 2 2
CONECT 2 1 1 3 4
CONECT 3 2 14 15 16
CONECT 4 2 5
CONECT 5 4 6 10 10
CONECT 6 5 7 7 17
CONECT 7 6 6 8 18
CONECT 8 7 9 9 19
CONECT 9 8 8 10 20
CONECT 10 9 5 5 11
CONECT 11 10 12 12 13
CONECT 12 11 11
CONECT 13 11 21
CONECT 14 3
CONECT 15 3
CONECT 16 3
CONECT 17 6
CONECT 18 7
CONECT 19 8
CONECT 20 9
CONECT 21 13
MASTER 0 0 0 0 0 0 0 0 21 0 21 0
END
~~~~~~~~~~~~~~~~~~~~~~~~~~
pdb 파일포멧은 단백질의 구조 정보를 다루는데 최적화된 파일 포멧이지만, 리간드도 포함할 수 있습니다.
HETATM 라 적인 부분이 주로 리간드 부분입니다.
UNL 부분은 단백질이나 펩타이드인 경우는 아미노산의 three code 이지만, 리간드일 때는 리간드 이름을 적습니다. U
NL은 unnamed ligand 인 듯...
CONECT는 bond 정보를 포함합니다. pdb 사이트에서 제공하는 pdb 파일은 connectivity만을 줄 뿐, bond order 정보가 빠져있습니다. 빠진 정보는 구조로부터 파싱을 하지만, 애석하게도 많은 툴들이 분자를 정확히 파싱해주지 못합니다. 대신 pdb 사이트에서 리간드의 정확한 bond 정보를 mol 파일이나 cif 파일로 제공하니까 그것을 참고해야 합니다.
하지만 몇몇 프로그램들은 자체적으로 bond order를 표기하는 방법을 사용합니다.
single일 때는 숫자를 1번, double일 때는 숫자를 2번, triple일 때는 숫자를 3번 반복합니다.
이것이 공식적 규칙이 아니다보니 프로그램마다 다르기도 한데, 숫자가 4번 반복될 경우는 obabel은 4중 결합, pymol은 aromatic을 의미합니다.
저는 가능하다면 obabel 스타일로 pdb conect에 bond order를 (aromatic 대신) Kekule form 으로 적는 것을 추천합니다. obabel이나 rdkit등은 aromatic form을 Kekulize하는데, 자주 실패하는 문제가 발생하기 때문입니다. 그럴경우 분자 구조가 잘못 식별됩니다. fixer 만드는데도 상당히 고생했습니다. fixer는 나중에 논문 쓸 때 업로드할 예정입니다.
사실 저는 이 pdb 파일을 그냥 사용하지 않습니다. pdb에서 atom name은 4개의 문자로 표현할 수 있는데, 여기선 원자 기호만을 사용했습니다. 이렇게 사용할 경우, pdbqt로 변환한 후에 다시 pdb로 변환할 때 문제가 생길 가능성이 높습니다. pdbqt 가 bond 정보를 날려버리고, 그걸 다시 pdb로 변환할 때 완벽하게 복원에 실패하기 때문입니다. 이것도 주로 Kekulize 부분이 문제입니다.
그래서 저는 다음처럼 atom_name 에 index를 함께 표기해줍니다. index 를 붙이는 기준은 사람마다 다릅니다.
COMPND UNNAMED
AUTHOR GENERATED BY OPEN BABEL 3.1.0
HETATM 1 O1 UNL 1 2.791 0.042 1.326 1.00 0.00 O
HETATM 2 C1 UNL 1 2.314 -0.175 0.222 1.00 0.00 C
HETATM 3 C2 UNL 1 0.844 -0.156 -0.080 1.00 0.00 C
HETATM 4 O2 UNL 1 3.021 -0.434 -0.960 1.00 0.00 O
HETATM 5 C3 UNL 1 4.434 -0.540 -0.991 1.00 0.00 C
HETATM 6 C4 UNL 1 5.146 -1.082 0.084 1.00 0.00 C
HETATM 7 C5 UNL 1 6.534 -1.205 0.014 1.00 0.00 C
HETATM 8 C6 UNL 1 7.215 -0.812 -1.136 1.00 0.00 C
HETATM 9 C7 UNL 1 6.508 -0.292 -2.220 1.00 0.00 C
HETATM 10 C8 UNL 1 5.111 -0.158 -2.162 1.00 0.00 C
HETATM 11 C9 UNL 1 4.448 0.413 -3.367 1.00 0.00 C
HETATM 12 O3 UNL 1 5.022 1.162 -4.140 1.00 0.00 O
HETATM 13 O4 UNL 1 3.188 -0.021 -3.578 1.00 0.00 O
HETATM 14 H1 UNL 1 0.286 -0.309 0.849 1.00 0.00 H
HETATM 15 H2 UNL 1 0.587 -0.972 -0.760 1.00 0.00 H
HETATM 16 H3 UNL 1 0.570 0.806 -0.516 1.00 0.00 H
HETATM 17 H4 UNL 1 4.645 -1.411 0.992 1.00 0.00 H
HETATM 18 H5 UNL 1 7.085 -1.607 0.861 1.00 0.00 H
HETATM 19 H6 UNL 1 8.297 -0.911 -1.189 1.00 0.00 H
HETATM 20 H7 UNL 1 7.051 0.007 -3.116 1.00 0.00 H
HETATM 21 H8 UNL 1 2.875 -0.657 -2.908 1.00 0.00 H
CONECT 생략
어떤 사람은 원자 종류 상관 없이 1, 2, 3, 4 순으로 붙이지만, 저는 원자별로 O1, O2, O3..., C1, C2, C3... 같이 붙여주었습니다. 수소 원자의 경우 붙어있는 heavy atom 과 연계해서 붙일 수도 있지만, 코드가 복잡해질 것 같아서 대충 했습니다.
이제 pdb파일을 pdbqt로 변환해줍시다. Autodock tools 의 prepare_ligand4.py를 사용합니다.
prepare_ligand4.py -l a.pdb -o a.pdbqt
REMARK 4 active torsions:
REMARK status: ('A' for Active; 'I' for Inactive)
REMARK 1 A between atoms: C1_2 and O2_4
REMARK 2 A between atoms: O2_4 and C3_5
REMARK 3 A between atoms: C8_10 and C9_11
REMARK 4 A between atoms: C9_11 and O4_13
ROOT
HETATM 1 C3 UNL 1 4.434 -0.540 -0.991 1.00 0.00 0.096 A
HETATM 2 C4 UNL 1 5.146 -1.082 0.084 1.00 0.00 0.038 A
HETATM 3 C5 UNL 1 6.534 -1.205 0.014 1.00 0.00 0.003 A
HETATM 4 C6 UNL 1 7.215 -0.812 -1.136 1.00 0.00 0.001 A
HETATM 5 C7 UNL 1 6.508 -0.292 -2.220 1.00 0.00 0.021 A
HETATM 6 C8 UNL 1 5.111 -0.158 -2.162 1.00 0.00 0.071 A
ENDROOT
BRANCH 1 7
HETATM 7 O2 UNL 1 3.021 -0.434 -0.960 1.00 0.00 -0.278 OA
BRANCH 7 8
HETATM 8 C1 UNL 1 2.314 -0.175 0.222 1.00 0.00 0.260 C
HETATM 9 O1 UNL 1 2.791 0.042 1.326 1.00 0.00 -0.265 OA
HETATM 10 C2 UNL 1 0.844 -0.156 -0.080 1.00 0.00 0.126 C
ENDBRANCH 7 8
ENDBRANCH 1 7
BRANCH 6 11
HETATM 11 C9 UNL 1 4.448 0.413 -3.367 1.00 0.00 0.179 C
HETATM 12 O3 UNL 1 5.022 1.162 -4.140 1.00 0.00 -0.646 OA
BRANCH 11 13
HETATM 13 O4 UNL 1 3.188 -0.021 -3.578 1.00 0.00 -0.773 OA
HETATM 14 H8 UNL 1 2.875 -0.657 -2.908 1.00 0.00 0.167 HD
ENDBRANCH 11 13
ENDBRANCH 6 11
TORSDOF 4
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Autodock 의 conformer 탐색 방식은 기본적으로 torsion search 입니다.
분자의 bond 중에서 rotatable을 찾고, 회전하지 않는 영역들은 하나의 (강체)조각으로 취급합니다.
탐색 과정에서 bond distance와, bond angle 는 변하지 않고, 강체들 사이의 torsion 만을 변화시킵니다.
pdbqt 파일은 이런 탐색에 적합하도록 변형된 파일포멧입니다.
REMARK 부분에 active torsions 이 기록되어 있습니다.
또한 에너지 계산식에도 rotatable bond 의 수가 포함되어있습니다.
fragment 강체들은 root와 branch 로 묶여있습니다. 여기서, 원래 pdb에서의 인덱스가 다 날아가버리고, bond 정보도 날아가버리는데, 이를 다시 pdb로 바뀌어도 원래 인덱스로 돌아오지 않습니다.
또한 N, O 등에 붙은 수소만을 남기고 C에 붙은 수소들은 제거됩니다. partial charge 를 명시하지 않습니다.
protein 같은 경우는 flexiable docking이 아닌 한에는 강체 취급이라 branch가 없습니다.
이런 특징들 때문에, pdbqt 파일을 pdb로 변환할 때 잘못 변환되는 일이 많습니다.
그래서 도킹 후 pdbqt를 pdb로 변환할 때는, 도킹 전 원본 pdb 파일을 남겨두고, 이 파일의 인덱스가 포함된 atom_name 을 이용해서, 원본 pdb를 참고해서 변환을 하는것이 좋습니다.
특히 파싱에 자주 실패하는 경우는, (질소가 포함된) 2개의 aromatic ring이 붙어있는 경우, 질소가 포함된 5-member aromatic ring 입니다.
코드는 다 만들어두긴 했는데, 아직 논문 작업중이라 공개하기 어렵네요.
분자동역학 하는 분들은 대체 무슨 툴을 쓰시는지 모르겠네요. 도킹 초기구조 필요할텐데...