Contents

RDKitによる低分子のコンホメーション探索¶

ETKTG法を使ったコンホメーション生成、エネルギー最適化¶

研究室内用のメモです。Jupyter Lab形式ファイル（sample.ipynb）をダウンロードして自分で動かしてみてください。

結論を先に書いておくと、今回の方法ではテスト分子のコンホメーションのサンプリングをうまくできませんでした。

今回テストで使用した9員環ラクタム分子（Indolactam-V）ではトランスアミド型がほぼ生成されず（2000個発生させて1個だけ）、それも実験的に分かっている安定なトランスアミド型コンホメーション（Sofa型）ではありませんでした。以前紹介したGROMACSのSimulated Annealingを使った方法だとシスアミド型もトランスアミド型も満遍なくサンプリングできているのでそちらを使った方がまだ安心のようです。

以下のサイトを参考にさせていただきました:

LibraryのImport¶

In [1]:

import os

import numpy as np
print("NumPy version:", np.__version__)

import rdkit # RDKit 2019.09.3
from rdkit import Chem 
from rdkit.Chem import AllChem, Draw, Descriptors, PandasTools, TorsionFingerprints
from rdkit.ML.Cluster import Butina
from rdkit.Chem.Draw import IPythonConsole
#IPythonConsole.ipython_3d=True

print("RDKit version:", rdkit.__version__)

import nglview # NGLView 2.7.1
print("NGLView version:", nglview.__version__)

import pandas as pd
print("Pandas version:", pd.__version__)

import pubchempy as pcp
print("PubChemPy version:", pcp.__version__)

NumPy version: 1.18.1

RDKit WARNING: [14:21:55] Enabling RDKit 2019.09.3 jupyter extensions

RDKit version: 2019.09.3

NGLView version: 2.7.1
Pandas version: 1.0.1
PubChemPy version: 1.0.4

SMILES stringの読み込み¶

Aplysiatoxin. Pubchem CID: 21672114.
Indolactam-V. Pubchem CID: 105000.

In [2]:

pubchemcid = 105000
smiles = pcp.get_properties('IsomericSMILES', pubchemcid)
print(smiles)
mol = Chem.MolFromSmiles(smiles[0]['IsomericSMILES'])
Draw.MolToImage(mol) # RDKitによる構造の表示

[{'CID': 105000, 'IsomericSMILES': 'CC(C)[C@H]1C(=O)N[C@@H](CC2=CNC3=C2C(=CC=C3)N1C)CO'}]

Out[2]:

コンホメーションを発生する関数、コンホマー数を定義する¶

ETKTG法を使用

Riniker, S.; Landrum, G. A. Better Informed Distance Geometry: Using What We Know To Improve Conformation Generation. J. Chem. Inf. Model. 2015, 55, 2562-2574.

In [3]:

ps = AllChem.ETKDGv2() # AllChem.ETKDG() がデフォルト
ps.pruneRmsThresh = 0.2 # 0 だと枝刈りしない
ps.numThreads = 0 #並列実行数。0は最大のスレッド数を使用。
#生成するコンホマー数; 例: n = 100
n = 2000

def generateconformations(m, n):
    m_h = Chem.AddHs(m) #水素原子の付加
    ids = AllChem.EmbedMultipleConfs(m_h, n, ps)
    AllChem.MMFFOptimizeMoleculeConfs(m_h, numThreads=0) #ETKDG法の場合は最適化されている
    prop = AllChem.MMFFGetMoleculeProperties(m_h)
    energies = []
    for cid in ids:
        ff = AllChem.MMFFGetMoleculeForceField(m_h, prop, confId=cid)
        energy = ff.CalcEnergy()
        energies.append((energy, cid))
    return m_h, list(ids), energies

# 入力した分子はm, 生成したコンホマーはconformersとして返す。

コンホメーションをn個発生¶

In [4]:

