部分構造で絞り込む話 - magattacaのブログ

Twitterで拝見した以下の処理の手順を参考に、指標を用いた粗い絞り込みまでを実施しました。

処理の手順

	処理	意図
①	分子量の範囲を設定	粗い絞り込み
②	部分構造で絞り込み	活性化合物を参照したLBVS（？）
③	粗いドッキング	タンパク質との相互作用を予測するSBVS（？）
④	ファーマコフォアを指定したより詳細なドッキング	SBVS ２回目（？）

指標による絞り込み基準

指標	分子量	LogP	水素結合供与体数	水素結合受容体数	回転可能結合数	極性表面積
Rule of 5	≦500	≦5	≦5	≦10
基準	>300	>3		>3	>3	>60

今回は指標による絞り込みの次の段階、部分構造での絞り込みを行いたいと思います。
参考にさせていただいている処理の方法では、オルト位に置換基の入ったビフェニルで絞り込んだとのことでしたが、いきなり置換基なども含めた処理を行うのは難しいので、まずはビフェニル構造での絞り込みを行いたいと思います。

活性化合物のデータセットでビフェニルを探す
ライブラリ化合物で部分構造による絞り込み
部分構造検索をKNIMEで
オルト位置換の認識に失敗
まとめ

活性化合物のデータセットでビフェニルを探す

そもそもなぜビフェニルなのか？　活性化合物のデータセット（共闘用シェアデータ）中の分子（マクロサイクル型を除く）が、ビフェニル構造をもつのか検証してみます。

「化学の新しいカタチ」さんのこちらの記事「 RDKitを用いた部分構造検索とMCSアルゴリズム」を参考にRDKitの部分構造検索メソッドを利用します。

まずは必要なものをimport・・・

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw

RDKitで部分構造検索を行う際、検索構造（query）はMolオブジェクトにする必要があるということなので、まずはビフェニル構造のMolオブジェクトをSMILESから作成します。

biphenyl = Chem.MolFromSmiles('c1ccc(cc1)c1ccccc1')
Draw.MolToImage(biphenyl)

f:id:magattaca:20190203191836p:plain

きちんとビフェニル構造のMolオブジェクトが作られていそうです。こちらをテンプレートとして用い、部分構造検索を行います。

#以前取り出した鎖状分子のSDFを使用する
chain_compounds_suppl = Chem.ForwardSDMolSupplier('chain_compounds.sdf')

mols = [mol for mol in chain_compounds_suppl if mol is not None]

HasSubstructMatchをもちいて、ビフェニル構造を持つ分子と持たない分子を区別し、それぞれ別々のリストに入れていきます。

#ビフェニルありのリスト
biphenyl_compounds = []
#ビフェニルなしのリスト
w_o_biphenyl = []

for mol in mols:
    if mol.HasSubstructMatch(biphenyl):
        biphenyl_compounds.append(mol)
    else:
        w_o_biphenyl.append(mol)

print('Number of compounds with biphenyl', len(biphenyl_compounds))
# Number of compounds with biphenyl 29
print('Number of compounds without biphenyl', len(w_o_biphenyl))
# Number of compounds without biphenyl 10

39個の分子のうち、29個にビフェニルが含まれました。確かに活性化合物は多くがビフェニル構造を持つようです。

最初の一つの分子を取り出して、ビフェニル構造の確認を行います。

GetSubstructMatchesでマッチした複数の原子インデックスを取得できるそうです。

test_mol = biphenyl_compounds[0]
matched_atoms = test_mol.GetSubstructMatches(biphenyl)
print(matched_atoms)
# ((23, 17, 18, 20, 21, 22, 24, 25, 26, 27, 28, 29),)

原子インデックスのタプルのタプルとして得られました。

原子インデックスを描画してみます。

RDKitのメーリングリストのディスカッションによると、オプションを「Draw.DrawingOptions.includeAtomNumbers=True」と設定することで、表示できるようです。少し書き方を変えて以下のようにしました。

Draw.MolToImage(test_mol, includeAtomNumbers = True)

