protein-ligand 결합 구조와 도킹 구조 사이의 RMSD 계산

 특정 타깃에 대해서, 그 타깃과 결합력 및 결합 구조를 예측하는 방법으로 도킹 시뮬레이션이 있습니다.

 아무래도 예측이다보니, docking이 얼마나 잘되는지 의문이 들 수가 있습니다.

 그래서 보통 제일 먼저, 특정 PDB 구조에서 receptor와 ligand를 분리한 후, 그 ligand를 receptor에 붙이는 도킹 시뮬레이션을 돌려봐서 얼마나 잘 붙는지 확인을 합니다. 그리고, 같은 단백질에 대한 PDB들에 대해서, 서로 receptor와 ligand를 교환해서 docking 시뮬레이션을 진행합니다. (cross-docking)

ligand를 원래 붙어있던 receptor 구조에 도킹을 하는 경우는 잘 되는 경우가 많습니다. 하지만 실험 구조에도 오류가 있을 수 있고, 크러시가 있는 경우는 잘 안 붙기도 합니다.

cross-docking은 잘 안되는 경우도 자주 나옵니다. 이것은 induced fit 문제와도 연관되어 있습니다. ligand 와의 결합에서 receptor의 구조에도 변화가 생기는데, ligand에 따라서 적절한 receptor의 구조가 다르기 때문입니다.

 

 어쨌건 이런 테스트를 할 때, pymol로 그려서 눈으로 보면 얼마나 비슷하게 생겼는지 바로 직관적으로 알 수 있지만, 좀 더 정량적으로 스코어를 사용해서 보고 싶다고 한다면 일단 RMSD (root mean square deviation를 스코어로 사용할 수 있습니다.

$$RMSD = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(r_{i}^{mol}-r_{i}^{ref})^{2}}$$

하지만 이것을 계산하는 것이 생각보다는 쉽지 않습니다. 왜냐하면, 분자에는 내부적으로 여러 대칭성이 있기 때문입니다.

서로 동일한 분자에 대해서 RMSD를 계산해주는 경우라면, 각 원자들을 서로에 대해 매칭 시켜주는 작업이 필요하고, 매칭 한 후에는 둘에 대한 절대적인 좌표의 차이를 계산해줍니다. 그런데, 이 매칭 하는 작업은 단순하지 않습니다. 아래의 예시처럼 동일한 분자에 대해서, 좌표의 index가 다를 수 있습니다.

이런 경우라면 {1:1, 2:2,3:3,4:4,5:5,6:6,7:7} 과 {1:1,2:7,3:6,4:5,5:4,6:3,7:2}의 두 가지 매칭이 가능합니다. 따라서 이 두가지 경우 모두에 대해서 RMSD를 계산하여야 하고, 그중에서 더 작은 값을 선택해야 합니다.

 

저는 이전에는 Yang Zhang Lab의 DockRMSD 를 사용했는데 https://zhanggroup.org/DockRMSD/, mol2 파일만 지원하거나, 에러가 생겼을 때 왜 나오는지 확인이 어려워서 사용이 좀 불편해서 그냥 rdkit을 이용해서 코드를 새로 만들었습니다.

 수소원자는 RMSD 계산에서 제외하였습니다.

dock_rmsd.py

#!/usr/bin/env python
import sys
import numpy as np
from rdkit import Chem


def main():
    if len(sys.argv) < 3:
        print('dock_rmsd.py mol_pdb ref_pdb')
        sys.exit()

    m_file = sys.argv[1]
    ref_file = sys.argv[2]
    m = Chem.MolFromPDBFile(m_file)
    m_ref = Chem.MolFromPDBFile(ref_file)

    smarts = Chem.MolToSmarts(m_ref)

    patt = Chem.MolFromSmarts(smarts)
    cs_list = m.GetSubstructMatches(patt, uniquify=False)
    cs_ref = m_ref.GetSubstructMatch(patt)

    num_atoms = len(cs_ref)
    conformers = m.GetConformers()
    conformer_ref = m_ref.GetConformer()
    num_conf = len(conformers)
    num_hatoms = m.GetNumHeavyAtoms()
    for i_conf in range(num_conf):
        conformer = conformers[i_conf]
        rmsd_min = 99.999
        for cs in cs_list:
            tmp = 0.0
            for i in range(num_atoms):
                idx = cs[i]
                idx_ref = cs_ref[i]
                atom = m.GetAtomWithIdx(idx)
                if atom.GetSymbol() == 'H':
                    continue
                p = conformer.GetAtomPosition(idx)
                p_ref = conformer_ref.GetAtomPosition(idx_ref)
                pos = np.array([p.x, p.y, p.z])
                pos_ref = np.array([p_ref.x, p_ref.y, p_ref.z])
                tmp += ((pos-pos_ref)**2).sum()
            rmsd = np.sqrt(tmp/num_hatoms)
            if rmsd < rmsd_min:
                rmsd_min = rmsd
        print(i_conf, rmsd_min)


if __name__ == '__main__':
    main()

 

m.GetSubstructMatches(patt, uniquify=False)

에서 uniquify=False 옵션이 위에서 언급한 index 문제를 해결해줍니다.