RDKitによる3次元構造の生成

化学

RDKitの分子Molオブジェクトを扱う」という記事では,RDKitの分子であるMolオブジェクトを用いて2次元構造の取得と描画方法について学んできました.今回はそれを3次元構造へと展開していきます.

ほぼ全ての分子は立体的な性質を持ちますから,分子を3次元で考えることは大切です.

必要なライブラリのインポート

まずはいつものように必要なライブラリーのimportをしておきましょう.

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole
import sys, py3Dmol
import pandas as pd
print(sys.version)
# 3.6.7 |Anaconda custom (64-bit)| (default, Oct 23 2018, 14:01:38) 
# [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]

print(rdBase.rdkitVersion)
# 2019.03.3

3次元構造の生成

AllChem.EmbedMolecule

Molオブジェクトの2次元構造はAllChem.Compute2DCoodinateで計算することができました.それでは3次元構造はというと,AllChem.EmbedMoleculeで行うことが可能です.この際,RDKitには以下に示す4つの方法が実装されています.2018.09バージョンからデフォルトの方法がETKDGに変更になっていることに注意してください.

方法 コメント
ディスタンス・ジオメトリー(DG)法 2018.03までのデフォルト方法
AllChem.KDG useBasicKnowledge=True
AllChem.ETDG useExpTorsionAnglePrefs=True
AllChem.ETKDG useBasicKnowledge=True
useExpTorsionAnglePrefs=True
2018.09以降のデフォルト方法
AllChem.ETKDGv2 新しいパラメータを用いたETKDGメソッド.python3.6以降で利用可能

立体配座生成の各アルゴリズムの説明

名前から想像できるように上記の3つは全てディスタンス・ジオメトリー(DG)法に改良を加えたものです.RDKitに実装されているDG法はStochasticな立体構造の生成方法としてはかなり精度がよいことがベンチマークによって示されていますが,芳香環が歪む傾向にあるなどといった明らかな難点も知られています.

そこでDG法で得た構造に対して,

  • ベンゼン環は平面
  • アルキンは直線

といった一般的な化学常識(basic knowledge)を加えたのがKDGになります.

またケンブリッジ結晶構造データベースの実験データから得た二面角の分布傾向(experimental torsion angle preferences)を考慮したものがETDGになります.そしてその両方を加えたのがETKDGです.

我々が普段使う分にはETKDGを使えばよいと思います.ETKDGについてさらに詳しく知りたい方は以下の論文を参照してください.

ディスタンス・ジオメトリー法で生成された分子は多くの場合,芳香環が歪んでいるなどあまり精度が高くよくありません.従ってDG法で得た構造は実際に分析を始める前に分子力学法(MM)で構造最適化する必要があります.分子力場としてはUFFMMFFの2つがRDKitでは提供されています.

ETKDGでも構造の補正を行いますが,DG法による構造生成とMMを用いた最適化にかかる時間を考慮すると計算コストはそこまで変わらないとも言われています.

立体配座生成の各アルゴリズムの精度比較

今回は「Platinum dataset」と呼ばれるハンブルク大の研究チームが高精度のX線構造から得た構造データセットを用いて,

  1. DG法で発生させた構造をUFFで最適化した構造
  2. ETDG法で発生させた構造
  3. ETKDGで発生させた構造
  4. ETKDGv2で発生させた構造

の4種類について,もとのX線構造と比較した精度の違いを調べてみます.

なおPlatinum datasetに関しての詳細は以下の論文を参照してください.

# 分子の読み込み
suppl = Chem.SDMolSupplier('./platinum_dataset_2017_01.sdf', removeHs=False)
mols = [x for x in suppl if x is not None]
len(mols) # 4548

# 分子を1度SMILESに変換して再構成
smiles = [Chem.MolToSmiles(m) for m in mols]
mols_from_sm = [Chem.MolFromSmiles(sm) for sm in smiles]

# それぞれの方法について4548分子の3次元構造を発生
def DG(mols):
    DG_uff = []
    for mol in mols:
        mh = Chem.AddHs(mol)
        AllChem.EmbedMolecule(mh, useBasicKnowledge=False, useExpTorsionAnglePrefs=False)
        if AllChem.UFFHasAllMoleculeParams(mh):
            AllChem.UFFOptimizeMolecule(mh)
            DG_uff.append(mh)
    return DG_uff

def ETDG(mols):
    ETDG_mols = []
    for mol in mols:
        mh = Chem.AddHs(mol)
        AllChem.EmbedMolecule(mh, AllChem.ETDG())
        ETDG_mols.append(mh)
    return ETDG_mols

def ETKDG(mols, version=1):
    ETKDG_mols = []
    for mol in mols:
        mh = Chem.AddHs(mol)
        if version == 1:
            p = AllChem.ETKDG()
        elif version == 2:
            p = AllChem.ETKDGv2()
        else:
            print('invalid input')
        AllChem.EmbedMolecule(mh, p)
        ETKDG_mols.append(mh)
    return ETKDG_mols

DG_m = DG(mols_from_sm)
ETDG_m = ETDG(mols_from_sm)
ETKDGv1_m = ETKDG(mols_from_sm)
ETKDGv2_m = ETKDG(mols_from_sm, 2)

