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

投稿日: 投稿者:
 

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法を使用

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)
 

クラスタリング

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)

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です

CAPTCHA


このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください。