ファーマコフォアを作ろうとする話

前回の記事でRDKitのファーマコフォアがどのような化学的構造を認識しているのか、その定義を眺めました。実際の化合物に適応するとどうなるのか、活性化合物のデータセット （共闘用シェアデータ ） から取り出した分子（マクロサイクル型を除いたもの）で試したいと思います。

1. 活性化合物で検証

1-1. フィーチャーの探索

まずは活性化合物群を読み込みます。

import os
from rdkit import RDConfig
from rdkit.Chem import AllChem
from rdkit import Chem

chain_suppl = Chem.SDMolSupplier('./chain_compounds.sdf', removeHs=False)
chain_mols = [x for x in chain_suppl if x is not None]

FDefファイルを読み込んでフィーチャーファクトリーを作成します。

fdef = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'))
print(fdef.GetFeatureFamilies())
#('Donor', 'Acceptor', 'NegIonizable', 'PosIonizable', 'ZnBinder', 'Aromatic', 'Hydrophobe', 'LumpedHydrophobe')

一つ目の化合物を取り出し、特性基（フィーチャー）の探索を行います。

test_m = chain_mols[0]
# そのままfeatureを検索
rawFeats = fdef.GetFeaturesForMol(test_m)
len(rawFeats)
#17

17個のフィーチャーが認識されました。各フィーチャーがどのフィーチャーファミリーのものか確認します。

feat_fam = []
for feat in rawFeats:
    feat_fam.append(feat.GetFamily())

Python標準ライブラリcollectionsを使って、上記のリストの各要素の個数を確認します。*1

import collections
c = collections.Counter(feat_fam)
print(c)
# Counter({'Hydrophobe': 5, 'Acceptor': 4, 'Aromatic': 3, 'Donor': 2, 'LumpedHydrophobe': 2, 'PosIonizable': 1})

用意されている８種類ファミリーのうち「'Hydrophobe', 'LumpedHydrophobe', 'ZnBinder'」の３つは特殊なように思えたので、より基本的な残りの５種類「'Donor', 'Acceptor', 'NegIonizable', 'PosIonizable', 'Aromatic'」を使うことにします。*2

keep = ('Donor','Acceptor','NegIonizable','PosIonizable','Aromatic')
featLists = []
for feat in rawFeats:
    # ファミリーを取得し、keepに含まれている時のみリストに追加
    if feat.GetFamily() in keep:
        featLists.append(feat)
len(featLists)
# 10

Hydrophobe（5個）とLumpedHydrophobe（2個）を除くと、基本的なフィーチャーは10個となりました。

1-2. フィーチャーの可視化

@iwatobipen先生の記事visualize-pharmacophore-in-rdkitを参考に、フィーチャーとして認識された構造を眺めます。

atom_ids_list = [] 
legend_list = []
for feat in featLists:
    # feat_family_list.append(feat.GetFamily())
    atom_ids_list.append(feat.GetAtomIds())
    # feat_type_list.append(feat.GetType())
    legend_list.append(feat.GetFamily() + ':' + feat.GetType())

from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw

Draw.MolsToGridImage([test_m]*10, molsPerRow=5, subImgSize=(200, 200), legends=legend_list, highlightAtomLists=atom_ids_list)

水素結合ドナー(Donor)としてはアミンやアミドのNH、アクセプターとしてはアミド結合のカルボニル酸素、エーテル結合の酸素原子が認識されています。また、ピリジンのNは塩基性があるので、プラス電荷を帯びる塩基性グループ（PosIonizable:BasicGroup）としても良さそうですが、デフォルトのフィーチャーファクトリーでは、水素結合アクセプターとして認識されるようです（下段左端）。

1-3. SMARTSの表現を再確認

前回フィーチャーのSMARTSを眺めましたが、ピリジンのNを認識していると思われる部分を再度確認しましょう。 AcceptorファミリーのフィーチャーAcceptor.SingleAtomAcceptor、「[n;+0;!X3;!\$([n;H1](cc)cc)]」あたりでしょうか?

SMARTSviewerを使ってみます。*3

「1」で囲まれた芳香族NH（ピロールのような構造？）ではない（「!:否定」）、芳香族Nと解釈できそうなので、ピリジンNを認識しそうです。

もっとかっこよく眺めたい！ということで、RDKitとPyMolの組み合わせを使ってみたいと思います。

1-4. PyMol x RDKit

@iwatobipen先生の記事Visualize pharmacophore with RDKitとGreg先生のノートShow_Ph4_Features_in_PyMOL.ipynbを参考にさせていただきました。

まず、ターミナルでpymolを xml-rpc（remote procedure call）サーバーモードで起動します。

# ターミナル
pymol -R

jupyter notebookに戻って以下を実行します。

from IPython import display
from rdkit.Chem import PyMol
IPythonConsole.ipython_useSVG=True

# viewerを作成
v = PyMol.MolViewer()

# viewerを一旦綺麗に掃除
v.DeleteAll()

# 先の構造を表示
v.ShowMol(test_m)

先に起動していたpymolを眺めると構造式が出ていました。グリグリ回してポーズを決めてから次のコマンドをjupyter notebookで実行するとノート上に静止画（PNG）が表示されました。

png=v.GetPNG()
display.display(png)

PyMol綺麗・・・

PyMolが無事動いたので、フィーチャーの可視化に進みます。

1-5. PyMolでフィーチャーを可視化

フィーチャーを表示させるための関数を作成します。さっぱりわからないのでGreg先生のipynbからコピペさせていただきます。

RDKitのFeatures.ShowFeatsモジュールを使っているようですが、CGOなど全くわからない・・・

# コピペ
from rdkit.Chem.Features import ShowFeats

def showMolFeats(mol,factory,viewer,name="mol",showOnly=True):
    featLabel = f'{name}-feats'
    dirLabel = f'{name}-feats-dirs'
    if(showOnly):
        ShowFeats._globalSphereCGO = []
        ShowFeats._globalArrowCGO = []
    # workaround for what is either a bug in the way pymol handles CGOs
    # or a gap in my understanding of how it should work
    viewer.server.resetCGO(featLabel)
    viewer.server.resetCGO(dirLabel)
    viewer.server.sphere((0, 0, 0), .01, (1, 0, 1), featLabel)
    viewer.server.cylinder((0, 0, 0), (.01, .01, .01), .01, (1, 0, 1), dirLabel)
    ShowFeats.ShowMolFeats(mol,factory,v,name=name,featLabel=featLabel,dirLabel=dirLabel,
                           useFeatDirs=False,showOnly=showOnly)
    viewer.server.renderCGO(ShowFeats._globalSphereCGO,featLabel,1)
    viewer.server.renderCGO(ShowFeats._globalArrowCGO,dirLabel,1)

フィーチャーを眺めたい分子（molオブジェクト）はtest_m、フィーチャーファクトリーはfdef、ビューワーはv、という名前にしているので、これらを引数にして上の関数を使います。

showMolFeats(test_m,fdef,v)

pymolに移動して回転してから、jupyter notebookに表示します。

png=v.GetPNG()
display.display(png)

格好いい！！

認識されたフィーチャーが球として表示されています。

先にMolsToGridImageで並べて見ていた時との違いは、黄土色の球の部位です。

表示されないようにFeatureFamilyから除いた「Hydrophobe」や「LumpedHydrophobe」のようです。

2. 共結晶構造をつかって検証

2-1. 共結晶構造のリガンドでフィーチャーを確認

かなりいい感じでフィーチャーを眺められるようになってきたので、PLIPを検証した記事で眺めた共結晶構造（PDB id: 5NIX）のリガンド（ID: 8YQ）でも、同じことをしてみようと思います。

まず、PDBからリガンドの構造をSDFでとってきて、Chain Aに含まれるリガンドの座標のみ残して保存しました（8YQ_noH.sdf）。

suppl_8YQ = Chem.SDMolSupplier('./8YQ_noH.sdf')
mol_8YQ = [x for x in suppl_8YQ if x is not None][0]

showMolFeats(mol_8YQ,fdef,v)
png=v.GetPNG()
display.display(png)

なかなかごつい感じになりました。ほとんどの原子がフィーチャーに認識されていて、単なるBall&Stick表示のようにも見えてしまいます。

このフィーチャーのなかから、PLIPで探した「複数の共結晶構造に共通する相互作用」に該当する「フィーチャー」を取り出して組み合わせればファーマコフォアのモデルができるはず！

