magattacaのブログ

日付以外誤報

LBVSしようとする話

さて、ファーマコフォアを用いた化合物の絞り込みも終わりました・・・と言いたい所ですが、またしても間違いに気づきました。ファーマコフォアサーチに用いたsdfは、オルト位置換基で絞り込む前の、ビフェニルで絞り込んだ段階のものでした！！！

順番は入れ替わってしまいますが、ファーマコフォアスクリーニング後の化合物群に対して部分構造（オルト位置換基）で絞り込みをかけても結果は変わらない、、、はず。。。ということで、改めて部分構造のフィルターをかけます。

from rdkit import Chem
from rdkit.Chem import AllChem

ortho_biphenyl= Chem.MolFromSmarts('c1ccc(c(*)c1)c1ccccc1')
P_Matched2_suppl = Chem.ForwardSDMolSupplier('./P_Matched2.sdf')
P_Matched2_mols = [mol for mol in P_Matched2_suppl if mol is not None]

print(len(P_Matched2_mols))
#997

#オルト置換ありのリスト
ortho_cpds = []
#オルト置換なしのリスト
w_o_ortho = []

for mol in P_Matched2_mols:
    if mol.HasSubstructMatch(ortho_biphenyl):
        ortho_cpds.append(mol)
    else:
        w_o_ortho.append(mol)

print('Number of ortho substituted compounds', len(ortho_cpds))
# Number of ortho substituted compounds 494
print('Number of ortho unsubstituted compounds', len(w_o_ortho))
# Number of ortho unsubstituted compounds 503

ファーマコフォアスクリーニングで残った化合物（997個）のうち、オルト位に置換基を有するものは約半数の494個でした。こちらを「P_Matched_ortho.sdf」という名前で保存しておきました。

さて、ようやく今回の本題に、、、

これまで既知の活性化合物の情報をもとに

  1. 部分構造（オルト置換ビフェニル）の有無
  2. ファーマコフォア（Aromatic x2、Acceptor x1）を満たすか否か

と、すこしずつ構造を抽象化（？）しながら数を絞り込んできたことになります。次のステップとして、さらに抽象化した基準で、スコアとして評価する、ということを目指したいと思います。

具体的には

  1. フィンガープリントの使用（複数種類）
  2. Tanimoto係数による評価

を行います。複数のフィンガープリントを使用するのは、前回の記事で参照した DeNA 望月氏の資料で単一のものを用いるよりも、統合した方が良い結果が出た、との記載があったからです。

1. フィンガープリントによるスコア化

1-1. スコア化手順

具体的な手順としては以下とします。

  1. 活性化合物群からファーマコフォア基準を満たすものを選ぶ
  2. フィンガープリントを計算
  3. 複数の活性化合物で共通するビットを重要なビットとして抜き出す
  4. 3.で抜き出したビットのベクトルと比較した類似度をスコアとする
  5. 複数のフィンフガープリントにおけるスコアを統合

1-2. step 1. 活性化合物の選別

まずは活性化合物群から、前回スクリーニングに用いたファーマコフォア基準を満たすものを選別します。すべての化合物を用いないのは、すでにライブラリ側の絞り込みをおこなっているため、条件を等しくした方がファーマコフォアに明示的に取り込めなかった重要な情報を反映させることができるのではないか、と考えたからです。

（ちょっと自分でも何を言っているかよく分からないですが、そんな気がしました。）

同じことの繰り返しなので詳細は省きますが、共闘用シェアデータから抜き出した、環状ペプチドを除く低分子 39個のうち、24個の分子で２段階のファーマコフォアスクリーニングの条件を満たしました。

この化合物群に、共結晶構造から抜き出したリガンド6個を加えて、合計30個をフィンガープリントを計算する対象の活性化合物群とします。(fp_molsというリストに格納しました)

1-3. step 2. フィンガープリントを計算

つづいてフィンガープリントを計算します。以下の６つを使います（Morganを2種類使用）。

  1. MACCS keys RDKit documentation / Python API
  2. RDKit RDKit documentation / Python API / Python API
  3. Morgan(Circular) RDKit documentation / Python API
    • ECFP-like (Morgan(radius=2))
    • FCFP-like (Morgan(radius=2, ,useFeatures=True))
  4. AtomPairs RDKit documentation / Python API
  5. Avalon Python API

これ以外に、Topological Torsions （ RDKit documentation / Python API )も候補にあったのですが、ビット数が多すぎて異音、発熱とともにjupyter notebookのカーネルが落ちたので省きました。

