オルト位置換基で絞り込む話

前回の記事で、SMARTS記法の扱い方がなんとなくわかりました。本題にもどって、ライブラリをオルト位に置換基の入ったビフェニルで絞り込みたいと思います。

その前に、現在残っている化合物数の確認です。

	Enamine_Premium	Enamie_Advenced	Enamine_HTS	UOS_HTS	total
指標でのフィルタリング後	4060	37431	414562	106948	563001
ビフェニル有	12	329	3182	697	4220 (0.7%)
ビフェニル無	4048	37102	411380	106251	558781 (99.3%)

合計4000個程度残っています。このうちオルト位に置換基が入ったものはいくつあるでしょうか？

オルト位置換基による絞り込み

PandasToolsでSDFをよみこむ

４種類のライブラリから抽出した、ビフェニルを部分構造として含む化合物のリストを連結（extend()を使った）し、SDFファイル「biphenyl_library.sdf」として保存しました。こちらをRDKitのPandasToolsをもちいてDataFrameでよみこみます。

PandasToolsの使い方に関しては化学の新しいカタチさんの記事「RDKitのPandasToolsでデータ分析を加速する」を参考にさせていただきました。

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem import PandasTools
import pandas as pd

biphenyl_df = PandasTools.LoadSDF('./biphenyl_library.sdf')
# 数を確認
len(biphenyl_df)
# 4220

こんな感じになりました。

Pandas上で部分構造検索

PandasToolsを使うと「>=（ge比較演算子）」で部分構造検索が行えるとのことですので試しに利用して見ます

ortho_biphenyl= Chem.MolFromSmarts('c1ccc(c(*)c1)c1ccccc1')
ortho_df = biphenyl_df[biphenyl_df['ROMol'] >= ortho_biphenyl]

len(ortho_df)
# 2353

オルト位に置換基を持つ化合物は約半数程度の2353個でした。

一致検索後のDataFrameでは、Molオブジェクトを格納したカラムにおいて、一致した部分構造が赤色でハイライトされて表示されています。

忘れずにSDFに書き出しておきます。SDFに残すプロパティを指定するため、DataFrameのカラム名をリストとして取得します。

property_names = list(biphenyl_df.columns)
print(property_names)
# ['ID', 'MW', 'MolLogP', 'NumHAcceptors', 'NumHDonors', 'NumRotatableBonds', 'ROMol', 'TPSA', 'original_id']

「ROMol」は不要なのでリストから除きます。

property_names.remove('ROMol')
print(property_names)
# ['ID', 'MW', 'MolLogP', 'NumHAcceptors', 'NumHDonors', 'NumRotatableBonds', 'TPSA', 'original_id']

出力します。

PandasTools.WriteSDF(ortho_df, 'ortho_biphenyl_library.sdf', properties = property_names)

以上で、オルト位に置換基をもつ部分構造での絞り込みが完了しました。

オルト位の置換基の構造的意義

なぜ「オルト位に置換基を持つ」ということを、絞り込みの基準として重要視するのか？ついでなので「オルト位に置換基を持つ」ということの意味について考えたいと思います。

・・・といっても立体障害でビフェニル骨格が捻れるというだけなのですが、、、、

3次元構造を眺める

オルト位の置換基の有無が３次元構造（コンフォメーション）にどのような影響を与えるか？ということを眺めて見ます。３次元構造の生成は「化学の新しいカタチ」さんの記事「RDKitによる3次元構造の生成」を参考にしました。

まずは構造を準備。無置換、オルト位にメチル基を1つもつもの、2つ持つものの3つを作りました。

biPh = Chem.MolFromSmiles('c1ccc(cc1)c1ccccc1')
ortho_Me_biPh = Chem.MolFromSmiles('c1ccc(c(C)c1)c1ccccc1')
ortho_diMe_biPh = Chem.MolFromSmiles('c1cc(C)c(c(C)c1)c1ccccc1')

biPh_list = [biPh, ortho_Me_biPh, ortho_diMe_biPh]
Draw.MolsToGridImage(biPh_list)

水素原子を付加し、ETKDG法で3次元構造を生成します。

biPh_3D_list = []

for mol in biPh_list:
    mol_H = AllChem.AddHs(mol)
    AllChem.EmbedMolecule(mol_H, AllChem.ETKDG())
    biPh_3D_list.append(mol_H)

生成した構造を眺めます。

from rdkit.Chem.Draw import IPythonConsole
import py3Dmol

v = py3Dmol.view(width=600, height=300, linked=False, viewergrid=(1,3))