# 結晶構造と発生させた構造のRMDSを取得
RMSD = []
for mol, dg, etdg, etkdg1, etkdg2 in zip(mols, DG_m, ETDG_m, ETKDGv1_m, ETKDGv2_m):
    DG_rms = AllChem.GetBestRMS(mol, dg)
    ETDG_rms = AllChem.GetBestRMS(mol, etdg)
    ETKDG1_rms = AllChem.GetBestRMS(mol, etkdg1)
    ETKDG2_rms = AllChem.GetBestRMS(mol, etkdg2)
    RMSD.append((DG_rms, ETDG_rms, ETKDG1_rms, ETKDG2_rms))
RMSD = pd.DataFrame(RMSD, columns=['DG+UFF', 'ETDG', 'ETKDG', 'ETKDGv2'])
RMSD.describe().round(2)
DG + UFF ETDG ETKDG ETKDGv2
count 4548.00 4548.00 4548.00 4548.00
mean 2.05 2.06 2.01 2.00
std 0.91 0.85 0.90 0.91
min 0.07 0.20 0.07 0.07
25% 1.40 1.46 1.39 1.37
50% 2.01 1.99 1.96 1.96
75% 2.62 2.60 2.60 2.58
max 5.96 5.69 5.89 6.00

上記の結果からRMDSの平均値はETKDGv2が最も小さいことがわかります.四分位点も総じて小さいですが,最大値はやや大きいことから,v2が最もよいとは簡単には結論付けられないかもしれません.DG+UFFと比べるとETKTGの2つは精度がよい方法と述べてよさそうです.

IPython上で%timeを用いて手許の計算機で実行時間(CPU times)を計測した結果が以下の表になりますが,DG+UFFとETKDG(及びv2)の実行速度は同程度と考えてよそうです.

method user sys total
DG + UFF 3min 9s 1.28 s 3min 11s
ETDG 2min 30s 518 ms 2min 31s
ETKDG 3min 2s 1.36 s 3min 3s
ETKDGv2 3min 5s 1.18 s 3min 6s

py3Dmolを用いた3次元構造の描画

続いて3次元構造を描画してみましょう.以前はPyMolを用いて描画するのが一般的だったと思いますが,py3Dmolの登場以降はJupyterノートブック上でこちらのライブラリを用いて描画する方が人気があるようです.

py3DMolをまだインストールしていない場合はpipでインストールしておきましょう.

pip install py3Dmol
py3Dmolを使って化学構造をJupyter上で美しく表示する」という記事でpy3Dmolの,「RDKitからPyMOLを利用する」という記事ではPyMOLの使い方を解説しています.参照してみてください.
py3Dmolを使って化学構造をJupyter上で美しく表示する
「RDKitによる3次元構造の生成」という記事でJupyterノートブック上で化合物の3次元構造を表示するにはpy3Dmolというライブラリを利用するのが流行っていると述べました.実際にpythonのケモインフォマティクス用ライブラリであるrdkitのIPythonConsole...
RDKitからPyMOLを利用する
本ブログではJupyter Notebook上で分子構造を描画するためのライブラリーとして「py3Dmolを使って化学構造をJupyter上で美しく表示する」という記事でpy3Dmolについて説明しました.一方で生命科学の分野で,分子構造を可視化する際に使われるポピュラーなソフト...

分子1つの描画

IPythonConsole.drawMol3D(mol, view)

1つの分子を描画するにはdrawMol3Dを使うのが簡単です.描画したい分子のみを指定した場合はpy3DmolのViewが作成されて,そこに分子が表示されます.

IPythonConsole.drawMol3D(mols[0])

複数分子を1つの枠に描画

py3Dmol.view.addModel()

下のコードはpy3Dmolのviewオブジェクトを1つ作成して,そこに描画する分子を追加しています.その際にviewオブジェクトにはMolBlock形式PDB形式で分子を渡す必要があるため,一度MolBlock形式に変換しています.IPythonConsoleを用いるとこの作業を内部でやってくれるわけです.

v = py3Dmol.view(width=400, height=400)
mols_show = [DG_m[30], ETDG_m[30], ETKDGv1_m[30], ETKDGv2_m[30]]
for mol in mols_show:
    mb = Chem.MolToMolBlock(mol)
    v.addModel(mb, 'sdf')
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick': {}})
v.zoomTo()
v.show()

複数分子を違う枠に分けて描画

py3Dmol.view(viewergrid)

今度ののコードではviewオブジェクトをグリッドを用いて分割し,そこに1分子ずつ描画しています.linked=Falseオプションを使うことで,グリッド毎に分子を回転させることができるようになります.今回は静止画だけしかお示ししていませんが,是非自分の手で表示させてグリグリと回転させてみてみてください.

v = py3Dmol.view(width=600, height=400, linked=False, viewergrid=(2,2))
for m, i in zip(mols[:4], [(0,0),(0,1),(1,0),(1,1)]):
    mb = Chem.MolToMolBlock(m)
    v.addModel(mb, 'sdf', viewer=i)
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick': {}})
v.zoomTo()
v.show()

終わりに

今回はRDKitに実装されている3次元構造生成のアルゴリズムとpy3Dmolを用いた立体構造の描画方法について学びました.今回はどれだけ結晶構造に近い構造を得られるかに着目しましたが,立体構造の生成では分子が取り得る多様なコンフォマーを生成するということも大事になります.

次回は1つの分子に対して複数のコンフォマーを発生させて,その構造分布を見ていきたいと思います.その過程でコンフォマーオブジェクトについて理解するとともに,今回簡単にだけ触れたRDKitにおける分子力学法の実装についても触れたいと思います.

>>次の記事:「RDKitによるコンフォマーの生成

コメント

タイトルとURLをコピーしました