計算にあたっては以下の記事、資料を参考とさせていただきました。

maccs_fps = [AllChem.GetMACCSKeysFingerprint(mol) for mol in fp_mols]
rdkit_fps = [AllChem.RDKFingerprint(mol) for mol in fp_mols]
ECFP_fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048) for mol in fp_mols]
FCFP_fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048, useFeatures=True) for mol in fp_mols]
Apair_fps = [AllChem.GetAtomPairFingerprint(mol) for mol in fp_mols]

from rdkit.Avalon import pyAvalonTools
avalon_fps = [pyAvalonTools.GetAvalonFP(mol) for mol in fp_mols]

1-4. step 3. 共通するビットの抜き出し

ついで各フィンガープリントから複数の活性化合物で共通するビットを抜き出します。

まずはフィンガープリントをnumpyのarrayに変換します。(上記、参考記事1 をご参照ください)
MACCS keysを使って、試し実験を行います。

import numpy as np

#関数作成（コピペ）
def fp2arr(fp):
    from rdkit import DataStructs
    arr = np.zeros((1,))
    DataStructs.ConvertToNumpyArray(fp, arr)
    return arr

#MACCS keysについてアレイを作成、行列化する
maccs_fpMtx = np.array([fp2arr(fp) for fp in maccs_fps])

この行列を使って、以前の記事でPLIFを作成した時と同様の可視化を行おうと思います。

import pandas as pd

#行列からDataFrameへ
maccs_df = pd.DataFrame(maccs_fpMtx)

#各列の合計値をもとめ、DataFrameに追加
sums = maccs_df.sum()
maccs_df = maccs_df.append(sums, ignore_index=True)

#Heatmapで可視化
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(20, 10))
sns.heatmap(maccs_df)

f:id:magattaca:20190324164344p:plain

少々見づらいですが、縦軸上から「0 ~ 29」が化合物(30個)、一番下の行「30」が各列の合計値です。横軸はMACCS keyの各ビット(１６７個)です。数に応じて着色されていて、一番下の行をみるとMACCS keyのビットの番号が大きいもので、特に多数の化合物でビットが立っている傾向があることがわかります。

予想以上にわかりにくいので、最後の行を抜き出して棒グラフにします。

plt.figure(figsize=(20, 10))
sums.plot(kind='bar', color='k', alpha=0.7)

f:id:magattaca:20190324165011p:plain

要約統計量を眺めて見ます。

sums.describe()
"""
count    167.000000
mean      10.377246
std       11.409474
min        0.000000
25%        0.000000
50%        3.000000
75%       21.000000
max       30.000000
dtype: float64
"""

各ビットは30個の化合物のうち、平均して10個でビットがたっているようです。10個以上の活性化合物でビットが立っていることを基準として、重要なビット、と捉えることにしてみます。

maccs_common_bit = []

for x in sums:
    if x >= 10:
        maccs_common_bit.append(1)
    else:
        maccs_common_bit.append(0)

print(sum(maccs_common_bit))
#75

MACCS keysフィンガープリントのうち、活性化合物群で重要と考えられるビットのベクトルが手に入りました。75 個 (167 個中)から構成されているようです。

上記処理を関数としますが、いちいちデータフレームにするのは意味がなさそうなので合計値と、ビットのアレイをそれぞれ作ることとします。まずはMACCS keysを例に手順を確認してみます。

maccs_size = len(fp2arr(maccs_fps[0]))
#合計値を格納するnumpyのarrayを用意
sum_arr = np.zeros(maccs_size)

#numpy行列の足し算
for i in maccs_fps:
    arr = fp2arr(i)
    sum_arr += arr

print(sum_arr.mean())
#10.377245508982035

先ほどと同じ平均値を返しました。DataFrameを介さずにすみそうです。
平均値の小数点以下を切り捨てた整数部分を閾値とするため、mathのfloorを使います

import math
maccs_threshold = math.floor(sum_arr.mean())
print(maccs_threshold)
#10

閾値を求めるところまで確認できたので関数とします。フィンガープリントのサイズと、閾値、各ビットで30個の化合物で立っている数、そして平均以上の数立っているビットのベクトルの４つを戻り値とします。

def common_bit(fps):
    fp_size = len(fp2arr(fps[0]))
    sumarr = np.zeros(fp_size)
    for i in fps:
        arr = fp2arr(i)
        sumarr += arr
    mean = sumarr.mean()
    threshold = math.floor(mean)
    
    common_bit = []
    for x in sumarr:
        if x >= threshold:
            common_bit.append(1)
        else:
            common_bit.append(0)
    return fp_size, threshold, sumarr, common_bit