%%time
m, ids, energies = generateconformations(mol,n)
print("Number of generated conformers:", len(ids))

Number of generated conformers: 209
CPU times: user 3min 16s, sys: 18.7 ms, total: 3min 16s
Wall time: 6.21 s

In [5]:

#コンホマーのアラインメント
AllChem.AlignMolConformers(m)

In [6]:

sort = np.argsort([energy[0] for energy in energies]) #arsortは値ではなく並び替えたインデックスを返す
energy_sorted = np.sort([energy[0] for energy in energies])

# エネルギー順にコンホマーを並べ替えた新しい分子を作る。
m_sorted = Chem.Mol(m)
m_sorted.RemoveAllConformers()
conf_ids = [conf.GetId() for conf in m.GetConformers()]
for i in sort:
    conf = m.GetConformer(conf_ids[i])
    m_sorted.AddConformer(conf, assignId=True)

クラスタリング¶

Schulz-Gasch, T.; Schärfer, C.; Guba, W.; Rarey, M. TFD: Torsion Fingerprints As a New Measure To Compare Small Molecule Conformations. J. Chem. Inf. Model. 2012, 52, 1499-1512

In [7]:

# 水素原子を削除
m_noh = Chem.RemoveHs(m_sorted)
# RMS行列（RMS matrix）
rmsmat = AllChem.GetConformerRMSMatrix(m_noh, prealigned=True)
# TFD行列（TFD matrix）
tfdmat = TorsionFingerprints.GetTFDMatrix(m_noh)
# クラスタリング
num = m_noh.GetNumConformers()
rms_clusters = Butina.ClusterData(rmsmat, num, 0.5, isDistData=True, reordering=True)
tfd_clusters = Butina.ClusterData(tfdmat, num, 0.04, isDistData=True, reordering=True)
print('Number of RMS clusters:', len(rms_clusters))
#for i in range(len(rms_clusters)):
#    print('RMS cluster_{}:'.format(i), rms_clusters[i])
print('\nNumber of TFD clusters:', len(tfd_clusters))

Number of RMS clusters: 47

Number of TFD clusters: 7

SDFファイルに保存¶

In [8]:

writer = Chem.SDWriter("./conformers.sdf")
for i in range(m_sorted.GetNumConformers()):
    writer.write(m_sorted, confId=i)
writer.close()

Pandasを使ってテーブルにする¶

In [9]:

df = PandasTools.LoadSDF('./conformers.sdf', removeHs=False) # 水素原子も読み込む
df['MMFF_Energy'] = [E for E in energy_sorted] # Dataframeに新たな列を追加

コンホメーションの可視化¶

NGLViewを使用して最安定コンホマーを表示

In [10]:

print([r[0] for r in rms_clusters])

[158, 156, 154, 109, 82, 73, 37, 192, 88, 57, 51, 160, 195, 141, 91, 56, 157, 100, 208, 203, 102, 26, 197, 186, 175, 159, 150, 131, 111, 75, 32, 3, 206, 140, 106, 74, 44, 199, 188, 187, 183, 172, 123, 76, 54, 53, 48]

In [11]:

confid = 0
view = nglview.show_rdkit(df.loc[confid, 'ROMol']) # MMFF力場での最安定コンホメーション
view

In [12]:

view.render_image()

QM計算による比較¶

手法： B97-3c: 2018年に発表されたDFT法の一種。

B97: 1997年に発表された汎関数
3c: “three-fold corrected”

基底関数形はdef2-mTZVPが使われる。

Intel CPUはHyper-Threading Technologyを利用して1コア/2スレッドで動作しますが、mpirunは “–use-hwthread-cpus” オプションを付けて実行しなければ、物理的なコア数までしか並列計算をしてくれない。物理コア数以内でも、qsubからジョブを送ると止まる。コマンドラインからだと大丈夫な時もあるが、止まる時は止まる。理由不明。

ORCA finished by error termination in SCF

