RDKitで3次元構造に挑戦! part 2 〜HALAVEN と Halichondrin Bの構造類似性比較〜
前回の記事でHALAVEN®️(eribulin mesylate)の元になった天然物 Halichondrin B の3次元構造を見ることができました。
今回は両者の比較を行って見たいと思います。
- 1. 3次元構造の比較
- 2. 2次元構造で共通部分構造を探索(MCS)
- 3. 2次元構造で類似性を比較(FingerPrint)
- 4. フィンガープリントへの各原子の寄与
- 5. フィンガープリントのビット配列を可視化
- まとめ
1. 3次元構造の比較
まずは前回の復習を兼ねてEribulinの3次元構造を生成します。
PubChemPyを用いた検索からEribulinは "CID: 11354606"ということが分かっていますので、こちらを利用したいと思います。
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 #Eribulin の3次元構造の生成 (参考記事①、②) #PubChem からEribulin (CID:11354606) のisomeric smilesを取得 df_eri = pcp.get_properties('IsomericSMILES', 11354606, 'cid', as_dataframe=True) eri_iso = df_eri.loc[11354606, 'IsomericSMILES'] eri_mol = Chem.MolFromSmiles(eri_iso) eri_molH = AllChem.AddHs(eri_mol) AllChem.EmbedMolecule(eri_molH, AllChem.ETKDG()) #Halicondrin B (CID : 548895) についても同様の処理(省略) #並べて表示 v = py3Dmol.view(width=1000, height=400, linked=False, viewergrid=(1,2)) eri_molH_mb = Chem.MolToMolBlock(eri_molH) v.addModel(eri_molH_mb, viewer=(0, 0)) hali_molH_mb = Chem.MolToMolBlock(hali_molH) v.addModel(hali_molH_mb, viewer=(0, 1)) v.setBackgroundColor('0xeeeeee') v.setStyle({'stick':{}}) v.zoomTo() v.show()
Fig. 1 よりEribulin はHalichondrin B の大員環状ケトン部位を基にしていることがわかります。
ただし完全に部分構造として抜き出されたものではなく、アナログ(構造類縁体)となっています。
ではどの部分構造がが両者で一致しているのでしょうか?
2. 2次元構造で共通部分構造を探索(MCS)
またしても「化学の新しいカタチ」さんの記事を参考にさせていただき、部分構造検索を実施したいと思います。(参考記事③)
使用するのは最大共通部分構造(MCS: Maximum Common Substructure)の探索[FindMCS]です。
では、早速・・・
#Moduleのimport from rdkit.Chem import rdFMCS #時間がかかるので水素原子を付加していない構造を使用し、探索終了の時間(timeout)も設定した mcs_01 = rdFMCS.FindMCS([eri_mol, hali_mol], timeout=60) #一致した部分構造の原子数と結合の数を出力 print('Num of Atoms:{}, Num of Bonds:{}'.format(mcs_01.numAtoms, mcs_01.numBonds))
探索時間を短くするため水素原子を付加しない状態で計算しました。(・・・我慢できない大人なので)
一致した原子の数と結合の数は下記のように出力されました。
「Num of Atoms:48, Num of Bonds:54」
構造としては下記のようになります。
#MCS探索で取得した部分構造のSmartsからMolオブジェクトを作成 mol_mcs_01 = Chem.MolFromSmarts(mcs_01.smartsString) #Eribulin、Halichondrin B それぞれから共通部分構造を取得 eri_match = eri_mol.GetSubstructMatch(mol_mcs_01) hali_match = hali_mol.GetSubstructMatch(mol_mcs_01) #それぞれで一致した原子番号の一覧を表示 print('Eribulin Matched Atoms:{}\n Halichondrin B Matched Atoms:{}'.format(eri_match, hali_match)) #マッチした部分構造をハイライトして描画 Draw.MolsToGridImage([eri_mol, hali_mol], molsPerRow=2, subImgSize=(400,400), highlightAtomLists=[eri_match, hali_match])
Fig 2. 構造式の赤くハイライトされた部分がMCSで共通部分構造として認識された構造となります。
大員環ケトン部位が確かに認識されていることがわかります。
また、大員環から伸びた部分についても、5員環と6員環の違いなどには関わらず、
原子としては一致として認識されていることがわかります。
もう少し格好良く描画したい・・・ということでSVG形式の描画というのに挑戦したいと思います。(参考記事④)
「svgファイル(Scalable Vector Graphics)はXMLをベースとした2次元ベクトル型の画像ファイルフォーマット」で「拡大・縮小しても画質が損なわれることが」ない、とのことです。(上記参考記事④より)
#コンテナとなるオブジェクトの作成 view = rdMolDraw2D.MolDraw2DSVG(300,300) #描画する分子の前処理 processed_eri_mol = rdMolDraw2D.PrepareMolForDrawing(eri_mol) #コンテナの描画設定 option = view.drawOptions() option.circleAtoms=False option.continuousHighlights=False #MCSで一致した部分構造をハイライトに設定 view.DrawMolecule(processed_eri_mol, highlightAtoms=eri_match) #コンテナをファイナライズ view.FinishDrawing() #コンテナに書き込んだデータを取り出す svg = view.GetDrawingText() #データを描画 SVG(svg.replace('svg:,',''))
・・・あまり先ほどと変わらない?
どうやらFig 2. では MolsToGridImage を用いましたが、こちらはChem.rdMolDraw2Dモジュールを使っているらしいので裏側で動いている仕組みが同じであるようです。
比較のためMolToImageで描画すると以下のようになります。
#比較用にMolToImageとMolsToGridImageでそれぞれ描画 Draw.MolToImage(eri_mol, highlightAtoms=eri_match) Draw.MolsToGridImage([eri_mol], molsPerRow=1, subImgSize=(400,400), highlightAtomLists=[eri_match])
確かにMolToImageよりもMolsToGrideImageの方が綺麗です。
ついでなのでもう少し、凝った描画に挑戦したいと思います。
#matplotlibによるRGBカラーコードの生成と色の指定 from matplotlib.colors import ColorConverter color_1 = {} color_2 = {} radius = {} for i, j in zip(eri_match, hali_match): color_1[i] = ColorConverter().to_rgb('skyblue') radius[i] = 0.5 color_2[j] = ColorConverter().to_rgb('khaki') radius[j] = 0.5 processed_eri_mol = rdMolDraw2D.PrepareMolForDrawing(eri_mol) processed_hali_mol = rdMolDraw2D.PrepareMolForDrawing(hali_mol) view2 = rdMolDraw2D.MolDraw2DSVG(1000, 300, 500, 300) #フォントサイズの設定 view2.SetFontSize(2*view2.FontSize()) option = view2.drawOptions() #多重結合の感覚を指定 option.multipleBondOffset = 0.2 #描画領域のマージンを設定 option.padding=0.13 #レジェンドのフォントサイズを設定 option.legendFontSize=18 view2.DrawMolecules([processed_eri_mol, processed_hali_mol], highlightAtoms=[eri_match, hali_match], highlightAtomColors=[color_1, color_2], highlightAtomRadii=[radius, radius], highlightBonds=[[],[]], legends=['Eribulin', 'Halichondrin B']) view2.FinishDrawing() svg = view2.GetDrawingText() SVG(svg.replace('svg:', ''))
原子ラベルのサイズなどを変更することで、かなり分かりやすい図になったように思います。
3. 2次元構造で類似性を比較(FingerPrint)
共通部分構造にもとづく比較を行ったので、次は新たな手法、
フィンガープリントを用いた類似性の検証に挑戦したいと思います。(参考記事⑤)
違いがわからないのでとりあえず全部計算! ・・・アホ丸出し
#比較する二つをまとめる mols = [eri_mol, hali_mol] #①MACCS Keys from rdkit import DataStructs maccs_fps = [AllChem.GetMACCSKeysFingerprint(mol) for mol in mols] maccs_Tanimoto = DataStructs.TanimotoSimilarity(maccs_fps[0], maccs_fps[1]) maccs_Dice = DataStructs.DiceSimilarity(maccs_fps[0], maccs_fps[1]) print('MACCS Keys\nTanimoto Similarity:{}, Dice Similarity:{}'.format(maccs_Tanimoto, maccs_Dice)) #②Topological フィンガープリント(RDKitフィンガープリント) from rdkit.Chem.Fingerprints import FingerprintMols rdkit_fps = [Chem.Fingerprints.FingerprintMols.FingerprintMol(mol) for mol in mols] rdkit_Tanimoto = DataStructs.TanimotoSimilarity(rdkit_fps[0], rdkit_fps[1]) rdkit_Dice = DataStructs.DiceSimilarity(rdkit_fps[0], rdkit_fps[1]) print('Topological Fingerprint\nTanimoto Similarity:{}, Dice Similarity:{}'.format(rdkit_Tanimoto, rdkit_Dice)) #③Morganフィンガープリント-1 ECFP類似 morgan_ECFP_fps = [AllChem.GetMorganFingerprint(mol, 2) for mol in mols] morgan_ECFP_Tanimoto = DataStructs.TanimotoSimilarity(morgan_ECFP_fps[0], morgan_ECFP_fps[1]) morgan_ECFP_Dice = DataStructs.DiceSimilarity(morgan_ECFP_fps[0], morgan_ECFP_fps[1]) print('Morgan Fingerprint(ECFP-like)\nTanimoto Similarity:{}, Dice Similarity:{}'.format(morgan_ECFP_Tanimoto, morgan_ECFP_Dice)) #③'Morganフィンガープリント-1 ECFP類似 (bit vector) morgan_ECFP_bit_fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048) for mol in mols] morgan_ECFP_bit_Tanimoto = DataStructs.TanimotoSimilarity(morgan_ECFP_bit_fps[0], morgan_ECFP_bit_fps[1]) morgan_ECFP_bit_Dice = DataStructs.DiceSimilarity(morgan_ECFP_bit_fps[0], morgan_ECFP_bit_fps[1]) print('Morgan Fingerprint(ECFP-like/bit vector)\nTanimoto Similarity:{}, Dice Similarity:{}'.format(morgan_ECFP_bit_Tanimoto, morgan_ECFP_bit_Dice)) #④Morganフィンガープリント-2 FCFP類似 morgan_FCFP_fps = [AllChem.GetMorganFingerprint(mol, 2, useFeatures=True) for mol in mols] morgan_FCFP_Tanimoto = DataStructs.TanimotoSimilarity(morgan_FCFP_fps[0], morgan_FCFP_fps[1]) morgan_FCFP_Dice = DataStructs.DiceSimilarity(morgan_FCFP_fps[0], morgan_FCFP_fps[1]) print('Morgan Fingerprint(FCFP-like)\nTanimoto Similarity:{}, Dice Similarity:{}'.format(morgan_FCFP_Tanimoto, morgan_FCFP_Dice)) #⑤Avalon フィンガープリント from rdkit.Avalon import pyAvalonTools avalon_fps = [pyAvalonTools.GetAvalonFP(mol) for mol in mols] avalon_Tanimoto = DataStructs.TanimotoSimilarity(avalon_fps[0], avalon_fps[1]) avalon_Dice = DataStructs.DiceSimilarity(avalon_fps[0], avalon_fps[1]) print('Avalon Fingerprint\nTanimoto Similarity:{}, Dice Similarity:{}'.format(avalon_Tanimoto, avalon_Dice)) #⑥AtomPair フィンガープリント #公式のDocumentationによるとatom-pairにはDice similarityを使うのが通常らしい from rdkit.Chem.AtomPairs import Pairs pair_fps = [Pairs.GetAtomPairFingerprint(mol) for mol in mols] pair_Tanimoto = DataStructs.TanimotoSimilarity(pair_fps[0], pair_fps[1]) pair_Dice = DataStructs.DiceSimilarity(pair_fps[0], pair_fps[1]) print('AtomPair Fingerprint\nTanimoto Similarity:{}, Dice Similarity:{}'.format(pair_Tanimoto, pair_Dice)) #⑥' AtomPair フィンガープリント(bit vector fingerprint) pair_bit_fps =[Pairs.GetAtomPairFingerprintAsBitVect(mol) for mol in mols] pair_bit_Tanimoto = DataStructs.TanimotoSimilarity(pair_bit_fps[0], pair_bit_fps[1]) pair_bit_Dice = DataStructs.DiceSimilarity(pair_bit_fps[0], pair_bit_fps[1]) print('AtomPair Fingerprint (bit vector)\nTanimoto Similarity:{}, Dice Similarity:{}'.format(pair_bit_Tanimoto, pair_bit_Dice)) #⑦Torsion フィンガープリント #bitが多すぎるのでこっちはAsBitVectは存在しないらしい from rdkit.Chem.AtomPairs import Torsions tts_fps = [Torsions.GetTopologicalTorsionFingerprintAsIntVect(mol) for mol in mols] tts_Tanimoto = DataStructs.TanimotoSimilarity(tts_fps[0], tts_fps[1]) tts_Dice = DataStructs.DiceSimilarity(tts_fps[0], tts_fps[1]) print('Torsion Fingerprint\nTanimoto Similarity:{}, Dice Similarity:{}'.format(tts_Tanimoto, tts_Dice)) #結果をPandasのデータフレームにまとめる import pandas as pd data = {"MACCS keys":{"Tanimoto": maccs_Tanimoto, "Dice":maccs_Dice}, "Topological FP":{"Tanimoto": rdkit_Tanimoto, "Dice":rdkit_Dice}, "Morgan FP(ECFP-like)":{"Tanimoto": morgan_ECFP_Tanimoto, "Dice":morgan_ECFP_Dice}, "Morgan FP(ECFP-like/bit vector)":{"Tanimoto": morgan_ECFP_bit_Tanimoto, "Dice":morgan_ECFP_bit_Dice}, "Morgan FP(FCFP-like)":{"Tanimoto": morgan_FCFP_Tanimoto, "Dice":morgan_FCFP_Dice}, "Avalon FP":{"Tanimoto": avalon_Tanimoto, "Dice":avalon_Dice}, "AtomPair FP":{"Tanimoto": pair_Tanimoto, "Dice":pair_Dice}, "AtomPair FP(bit vector)":{"Tanimoto": pair_bit_Tanimoto, "Dice":pair_bit_Dice}, "Torsion FP":{"Tanimoto": tts_Tanimoto, "Dice":tts_Dice}, } FP_df = pd.DataFrame(data).T
計算結果は以下のような感じ・・・・
Tanimoto係数とJaccard係数は(ほぼ)同じ・・・?
今回計算した7つのフィンガープリント全てで 「Dice 係数 > Tanimoto係数」となっています。
定義上、(Simpson ≧) Dice ≧ Tanimoto(Jaccard) となるようです。
EribulinとHalichondrin B の Fingerprintをそれぞれ a, b としてみると、
・Eribulinのほぼ全体が共通部分に含まれるので 「a ≒ c」 より 「Tanimoto係数 ≒ c / b」
・EribulinはHalichondrin B の(ほぼ)部分構造なので 「a < b」として 「Dice係数 > {2c / (b + b) } = c/b 」
と考えてみても、まあ、なるほどそうなるかなあ・・・という感じです。
ところで構造式のみをみると、共通部分構造(c)は8割くらいの構造を占めているので、
係数は全て0.8程度にはなるのではないかと思いましたが、意外と0.5前後の値もあります。
特にAtomPairやMorgan フィンガープリントはそこまで大きな値ではないのが気になります。
これらのフィンガープリントでは構造のどのあたりを考慮して計算されているのでしょうか?
各原子の寄与を可視化する機能がある、とのことなので試して見たいと思います。
4. フィンガープリントへの各原子の寄与
では早速AtomPairから・・・参照分子をEribulin、描画する分子をHalicondrin B として可視化します。
from rdkit.Chem.Draw import SimilarityMaps #AtomPair(規格化) fig1, maxweight = SimilarityMaps.GetSimilarityMapForFingerprint(eri_mol, hali_mol, SimilarityMaps.GetAPFingerprint) #AtomPair(ウェイトの絶対値を利用) weight = SimilarityMaps.GetAtomicWeightsForFingerprint(eri_mol, hali_mol, SimilarityMaps.GetAPFingerprint, metric=DataStructs.DiceSimilarity) fig2 = SimilarityMaps.GetSimilarityMapFromWeights(hali_mol, weight, size=(400, 400))
規格化も絶対値もほぼ同一のようにみえるので、sizeの設定可能な絶対値を利用して行きたいと思います。
次にTorsionとMorgan Fingerprint、Topological Fingerprint(RDKit Fingerprint)を見て見たいと思います。
#Torsionフィンガープリント weight = SimilarityMaps.GetAtomicWeightsForFingerprint(eri_mol, hali_mol, SimilarityMaps.GetTTFingerprint, metric=DataStructs.DiceSimilarity) fig3 = SimilarityMaps.GetSimilarityMapFromWeights(hali_mol, weight, size=(400, 400)) #Morganフィンガープリント(デフォルトは[radius=2, fpType='bv', useFeatures=False]なので ECFP-like/bit vector に対応) weight = SimilarityMaps.GetAtomicWeightsForFingerprint(eri_mol, hali_mol, SimilarityMaps.GetMorganFingerprint, metric=DataStructs.DiceSimilarity) fig4 = SimilarityMaps.GetSimilarityMapFromWeights(hali_mol, weight, size=(400, 400)) #Topological フィンガープリント (RDKitFingerprint) weight = SimilarityMaps.GetAtomicWeightsForFingerprint(eri_mol, hali_mol, SimilarityMaps.GetRDKFingerprint, metric=DataStructs.DiceSimilarity) fig5 = SimilarityMaps.GetSimilarityMapFromWeights(hali_mol, weight, size=(400, 400))
下に行くほど分子全体が緑色となっています。
Dice係数はMorgan Fingerprint よりもTorsion Fingerprintのほうが大きいですが、前者の方が全体的に緑色となっています。
また、RDKit フィンガープリントではHalichondrin B のうち Eribulinと重ならない大員環ケトン外部も緑色となっています。
これらは単純に構造に落とし込めるわけではなさそうです。
そこで最後にビット配列の可視化を行っておきたいと思います。
5. フィンガープリントのビット配列を可視化
それではMorganフィンガープリントおよび RDKit フィンガープリントのビット配列を眺めて見たいと思います。(参考記事⑥)
#Morgan フィンガープリント (Eribulin) bitI_morgan_eri = {} fp_morgan_eri = AllChem.GetMorganFingerprintAsBitVect(eri_mol, 2, bitInfo=bitI_morgan_eri) print(len(bitI_morgan_eri)) ###80 morgan_tuples_eri = ((eri_mol, bit, bitI_morgan_eri) for bit in list(bitI_morgan_eri.keys())[:12]) Draw.DrawMorganBits(morgan_tuples_eri, molsPerRow=6, legends=['bit: '+str(x) for x in list(bitI_morgan_eri.keys())[:12]]) #Morgan フィンガープリント (Halicondrin B) bitI_morgan_hali = {} fp_morgan_hali = AllChem.GetMorganFingerprintAsBitVect(hali_mol, 2, bitInfo=bitI_morgan_hali) print(len(bitI_morgan_hali)) ###93 morgan_tuples_hali = ((hali_mol, bit, bitI_morgan_hali) for bit in list(bitI_morgan_hali.keys())[:12]) Draw.DrawMorganBits(morgan_tuples_hali, molsPerRow=6, legends=['bit: '+str(x) for x in list(bitI_morgan_hali.keys())[:12]])
Eribulin でが全部で80個、Halichondrin B では全部で93個のMorgan フィンガープリントがONとなっています。
最初の12個の0でないビットについて両者を比較すると下図のようになります。
それぞれ枠で囲んだ部分が差異となります。
同様にRDKitフィンガープリントでも実施します。
#RDKit フィンガープリント (Eribulin) bitI_RDK_eri = {} fp_RDK_eri = Chem.RDKFingerprint(eri_mol, bitInfo=bitI_RDK_eri) print(len(bitI_RDK_eri)) ###779 RDK_tuples_eri = ((eri_mol, bit, bitI_RDK_eri) for bit in list(bitI_RDK_eri.keys())[:12]) Draw.DrawRDKitBits(RDK_tuples_eri, molsPerRow=6, legends=['bit: '+str(x) for x in list(bitI_RDK_eri.keys())[:12]]) #RDKit フィンガープリント (Halichondrin B) bitI_RDK_hali = {} fp_RDK_hali = Chem.RDKFingerprint(hali_mol, bitInfo=bitI_RDK_hali) print(len(bitI_RDK_hali)) ###753 RDK_tuples_hali = ((hali_mol, bit, bitI_RDK_hali) for bit in list(bitI_RDK_hali.keys())[:12]) Draw.DrawRDKitBits(RDK_tuples_hali, molsPerRow=6, legends=['bit: '+str(x) for x in list(bitI_RDK_hali.keys())[:12]])
Eribulin でが全部で779個、Halichondrin B では全部で753個のMorgan フィンガープリントがONとなっています。
興味深いことにRDKit フィンガープリントではHalichondrin Bの方がONになっているフィンガープリントの数が少なくなっています。
最初の12個の0でないビットについて両者を比較すると下図のようになります。
エステル部位が差異となります。
まとめ
以上、今回はEribulinとHalichondrin B について、その構造の比較、類似性の検証をRDKitの様々な指標計算、描画をもちいて行いました。
単純に構造式を目でみて判断するのとは異なった視点で検証ができ非常に興味深い結果でした。
コピペで手一杯で、裏側でどのような計算が行われているかについては全く踏み込めませんでしたが、これはこれで良しとしていきたい!!
いつか勉強する・・・多分・・・
追記
Twitter で より立体を反映した2次元構造の描写方法(rdDepictor.SetPreferCoordGen )を教えていただきました。
Nice work :) 話の本筋から逸れますが rdDepictor.SetPreferCoordGen
— yamasaKit_ (@yamasaKit_) 2018年11月25日
を使うともう少し綺麗な描写になりますよ。それでも少し変ですがw pic.twitter.com/ZJh4418gSU
@yamasaKit_ 先生ありがとうございました!
忘れてしまいそうなのでコードを転記させていただきます。
from rdkit import Chem, rdBase from rdkit.Chem import rdDepictor from rdkit.Chem.Draw import IPythonConsole # True or False で設定変更 rdDepictor.SetPreferCoordGen(True) eri_mol # isomeric smilesから生成したオブジェクトを使用
元々の描画(右)では、大きな丸となっており環員数もわかりづらくなっていましたが、
「rdDepictor.SetPreferCoordGen」(左)を用いると、エッジがはっきりしてよりマニュアルの描画に近い印象になりました。
こちらも使っていこう
rdkit.Chem.rdDepictor module — The RDKit 2018.09.1 documentation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
上記記事の作成にあたって「化学の新しいカタチ」の以下の記事を参考にさせていただきました。
また、本記事の中に誤りがある場合、当ブログ作成者の責任です。
参考記事① : 化合物データベースPubChemをpythonで使いこなす
参考記事②: RDKitによる3次元構造の生成
参考記事③: RDKitを用いた部分構造検索とMCSアルゴリズム
参考記事④: RDKitでの構造式描画を詳しく解説
参考記事⑤: RDKitでフィンガープリントを使った分子類似性の判定
参考記事⑥: RDKitでフィンガープリントの可視化
RDKit documentation① 構造の描画
rdkit.Chem.Draw package — The RDKit 2018.09.1 documentation
RDKit documentation② MACCS Keys
http://www.rdkit.org/docs/GettingStartedInPython.html#maccs-keys
RDKit documentation③ Topological Fingerprints
http://www.rdkit.org/docs/GettingStartedInPython.html#topological-fingerprints
RDKit documentation④ Morgan Fingerprints
http://www.rdkit.org/docs/GettingStartedInPython.html#morgan-fingerprints-circular-fingerprints
RDKit documentation⑤Avalon フィンガープリント
http://rdkit.org/docs/source/rdkit.Avalon.pyAvalonTools.html
RDKit documentation⑥Similarity Maps using Fingerprints
https://www.rdkit.org/docs/GettingStartedInPython.html#generating-similarity-maps-using-fingerprints
集合類似度(係数)に関して参考にしたページ
https://www.conectron-networks.com/conectron/52
https://mieruca-ai.com/ai/jaccard_dice_simpson/
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~