for m, i in zip(biPh_3D_list, [(0,0), (0,1), (0,2)]):
    mb = Chem.MolToMolBlock(m)
    v.addModel(mb, 'sdf', viewer=i)

v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

無置換、モノ置換、ジ置換の順にねじれ角度が大きくなっているのがわかります。

もっとわかりやすくするために、重ねて見ます。

アラインメントをとって重ねる

RDKitブログの記事を参考にして３つとも重ねて見ます。

rdMolAlign.AlignMolを使って部分構造のアラインメントを取り３次元構造を重ね合わせます。重ねたい部分構造（template）を明示的に与えるため引数atomMapに無置換のベンゼン環に相当するatom numberを与えます。

こちらや、こちらを参考にエラーに対処していたらやたらとlist()を使うことになってしまった・・・

from rdkit.Chem import rdMolAlign

# 無置換フェニルのテンプレートをSmartsから作成
template = Chem.MolFromSmarts('[cH1]1[cH1][cH1][cH1][cH1]c1-c')

# テンプレートに一致する無置換ビフェニルの原子を基準にする
# タプルのタプルで得られる　
template_atoms = biPh_3D_list[0].GetSubstructMatches(template)

# アラインメントを取る
for m in biPh_3D_list:
    matched_atoms = m.GetSubstructMatches(template)
    rdMolAlign.AlignMol(m, biPh_3D_list[0], atomMap = list(zip(list(matched_atoms[0]), list(template_atoms[0]))))  
# 描く
v = py3Dmol.view(width=300, height=300)

for m in biPh_3D_list:
    mb = Chem.MolToMolBlock(m)
    v.addModel(mb, 'sdf')

v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()

静止画ではわかりにくいですが、無置換、モノ置換、ジ置換の順に二面角（Torsion Angle）が大きくなっているのがわかります。

二面角を計算してみる

それぞれ一つずつコンフォマーを取得して二面角を計算して見ます。

まず二面角を計算するために、目的の部分のatom numberを確認します。無置換ビフェニルの場合・・・

Draw.MolToImage(biPh_3D_list[0], includeAtomNumbers = True)

atom numberがわかったので目的の４原子を指定して、２面角を求めます。

from rdkit.Chem import rdMolTransforms
conf=biPh_3D_list[0].GetConformer(0)
rdMolTransforms.GetDihedralDeg(conf, 2,3,6,7)
# 148.88976772600805

モノ置換、ジ置換も同様に計算し、0° - 90°の範囲に値を修正すると以下のようになります。

	無置換	モノ置換	ジ置換
二面角	31°	50°	81°

順に二面角が大きくなり、特にジ置換ではほぼ直角近くになっています。

フィンガープリントではどう認識される？

上記の二面角の差異がフィンガープリントとしてはどのように認識できるか、少し見てみたいと思います。Atom PairsとTopological Torsionsでそれぞれ類似性を比較した際にどのようになるか試してみます。

この二つを選んだ理由は名前の雰囲気です・・・（計算方法を理解していません。すみません）。

from rdkit import DataStructs

from rdkit.Chem.AtomPairs import Pairs
pair_fps = [Pairs.GetAtomPairFingerprint(mol) for mol in biPh_3D_list]

from rdkit.Chem.AtomPairs import Torsions
tts_fps = [Torsions.GetTopologicalTorsionFingerprintAsIntVect(mol) for mol in biPh_3D_list]

# Tanimoto係数、Dice係数をそれぞれの組み合わせで計算
pair_Dice_0_1 = DataStructs.DiceSimilarity(pair_fps[0], pair_fps[1])
pair_Tani_0_1 = DataStructs.TanimotoSimilarity(pair_fps[0], pair_fps[1])

# 残りの組み合わせ、TorsionFingerprintについても同じ計算をしたので省略

結果をまとめると以下になりました。

	無置換-モノ置換	無置換-ジ置換	モノ置換-ジ置換
AtomPair_Tanimoto係数	0.76	0.59	0.77
AtomPair_Dice係数	0.86	0.75	0.87
Torsion_Tanimoto係数	0.87	0.77	0.88
Torsion_Dice係数	0.93	0.86	0.94

Torsionの方が差異が大きくなると面白いと思っていたのですが、今回の結果ではAtom Pairの方が置換の数の差を反映しているようにも見えます。いずれにせよ、メチルの数が１つ増えるか２つ増えるかで類似度が下がっていくことに代わりわなさそうです。

オルト位置換を創薬の話に繋げようとする