f:id:magattaca:20190203192126p:plain

構造式の上の方に、認識された原子インデックスがあるのがわかります。

こちらをハイライトしてみます。

#タプルのタプルなのでインデックス[0]が必要
Draw.MolToImage(test_mol, highlightAtoms = matched_atoms[0])

f:id:magattaca:20190203192244p:plain

順番に眺めます。

from ipywidgets import interact,fixed,IntSlider
import ipywidgets

def biphenyl_viewer(idx):
    mol = biphenyl_compounds[idx]
    matched_atoms = mol.GetSubstructMatches(biphenyl)
    return(Draw.MolToImage(mol, highlightAtoms = matched_atoms[0]))

interact(biphenyl_viewer, idx=ipywidgets.IntSlider(min=0,max=len(biphenyl_compounds)-1, step=1));

f:id:magattaca:20190203192518g:plain

ビフェニル構造は末端近くにあるものが多いようですが、いくつかは分子の中央部分にもあります。共通構造といっても分子の中で占める位置が異なるものもあるようです。
中央部分にビフェニルを含むものは、ビフェニルを中心とした対称構造のようにもみえるので、結合様式等（？）の知見を反映してデザインされた化合物なのかもしれません。

ライブラリ化合物で部分構造による絞り込み

それでは、ライブラリ化合物の中からビフェニル構造を持つ化合物を取り出してみたいと思います。

まずは指標を計算済みのEnamine_Premium_collectionについて、指標でのフィルタリングを実施します。

EPc_suppl = Chem.ForwardSDMolSupplier('./Enamine_Premium_collection.sdf')

EPc_mols_filt = []

for x in EPc_suppl:
    if x is not None:
        mol = x
    MW_value = mol.GetDoubleProp('MW')
    MolLogP_value = mol.GetDoubleProp('MolLogP')
    NumHDonors_value = mol.GetDoubleProp('NumHDonors')
    NumHAcceptors_value = mol.GetDoubleProp('NumHAcceptors')
    NumRotatableBonds_value = mol.GetDoubleProp('NumRotatableBonds')
    TPSA_value = mol.GetDoubleProp('TPSA')
    
    #フィルタリング基準を満たす場合のみリストに追加する
    if MW_value > 300 and MW_value <= 500 \
    and MolLogP_value > 3 and MolLogP_value <=5\
    and NumHDonors_value <= 5 \
    and NumHAcceptors_value > 3 and NumHAcceptors_value <= 10 \
    and NumRotatableBonds_value > 3 \
    and TPSA_value > 60 :
        filtered_mols.append(mol)

フィルタリングを満たした化合物のリストを、活性化合物群に対して行なったのと同様にビフェニルの有無で別々のリストにします。

biphenyl = Chem.MolFromSmiles('c1ccc(cc1)c1ccccc1')
biphenyl_compounds = []
w_o_biphenyl = []
<200b>
for mol in filtered_mols:
    if mol.HasSubstructMatch(biphenyl):
        biphenyl_compounds.append(mol)
    else:
        w_o_biphenyl.append(mol)

print(len(biphenyl_compounds))  
#12
print(len(w_o_biphenyl))
#4048

Enamine_Premium_collectionで条件を満たすものは12個となりました。構造を眺めてみます。

Draw.MolsToGridImage(biphenyl_compounds, molsPerRow=4, subImgSize=(200,200))

f:id:magattaca:20190203193038p:plain

なかなかPremium感のあるいい感じの化合物が残っている気がします。

他のライブラリに適用した結果とあわせてまとめます。

	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%)

全体を合わせて 4220 個の分子が残りました。まだまだ多いですが、かなり数を絞り込めてきているように思います。

部分構造検索をKNIMEで

Twitter話題のKNIME・・・、@sumatat_ さんのブログ非プログラマーのためのインフォマティクス入門。（仮）のコンテンツが素晴らしいのでKNIMEブームに便乗することにしました。

（ブログタイトルがまさにこれを求めていた！！って感じです）

今回行った流れは、

