RDKitで３次元構造に挑戦！part 1 〜HALAVEN (Eribulin mesylate)を題材に〜

以前の記事で少し取り上げたHalaven (Eliburin)・・・

エリブリン - Wikipedia

私は構造式をみても全く３次元的な構造がイメージできません。

ところで、ケモインフォマティクス について調べていると度々目にする"RDKit"。

どうやらタダで３Dの描画もできる様です。
さらには日本語で使い方を解説してくださっている素晴らしいブログ（↓）もあります。
future-chem.com

これは入門するしかない！　
ということでRDKit に入門（コピペ）することにしました。。
備忘録を兼ねてコードを貼りますが主に上記ページの引用です。

Pythonも初心者なので怪しい記述がたくさん出てきます。ご容赦ください。

最終的には

　1. Eribulin と 元になった天然物 Halichondrin B の比較

　2. RDKitで出した構造と結晶構造の比較

をしてみたいのですが、とりあえずはHalichondrin B の描画に挑戦しました。
とりあえず、分かったことは

　1. SMILES にもCanonical やら Isomeric やら色々ある

　2. 立体構造の描画にプロトンを無視してはいけない

です。
 

(1) PubChemからの情報取得 (参考記事①)

まず、PubChemから化合物の構造に関する情報を取得します。

Halichondrin B と Eribulin の情報をそれぞれpandasのdataframeとして取得しました。

#PubChemPyをpcpとしてimport
import pubchempy as pcp

#列が多すぎるので全列表示するようにpandasを設定する
import pandas as pd
pd.set_option('display.max_columns', None) 

#Halichondrin B とEribulinを名前検索で情報取得
df_hali = pcp.get_compounds('Halichondrin B', 'name', as_dataframe=True)
df_eri = pcp.get_compounds('Eribulin', 'name', as_dataframe=True)

取得される情報は右に記載されていました。API documentation — PubChemPy 1.0.4 documentation

cid	PubChem Compound Identifier	elements	含まれる元素記号のリスト
atoms	原子のリスト	bonds	結合のリスト
synonyms	別名	sids	PubChem Substance Identifier
aids	PubChem Assay Identifier	charge	形式電荷
molecular_formula	分子式	molecular_weight	分子量
canoical_smiles	Canoniacal SMILES	isomeric_smiles	Isomeric SMILES
inchi	InChI string	inchikey	InChIKey
iupac_name	IUPAC name	xlogp	XLogP
exact_mass	Exact mass	monoisotopic_mass	Monoisotopic mass
tpsa	Topological Polar Surface Area	complexity	Complexity
h_bond_donor_count	水素結合ドナー数	h_bond_acceptor_count	水素結合アクセプター数
rotatable_bond_count	Rotatable bond count	fingerprint	hex-encoded fingerprint
cactvs_fingerprint	PubChem CACTVS fingerprint	heavy_atom_count	Heavy atom count
isotope_atom_count	isotope atom count	atom_stereo_count	Atom stereocenter count
defined_atom_stereo_count	Defined atom stereocenter count	undefined_atom_stereo_count	Undefined atom stereocenter count
bond_stereo_count	Bond stereocenter count	defined_bond_stereo_count	Defined bond stereocenter count
undefined_bond_stereo_count	Undefined bond stereocenter count	covalent_unit_count	Covalently-bonded unit count

"英語のまま = 意味がわからないからコピペした"

Halichondrin B はエントリーが7つでてきますが、どれも同じようなのでWikipediaにものっている”CID(Compound ID) : 5488895" を使用することにします。

#どうせWikipediaのCIDつかうならこうすりゃよかった。 
pcp.get_compounds(5488895, 'cid', as_dataframe=True)

Eribulin は" CID: 11354606" のみでした。
Wikipedia にはメシル酸塩 "CID: 17755248" が記載されていますが、こちらは出てきませんでした。
名前検索は部分一致ではなく、全体一致なのでしょうか？

(2) Canonical SMILES からの２次元描画 ( 参考記事② )

次にHalichondrin B の Canonical SMILESから2次元構造を描画してみます。

#必要なものをimport
from rdkit import Chem
from rdkit.Chem import AllChem, Draw

#取得したDataFrameからCanonical SMILESの情報を取り出す。
hali_smi_cano = df_hali.loc[5488895, "canonical_smiles"]

#smilesからMolオブジェクトを作成
hali_mol_cano = Chem.MolFromSmiles(hali_smi_cano)

#2次元座標の計算
AllChem.Compute2DCoords(hali_mol_cano)

#2次元描画
Chem.Draw.MolToImage(hali_mol_cano)

できました！！

なるほど・・・そうきましたか。
流石にこういう感じ（↓）を自動生成は無理ですよね。
ハリコンドリンB - Wikipedia

上記は水素原子が省かれているのでつけてみます。

#プロトンをつけて描画
hali_mol_canoH = AllChem.AddHs(hali_mol_cano)
AllChem.Compute2DCoords(hali_mol_canoH)
Chem.Draw.MolToImage(hali_mol_canoH)

ぞわぞわする・・・サクサク行こう。

(3) ３次元構造を生成・描画 (参考記事③)

では、本題の３次元構造を作ってみましょう。
まずは見やすそうなので、そのまま"水素原子なし"でやって見ます。

#必要なものをimport
from rdkit.Chem.Draw import IPythonConsole
import py3Dmol

#3次元構造を生成
AllChem.EmbedMolecule(hali_mol_cano, AllChem.ETKDG())

#Jupyter notebook上で描画
IPythonConsole.drawMol3D(hali_mol_cano)

おおおお！格好良い！これだ、これが見たかった！
Jupyter notebookでは回転とズームもできます！（上は画像なので無理です）

"水素原子あり"も作って、二つ並べて見ましょう。

#プロトン付きでも３次元構造を生成
AllChem.EmbedMolecule(hali_mol_canoH, AllChem.ETKDG())

#並べて表示するためにpy3Dmolを使用
#まずはviewオブジェクトを生成
v = py3Dmol.view(width=1000, height=400, linked=False, viewergrid=(1,2))

# "プロトンあり"と"なし"をそれぞれMolBlock形式としてグリッドに渡す。
hali_mol_cano_mb = Chem.MolToMolBlock(hali_mol_cano)
v.addModel(hali_mol_cano_mb, viewer=(0, 0))
hali_mol_canoH_mb = Chem.MolToMolBlock(hali_mol_canoH)
v.addModel(hali_mol_canoH_mb, viewer=(0, 1))

#格好良くする設定
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

" linked=False " と設定しているので左と右それぞれの別々に動かせます。
うひょーっ・・・ん？ なんか変。

なんと、水素あり、水素なしで立体が変わってしまっています。
水素付加の仕方に問題があるのでしょうか？？

いろいろ調べてみると気になる記述が・・・SMILESにも色々種類があるそうです。

・generic SMILES: 原子と結合のみを記述
・isomeric SMILES: 同位体や不斉中心についての記述を含む
・canonical SMILES：generic SMILESをある定義に従って一義的に作成したもの言う
ESI友の会 - 2. SMILESとメタボロミクス

Dr Tomさんも別の記事でSMILES記法について解説してくださっていました。 （参考記事④）

今回、立体構造を生成するにあたって使用したのはcanonical SMILESですが、こちらには不斉中心の情報が含まれていないとのこと。
立体の情報をもたない線形表記法から図を作成したことに問題がありそうです。

ならば、isomeric SMILESでやって見ましょう。

(4) isomeric SMILES からの描画

PubChemから取得した情報にisomeric SMILESもあったので、Halichondrin B の構造をisomeric SMILESからおこしたいと思います。

まずは2次元。

#取得したDataFrameからisomeric SMILESの情報を取り出す。
hali_smi_iso = df_hali.loc[5488895, "isomeric_smiles"]

#smilesからMolオブジェクトを作成し、座標の計算と描画
hali_mol_iso = Chem.MolFromSmiles(hali_smi_iso)
AllChem.Compute2DCoords(hali_mol_iso)
Chem.Draw.MolToImage(hali_mol_iso)