フィンガープリント6種に対して適用します。

fingerprints = [maccs_fps, rdkit_fps, ECFP_fps, FCFP_fps, Apair_fps, avalon_fps]
fp_names = ['maccs_fps', 'rdkit_fps', 'ECFP_fps', 'FCFP_fps', 'Apair_fps',  'avalon_fps']

#各返り値を格納するためのリストを用意
fp_size_list = []
threshold_list = []
sumarr_list = []
common_bit_list = []

#計算実行
for fps in fingerprints:
    sz, th, ar, cb = common_bit(fps)
    fp_size_list.append(sz)
    threshold_list.append(th)
    sumarr_list.append(ar)
    common_bit_list.append(cb)

各フィンガープリントのサイズと閾値を確認してみます。

for name, size, threshold in zip(fp_names, fp_size_list, threshold_list):
    print(name, ':', size, ':', threshold)

"""
maccs_fps : 167 : 10
rdkit_fps : 2048 : 14
ECFP_fps : 2048 : 0
FCFP_fps : 2048 : 0
Apair_fps : 8388608 : 0
avalon_fps : 512 : 10
"""

ECFP、FCFP、Atom_pairでは閾値が0となってしまっています。ビットベクトルのサイズの大きさの割に立っているビットが少なく（スパース？）、平均値がとても小さな値となってしまっているのでしょうか？？？

これでは閾値として不適合なので、フィンガープリント全体ではなく、合計値がゼロではないビットのみについての平均値を代わりに用いたいと思います。

MACCS keysで実験してみます。

#np.nonzero()は戻り値がインデックスであることに注意
non_zero_sum = sum_arr[np.nonzero(sum_arr)]
print(non_zero_sum.mean())
#16.349056603773583

今度は閾値が16という結果となりました。MACCS keyにおいて非ゼロのビットのみを考慮すると、30個の化合物群のうち、平均して半数で同一のビットが立っているということになります。（この解釈で正しいでしょうか？自信がありません。）

上記の新しい閾値をとりいれた関数に書き換えます。

def common_bit2(fps):
    fp_size = len(fp2arr(fps[0]))
    sumarr = np.zeros(fp_size)
    for i in fps:
        arr = fp2arr(i)
        sumarr += arr
    NonZeroArr = sumarr[np.nonzero(sumarr)]
    NonZeroMean = NonZeroArr.mean()
    threshold = math.floor(NonZeroMean)
    
    common_bit = []
    for x in sumarr:
        if x >= threshold:
            common_bit.append(1)
        else:
            common_bit.append(0)
    return fp_size, threshold, sumarr, common_bit

実行の部分は同じなので省略し、結果のみ記載します。

for name, size, threshold in zip(fp_names, fp_size_list2, threshold_list2):
    print(name, ':', size, ':', threshold)
"""
maccs_fps : 167 : 16
rdkit_fps : 2048 : 14
ECFP_fps : 2048 : 4
FCFP_fps : 2048 : 5
Apair_fps : 8388608 : 14
avalon_fps : 512 : 11
"""

全て閾値が0ではなくなりました！スパース性を解決した！
レイドバトルは非ゼロサムゲーム！（知らんけど）

くどいですが再度MACCS keysの値を確認します。

print('MACCS NonZero: ', len(non_zero_sum))
#MACCS NonZero:  106
print('MACCS base bits: ', sum(common_bit_list2[0]))
#MACCS base bits:  61

MACCS keyに関しては「非ゼロのビットが106個で、このうち61個で閾値よりを満たす」、となりました。先の棒グラフを見てもそこまでおかしな値とはなっていなさそうです。

取り急ぎ、手に入れたビットベクトルを基準としてもちいることとし、辞書に格納しておきます。

base_BitVect_dict={}

for k, v in zip(fp_names, common_bit_list2):
    base_BitVect_dict[k] = v

1-5. step 4. 類似度の指標を作成

基準となるビットベクトルが用意できたので、このベクトルに対してスクリーニング対象の化合物のビットベクトルの類似度を計算したいと思います。

タニモト係数を使いたいのですが、numpyのarrayとしており、RDKitのオブジェクトではなくなってしまったのでDataStructs.TanimotoSimilarity()を使えるかわかりません。

