RDKitによるコンフォマーの生成

RDKitによる3次元構造の生成」という記事ではSMILESなどから作成したMolオブジェクトから,どのように3次元構造を作り出すかを学びました.その際にRDKitに実装されているアルゴリズムはランダム(stochastic)な手法であるディスタンス・ジオメトリー法を基にしていることにも触れました.今回は1つの分子から複数の構造を生成することにより,分子が取りうる色々な配座について見ていきたいと思います.
まずは必要なライブラリなどのimportをしましょう.今回は「Platinum dataset」の最初の分子を扱うことにします.なお適切な3次元構造の生成には水素原子の存在が重要となりますので,SMILESなどから作成したMolオブジェクトを使用するさいにはChem.AddHsで水素を付加してから作業を行いましょう.

# ライブラリのimport
from rdkit import rdBase, Chem
print(rdBase.rdkitVersion) # 2017.09.1
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
import py3Dmol
import pandas as pd
# 分子の読み込み
suppl = Chem.SDMolSupplier('./platinum_dataset_2017_01.sdf', removeHs=False)
mols = [x for x in suppl if x is not None]
len(mols) # 4548
mol = mols[0]

コンフォマーの生成

AllChem.EmbedMultipleConfs(mol, numConfs, pruneRmsThresh)

様々な引数によって細かい設定が可能ですが,代表的なものは上にあげた3つになります.すなわち,

  1. どの分子の3次元構造を生成するか(mol)
  2. いくつのコンフォマーを生成するか(numConfs)
  3. どの程度類似した構造を排除するか(pruneRmsThresh)

です.注意点としては生成したコンフォマーの枝刈り(プルーニング)は指定数のコンフォマーを生成し終わった後に実行されるため,例えば100個のコンフォマーを生成したものの閾値が大きかったために実際には10個のコンフォマーしか得られなかったということがありえます.EmbedMultipleConfsは実際に得られたコンファオマーのID数をリストで返しますので,そのリストを調べることでいくつのコンフォマーが得られたかがわかります.下のコードでは1000個のコンフォマーを発生させて,pruneRmsThreshの大きさによってどの程度枝刈りされるかを見てみます.その際randomSeedを設定することによって生成するコンフォマーは同じものが得られるようにしています.

def conf_generation(mol, numConfs, rms):
    sm = Chem.MolToSmiles(mol)
    m = Chem.MolFromSmiles(sm)
    m_h = Chem.AddHs(m)
    num_of_confs = []
    for i in rms:
        confids = AllChem.EmbedMultipleConfs(m_h, numConfs=numConfs,
                                             randomSeed=1234, pruneRmsThresh=i,
                                             numThreads=0)
        num_of_confs.append((i,len(confids)))
    return num_of_confs

rms = [0.5, 1.0, 1.5, 2.0] # 4レベルの閾値設定
num_of_confs = conf_generation(mol, 1000, rms)
threshold conformers
0.5 1000
1.0 1000
1.5 987
2.0 259

今回の分子は回転可能な結合数が多めなので,閾値が1.5程度でもほとんどのコンフォマーが残っていますが,もっと自由度の低い分子ではこの限りではありません.

分子力学法による構造最適化

ディスタンス・ジオメトリー法によって生成させた構造は,芳香環が平面構造でないなどの問題を含むために構造最適化を行う必要があるということを以前学びました.今回はRDKitにおける分子力学法の実装について詳しくみていきたいとおもいます.大きくわけると,1) AllChemに実装されているメソッドによって直接最適化を行う方法と,2) 最初に力場(フォースフィールド)オブジェクトを作成してから構造最適化を行う方法,の2つに分けられます.構造最適化のみが行いたい場合は前者を,エネルギー計算なども行いたい場合は後者を用いるのがよいでしょう.これはUFF・MMFFともに共通です.

メソッドによる直接最適化

AllChem.’XXX‘OptimizeMoleculeConfs(mol, numThreads, maxIters)
‘XXX’ = UFF or MMFF

下のコードでは枝刈りの閾値を2に設定して得たコンフォマー間のRMSを取得しています.十分にばらけた構造が得られていることがわかります.

def mm_opt(a, ff):
    m = Chem.MolFromSmiles(Chem.MolToSmiles(a))
    m_h = Chem.AddHs(m)
    cids = AllChem.EmbedMultipleConfs(m_h, numConfs=100, randomSeed=1234,
                                     pruneRmsThresh=2, numThreads=0)
    if ff == 'uff':
        AllChem.UFFOptimizeMoleculeConfs(m_h, numThreads=0)
    if ff == 'mmff':
        AllChem.MMFFOptimizeMoleculeConfs(m_h, numThreads=0)
    rmsd = []
    for cid in cids:
        rmsd.append(AllChem.GetConformerRMS(m_h, 0, cid))
    return rmsd

