AutoDock Vina 및 변종 사용법

(추가) AutoDock Vina에서 사용하는 pdbqt 확장자에 대한 설명은 다음 글을 참고하세요.

novelism.tistory.com/32

(추가) GUI를 이용한 도킹이나 실행 예제는 

https://novelism.tistory.com/259

https://novelism.tistory.com/260

https://novelism.tistory.com/261

 

를 참고하세요.

 

 

protein ligand docking 이란? 

 도킹 시뮬레이션의 1차적인 목표는 단백질-리간드의 결합 구조의 탐색 입니다. 

 결합 구조를 바꿔가면서 에너지 (혹은 스코어)를 계산해서 에너지가 최소화 되는 구조를 탐색하는 것입니다. 
에너지 모델을 사용하니 결합 에너지 (다른 표현으로 binding free energy, 결합력, affinity, binding scoring) 예측도 해줍니다. 

이론적으로는 프리에너지가 미니멈이 되는 것이 가장 적절한 구조이고, 에너지 모델이 완벽하다면 바인딩 프리 에너지의 예측과 구조 탐색이 별개가 아니겠지만, 현실적인 이유로 구조 예측과 에너지 시뮬레이션에 차이가 있습니다. 

 일단 정확한 에너지 모델은 계산량이 너무 높습니다. 구조 탐색에선 효율성을 중시하기 때문에 단순화한 이펙티브 에너지 모델을 사용하거나, 비 물리적인 에너지 모델들을 사용하기도 합니다. 하지만 너무 심하게 비현실적인 에너지 모델로 좋은 구조를 탐색할 수 없겠죠. 그래서 적당히 빠르면서 크게 문제 없는 에너지 모델을 사용하고, 도킹을 스코어링용으로도 사용할 수 있습니다. 

 

 

Vina 및 변종

AutoDock Vina 는 무료, 자유 소프트웨어이고 쉽게 사용할 수 있는 도킹 시뮬레이션 프로그램입니다. 자유 소프트웨어이다보니 변종도 여러가지 있습니다. 

 

 Quick Vina 2 (qvina2)

이것은 Vina와 사용법은 동일한데, 최적화를 (논문을 안읽어서 정확히 뭘 했는지 모르겠지만) 해서 Vina보다 빠릅니다. 특별히 탐색이 잘못되거나 하는 문제도 없는 것 같습니다. Vina 쓸거면 그냥 Quick Vina 2 를 사용하는것을 추천합니다.

 단, conda 에서 설치하는 qvina 는 이상하게 속도가 빠르지 않습니다. (vina와 동일한 속도)

github.com/QVina/qvina 에서 직접 다운로드 받아야 합니다.

컴파일 된 파일 (qvina02 )이 bin 안에 있으므로 직접 안해도 됩니다. 
백만개 이상 분자에 대해 docking을 돌리고 싶다면 qvina02 를 추천합니다. 

 

smina

사용이 편리하도록 변형한 버전입니다. vina 사용할 때 귀찮은 점이 파라미터가 너무 길고, box_size 집어넣기 귀찮다는건데, smina는 파라미터 이름을 단축했고 (예를들어 --receptor 대신 -r 사용가능) ref ligand 로 box_size를 계산할 수 있습니다. 대신 속도는 qvina에 비해 느립니다. 다양한 커스터마이즈가 가능합니다. 에너지 텀을 추가할 수도 있습니다.  openbabel 라이브러리를 활용해서 pdb to pdbqt 변환을 자동으로 해주기 때문에 입력, 출력 파일에 pdb 파일도 사용이 가능합니다. 또한 vina는 multi conformer 파일 입력이 불가능한데 (vina_split으로 파일을 쪼개야 함...) smina는 pdb 파일 한정으로 multi conformer 를 읽을 수 있습니다. pdbqt파일일 때는 안되네요. (--score_only 옵션 사용시 유용합니다. )

 

소스를 조금 고치면 더 재미있는 작업도 가능합니다. 비물리적인 방법이지만, 금속이온과의 결합이나 공유결합 같은 것도 사용할 수 있습니다. 

그냥 사용할것이라면 conda 로 다운받으면 되고, 소스를 고치거나 컴파일을 하고 싶다면 

sourceforge.net/projects/smina/ 에서 다운받으면 됩니다. git-hub 버전도 있지만, sourceforge 버전이 더 최신이고 cmake를 지원합니다. soruceforge 버전은 openbabel3 를 사용합니다. github 버전은 openbabel 2를 사용합니다.  대신 sourceforge 버전에 또 문제가 있는데, multi conformer 를 읽을 수 없습니다.

 

준비물 

receptor

일단 receptor 구조를 알아야 하고, receptor의 결합 위치 (pocket)를 알아야 합니다. 모르면 훨씬 광역적으로 탐색을 할 수 있지만 (quick vina W) 잘 안될 가능성이 있습니다. 분자가 단백질의 어느 위치에 결합하는가에 따라 역할이 달라지니 가능하면 결합 위치를 먼저 선택하는편이 좋습니다. 그리고 결합 전후로 단백질의 구조가 변하는데 (induced fit), 단백질은 unbound form보다 bound form이 좋습니다. 리간드의 종류에 따라서 induced fit의 형태도 달라진다면 훨씬 어려운 문제가 되겠네요. docking에선 ligand가 존재할 수 있는 영역을 box (center, size) 파라미터로 받습니다. 

 