2-2. 共結晶構造の相互作用をPyMolで眺める

まずは、取り出した相互作用形式を再度確認します。

「PDB id:5NIX(Ligand:8YQ)」の場合・・・

	残基番号	アミノ酸	距離	相互作用		残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)					8YQ(chain C/D)	54C	ILE	3.89	Hydrophobic
	56B	TYR	3.82	Hydrophobic		56C	TYR	3.42/3.62	Hydrophobic
	66B	GLN	3.72	Hydrophobic
	115A/B	MET	3.70/3.80	Hydrophobic		115C/D	MET	3.85/3.82	Hydrophobic
	121A/B	ALA	3.69/3.75	Hydrophobic		121C/D	ALA	3.96/3.57	Hydrophobic
	123A	TYR	3.80/3.52/3.85	Hydrophobic		123D	TYR	3.72	Hydrophobic
	124A	LYS	5.03/3.18	π-cation/Salt Bridges		124D	LYS	3.70	Water Bridges/Salt Bridges

このうちChain A/B の残基（左半分）についてpymol上でハイライトして眺めて見ます。

pymolの設定などに関しては或る化みす途さんのブログ話題のXofluzaの変異をPymolで見てみるを参考にさせてきただきました。

「5nix.pdb」を表示させた後、

1. set cartoon_transparency, 0.8
2. リガンドを中心に持ってくる（Command+左クリック）
3. 水を消す（object panelのallでH(hide)からwatersを選択）
4. 見たい残基を選択（sequenceを表示「右下Sをクリック」して順番に選択）
5. 残基をLine表示（object panelの(sele)でS(show)からas→wire→lines）

とします。スナップショットを保存するには、右上の「Draw/ray」というボタンから出力したいサイズを選んでDraw→save image to fileとしてやれば良いそうです。以下のようになりました。

2-3. ファーマコフォア作成に向けた偽アトムの設定

ファーマコフォアを作るには、フィーチャー間の距離や角度といった情報をつかって定義するようです。 フィーチャーが原子一つに帰着できるなら良いですが、芳香環みたいな複数の原子からなるグループ（官能基）の場合は、どの点を基準にして距離計測を行うかで結果が変わりそうです。

先に描画した図では芳香環（Aromatic フィーチャー）の真ん中に球 が表示されていました。同じような点（偽原子）を作って見たいと思います。

PyMolWikiを参考にPseudoatomというのを作って見ました。

1. 右下緑文字のSelectingをクリックしてAtomsに変更
2. 芳香環の各原子を選択
3. pseudoatom arom1, sele

と打ち込みました。これで真ん中の芳香環の中に、新しい点(arom1という名前のobjectにしました)ができました

同様にして他の環にもpseudoatomを作っていきました。

2-4. 相互作用を順番に

~ Tyr / Gln ~

だいたい準備ができた感じがするので一つ一つの相互作用部位を見ていきます。

	残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)	56B	TYR	3.82	Hydrophobic
	66B	GLN	3.72	Hydrophobic

距離（点線）はWizardのMeasurementを選択してから、atomをクリックしていくと表示されました。 Chain Bの残基56 Tyrと66 GLNは芳香環と疎水性相互作用をしているようです。

~ Ala / Met ~

	残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)	115A/B	MET	3.70/3.80	Hydrophobic
	121A/B	ALA	3.69/3.75	Hydrophobic

Chain Aの残基115 METとChain Bの残基121 ALA、Chain Bの残基115 METとChain Aの残基121 ALAがそれぞれ組み合わせとなって別々の芳香環と疎水性相互作用をしているようです。

~ Tyr / Lys ~

	残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)	123A	TYR	3.80/3.52/3.85	Hydrophobic
	124A	LYS	5.03/3.18	π-cation/Salt Bridges

Chain Aの残基123 TYR（黄色）は芳香環と疎水性相互作用、Chain Aの残基124 LYS（オレンジ色）は芳香環との相互作用に加えて、末端カルボキシ基とも相互作用が見られるようです。

3. 共結晶構造をもとにファーマコフォアを作成

3-1. フィーチャーの選択と相対配置

タンパク質残基との相互作用を眺めることで、リガンド側で大事そうな部位がわかってきました。