uff_rmsd = mm_opt(mol, 'uff')
mmff_rmsd = mm_opt(mol, 'mmff')
df = pd.DataFrame({'uff': uff_rmsd,
                  'mmff': mmff_rmsd})
df.describe()
uff mmff
count 62.000000 62.000000
mean 3.075043 3.021433
std 0.566500 0.511612
min 0.000000 0.000000
25% 2.832945 2.836621
50% 3.106862 3.115341
75% 3.407391 3.289519
max 4.197447 3.709106

フォースフィールドオブジェクトの生成

AllChem.UFFGetMoleculeForceField(mol, confId)
AllChem.MMFFGetMoleculeForceField(mol, Properties, confId)

UFFとMMFFの違いは力場パラメータが引数に必要かどうかになります.パラメータはMMFFGetMoleculePropertiesによって得られます.このようにして作成したフォースフィールドオブジェクトにはCalcEnergyMinimizeという2つのメソッドが実装されていて,それぞれエネルギー計算構造最適化に用います.エネルギーはkcal/mol単位で表示されます.ここでは発生させたコンフォマーを最適化したのちに,エネルギー計算を行い,安定な5つのコンフォマーを取り出してそのエネルギー差を見てみます.

def opt_sp_mm(a, ff):
    m = Chem.MolFromSmiles(Chem.MolToSmiles(a))
    m_h = Chem.AddHs(m)
    cids = AllChem.EmbedMultipleConfs(m_h, numConfs=100, randomSeed=1234,
                                       pruneRmsThresh=2, numThreads=0)
    energy = []
    if ff == 'uff':
        for cid in cids:
            uff = AllChem.UFFGetMoleculeForceField(m_h, confId=cid)
            uff.Minimize()
            energy.append((uff.CalcEnergy(), cid))
    if ff == 'mmff':
        prop = AllChem.MMFFGetMoleculeProperties(m_h)
        for cid in cids:
            mmff = AllChem.MMFFGetMoleculeForceField(m_h, prop, confId=cid)
            mmff.Minimize()
            energy.append((mmff.CalcEnergy(), cid))
    energy.sort()
    return [(i-energy[0][0],j) for i,j in energy]
uff_e = opt_sp_mm(mol, 'uff')
mmff_e = opt_sp_mm(mol, 'mmff')

結果をまとめたものが下の表になりますが,安定なコンフォマーのIDが計算手法によって随分異なるのが気になります.コンフォマー間のエネルギー差も小さいので,本当に異なる構造なのか次のセクションで可視化してみることにしましょう.

力場 エネルギー差(kcal/mol) コンフォマーID
UFF 0.0 26
0.73 25
1.79 16
1.94 17
2.23 48
MMFF 0.0 46
0.10 9
0.19 1
0.89 56
1.67 29

コンフォマーの整列

AllChem.AlignMolConformers(mol, atomIds, RMSlist)

atomIdsをテンプレートとして用いることでどの構造を中心に整列するかを指定できます.RMSlistに空のリストを指定することによってコンフォマー間のRMSリストの一括取得が可能です.または,
AllChem.GetConformerRMS(mol, confId1, confId2)
を用いることで任意のコンフォマー間のRMSを直接計算することも可能です.

コンフォマーの座標を整列させた後に描画することで,どのあたりの構造が異なるのかが目で見て確認することが可能になるかもしれません.今回はベンゼン環部位をテンプレートとして構造を整列させてみます.その際にベンゼン環部位を形成するatomIdが必要になりますが,Mol.GetSubstructMatchを用いることで取得できます.

def get_confs(a, ff):
    m = Chem.MolFromSmiles(Chem.MolToSmiles(a))
    m_h = Chem.AddHs(m)
    cids = AllChem.EmbedMultipleConfs(m_h, numConfs=100, randomSeed=1234,
                                       pruneRmsThresh=2, numThreads=0)
    if ff == 'uff':
        for cid in cids:
            uff = AllChem.UFFGetMoleculeForceField(m_h, confId=cid)
            uff.Minimize()
    if ff == 'mmff':
        prop = AllChem.MMFFGetMoleculeProperties(m_h)
        for cid in cids:
            mmff = AllChem.MMFFGetMoleculeForceField(m_h, prop, confId=cid)
            mmff.Minimize()
    return Chem.RemoveHs(m_h)