receptor는 아마도 pdb 파일일테니, autodock tools 의 prepare_receptor4.py 로 pdbqt 파일로 변환합니다. 

autodock tools는 기본 python2 버전인데, github에  python 3버전도 있습니다. 
github.com/jaimergp/autodocktools-prepare-py3k

github.com/Valdes-Tresanco-MS/AutoDockTools_py3 
저는 위에거만 써봤습니다. 살짝 문제가 있어서 코드 한줄 고치긴 하는데...

아래거는 잘 되는지 모르겠네요.  

옵션은 다음과 같습니다. preperation도 포함되었지만, 이 과정은 다른 프로그램으로도 할 수도 있습니다. 

 수소가 없을 시 수소를 붙임 (-A hydrogens )

 극성 수소 이외 제거, lone pairs 제거, water 제거, altLoc 제거 (-U nphs_lps_waters_deleteAltB )
 prepare_receptor4.py -r receptor.pdb -A hydrogens -U nphs_lps_waters_deleteAltB -o receptor.pdbqt 

혹은 preperation이 이미 잘 된 파일이라면,
obabel receptor.pdb -O receptor.pdbqt
로도 변환 가능합니다. 

 ligand

 ligand가 conformer로 제공되어 있으면 거기서 출발해도 되지만 보통 분자 SMILES에서 출발할 것입니다.

 분자 SMILES가 주어지면 이성질체 (isomer) 정보가 충분한지 확인합니다. ( [C@H], [C@@H] 같은 것들...) 부족하다면 경우에 따라 전부 만들어서 시뮬레이션하는게 원칙입니다. 왜냐하면 isomer는 실제로 서로 다른 분자이고, 한 폼이 단백질에 결합할지라도, 다른 폼은 결합하지 않을 수 있기 때문입니다. 
/=/ 와 /=\ 는 trans, cis 를 구분하는 표기입니다. (이 정보가 빠졌을 정도면 좀 심각한게 아닌지...?) 
일단 저는 지금은 편의상 이성질체까지 따로 만들지는 않고 있습니다.

 

다음 과정은 openbabel 3 버전 이상을 사용하여 진행하였습니다. 
그 다음은 적절한 protonation state를 찾아야 합니다. 순서는 다음과 같습니다. 

 먼저 분자 SMILES에서 분자를 중성화합니다.  (nitro group처럼 안떨어지는 경우도 있습니다.) 

"SMILES"에 분자 SMILES 구조를 대입합니다. 

 

obabel -:"SMILES" -osmi --neutralize   (openbabel 3부터 지원)

다음은 pH 값에 적절한 protonation state를 찾습니다. 

obabel -:"SMILES" -osmi -p 7.4           ( pH=7.4인 경우) 
원래는 protonation state는 확률적으로 존재하는 것이니, 가능한 상태들에 대해서 전부 docking을 하는게 좋겠지만, docking 의 경우 어떻게 돌려도 큰 차이가 나지 않으므로 그냥 대표값 하나만 해도 무방합니다.

 

여담: obabel에 수소를 붙이는 옵션은 두가지가 있는데, 하나는 -p pH 이고, 다른 하나는 -h 입니다. -h의 경우는 3D conformer 입력에 대해서 implicit 인 수소를 explicit 하게 붙이라는 것이고, -p는 pH 상태를 고려하라는 의미입니다.  obabel에 버그가 많은데, conformer 파일의 경우 수소가 explicit하게 붙어있지 않으면 문제가 발생하는 경우가 많습니다. 가능한 implicit은 SMILES 단계에서만 사용하고 conformer 로 전환된 후부터는 꼭 explicit 모델로 가는게 좋습니다. 또한 protonation state도 SMILES 단계에서 끝내는게 좋습니다. 

 

 다음은 3D conformer 를 만들어야 합니다. 

obabel -:"SMILES" --gen3D -O ligand.pdb 

이러면 ligand.pdb 파일이 생성됩니다. (explicit hydrogen model)

이를 pdbqt 로 바꾸어야 합니다. vina에선 수소 중에 C에 붙은 수소는 제외하고 극성 (N과 O에 붙은 수소) 만을 취급합니다. 원자가 수소결합 도너인가 아닌가를 판단할 때 수소를 체크합니다. C는 당연히 도너가 아니므로 무시하는 듯...

옵션: 극성 수소 이외 제거, lone pairs 제거 -U nphs_lps 

 prepare_ligand4.py -l ligand.pdb -U nphs_lps -o ligand_0.pdbqt

혹은 obabel ligand.pdb -O ligand.pdbqt 

 

 

여기까지 되면 이제 준비가 되었습니다. 

 

docking 사용법

smina나 qvina는 vina의 옵션을 그대로 사용 가능하므로 vina 기준으로 설명합니다. 