SAR news No.19の記事を参考に、aromatic ２つとacceptor 1つの３点のフィーチャーをファーマコフォアとして選んで見ました。

（大事な点が３つあります！っていうと格好いいから３にした）

距離をPymol上で確認して見ます。

1. Mouse Modeを3-Button Editingに変更
2. ２点を選ぶと距離、３点では角度、４点で２面角が表示される

結果をまとめると以下のようです。

便宜的にAcceptorはカルボキシ基の炭素原子との距離を測りました。 酸素原子を使った場合の距離はそれぞれ、5.1Åが[5.0, 5.6] 、10.4Åが[9.8, 11.2]となります。

着目するフィーチャー３点と、その位置関係が決まりました。このファーマコフォアを定義するため、フィーチャーの座標（位置情報）を確認します。

3-2. 共結晶構造におけるフィーチャーの座標の取得

pymolでshowfeatを使う際に、--writeFeatsのオプションをつかうことで、フィーチャーの座標を出力できるようになるそうです。@iwatobipen先生の記事Visualize pharmacophore with RDKitのコードをまたしてもそのままコピペさせていただきました。

(Greg先生の先の関数「showMolFeats」のどこかに--writeFeatsを書き込めばよかったのかもしれませんが私の能力ではわかりませんでした・・・)

import subprocess
import pprint
showfeatpath = os.path.join(RDConfig.RDCodeDir, 'Chem/Features/ShowFeats.py')

v = PyMol.MolViewer()
v.DeleteAll()
process = subprocess.Popen(['python', showfeatpath, '--writeFeats','./8YQ_noH.sdf'], stdout=subprocess.PIPE)
stdout = process.communicate()[0]

res = stdout.decode('utf-8').replace('\t', ' ').split('\n')
pprint.pprint(res)

下の図のように、座標が出力されます。（一部抜粋。全フィーチャーが並んでいる。）

ここから必要なフィーチャーの座標を抜き出します。フィーチャーを間違えないようにAtom_Idを確認しておきます。

Draw.MolToImage(mol_8YQ, includeAtomNumbers = True)

Atomid と Familyから該当の座標を取り出してきます。Draw.MolToImageで表示させたAtomNumberと、ShowFeatのAtom Idは何故か１ずつ番号がずれていました。理由はよくわかりません・・・とりあえず先に進みます。

芳香環２つの座標は以下・・・

• 'Aromatic -7.625 10.407 -21.653 1.0 # 1 10 11 36 35 34'
• 'Aromatic -1.851 11.051 -24.019 1.0 # 6 8 37 38 39 19'

アクセプター２点（カルボキシル基酸素原子）の間が、ちょうど'NegIonizable'として認識されているのでこちらの座標を使いたいと思います。

• 'NegIonizable -0.100 13.586 -28.137 1.0 # 45 28 20'
• 'Acceptor 0.862 14.067 -27.897 1.0 # 20', 'Acceptor -1.171 13.517 -28.363 1.0 # 28'

フィーチャーの座標位置とその位置関係が決まりました。いよいよRDKitを使ってファーマコフォアを設定したいと思います。

4. RDKitでファーマコフォアを設定

4-1. ファーマコフォア設定方法の流れ

SAR News No.19の吉森氏による記事「Chemoinformatics Toolkits を用いた創薬システム開発における ラピッドプロトタイピング」での流れを確認します。

1. 座標を使って特性基を定義
2. 距離情報を使ってファーマコフォアを設定
3. ファーマコフォアサーチ
   1. 対象のフィーチャーをそもそも含んでいるか？
   2. フィーチャー間の距離が条件を満たすか？
   3. 距離情報を拘束条件にして3D構造を発生させる。

以上の流れを辿って見ます。テストとして、記事冒頭に用共闘用シェアデータから取り出した分子を用います。

4-2. ファーマコフォアの設定を実践

フィーチャーの定義にはChemicalFeaturesモジュールを使うようです。

ChemicalFeatures.FreeChemicalFeature(Feature Family, Feature Type, 位置) とすることで、フィーチャーを定義できます（Typeは省略可能）。 位置の指定にはGeometryのPoint3Dをもちいます。

from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.Pharm3D import Pharmacophore
from rdkit import Geometry