# 前項で得た安定構造のID
uff_confIds = [26, 25, 16, 17, 48]
mmff_confIds = [46,9,1,56,29]
# コンフォマーの取得
uff = get_confs(mol, 'uff')
mmff = get_confs(mol, 'mmff')
# テンプレートとの一致部位の取得
core = uff.GetSubstructMatch(Chem.MolFromSmiles('C1=CC=CC=C1'))
# テンプレートを用いてコンフォマーの整列
AllChem.AlignMolConformers(uff, atomIds=core)
AllChem.AlignMolConformers(mmff, atomIds=core)
# 描画
v = py3Dmol.view(width=600, height=600)
for cid in uff_confIds[:3]:
    IPythonConsole.addMolToView(uff_H, confId=cid, view=v)
v.setBackgroundColor('0xeeeeee')
v.zoomTo()
v.show()

まずUFFで安定化した上位3つの構造です.

続いてMMFFで安定化した上位3つの構造です.

最後に最安定構造同士を重ね合わせたものです.

目で見るとすぐにわかりますが,構造が全然違います.結論としては,今回のように回転可能な結合数が多数ある分子を扱う場合には再安定構造を見つけるのは非常に困難です.これはより高レベルの計算手法である分子軌道法やDFTを用いたとしても変わりません.特に今回はわざわざRMSの閾値を大きめに取ったため,なおさらコンフォマー間の構造の違いが大きくなっています.実際には今回のようなフレキシブルな分子は,まずまずエネルギーが低いコンフォマー間の集団(アンサンブル)として振る舞うことから,そのエネルギー差をふまえたボルツマン分布を考慮して分子の性質を考えていく必要があります.

コンフォマーオブジェクトと原子の座標

Mol.GetNumConformers()
Mol.GetConformer(id)
Mol.GetConformers()

以前にRDKitにおける原子であるAtomオブジェクトを扱った際に,プロパティとして原子の座標がなかったことを不思議に思った人がいたかもしれません.今回見てきたようにMolオブジェクトはいくつものコンフォマーを含む箱のようなものであって,原子の座標を確定させるためには原子番号だけではなく,どのコンフォマーに関する座標かを指定するコンフォマーIDが必要になります.
しかし,これまでにもMolBlockへと分子を書き出す際に原子の座標が表示されていました.実はMolToMolBlockには引数にデフォルトではconfId=-1が設定されているため,特にコンフォマーオブジェクトの存在を意識することなく原子の座標を取得できたのです.従って別のコンフォマーの座標を書き出したい場合には,引数に任意のIDを設定することでMolBlockへと書き出すことが可能です.

コンフォマーオブジェクトの取り扱い

Conformer.GetId()
Conformer.GetAtomPosition(atomId)
Conformer.GetPositions()
Conformer.GetOwningMol()

Molオブジェクトに存在するコンフォマーにアクセスするには2つの方法があり,IDを指定する方法全てのコンフォマーを一気に取得する方法があります.いずれの場合にもコンフォマーオブジェクトが得られます.コンフォマーオブジェクトには原子座標を取得するメソッドや,自身の入れ物であるMolオブジェクトを取得するメソッドなどが設定されています.ここでは簡単に最初の10原子について元素記号と(x,y,z)の座標を表示してみます.シンプルですが色々と他のソフトとの連携が思い描けると思います.

conf = mol.GetConformer(0)
for i,(x,y,z) in zip(range(10),conf.GetPositions()[:10]):
    atom = mol.GetAtomWithIdx(i)
    print('{}\tx: {:.2f}\ty: {:.2f}\tz: {:.2f}'.format(atom.GetSymbol(),x,y,z))
C  x: 46.82    y: 28.74    z: 39.61
C   x: 47.36    y: 28.03    z: 38.34
C   x: 47.59    y: 29.05    z: 37.23
C   x: 48.41    y: 28.49    z: 36.08
C   x: 47.16    y: 31.04    z: 40.17
C   x: 46.62    y: 27.70    z: 40.71
C   x: 45.03    y: 26.13    z: 41.65
C   x: 44.11    y: 26.45    z: 42.83
C   x: 49.26    y: 34.15    z: 39.31
C   x: 48.27    y: 33.01    z: 39.17

終わりに

今回の記事ではMolオブジェクトからのコンフォマー生成の方法から始まり,RDKitにおける分子力学法の取り扱い方,コンフォマーオブジェクトから原子の座標データを取り出す方法などを学びました.次回はMolオブジェクトを用いて算出できる記述子について見ていくことで,分子の性質をどのように表現するかを学んでいきます.化学における各記述子は機械学習モデルの特徴量に相当するため,モデリングや分析といった仕事に直結する内容と言えるでしょう.

>>次の記事:「RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ

シェアする

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

フォローする