活性化合物群をビフェニルの有無で分ける
フィンガープリントを計算し、PCAで次元圧縮したうえでプロット
インタラクティブに遊ぶ

です。

以下の記事を参考にさせていただきました

@Mochimasa さんの記事化合物をベクトルにして比較しプロットする
@sumatat_ さんの記事 KNIMEで化合物をクラスタリング＆可視化してみよう
@iwatobipen先生の記事 Make interactive plot with Knime
@PKさんの記事【KNIME】CellProfilerの画像解析結果を可視化する

早速ワークフロー全体です。２通りの可視化を行うため、プロットのためののメタノード（灰色）が２つあります。

f:id:magattaca:20190203193203p:plain

KNIMEの視認性の高さはすごいですね！何となくやっていることがわかります。 Jupyter NoteBookは３日経つと自分の操作の意味すらわからなくなります（私だけ？）が、これなら思い出しやすそうです。

一応ワークフローの説明をすると

SDFの読み込み(SDF Reader)
SDFの構造ををRDKitのMolオブジェクトに変換（RDKit To Molecule）
ビフェニル構造の有無で２つに分割（RDKit Substructure Filter）
分割したテーブルにクラス名のカラムを追加（Constant Value Column）
２つのテーブルを一つに再結合（Concatenate）
フィンガープリント（Morgan/ 1024ビット/ 半径2）を計算（RDKit Fingerprint）
フィンガープリントのビットベクトルをPCAに渡すため整数値のカラムに分割（Expand Bit Vector）
PCA（PCA）
分子の構造を描画のためSVGに変換（Renderer to Image）
クラスに応じた色分け（Color Manager）
プロット

プロットでは、ビフェニルの有無での色分けと、各化合物の出元（特許の出願人（applicant））での色分けの２通りを行っています。前者ではTile View (KNIME 3.7.1 ではCard Viewの名称が変更となったそうです)、後者ではTable Viewを使用しました。

メタノードの中身はこちら

f:id:magattaca:20190203193246p:plain f:id:magattaca:20190203193252p:plain

実行した結果はこのようになります。まずはビフェニルの有無で色分けした場合・・・

f:id:magattaca:20190203193456g:plain

ビフェニルを含まない構造（赤色）の多くは、プロットの右下に集まっています。これらは全てapplicantがAurigeneとなっています。他のapplicantとは系統の異なる化合物を報告しているようです。また、ビフェニルを含む構造（緑色）でも左下には中国医学科学院の化合物群、中央上部にはIncyteの化合物群が集まっているように見えます。興味深いことに、Morganフィンガープリントにもとづくプロットで化合物の出所がうまく分かれているようです。

applicantによる色分けを眺めてみます。

f:id:magattaca:20190203201707g:plain

先のプロットで想定した通り、applicantごとにプロット上でうまく分かれていそうです。

見よう見まねでも何となく使えてしまうとはKNIMEの完成度高い！

オルト位置換の認識に失敗

元々の部分構造フィルタリングの目標は、オルト位に置換基の入ったビフェニルだったので、SMILESのワイルドカード（*）をつかって、以下のコードを実施しました。

ortho_biphenyl = Chem.MolFromSmiles('c1ccc(c(*)c1)c1ccccc1')
Draw.MolToImage(ortho_biphenyl)

f:id:magattaca:20190203200615p:plain

オルト位置換（なんでも良い）が再現できていそうです。以下のテスト分子で試します。

f:id:magattaca:20190203200739p:plain

test_mol.HasSubstructMatch(ortho_biphenyl)
# Flase

・・・ダメでした。置換基あるのに・・・

まとめ

今回、「ビフェニル構造を部分構造としてもつか否か」を基準としてライブラリのさらなる絞り込みを実施しました。その結果、約4000 個にまで化合物数を絞り込むことができました。また、KNIMEをもちいた可視化を実施し、フィンガープリントを持ちいることで化合物群を出所（出願者）に応じてうまく識別できる可能性があることがわかりました。

残念ながら、オルト位に置換基の入ったビフェニル、というフィルタリングはうまくできませんでした。私、応用力無さすぎ・・・