feat_1=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-7.625, 10.407, -21.653))
feat_2=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-1.851, 11.051, -24.019))
feat_3=ChemicalFeatures.FreeChemicalFeature('Acceptor',Geometry.Point3D(-0.100, 13.586, -28.137)) 
feats = [feat_1,feat_2,feat_3]

Aromatic ２つとAcceptor１つが定義できました。

続いてPharm3D.Pharmacophoreモジュールを使って距離情報を設定します。先に定義したフィーチャーのリストでファーマコフォアを設定した後、フィーチャー間の距離のとりうる範囲を下限値（setLowerBound）、上限値（setUpperBound）という形で与えます。

SAR Newsの記事では距離の許容範囲の値として、下限値に0.5（d_lower）、上限値に1.5（d_upper）が用いられていました。

pcophore = Pharmacophore.Pharmacophore(feats) # ファーマコフォアの設定
d_upper = 1.5  # 特性基間の距離の許容範囲(上限値)
d_lower = 0.5 # 特性基間の距離の許容範囲(下限値)

# feat_1とfeat_2の距離 6.3Åを基準に、下限と上限を設定
pcophore.setLowerBound(0,1, 6.3-d_lower)
pcophore.setUpperBound(0,1, 6.3+d_upper)

# Acceptor(feat_3)の代表点の選び方によってfeat_2との距離は[5.0~5.6]の値をとる
pcophore.setLowerBound(1,2, 5.0-d_lower)
pcophore.setUpperBound(1,2, 5.6+d_upper)

# 同様にfeat1とfeat_3は[9.8~11.2]の値をとる
pcophore.setLowerBound(0,2, 9.8-d_lower)
pcophore.setUpperBound(0,2, 11.2+d_upper)

print(pcophore)

以上でファーマコフォアの設定が終わりました。printで出力したPharmacophoreを見てみると、距離行列となっており、対角線右上三角に上限値、左下三角に下限値の情報を持つようです。

4-3. ファーマコフォアサーチ ~step 1. フィーチャーの有無~

それではテスト分子を使ってファーマコフォアサーチを行います。３次元構造を扱うので水素原子を付加しておきます。

test_mH = AllChem.AddHs(test_m)
AllChem.Compute2DCoords(test_mH)

まずはフィーチャーをそもそも持っているか？

Parm3D.EmbedLibモジュールのMatchPharmacophoreToMolを使います。

EmbedLib.MatchPharmacophoreToMol(molオブジェクト,フィーチャーファクトリー, ファーマコフォア)と、することで対象のmolオブジェクトにファーマコフォアのマッピングを行い、可能なマッピングのリストを作成します。

戻り値は2要素のタプルで、

1. 全てのフィーチャーが見つけられたか否かのブール値
2. フィーチャーの並びのリスト

となっています。

テスト分子をtest_mH、フィーチャーファクトリーをfdef、ファーマコフォアをpcophoreとして実行します。

from rdkit.Chem.Pharm3D import EmbedLib

match, mList = EmbedLib.MatchPharmacophoreToMol(test_mH, fdef, pcophore)
print(match)
# True

どうやら全てのフィーチャーのマッチングはOKだったみたいです。確認のためリストからフィーチャーの情報を取り出して見ます。

print(len(mList))
#3
print(len(mList[2]))
#4
print(mList[2][0].GetFamily())
#Acceptor
print(mList[2][0].GetType())
#SingleAtomAcceptor
print(mList[2][2].GetAtomIds())
#(15,)

心配なのでAcceptorについて描画して見ます。

highlight_list = [mList[2][x].GetAtomIds() for x in range(len(mList[2]))]
Draw.MolsToGridImage([test_m]*4, molsPerRow=4, subImgSize=(200, 200), highlightAtomLists=highlight_list)

うまく機能していそうです。

4-4. ファーマコフォアサーチ ~step 2. 距離条件~

ついで距離条件の検証を行います。

まず、rdDistGeomモジュールのGetMoleculeBoundsMatrixを使って分子の距離行列を取得します。

Parm3D.EmbedLibモジュールのGetAllPharmacophoreMatchesを使って、先に取得した分子の距離行列が、ファーマコフォアで定義されている距離の条件を満たすか否かを判定します。