Tanimoto類似度の説目を見る限り、ビットが立っている要素について積集合を和集合で割れば良いということのようですから、比較対象のビットベクトル同士を足し合わせて、「2の要素数」を「1 or 2の要素数」で割ってしまえば良さそうです。

まずはファーマコフォアで絞り込んだスクリーニング対象の化合物（オルト位置換も考慮済みのもの）をSDFから読み込み、scr_molsというリストに格納しました。

一番最初の分子をテスト化合物として取り出し、各フィンガープリントを計算したのち、アレイに変換し、辞書型として格納します。

test_mol = scr_mols[0]

#関数を作成
def array_dict(mol):
    maccsFP_arr = fp2arr(AllChem.GetMACCSKeysFingerprint(mol))
    rdkitFP_arr = fp2arr(AllChem.RDKFingerprint(mol))
    ECFP_arr = fp2arr(AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048))
    FCFP_arr = fp2arr(AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048, useFeatures=True))
    Apair_arr = fp2arr(AllChem.GetAtomPairFingerprint(mol))    
    avalon_arr = fp2arr(pyAvalonTools.GetAvalonFP(mol))
    
    tmp_dict = {}
    fp_names = ['maccs_fps', 'rdkit_fps', 'ECFP_fps', 'FCFP_fps', 'Apair_fps',  'avalon_fps']
    fp_arrs = [maccsFP_arr, rdkitFP_arr, ECFP_arr, FCFP_arr, Apair_arr, avalon_arr]
    for k, v in zip(fp_names, fp_arrs):
        tmp_dict[k] = v
    
    return tmp_dict

#テスト分子で実行
test_array_dict = array_dict(test_mol)

MACCS key に関して、テスト分子と、先に作成した基準ビットベクトルでTanimoto係数を計算してみます。

#それぞれのビットベクトル（array）を読み込む
base_maccsAR = base_BitVect_dict['maccs_fps']
test_maccsAR = test_array_dict['maccs_fps']

#array同士を足し算
maccsAR_sum = base_maccsAR + test_maccsAR

#各値の数を確認
test_count_1 = np.sum(maccsAR_sum == 1)
print('1 :', test_count_1)
# 1 : 32
test_count_2 = np.sum(maccsAR_sum == 2)
print('2 :', test_count_2)
# 2 : 36

#Tanimoto類似度を計算
test_Tanimoto = test_count_2 / (test_count_1 + test_count_2)
print('Tamimoto :', test_Tanimoto)
#Tamimoto : 0.5294117647058824

うまくいったようです。関数化します。

def TaniSimCal(arr1, arr2):
    arr_sum = arr1 + arr2
    count1 = np.sum(arr_sum == 1)
    count2 = np.sum(arr_sum == 2)
    Tanimoto = count2 / (count1 + count2)
    return Tanimoto

スクリーニング対象の化合物について、6種類のフィンガープリントのそれぞれのTanimoto類似度を計算します。

Scores = [MACCS_Scores, RDkit_Scores, ECFP_Scores, FCFP_Scores, Apair_Scores, Avalon_Scores]

MACCS_Scores = []
RDkit_Scores = []
ECFP_Scores = []
FCFP_Scores = []
Apair_Scores = []
Avalon_Scores = []

for mol in scr_mols:
    #スクリーニング対象のFPのアレイを作成（辞書）
    fp_arr_dict = array_dict(mol)
    
    #基準の各ビットベクトルに対するTanimoto類似度を計算
    for fp, sc in zip(fp_names, Scores):
        base_arr = base_BitVect_dict[fp]
        scr_arr = fp_arr_dict[fp]
        
        TaniScore = TaniSimCal(scr_arr, base_arr)
        
        #各フィンガープリントのスコアリストに追加
        sc.append(TaniScore)

スクリーニング対象４９４化合物の、６種類のフィンガープリントについて、活性化合物群から導いた基準ビットベクトルに対する類似度をプロットしてみました。

for i in range(6):
    plt.subplot(2, 3, i+1)
    plt.hist(Scores[i], bins=20, range=(0, 1))
    plt.title(fp_names[i])

f:id:magattaca:20190324173142p:plain

1-6. step 5. 統合したスコアの作成

上記で求めた値を統合して、各化合物の統合スコアとしたいのですが、このままではフィンガープリントごとにスコアの分布のばらつきが大きいです。Tanimoto類似度にこのような処理をするのが適切かわかりませんが、各フィンガープリントにおいてスコアを標準化（平均 0, 分散 1）したのち、和をとった値を統合したスコアとしたいと思います*4。