オルト位に置換基を加えることで分子がねじれ、立体的な構造が変化するということがわかりました。では、このねじれを医薬品に活用した事例は？？？ということでFasiglifam(TAK-875)の構造を眺めて見たいと思います。２型糖尿病を適応として開発が行われていた武田薬品の化合物で、GPCRの一つGPR40の作動薬です。残念ながら臨床試験第三相で肝臓における安全性の懸念から開発中止となり、上市にはいたらなかったそうです。

Fasiglifam(TAK-875)の構造

なぜこちらを取り上げるのか？　構造を見ていただくとお分かりいただけると思います。 印象的なビフェニル、オルト-ジ置換構造を持っています。

Fasiglifam_mol = Chem.MolFromSmiles('CC1=CC(=CC(=C1C2=CC(=CC=C2)COC3=CC4=C(C=C3)[C@@H](CO4)CC(=O)O)C)OCCCS(=O)(=O)C')

from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

highlights = Fasiglifam_mol.GetSubstructMatch(ortho_diMe_biPh)
view = rdMolDraw2D.MolDraw2DSVG(500,500)
pre_mol = rdMolDraw2D.PrepareMolForDrawing(Fasiglifam_mol)
view.DrawMolecule(pre_mol, highlightAtoms=highlights)
view.FinishDrawing()
svg = view.GetDrawingText()

SVG(svg.replace('svg:,',''))

Fasiglifamを見出した構造活性相関の研究に関していくつかの論文が出されています。ずいぶんと以前に読んだ論文のため引用しようにも忘れてしまいました（おそらくJ.Med.ChemやACS Med. Chem. Lett.だったような・・・）。その中でビフェニル構造に置換基を加えることで２面角が変化すること、その構造の変化が活性に大きな影響を与えたことが議論されていたように思います。（いろいろぼんやりとした記憶で不確かですみません・・・）

では立体的にはどのような構造なのか？眺めて見ます。PubChemのID「cid: 24857286」をつかってpy3Dmolで描画して見ます。

化学の新しいカタチさんの記事「py3Dmolを使って化学構造をJupyter上で美しく表示する」を参考にさせていただきました。

view = py3Dmol.view(width= 800, height=300, viewergrid=(1,3), query='cid:24857286')
view.setStyle({'cartoon': {}}, viewer=(0,0))
view.setStyle({'stick': {}}, viewer=(0,1))
view.setStyle({'sphere': {}}, viewer=(0,2))
view.show()

ビフェニル構造がほぼ直角にねじれています。

受容体との共結晶構造

FasiglifamとGPR40との共結晶構造が解かれています（PDB id: 4PHU*1）。ついでに眺めて見ます。 *2*3

view = py3Dmol.view(width=680, height=480, query='pdb:4PHU', viewergrid=(1,1))
view.setStyle({'cartoon': {'color': 'spectrum'}}, viewer=(0,0))
view.setStyle({'resn': '2YB'}, {'sphere': {'colorsheme': 'default'}}, viewer=(0,0))
view.addSurface(py3Dmol.SES, {'opacity': 0.6, 'color': '#FADBD8'}, viewer=(0,0))
view.show()

とても興味深いことに、Fasiglifamはその構造の半分程度がGPR40のヘリックスの合間から外に向かって飛び出しています。

上の図では少しわかりにくいのでRCSD PDBのビューワーでリガンドと受容体の相互作用部位を拡大して眺めます。

こんな感じ・・・

ちょうど直角にねじれたビフェニル構造の部分が受容体内-外の境目に来ています。最近ではGPCRの構造が多数解かれるようになってきますが、この構造が公開された当時（2014年）はまだそこまで多くなく、さらにリガンドの結合様式も見たことがないものでこんな結合の仕方あるのか！！！と、とても驚いたのを覚えています。

まとめ

以上、今回はオルト位に置換基を含む化合物のライブラリからの抽出と、その構造的意義についてざっと眺めて見ました。オルト位に置換基を加えることでライブラリ化合物は2400弱にまで絞り込むことができました。後半、古い記憶を辿ってFasiglifamの話などをし始めたのは、ライブラリスクリーニングの方法に行き詰まったからです・・・困難から逃走するタイプ・・・

*1:Crystal structure of Human GPR40 bound to allosteric agonist TAK-875 DOI: 10.2210/pdb4PHU/pdb

*2:High-resolution structure of the human GPR40 receptor bound to allosteric agonist TAK-875. (2014) Nature 513: 124-127

*3:AS Rose et al. (2018) NGL viewer: web-based molecular graphics for large complexes. Bioinformatics dio:10.1093/bioinformatics/bty419