from rdkit.Chem import rdDistGeom
# 距離行列の取得
bounds = rdDistGeom.GetMoleculeBoundsMatrix(test_mH)
#ファーマコフォアのマッチング 
pList =EmbedLib.GetAllPharmacophoreMatches(mList,bounds,pcophore) 
print(len(pList))
# 4

ファーマコフォア（３つのフィーチャー）の条件を満たす組み合わせは４つあったようです。一つ目の組み合わせのフィーチャーを構成する原子のIDを取得して見ます。

AtomIds_list = []

for i in range(len(pList[0])):
    p = pList[0][i]
    print(p.GetFamily(), ':', p.GetAtomIds())
    AtomIds_list.append(p.GetAtomIds())

Draw.MolsToGridImage([test_m]*3, subImgSize=(200, 200), highlightAtomLists=AtomIds_list)

他の組み合わせでは以下になりました。

想定していたのと異なり、一番遠くの芳香環を使った組み合わせとなっていました。
ところで、このマッチングは距離を条件につかっているのですが、検索対象としているmolオブジェクトは３次元構造を生成させていません（pymolで眺めた時もぺちゃんこだった）。

SAR Newsの記事の該当部分を引用すると

一般的なファーマコフォアサーチにおいては、分子の3D構造生成後、ファーマコフォアのサーチを行うが、RDKitにおいては、事前に分子の3D構造を生成させるのではなく、ファーマコフォアを満たす3D構造が生成できるか否かを判定基準として、ファーマコフォアのサーチを行うことができる。

とのことだそうです。それでは 距離条件をみたすような３次元構造が本当に作れるのかどうか、検証したいと思います。

4-5. ファーマコフォアサーチ ~step 3. 3D構造の発生~

距離情報を拘束条件にして3D構造を発生させます。まずはマッチしたフィーチャーの原子IDのリストを作成します。コピペ・・・

num_match = len(pList)
phMatches = []
for i in range(num_match): 
    num_feature = len(pList[i])
    phMatch = []
    for j in range(num_feature):
        phMatch.append(pList[i][j].GetAtomIds())
    phMatches.append(phMatch)

距離情報を拘束条件にした3D構造の発生にはPharm3D.EmbedLibのEmbedPharmacophoreを使います。

引数がたくさんありますが、対象のmolオブジェクト(mol)に対して、ファーマコフォアのフィーチャーを構成する原子のid(atomMatch)をつかって、ファーマコフォア(pcophore)の条件を満たすように3次元構造を生成(embedding)することができるようです。 生成させる構造の最大数をcountで設定します。（silentはよくわかりません・・・）

戻り値は３要素のタプルで、

1. ファーマコフォアに合うように調整された拘束条件の距離行列
2. ３次元構造（コンフォマー１つを有する分子）のリスト
3. ３次元構造の生成に失敗した回数

となっています。

生成させたのち、ETKDG法を使って３次元構造の最適化を実施します。

confs = []

for phMatch in phMatches:
    bm,embeds,nFail =EmbedLib.EmbedPharmacophore(test_mH, phMatch,
                                                 pcophore,count=20,
                                                 silent=1)
    print("Generated embeds: ",len(embeds))
    print("Failed attempts: ",nFail)
    for embed in embeds:
        AllChem.EmbedMolecule(embed, AllChem.ETKDG())
        confs.append(embed)

フィーチャーの組み合わせによって、構造生成の成否回数が異なっているようです。

SAR newsの記事ではUniversal Force Field法を使って構造最適化を行うと記載してあったのですが、以下のようなエラーメッセージが出てしまいました。

RuntimeError                              Traceback (most recent call last)
<ipython-input-254-6895836e8c95> in <module>()
      8     print("Failed attempts: ",nFail)
      9     for embed in embeds:
---> 10         AllChem.UFFOptimizeMolecule(embed)
     11         confs.append(embed)

RuntimeError: Invariant Violation
    bad direction in linearSearch
    Violation occurred on line 234 in file Code/Numerics/Optimizer/BFGSOpt.h
    Failed Expression: status >= 0
    RDKIT: 2018.09.1
    BOOST: 1_65_1

よくわからない・・・。飛ばそう・・・

念のためETKDGで最適化した３次元構造を一つ眺めて見ます。

from rdkit.Chem.Draw import IPythonConsole
import py3Dmol
IPythonConsole.drawMol3D(confs[0])

いい感じにできていそうです。

5. ファーマコフォアを基準に共結晶構造に重ね合わせ

活性化合物に関して、ファーマコフォアの基準を満たす３次元構造の取得までが確認できました。この3D構造は、共結晶構造のリガンド構造と類似しているはず・・・なので、共結晶構造に重ね合わせてタンパク質と一緒に描いて見たいと思います。

5-1. PDB座標にアラインメント

ファーマコフォア作成に用いたリガンドの座標はPDBの共結晶構造を元にしているので、この座標をテンプレートとしてアライメントをとります。

Numerics.rdAlignmentモジュールのGetAlignmentTransformを使って、参照分子の座標(refPoints)と目的の分子の座標(probePoints)とで、RMSDが最小となる最適なアラインメントを計算します。

rdkit.Numerics.rdAlignment.GetAlignmentTransform(refPoints, probePoints)の戻り値として、SSD値（Sum of Squared Difference）と、4x4変換行列（transform matrix）のアレイを要素として持つ、2-タプルが得られます。

得た変換をAllChem.TransformMolを使って、適応することで、アラインメントをとった座標を持つmolオブジェクトを得ます。*4

# 生成した構造のフィーチャーを取得
match, mList_0 = EmbedLib.MatchPharmacophoreToMol(confs[0], fdef, pcophore)
bounds_0 = rdDistGeom.GetMoleculeBoundsMatrix(confs[0])
pList_0 =EmbedLib.GetAllPharmacophoreMatches(mList_0,bounds_0,pcophore) 

# ファーマコフォアで設定したフィーチャーのリストを参照(ref)として使用する
refFeats = feats

# pList[0]で認識されているフィーチャーを
#ファーマコフォアのフィーチャーの定義の順番に並べ替えたリストにする
probFeats = [pList_0[0][1], pList_0[0][0], pList_0[0][2]]

# ref、prob、それぞれのフィーチャーの座標を取得し、リストにする
probPts = [list(x.GetPos()) for x in probFeats]
refPts = [list(x.GetPos()) for x in refFeats]

# ２つの座標をつかって変換行列を取得
ssd,tform = rdAlignment.GetAlignmentTransform(refPts, probPts)
# 変換行列を使って、生成した構造の座標を変換
AllChem.TransformMol(confs[0], tform)

#変換したmolオブジェクトをmolファイルとして出力
Chem.MolToMolFile(confs[0], 'conf.mol')

出力したmolファイルと、元にした共結晶構造 5NIX.pdb を同時にpymolで表示させました。

もともとのリガンドの色がシアンで、重ね合わせたリガンドが緑色です。 なかなかいい感じで重なっているようにみえます。今回ファーマコフォアとして指定しなかった末端の芳香環についてはズレが大きいように見える一方で、ファーマコフォアマッチングした部分は重なりがよく見えます。うまくいった！！！　・・・のか？？？

まとめ

以上、今回はファーマコフォアの作成とテスト分子を使った動作検証を行って見ました。ファーマコフォアの基準を満たす３次元構造の取得までが確認できましたので、この取得ができるか否かまでを、ファーマコフォアを満たすことができるか否かという判定基準として活用していけば多分いい感じでスクリーニングになるはず。。。

あまりにもわからないことが多すぎて先生方の記事からひたすらコピペを繰り返しましたが、それでも正しく使えているのか自信がまったくありません。特にBounds Matrixが理解できなかった・・・（3次元に関する条件にしては行列の各要素が3個でもないし・・・）。

色々と間違いがあると思うのでご指摘いただければ幸いです。

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

*1:参考: PythonのCounterでリストの各要素の出現個数をカウント

*2:参考: RDKitブログ記事Using Feature Maps

*3:citatition:

• SMARTS viewer https://smartsview.zbh.uni-hamburg.de
• 1. Schomburg, H.-C. Ehrlich, K. Stierand, M.Rarey; From Structure Diagrams to Visual Chemical Patterns, J. Chem. Inf. Model., 2010, 50 (9), pp 1529-1535

     *4:参考: RDKitのメーリングリスト