import statistics

Scores_STD = []
for l in Scores:
    l_mean = statistics.mean(l)
    l_stdev = statistics.stdev(l)
    s = [(i-l_mean) / l_stdev for i in l]
    Scores_STD.append(s)

標準化したスコアを先と同様にプロットします。

for i in range(6):
    plt.subplot(2, 3, i+1)
    plt.hist(Scores_STD[i], bins=20)
    plt.title(fp_names[i])

f:id:magattaca:20190324173742p:plain

標準化が完了したので、各スコアを足しあわせて統合スコアとします。なんとなく0以上の値にしたいので最小値分シフトさせます。

Integrated_SC_arr = np.zeros(494)
for i in Scores_STD:
    Integrated_SC_arr += np.array(i)

#最小値を求める
int_min = min(Integrated_SC_arr)
#最小値の絶対値分シフト
Integrated_SC = [i + abs(int_min) for i in Integrated_SC_arr]

#グラフ化
plt.hist(Integrated_SC, bins=20)

f:id:magattaca:20190324174003p:plain

スコア化が完了したので、SetDoublePropを使って、Molオブジェクトにプロパティ（Similarity_Score）として持たせ、SDFで出力しました（SCR_compounds_SimScore.sdf）。

2. 上位化合物の検証

スコアの付与とデータの書き出しが完了しました。今回設定したスコアで、上位となった化合物はどのようなものか、確認して見たいと思います。PandasToolsでDataFrameに読み込みます。

from rdkit.Chem import PandasTools

SIM_df = PandasTools.LoadSDF('./SCR_compounds_SimScore.sdf')

スコアの順に並べ替えたいので、typeをstrからfloatに変換しておきます。

SIM_df['Similarity_Score'] = SIM_df['Similarity_Score'].astype(float)

新しく追加したスコアをつかって降順にソートします。

SIM_df_st = SIM_df.sort_values('Similarity_Score', ascending=False)

スコアの良いものから１０個取り出してスコアととともに描画します。

from rdkit.Chem import Draw

#Molオブジェクトの取り出し
Top_10 = SIM_df_st.iloc[:10, 8]
#スコアの取り出しとstrへの変換
Top_10_scores = SIM_df_st['Similarity_Score'][:10]
legends = [str(i) for i in Top_10_scores]

Draw.MolsToGridImage(Top_10, legends = legends)

f:id:magattaca:20190324174512p:plain

かなりそれらしい化合物が選択できてきているのではないでしょうか？

少し想定外だったのは、オルト位に置換基の入ったビフェニル構造でスクリーニングした際に、置換基の大きさ（原子数？）をしていなかったためか、オルト位からそのまま分子が伸長したような化合物が上位にヒットしてきているということです。

例えばこれまで何度かみてきた下の構造のように、オルト位置換基をもちつつ、他の位置（メタ位など）からさらに分子が伸長している、というのがぼんやりとした理想だったのですが、、、

f:id:magattaca:20190324174635p:plain

3. まとめ

今回は、リガンドをベースにしたスクリーニングを行うにあたって、これまでよりも抽象化し、フィンガープリントを使って見ました。また、使用にあたって

  • 複数の活性化合物の情報の取り組み
  • 複数の種類のフィンガープリントの組み合わせ

というのを試みて見ました。
上位ランク化合物は、なんとなく頭の中で「こういう化合物がひっかかったらいいな」と思っていた構造とはことなりました。適切な情報にまで落とし込むことと、それを表現することがまだまでできていないようです。

ひとまずLBVSっぽいということにしたい！　

ところで、創薬レイドバトルの候補化合物絞り込みトップ500と、トップ10の提出みたいなのですが、すでに500個ない。 どうしよう・・・

それっぽいことをやって見たいというだけで、スコアの標準化や統合（アンサンブル？）などかなり適当なことをしています。間違っている点等あればご指摘いただければ幸いです。

*1:化合物をベクトルにして比較しプロットする https://qiita.com/Mochimasa/items/f1b60246ece7da46f6a9

*2: AI創薬のためのケモインフォマティクス入門 https://github.com/Mishima-syk/py4chemoinformatics/blob/

*3: 化学の新しいカタチ「RDKitでフィンガープリントを使った分子類似性の判定」https://future-chem.com/rdkit-fingerprint/

*4: Pythonで正規化・標準化(リスト、NumPy配列、pandas.DataFrame) https://note.nkmk.me/python-list-ndarray-dataframe-normalize-standardize/