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

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

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

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

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

import pandas as pd
from rdkit import rdBase, Chem
print(rdBase.rdkitVersion)
# 2017.03.3
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole
import py3Dmol

3次元構造の生成

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

方法 コメント
ディスタンス・ジオメトリー(DG)法 2017.09までのデフォルト方法
AllChem.KDG useBasicKnowledge=True
AllChem.ETDG useExpTorsionAnglePrefs=True
AllChem.ETKDG 2018.03以降のデフォルト方法

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

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

そこでDG法で得た構造に対してベンゼン環は平面,アルキンは直線といった一般的な化学常識(basic Knowledge)を加えたのがKDGになります.またケンブリッジ結晶構造データベースの実験データから得た二面角の分布傾向(Experimental Torsion angle preferences)を考慮したものがETDGになります.そしてその両方を加えたのがETKDGです.

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

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

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

なお分子力場としてはUFFMMFFの2つがRDKitでは提供されています.

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

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

  1. DG法の後にUFFで最適化した構造
  2. MMFF94で最適化した構造
  3. ETKDGで発生させた構造

の3種類の方法で得た構造を,もとのX線構造と比較してRMSDを調べてみます.

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

# 分子の読み込み
suppl = Chem.SDMolSupplier('./SDF/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(x) for x in mols[:10]]
mol_from_sm = [Chem.MolFromSmiles(sm) for sm in smiles]

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

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

def ETKDG(mols):
    ETKDG = []
    for m in mols:
        mh = AllChem.AddHs(m)
        AllChem.EmbedMolecule(mh, AllChem.ETKDG())
        ETKDG.append(mh)
    return ETKDG

DG_m = DG(mol_from_sm)
ETDG_m = ETDG(mol_from_sm)
ETKDG_m = ETKDG(mol_from_sm)

# 結晶構造と発生させた構造のRMDSを取得
RMSD = []
for i,j,k,n in zip(DG_m, ETDG_m, ETKDG_m, range(len(smiles))):
    a = AllChem.GetBestRMS(mols[n], i)
    b = AllChem.GetBestRMS(mols[n], j)
    c = AllChem.GetBestRMS(mols[n], k)
    RMSD.append((a,b,c))
df_RMSD = pd.DataFrame(RMSD, columns=['DG','ETDG','ETKDG'])
df_RMSD.describe()
DG + UFF ETDG ETKDG
count 4548.000000 4548.000000 4548.000000
mean 2.101152 2.100846 2.044582
std 0.925544 0.857655 0.910059
min 0.048141 0.219115 0.074690
25% 1.461596 1.477833 1.420255
50% 2.075114 2.042311 2.023980
75% 2.725340 2.657036 2.626488
max 6.054957 5.736471 6.214207


上記の結果からRMDSの平均値meanはETKDGが最も小さいことがわかります.四分位点も総じて小さいことから,一番精度のよい方法と述べてよさそうです.

因みにIPython上で%timeを用いて実行時間(CPU times)を計測した結果が以下の表になりますが,DG+UFFとETKDGの実行速度は同程度と考えてよそうです.実行にはMacBookPro 13′ 2017を用いていました.

method user sys total
DG + UFF 3min 13s 1.59 s 3min 15s
ETDG 2min 41s 1.23 s 2min 43s
ETKDG 3min 9s 7.48 s 3min 17s

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

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

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

pip install py3Dmol
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)
for i in mols[:4]:
    mb = Chem.MolToMolBlock(i)
    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によるコンフォマーの生成

シェアする

  • このエントリーをはてなブックマークに追加

フォローする