qvina02 --receptor receptor.pdbqt --ligand ligand.pdbqt --center_x 10 --center_y 30 --center_z 10 --size_x 20 --size_y 20 --size_z 20 --num_modes 10 --exhaustiveness 10 --cpu 4  --out ligand_out.pdbqt
설명:

qvina02 : 사용할 프로그램 이름, qvina02, vina, smina 등..

--receptor: receptor 파일명

--ligand: ligand 파일명

--center_?: ligand box center의 위치 x, y, z 좌표

--size_?: ligand box size x, y, z 방향 적당히 20 근처로 잡았지만, 포켓 구조를 참고해서 옵티멀하게 잡는게 좋음

--num_modes : 아웃풋으로 출력할 최대 구조의 수, 탐색 옵션과는 무관하게 아웃풋만을 위한 옵션

--exhaustiveness: 탐색을 얼마나 할지 정하는 파라미터, 기본적으로 1단위이고, (아마도) 몽테카를로 샘플링을 하는데 그것을 몇번 평행하게 진행할지 정함 10이면 10단위의 몽테카를로 샘플링을 함, 정수값만 가능 (1단위에서 무엇을 얼마나 하는지는 모르겠음...)
--cpu: 사용할 코어의 수, 코어마다 별개의 몽테카를로 샘플링에 투입됨, 즉 exhaustiveness 의 수보다 높은 cpu 수는 사용하지 않음 

--out : output ligand 파일명

 

후처리

docking의 목적에 따라 후처리가 다름니다. 

만약, scoring 용으로만 쓴다면 그냥 결과에서 나온 에너지 값만 있으면 되는데,

docking 구조를 활용해서 뭔가 하려면 별개의 후처리 과정이 필요합니다. 

특히 pdbqt 파일은 완전한 explicit hydrogen model 이 아니고, 극성 수소만을 포함하기에 obabel로 읽을 때 에러가 나올 가능성이 높습니다. (부족한 수소를 죄다 라디칼 취급함, formal charge 가 표기 안되서 문제 발생 )

이때문에 pdb로 변환할 때 주의해야 합니다.

 

일단 obabel로 pdbqt를 pdb로 변환합니다. 

obabel ligand_out.pdbqt -O ligand_out.pdb
이렇게 변환해두면 ligand_out.pdb 파일은 수소가 explicit도 아니고 implicit도 아닌 애매한 모델이 됩니다. 
(부족한 수소를 죄다 라디칼 취급하고, formal charge 가 표기 안되서 문제 발생 합니다. )

그래서 여기에 다시 수소를 붙여어야 하지만, 그냥 붙이면 안됩니다. 

(obabel ligand_out.pdb -h -O ligand_out2.pdb 이렇게 바로 하면 안됩니다. )
일단 obabel 이 pdb 파일을 파싱할 때, ATOM의 charge와 CONECT와 수소의 수를 참고하는데, 이게 망가져있습니다.

예를들면 N은 3개의 결합을 가지고, O는 2개의 결합을 가져야 하는데, + 혹은 - 차지를 가지고 있으면  결합의 수가 다를 수 있음. 하지만 차지 정보가 없기 때문에 -h 옵션으로 수소를 붙여버리면 엉뚱한 수의 수소를 붙여버리게 됩니다.

혹은 수소를 다 떼고 다시 pH 계산을 하는데, 이것도 안됩니다. 

(obabel ligand_out.pdb -d -O ligand_out2.pdb     

obabel ligand_out2.pdb -p7.4 -O ligand_out3.pdb
이렇게 하면 안됩니다. 

)

 그리고 이미 N과 O에는 적절한 수의 수소가 붙어있는 상태이니 더 붙이거나 뗄 필요가 없습니다.

기본적으로 원자는 N, O, C 만 생각하면 됩니다.... (P, S 같은거 별로 생각할거 없고...)

여기서부터는 직접 코드 작성이 필요합니다. (한두개 파일이면 손으로 고쳐도 되고... )

ligand_out.pdb 파일을 열고

N과 O 원자들에 대해서 CONECT의 수를 조사합니다.

CONECT에서 단일 결합이면 숫자가 1개만, 이중결합 이상이면 2개 이상 찍혀있습니다.

CONECT 1 2 2 3

CONECT 1 4

이렇게 적혀있다면, 1번 원자와 2번원자는 2중결합, 1번 원자와 3번 원자는 단일결합, 1 번 원자와 4번 원자는 단일결합을 한 것입니다. (줄이 길어지면 줄바뀜 가능)

이 정보로부터 N에 결합의 수가 (2중결합은 2로 카운트, 3중 결합은 3으로 카운트) 3이면 charge는 0, 4이면 charge는 1+입니다.

이것에 맞게 ATOM 이나 HETATM 줄에서 차지 정보를 추가해줍니다. 

O의 경우 기본 결합은 2이니, 결합수가 2이면 charge는 0, 결합 수가 1이면 charge는 1- 임 (pdb 포멧이..-1이던가 1-던가... 잘 기억이 안나네요.)

이렇게 고쳐서 ligand_out2.pdb 파일에 저장합니다.

그리고 obabel ligand_out2.pdb -h -O ligand_out3.pdb 을 실행해주면 됩니다.