…. aborting the run

In [13]:

def GenOrcaInput(i, method='B97-3c', thread=4):
    conf = df.at[i, 'ROMol']
    molBlock = Chem.MolToMolBlock(conf)
    # 原子数を取得
    num = conf.GetNumAtoms()
    # 全ての行をリストとして保存する
    lst_molBlock = molBlock.split("\n")
    # 4行目から（4+原子数）行目までの各行を抜き出す
    coordinates_part = lst_molBlock[4 : 4 + num]
    # 原子の種類と各座標をタプルとして保存
    coordinates = [(atoms[31],atoms[0:30]) for atoms in coordinates_part]
    os.makedirs('conf{}'.format(i), exist_ok=True)
    with open('conf{}/conf{}_{}.inp'.format(i, i, method.replace(' ', '_')), 'w') as f:
        print('#', file=f)
        print('# conformer{}'.format(i), file=f)
        print('#', file=f)
        print('! {}'.format(method), file=f)
        print('! OPT FREQ', file=f) # 構造最適化と振動計算
        print('%pal', file=f) # 並列計算
        print('nprocs {}'.format(thread), file=f) # スレッド数
        print('end', file=f)
        print('* xyz 0 1', file=f)
        for atom, xyz in coordinates:
            print('{}\t{}'.format(atom, xyz), file=f)
        print('*', file=f)

In [14]:

num_calc = 4 #計算する上位コンホマー数
method = "B97-3c" # 汎関数 (functional), 基底関数系 (basis set) 
#method="B3LYP def2-TZVP"
#method="wB97X-D3 def2-TZVP"

for i in range(num_calc):
    GenOrcaInput(i, method=method, thread=8)

qsub用job scriptの作成¶

In [15]:

file=open("orca.sh","w")
file.write("""#!/bin/bash
#PBS -l nodes=1:ppn=8

cd $PBS_O_WORKDIR

source /opt/orca/mpivars.sh

cd $PBS_O_WORKDIR/conf0/
/opt/orca/orca_4_2_1_linux_x86-64_openmpi314/orca conf0_B97-3c.inp > conf0_B97-3c.out
""")
file.close()

In [21]:

#out=!{"/opt/orca/orca_4_2_1_linux_x86-64_openmpi314/orca conf0/conf0_B97-3c.inp > conf0/conf0_B97-3c.out"}
#out=!{"qsub orca.sh"}

In [22]:

#grepout = !{"grep 'Final\sGibbs\sfree\senergy' conf0/conf0.out"}
#print(grepout)

In [23]:

# out=!{"/opt/orca/orca_4_2_1_linux_x86-64_openmpi314/orca conf1.inp > conf1.out"}
# 15分弱かかった

In [24]:

#grepout = !{"grep 'Final\sGibbs\sfree\senergy' conf1.out"}
#print(grepout)

RDKitを使った低分子のコンホメーション探索

RDKitによる低分子のコンホメーション探索¶

ETKTG法を使ったコンホメーション生成、エネルギー最適化¶

LibraryのImport¶

SMILES stringの読み込み¶

コンホメーションを発生する関数、コンホマー数を定義する¶

コンホメーションをn個発生¶

クラスタリング¶

SDFファイルに保存¶

Pandasを使ってテーブルにする¶

コンホメーションの可視化¶

QM計算による比較¶

qsub用job scriptの作成¶

コメントを残すコメントをキャンセル

RDKitによる低分子のコンホメーション探索¶

ETKTG法を使ったコンホメーション生成、エネルギー最適化¶

LibraryのImport¶

SMILES stringの読み込み¶

コンホメーションを発生する関数、コンホマー数を定義する¶

コンホメーションをn個発生¶

クラスタリング¶

SDFファイルに保存¶

Pandasを使ってテーブルにする¶

コンホメーションの可視化¶

QM計算による比較¶

qsub用job scriptの作成¶

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル