magattacaのブログ

日付以外誤報

CXSMILESでSMILESに情報を追加する話

マクロ分子の表現HELMについて調べているうちに見慣れないSMILES表現に出会いました。Pistoia Alliance HPのHELMの解説でMonomer (Alanine) の表現として記載されていたものです。

f:id:magattaca:20201123113248p:plain

色々と調べた結果おそらくCXSMILESなのかな？という結論に達しました。

というわけで、CXSMILESについての記事です。ついでにフェロセンをどう表すか問題についても遊んでみました。

CXSMILES

CXSMILES(ChemAxon Extended SMILES)はChemAxon社によるSMILESの拡張表現です（そのまんま）。

簡単にいうとSMILES文字列のあとに分子の特徴を追加で記載する表現形式です。基本的な形式は以下の通り

SMILES_String |<feature1>,<feature2>,...|

SMILESではスペースやタブで区切って書かれた情報をコメントとして無視することになっているそうです。これを利用して、区切りの後にバーティカルバー「|」で挟む形で特徴を記載し、＋αの情報を追加する、というのが基本的な内容です。

記載できる特徴としては原子の座標、アトムラベル、配位結合・水素結合、立体配置といったものがあります。

ChemAxson社のこちらのページ (ChemAxson Extended SMILES and SMARTS - CXSMILES and CXSMARTS )に追加できる特徴やその記載方法が書かれているのでご参照ください。*1

RDKitとCXSMILES

CXSMILESはRDKitでも取り扱うことができます。The RDKit BookのCXSMILES extentions項の記載をChemAxonの説明と比較しつつ表にしてみました。*2

構文解析（parse）できる特徴 記載方法 MolToCXSMILESで書き出せるか？
原子座標
atomic coordinate 		(x,y,z);(x,y,z); できる
原子価
atomic values		 $_AV: できる
原子のラベル/エイリアス $で挟み、各ラベルは;で区切る できる
原子のプロパティ atomprop:0.key1.value1:0.key2.value2:1.key3.value3
アトムインデックス、キー、バリューを１セットにコロンで区切る		 できる
配位結合 C:結合の最初の原子のインデックス.結合のインデックス,
RDKitでは２重結合に翻訳される
ラジカル 1価のラジカル中心は^1:の後、
2価のラジカル中心は^1:の後...		 できる
立体化学表現の拡張
enhanced stereo		 絶対立体配置 a:<atom index>,...
OR o:<atom index>,...
AND &:<atom index>,... 		できる
Link node LN:<atom index>:<minimum repetition>.<maximum repetition>.<first outer atom index>,<second outer atom index>,...
多重中心
multi-center attachments		 m:<multicenter atom index>:<atom index>.<atom index>...
環結合数の指定
rinb bond count specification 		rb
非水素置換数の指定 s
不飽和の指定 u

RDKitで原子のラベル/エイリアスとして認識されるものは_AP、star_e、Q_e、QH_p、AH_P、X_p、XH_p、M_p、MH_p、*です。

enhanced stereoはRDKitではStereoGroupsに変換されます。ABS、AND、ORの説明と利用方法はThe RDKit BookのSuppoer for Enhanced Stereochemistryの項に記載があります。

Link Nodeは実例があまりなく何を表すのかわかりませんでした。ご存知の方教えていただけると嬉しいです。

Pistoia Alliance アラニンモノマーの解釈

CXSMILESの記法だけを並べてもよくわからないので、実例として冒頭のアラニンを見直してみます。

f:id:magattaca:20201123111214p:plain

まず、SMILES文字列中のアスタリスク(*）は「水素原子以外の全ての原子」を表すワイルドカードなのでSMARTSにあたると思います。*3

スペースを挟んで追記されている特徴は以下の通りです。

  • 「r」 相対的立体配置(relatice stereoconfiguration)を表します。
  • 「$;;;_R1;;_R2;$」 原子のラベル/エイリアスを表します。

原子一つずつをセミコロンで区切ってあらわすので、SMILESにおける４番目と６番目の（水素以外の）原子、つまり２つのワイルドカード(*)に対して「_R1、_R２」というラベルを与えています。

このCXSMILES表記はHELMモノマーのアラニンを表すので、「ペプチド結合の相手に何が来ても良いようにSMILESを表現する」という要件を満たすように拡張表記が活用されている、という様子ですね。

Atom label/aliaseについてはRDKitでも取り扱えるとのことですので、早速読み込んでみましょう。

from rdkit import rdBase, Chem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
print('rdkit version: ', rdBase.rdkitVersion)
# rdkit version:  2020.03.2

alanine_monomer = Chem.MolFromSmiles('C[C@H](N[*])C([*])=O |r,$;;;_R1;;_R2;$|')
alanine_monomer

f:id:magattaca:20201123111331p:plain

無事、_R1や_R2といったラベルを含めたMolオブジェクトを生成することができました。原子の情報を確認してみます。*4

# 原子の取得
atoms = [a for a in alanine_monomer.GetAtoms()]

print('atom index 3: Symbol: ', atoms[3].GetSymbol())
print('atom index 3: Properties: ', atoms[3].GetPropsAsDict())

# atom index 3: Symbol:  *
#  atom index 3: Properties:  {'dummyLabel': '*', 'atomLabel': '_R1', '__computedProps': <rdkit.rdBase._vectNSt3__112basic_stringIcNS_11char_traitsIcEENS_9allocatorIcEEEE object at 0x11ea3ae10>, '_CIPRank': 1}

ワイルドカードにはdummy labelというKeyが、_R1にはatomLabelというKeyが割り当てられているようです。

セミコンマによる区切りがわかりづらいので数字をそのまま割り当ててみます。

alanine_fig = Chem.MolFromSmiles('C[C@H](N[*])C([*])=O |$0;1;2;3;4;5;6$|')
alanine_fig

f:id:magattaca:20201123111919p:plain

なるほど。

出力と読み込み

CXSMILESの動作が分かってきたので、もう一例CXSMILESの読み書きを行ってみます。

例としてRDKit Blog「 New drawing options in the 2020.03 release」で取り上げられているジクロフェナクをお借りします。

from rdkit.Chem.Draw import rdMolDraw2D
from rdkit.Chem import rdDepictor
rdDepictor.SetPreferCoordGen(True)
from rdkit.Chem.Draw import IPythonConsole
from IPython.display import SVG

diclofenac = Chem.MolFromSmiles('O=C(O)Cc1ccccc1Nc1c(Cl)cccc1Cl')

cp = Chem.Mol(diclofenac)
cp.GetAtomWithIdx(2).SetProp("atomNote","pKa=4.0")
d2d = rdMolDraw2D.MolDraw2DSVG(350,300)
d2d.drawOptions().setHighlightColour((0.8,0.8,0.8))
d2d.DrawMolecule(cp,highlightAtoms=[2])
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20201123112003p:plain

カルボキシル基の酸素原子にpKaのプロパティをセットしています。

これをCXSmilesに変換してみます。

cp_cxsmi = Chem.MolToCXSmiles(cp)
print(cp_cxsmi)
# 'O=C(O)Cc1ccccc1Nc1c(Cl)cccc1Cl |atomProp:2.atomNote.pKa=4.0|'

SMILES文字列に続いてatomPropが加わっています。Atom Index 2の原子に対して、atomNoteというKeyにpKa=4.0というValueがセットされています。

このCXSMILESから再度Molオブジェクトを構築してみます。

cp_reconst = Chem.MolFromSmiles(cp_cxsmi)
cp_reconst

f:id:magattaca:20201123112058p:plain

再構築したMolオブジェクトでも原子のプロパティがうまく認識されているようです。

線形表記にプロパティを追加でき、再構築しても情報が再現できるというのは多数の分子を取り扱う際にデータ量の節約になりそうです。

フェロセン問題

つづいて話題をかえて「フェロセンどうやって表記するの？」問題に取り組みます。

唐突感がありますが、フェロセンはどうもSMILESでは表しづらい分子だそうです。

フェロセンのSMILES表記

SlideShareに公開してくださっていた資料　「第11回分子化学討論会　PubChemQCプロジェクト: 分子データベース構築と機械学習による電子構造の推定」のスライド p13には「SMILESで表現できない例」というタイトルでFerroceneが例示されています。

スライドには表現例として以下の２つが記載されています。２つめはPubChem（Ferrocene） の表記と同じです。

  1. C12C3C4C5[Fe]23451234C5C1C2C3C45
  2. [CH-]1C=CC=C1.[CH-]1C=CC=C1.[Fe+2]

それぞれRDKitでMolオブジェクトに起こしてみましょう。まずは１つ目。。。

Fer_1 = Chem.MolFromSmiles('C12C3C4C5C1[Fe]23451234C5C1C2C3C45')
Fer_1

f:id:magattaca:20201123112124p:plain!

・・・・スターウォーズにこういう飛行機あった気がする。

for b in Fer_1.GetBonds():
    print('{}: {}'.format(b.GetIdx(), b.GetBondType()))

# 0: SINGLE
# 1: SINGLE
# 2: SINGLE
# 3: SINGLE
# 4: SINGLE
# 5: SINGLE
# 6: SINGLE
# 7: SINGLE
# 8: SINGLE
# 9: SINGLE
# 10: SINGLE
# 11: SINGLE
# 12: SINGLE
# 13: SINGLE
# 14: SINGLE
# 15: SINGLE
# 16: SINGLE
# 17: SINGLE
# 18: SINGLE
# 19: SINGLE

こちらの表現ではFe(II)イオンとシクロペンタジエニルアニオンの結合が明示されている一方で、20本全ての結合が単結合(SINGLE）となっており、芳香属性を表すことができていません。

つづいて２つ目の表現です。

Fer_2 = Chem.MolFromSmiles('[CH-]1C=CC=C1.[CH-]1C=CC=C1.[Fe+2]')

d2d = rdMolDraw2D.MolDraw2DSVG(200,200)
d2d.drawOptions().addBondIndices=True
d2d.DrawMolecule(Fer_2)
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20201123112309p:plain

for b in Fer_2.GetBonds():
    print('{}: {}'.format(b.GetIdx(), b.GetBondType()))

# 0: AROMATIC
# 1: AROMATIC
# 2: AROMATIC
# 3: AROMATIC
# 4: AROMATIC
# 5: AROMATIC
# 6: AROMATIC
# 7: AROMATIC
# 8: AROMATIC
# 9: AROMATIC

分かりやすさのためい結合のインデックスも表記しています。描画ではケクレ構造となっていますが、2つめの表現ではシクロペンタジエニルアニオンの結合がすべて芳香属性として認識されています。

一方で、こちらではFe(II)イオンとの間の結合が表現できていません。

SMILESの拡張表現を使えばうまくフェロセンを表現できるのか？検証してみましょう。

CXSMILESの配位結合表記

ChemAxonのCXSMILES解説には、「CXSMILESの一部として現れるSMILESが、必ずしも標準的なSMILESに一致しない」という事例としてフェロセンが挙げられています。

  • plain SMILES: [Fe].c1cccc1.c1cccc1
  • CXSMILES: c12c3c4c5c1[Fe]23451234c5c1c2c3c45 |C:4.5,0.6,1.7,2.8,3.9,7.12,6.10,9.16,10.18,8.14|

但しこれらはどちらもRDKitではエラーが出てしまい扱うことができません。RDKitではシクロペンタジエニルアニオン部位を芳香属性のc1cccc1であらわすとケクレ化ができないとなってしまいます。sanitize=Falseの設定をすることでオブジェクを作ることはできますがDrawでエラーがいっぱい出てきます。

cpa = Chem.MolFromSmiles('c1cccc1', sanitize=False)
cpa
# エラーメッセージは省略

f:id:magattaca:20201123112426p:plain

話を簡単にするために、ここではCXSMILESの配位結合の表現の利用箇所だけを検証したいと思います。

配位結合の表現は「C:結合の最初の原子のインデックス:結合のインデックス」ですので、Feに結合するBondのインデックスと結合に関わる原子のインデックスを取得します。原子のインデックスは配位結合の出発点、FeではなくCとなるようにします。

Fe_connected_bonds = []
C_idx_list = []

for b in Fer_1.GetBonds():
    bond_idx = b.GetIdx()
    
    begin_atom = b.GetBeginAtom()
    begin_symbol = begin_atom.GetSymbol()
    begin_idx = begin_atom.GetIdx()
    
    end_atom = b.GetEndAtom()
    end_symbol = end_atom.GetSymbol()
    end_idx = end_atom.GetIdx()
    
    if begin_symbol == 'Fe':
        Fe_connected_bonds.append(bond_idx)
        C_idx_list.append(end_idx)
        
    elif end_symbol == 'Fe':
        Fe_connected_bonds.append(bond_idx)
        C_idx_list.append(begin_idx)        

print(Fe_connected_bonds)
#    [4, 5, 11, 12, 13, 14, 15, 16, 17, 18]

print(C_idx_list)
#   [4, 6, 7, 0, 8, 1, 9, 2, 10, 3]

追記すべき配位結合の表現は以下となります。

dative_bonds = 'C:'
n = len(Fe_connected_bonds)

for i in range(n-1):
    dative_bonds = dative_bonds + str(C_idx_list[i]) + '.' + str(Fe_connected_bonds[i]) + ','

dative_bonds = dative_bonds + str(C_idx_list[n-1]) + '.' + str(Fe_connected_bonds[n-1])

print(dative_bonds)
#   C:4.4,6.5,7.11,0.12,8.13,1.14,9.15,2.16,10.17,3.18

従ってCXSMILESはこんな感じになりそうです。

Fer_1_CXSMILES = 'C12C3C4C5C1[Fe]23451234C5C1C2C3C45 |C:4.4,6.5,7.11,0.12,8.13,1.14,9.15,2.16,10.17,3.18|'
Fer_1_CX = Chem.MolFromSmiles(Fer_1_CXSMILES)
Fer_1_CX

f:id:magattaca:20201123112600p:plain

できました！　あいかわらずシクロペンタジエニルアニオンは単結合のままですが配位結合はそれとして認識できているようです！

念のため確認...

for b in Fer_1_CX.GetBonds():
    print('{}: {}'.format(b.GetIdx(), b.GetBondType()))

# 0: SINGLE
# 1: SINGLE
# 2: SINGLE
# 3: SINGLE
# 4: DATIVE
# 5: DATIVE
# 6: SINGLE
# 7: SINGLE
# 8: SINGLE
# 9: SINGLE
# 10: SINGLE
# 11: DATIVE
# 12: DATIVE
# 13: DATIVE
# 14: DATIVE
# 15: DATIVE
# 16: DATIVE
# 17: DATIVE
# 18: DATIVE
# 19: SINGLE

CXSMILESから起こした構造は配位結合をDative bondとして認識することができているようです。

RDKitで配位結合を形成する

次に、RDKitで配位結合を形成し、それをCXSMILESに変換するという方法を試してみます。

シクロペンタジエニルアニオンとFeの間に結合の無かったオブジェクトに対して配位結合を追加していきます。 まずは原子のインデックスを確認します。

d2d = rdMolDraw2D.MolDraw2DSVG(200,200)
d2d.drawOptions().addAtomIndices=True
d2d.DrawMolecule(Fer_2)
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20201123112712p:plain

Fe (Atom Index 10)に対して、シクロペンタジエニルアニオン (Atom Index: 0 ~ 9)から配位結合を形成すれば良いようです。

結合を編集するために、RWMolにオブジェクトを変換します*5

Fer_2_RW = Chem.RWMol(Fer_2)
type(Fer_2_RW)
#  rdkit.Chem.rdchem.RWMol

for i in [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]:
    Fer_2_RW.AddBond(i,10, Chem.BondType.DATIVE)

Fer_2_RW

f:id:magattaca:20201123112754p:plain

for b in Fer_2_RW.GetBonds():
    print(b.GetIdx(), b.GetBeginAtomIdx(),  b.GetEndAtomIdx(), b.GetBondType())
# 0 0 1 AROMATIC
# 1 1 2 AROMATIC
# 2 2 3 AROMATIC
# 3 3 4 AROMATIC
# 4 5 6 AROMATIC
# 5 6 7 AROMATIC
# 6 7 8 AROMATIC
# 7 8 9 AROMATIC
# 8 4 0 AROMATIC
# 9 9 5 AROMATIC
# 10 0 10 DATIVE
# 11 1 10 DATIVE
# 12 2 10 DATIVE
# 13 3 10 DATIVE
# 14 4 10 DATIVE
# 15 5 10 DATIVE
# 16 6 10 DATIVE
# 17 7 10 DATIVE
# 18 8 10 DATIVE
# 19 9 10 DATIVE

シクロペンタジエニルアニオンは全てAROMATIC、Feとの間はDATIVE結合となったオブジェクトを作成することができました！！

これをSMILESおよびCXSMILESに変換してみましょう。まずはSMILES。

Fer_2_smi = Chem.MolToSmiles(Fer_2_RW)
print(Fer_2_smi)
# c12->[Fe+2]3456789(<-c1c->3[cH-]->4c->52)<-c1c->6c->7[cH-]->8c->91

The RDKit Book SMARTS Support and Extensionsの項の記載によると、RDKitではSMILESの拡張として配位結合を<-、->で表せるようになっているそうです。

この配位結合の表現と芳香属性炭素(c)の表現が組み合わさったSMILESになっているようです。

つづいてCXSMILESです。

Fer_2_CXSMILES = Chem.MolToCXSmiles(Fer_2_RW)
print(Fer_2_CXSMILES)
#  c12->[Fe+2]3456789(<-c1c->3[cH-]->4c->52)<-c1c->6c->7[cH-]->8c->91 |^2:1|

先にみたように、配位結合はMolToCXSmilesで書き出せるプロパティの中に入っていません。
ですのでSMILES文字列のあとの特徴には配位結合の情報が含まれていないようです。

追記されている特徴^2:1はラジカルを表す表現のようです。ここではAtom Index 1のFeが２価のラジカルとして認識されています。確認のためCXSMILESからMolオブジェクトを作成し、Atom Indexと共に描画してみましょう。

mol = Chem.MolFromSmiles(Fer_2_CXSMILES)

d2d = rdMolDraw2D.MolDraw2DSVG(200,200)
d2d.drawOptions().addAtomIndices=True
d2d.DrawMolecule(mol)
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20201123112948p:plain

Atom Index 1のFe2+に対してシクロペンタジエニルアニオンから配位結合がでていることが無事確認できました。

CookBookのやり方

と、ここまで書いてきてRDKit CookBookにOrganometallics with Dative Bondsという項目があることに気づきました。そのまんまやん！

CookBookで書かれている方法を参考にGetNeighbors()をつかって配位結合を設定してみたいと思います。

 Fer_3 = Chem.MolFromSmiles('c12c3c4c5c1[Fe]23451234c5c1c2c3c45', sanitize=False)

# RWMolに変換
rw_mol = Chem.RWMol(Fer_3)

# Fe(原子番号26)のAtomを取得
for atm in rw_mol.GetAtoms():
    if atm.GetAtomicNum() == 26:
        metal = atm

# Feの近接原子を取得
nbrs = metal.GetNeighbors()

# Feとの結合をDativeBondに書き換える
for nbr in nbrs:
    rw_mol.RemoveBond(nbr.GetIdx(), metal.GetIdx())
    rw_mol.AddBond(nbr.GetIdx(), metal.GetIdx(), Chem.BondType.DATIVE)

rw_mol
# エラーは省略

f:id:magattaca:20201123113031p:plain

print(Chem.MolToSmiles(rw_mol))
# c12->[Fe]3456789(<-c1c->3c->4c->52)<-c1c->6c->7c->8c->91

先の２例目と異なりシクロペンタジエニルアニオンのアニオンを設定していないので、Kekule化できないエラーはそのままですが配位結合は設定できているようです。

気づかぬ間にRDKit CookBook内容がすごい更新されてました。はじめからこっち見ておけばよかった。。。

まとめ

以上、今回はSMILESの拡張表現の一つ、CXSMILESとその応用としてFerroceneをどのようにSMIELSで表せるか、という問題にアプローチしてみました。

特徴の追記ができることで、線形表現のままより複雑な情報が表せるというのはとても便利そうです。

一方で、副次的ではありますがChemAxon社のページで記載されていた小文字のｃを用いたフェロセンの表記（[Fe].c1cccc1.c1cccc1）がRDKitではエラーが出てしまうこと、ソフトウェアによるSMILES表記の取り扱いの違いがわかった点で興味深いと思いました。

適当なことをあれこれ書いているので間違い等ご指摘いただければ幸いです。

*1:以下の記載はChemAxson社ページの記載を引用しています

*2:RDKit Bookの内容をCC BY 4.0のライセンスのもとに利用させて頂いています。

*3:ChemAxsonのドキュメントではSMILESに記載されていますがおそらくSMARTの方が正しいと思います。「化学の新しいカタチ」さんの記事 SMILES記法は化学構造の線形表記法がとても参考になります。

*4:参考：化学の新しいカタチ「RDKitの分子Molオブジェクトを扱う」

*5:RDKitドキュメント Getting Started In Python