おお！立体を認識してそうです。

ちなみにそれぞれ線形表記は以下のようになります。＠で不斉点を表している模様・・・

canonical SMILES	CC1CC2CCC3C(=C)CC(O3)CCC45CC6C(O4)C7C(O6)C(O5)C8C(O7)CCC(O8)CC(=O)OC9C(C3C(CC4C(O3)CC3(O4)CC4C(O3)C(CC3(O4)CC(C4C(O3)CC(O4)C(CC(CO)O)O)C)C)OC9CC(C1=C)O2)C
isomeric SMILES	C[C@@H]1C[C@@H]2CC[C@H]3C(=C)C[C@@H](O3)CC[C@]45C[C@@H]6[C@H](O4)[C@H]7[C@@H](O6)[C@@H](O5)[C@@H]8[C@@H](O7)CC[C@@H](O8)CC(=O)O[C@@H]9[C@H]([C@H]3[C@H](C[C@@H]4[C@H](O3)C[C@@]3(O4)C[C@H]4[C@@H](O3)[C@H](C[C@]3(O4)C[C@@H]([C@H]4[C@@H](O3)C[C@H](O4)[C@H](C[C@H](CO)O)O)C)C)O[C@H]9C[C@H](C1=C)O2)C

それでは、３次元構造をつくりましょう。
今回も"水素なし"と"水素あり"でやってみます。

#水素原子なしの３次元構造生成
AllChem.EmbedMolecule(hali_mol_iso, AllChem.ETKDG())

#水素原子を付加し３次元構造を生成
hali_mol_isoH = AllChem.AddHs(hali_mol_iso)
AllChem.EmbedMolecule(hali_mol_isoH, AllChem.ETKDG())

#並べて表示
v = py3Dmol.view(width=1000, height=400, linked=False, viewergrid=(1,2))
hali_mol_iso_mb = Chem.MolToMolBlock(hali_mol_iso)
v.addModel(hali_mol_iso_mb, viewer=(0, 0))
hali_mol_isoH_mb = Chem.MolToMolBlock(hali_mol_isoH)
v.addModel(hali_mol_isoH_mb, viewer=(0, 1))
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

水素を付加していない構造（左）は平面になってしまいました。
無駄に回転可能です・・・これはこれで嫌いじゃない！

水素付加したもの（右）は望み通り3次元構造となりました。

以上をまとめると、不斉中心が重要な医薬品に関しては
・立体情報を含むisomeric SMILESを使用
・３次元構造は水素を付加してから生成
と、するのが良さそうです。

次回からは関数っぽくして使おう・・・

import pubchempy as pcp
import pandas as pd
from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole
import py3Dmol

def rittai(PubChemCID): #PubChemのCIDを別途用意して代入 
    #CIDからDataFrameとして化合物の情報を得る
    CID_df = pcp.get_compounds(PubChemCID, 'cid', as_dataframe=True)
    
    #DataFrameからisomeric smilesを抜き出す
    iso_smi = CID_df.loc[PubChemCID, 'isomeric_smiles']
    
    #Molオブジェクトに変換
    iso_mol = Chem.MolFromSmiles(iso_smi)
    
    #水素原子を付加
    iso_molH = AllChem.AddHs(iso_mol)
    
    #3次元構造を生成
    AllChem.EmbedMolecule(iso_molH, AllChem.ETKDG())
    
    #Jupyter notebookで描画
    IPythonConsole.drawMol3D(iso_molH)

それにしても初心者でもコピペで使えるとはDr Tom先生の記事クオリティ高すぎです。

創薬の基本は鍵と鍵穴の関係

*********************************************************************

上記記事の作成にあたって「化学の新しいカタチ」の以下の記事を参考にさせていただきました。
また、本記事の中に誤りがある場合、当ブログ作成者の責任です。
参考記事① : 化合物データベースPubChemをpythonで使いこなす
参考記事② : RDKitの分子Molオブジェクトを扱う
参考記事③: RDKitによる3次元構造の生成
参考記事④: SMILES記法は化学構造の線形表記法

********************************************************************