Lipinski先生に従うことにした話
以前の記事 でライブラリの指標の計算を行いました。その結果、さらに前の記事 で活性化合物群から求めた指標の閾値では、全く化合物の数が絞れないことがわかりました。
やはり、素人が付け焼き刃で基準値を考えるのは無理がある・・・ここは偉大な先人の後についていくしかない・・・ということでまずは、以下の２つのフィルタリングで化合物数が絞り込めるか試してみたいと思います。
- LipinskiのRule of 5を満たすもの
- フラグメント指標 Rule of Threeを元にしたフラグメントの削除
1.のLipinskiのRule of 5 に加えて、2.の基準を加えた理由は、
「今回の目標はフラグメントレベルではなくある程度の活性を達成できるような化合物を見出すことなので、むしろフラグメント以下の分子は除いてしまった方が良い」
ということです。
LipinskiのRule of 5の適用
Lipinskiの法則に関しては「化学の新しいカタチ」さんのこちらの記事（ RDKitにおける記述子の扱い方をリピンスキーの法則を通して学ぶ ）を参考にさせていただきました。
| 指標 | 分子量 | LogP | 水素結合供与体数 | 水素結合受容体数 |
|---|---|---|---|---|
| Rule of 5 | ≦500 | ≦5 | ≦5 | ≦10 |
指標の計算値のみを集めたcsvファイルを出力しておいたのでこちらを利用して、Rule of 5を満たす化合物の数を取得します。
まずはPandasのDataFrameで読み込みます。
df_EAc = pd.read_csv('./Enamine_Advanced_collection_desc.csv')
Lipinskiの法則を判断する関数を作成し、DataFrameにapplyで適用しようとしたのですが、どうしてもうまくできませんでした。また、mapでは遅いという記事を見かけたので、こちら（pandasで複数カラムを参照して高速に1行1行値を調整する際のメモ）を参考に処理を行いました。
PandasのDataFrameにapplyで関数を適用すると、Seriesとして１行ずつ処理することになるようですが、こちらのSeriesへのアクセスが遅いのが問題とのことです。解決策として、Pythonの辞書へのアクセスが高速であることを利用すれば良い、とのことでした。
# Lipinskiの判別に用いる指標を辞書として取り出す。(to_dict関数) MW_dict = df_EAc['MW'].to_dict() MolLogP_dict = df_EAc['MolLogP'].to_dict() NumHDonors_dict = df_EAc['NumHDonors'].to_dict() NumHAcceptors_dict = df_EAc['NumHAcceptors'].to_dict() # 辞書のKeyをDataFrameのindexの値としているので、indexの値を格納するカラムをDataFrameに追加する。 df_EAc['index_val'] = df_EAc.index # Lipinskiのルールを判別する関数を作成 def lipinski(index): # indexの分子の指標を取り出す MW = MW_dict[index] MolLogP = MolLogP_dict[index] NumHDonors = NumHDonors_dict[index] NumHAcceptors = NumHAcceptors_dict[index] # Lipinskiのルールに合致するならTrue、しないならFalseとする if MW <= 500 \ and MolLogP <=5\ and NumHDonors <=5 \ and NumHAcceptors <=10: return True else: return False # 上記関数を適用し、新しいカラム(Lipinski)に当てはめる df_EAc['Lipinski'] = df_EAc['index_val'].apply(lipinski)
以上で、Lipinskiのルールの判別が完了しました。あとはTrueの数を数えれば合致した分子の数がわかります。
Trueは1、Falseは0なので、そのままsum()を実行することで条件を満たす要素の数が得られるそうです。（参考記事）
EAc_Lipinski_True = df_EAc['Lipinski'].sum() print(EAc_Lipinski_True)
483858となりました。Enamine_Advanced_collectionは元々が486322だったので、2千個程度Lipinskiの法則を満たさないものが含まれているようです。
他のライブラリの計算結果と合わせると下記の通りです。
| Lipinski | Enamine_Premium (分子量300以上) |
Enamine_Advenced | Enamine_HTS | UOS_HTS |
|---|---|---|---|---|
| 総数 | 109602 | 486322 | 1921489 | 516664 |
| True | 109590 | 483858 | 1843248 | 437821 |
| False | 12 | 2464 | 78241 | 78843 |
Enamine_Premium_collectionのみFalseとなった数が12と少ないのが興味深い結果です。Premiumだけに何らかのプレミアムな基準で選ばれた優良な化合物たちなのでしょうか？？
Enamine_Premium_collectionだけでも11万個の化合物があるので、こちらのみから絞り込むのでも良い気がしてきました・・・
フラグメントライクな化合物
Lipinskiの基準は、ある閾値以下の化合物を選抜するものです。前回参照したSAR News Np.19の記事には、フラグメント指標「Rule of Three」が記載されていました。
| 指標 | 分子量 | LogP | 水素結合供与体数 | 水素結合受容体数 | 回転可能結合数 | 極性表面積 |
|---|---|---|---|---|---|---|
| 基準 | ≦300 | ≦3 | ≦3 | ≦3 | ≦3 | ≦60 |
こちらを頼りにしてフラグメントライクな分子を取り除こうと思います。
以下に、一つずつの指標でフィルタリングした場合と、全てを満たす場合の数をまとめました。
| 指標 | 総分子数 | 分子量 | LogP | 水素結合供与体数 | 水素結合受容体数 | 回転可能結合数 | 極性表面積 | すべて満たすもの |
|---|---|---|---|---|---|---|---|---|
| 基準 | >300 | >3 | >3 | >3 | >3 | >60 | ||
| Enamine_Premium (分子量300以上) |
109602 | 109602 (100%) |
13327 (12%) |
49 (0.04%) |
99362 (91%) |
78633 (72%) |
84269 (77%) |
0 |
| Enamine_Advenced | 486322 | 305230 (63%) |
129679 (27%) |
902 (0.2%) |
325234 (67%) |
305416 (63%) |
277665 (57%) |
64 (0.01%) |
| Enamine_HTS | 1921489 | 1609119 (84%) |
904582 (47%) |
6847 (0.4%) |
1492627 (78%) |
1547785 (81%) |
1319649 (69%) |
87 (0.05%) |
| UOS_HTS | 516664 | 439751 (85%) |
284727 (55%) |
2461 (0.5%) |
400416 (78%) |
433231 (84%) |
351715 (68%) |
746 (0.14%) |
6つの指標すべてでフラグメント指標 Rule of Threeよりも大きいものの和集合、という基準（一番右の列）にしてしまうと、ほとんどの分子が除外されてしまいます。
特にプレミアム感のある期待のライブラリEnamine_Premium_collection がすべてなくなってしまいます。これではちょっとやりすぎ感があります。
上記の表を見ると、特に水素結合供与体数（>3）が削減率が高く、ついで LogP（>3） で削られているものも多そうです。
計算方法が間違っているといけないので、以下にコードを転記しておきます。（かなり冗長です）
#一つずつ計算する場合 df_MW_300 = df_EAc['MW'] > 300 df_MolLogP_3 = df_EAc['MolLogP'] > 3 df_NumHD_3 = df_EAc['NumHDonors'] > 3 df_NumHA_3 = df_EAc['NumHAcceptors'] > 3 df_NumRB_3 = df_EAc['NumRotatableBonds'] > 3 df_TPSA_60 = df_EAc['TPSA'] > 60 print(df_MW_300.sum()) print((df_MW_300.sum() / len(df_EAc))*100) print(df_MolLogP_3.sum()) print((df_MolLogP_3.sum() / len(df_EAc))*100) print(df_NumHD_3.sum()) print((df_NumHD_3.sum() / len(df_EAc))*100) print(df_NumHA_3.sum()) print((df_NumHA_3.sum() / len(df_EAc))*100) print(df_NumRB_3.sum()) print((df_NumRB_3.sum() / len(df_EAc))*100) print(df_TPSA_60.sum()) print((df_TPSA_60.sum() / len(df_EAc))*100) # すべての基準の和集合 df_all =((df_EAc['MW'] > 300) & \ (df_EAc['MolLogP'] > 3) & \ (df_EAc['NumHDonors'] > 3) & \ (df_EAc['NumHAcceptors'] > 3) & \ (df_EAc['NumRotatableBonds'] > 3) & \ (df_EAc['TPSA'] > 60)) print(df_all.sum()) print((df_all.sum() / len(df_EAc)) * 100)
PCAを用いた次元圧縮
フラグメントライブラリの指標をすべて逆転して用いてしまうと、化合物数を減らしすぎてしまう・・・でもどの基準を残せば良いかわからない・・・。
指標を選択する基準はないか？？？ということでPCAを用いた次元圧縮を行ってみたいと思います。 （たぶん使い方間違ってる）
まずはすべてのDataframeを統合します。Enamine_Advanced_collectionのMolLogPの計算値には欠損値(NaN)があるのでこれは除いておきます。
# NaNを含む列の削除 df_EAc_w_o_NaN = df_EAc.dropna() # 確認（isnullでTrueとなる数をカウントして合計が0ならNaNは無い） print(df_EAc_w_o_NaN.isnull().values.sum() == 0) # True # 4つのライブラリを統合 df_all = pd.concat([df_EPc, df_EAc_w_o_NaN, df_EHc, df_UH]) # 化合物総数 print(len(df_all)) #3034076
全部で約 300万個あります。
これにPCAを行ってみます。
from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA import matplotlib.pyplot as plt sc = StandardScaler() # 指標の値をPandasからnumpyのndarrayに変換する (values) X = df_all[descriptors].values # 標準化 X_std = sc.fit_transform(X) # PCA(2成分) pca = PCA(n_components =2) X_pca = pca.fit_transform(X_std) #可視化 plt.figure() plt.scatter(X_pca[:, 0], X_std[:, 1]) plt.xlabel('PC1') plt.ylabel('PC2') plt.show() #次元ごとの寄与率 print(pca.explained_variance_ratio_)
コピペしてみたが、何がおきているかよくわからない！！
とりあえず第１主成分と、第２主成分をあわせても累積寄与率（？）は７割以下みたいです。
因子負荷量の計算
よくわからないが理解するのを待っていたらいつまでかかるかわからない！突き進むしかない！
各主成分に対してどの指標がどの程度相関しているかを眺めるため、因子負荷量を計算します。
import numpy as np # 因子負荷量を計算 pca_components = pca.components_*np.c_[np.sqrt(pca.explained_variance_)] # 眺めるための新しいDataFrameを作成 column_names = descriptors PCA_df = pd.DataFrame([pca_components[0], pca_components[1]], columns=column_names)
以下のようになりました。
ついでにグラフ化しておきます。（グラフでは絶対値にしています）
labels = ['MW', 'MolLogP', 'NumHD', 'NumHA', 'NumRB', 'TPSA'] # NumPyのfabsを使って絶対値に修正 PCA_df_abs = np.fabs(PCA_df) height_0=PCA_df_abs.loc[0] height_1=PCA_df_abs.loc[1] left= np.arange(len(height_1)) width=0.3 plt.bar(left, height_0, color='r', width =width) plt.bar(left+width, height_1, color='b', width =width) plt.xticks(left, labels) plt.show()
第１主成分（赤色）は、分子量(MW)、水素結合受容体数(NumHA)、回転可能結合数（NumRB）、極性表面積（TPSA）との相関が、第２主成分（青色）はLogP(MolLogP)との相関が強そうです。
水素結合供与体数(NumHD)は第１主成分よりも第２主成分との相関が強そうですが、MolLogPと比べると見劣りがします。
もう、NumHDを削ってしまっても良さそうな気がしてきました。（・・・無理やりこじつけた）
# NumHDonor以外の基準の和集合 df_all_w_o_NumHD =((df_all['MW'] > 300) & \ (df_all['MolLogP'] > 3) & \ (df_all['NumHAcceptors'] > 3) & \ (df_all['NumRotatableBonds'] > 3) & \ (df_all['TPSA'] > 60)) print(df_all_w_o_NumHD.sum()) print((df_all_w_o_NumHD.sum() / len(df_all)) * 100)
フラグメント指標から水素結合供与体数を除いた残りの５つを使って、５つ全てを満たす化合物数を求めると645704個（21%）となりました。
これならまだマシそうなので、こちらを用いたいと思います。
二つの基準を満たす分子の数
以上、見てきた内容をまとめると適用する基準は下記となります。
| 指標 | 分子量 | LogP | 水素結合供与体数 | 水素結合受容体数 | 回転可能結合数 | 極性表面積 |
|---|---|---|---|---|---|---|
| Lipinskiより | ≦500 | ≦5 | ≦5 | ≦10 | ||
| フラグメント指標より | >300 | >3 | >3 | >3 | >60 |
２つの基準を満たす分子はどの程度あるでしょうか？
df_all_both = ((df_all['MW'] > 300) & (df_all['MW'] <= 500) & \ (df_all['MolLogP'] > 3) & (df_all['MolLogP'] <= 5) &\ (df_all['NumHDonors'] <= 5) & \ (df_all['NumHAcceptors'] > 3) & (df_all['NumHAcceptors'] <= 10) &\ # 5以下になっていたので10以下に修正（01/28） (df_all['NumRotatableBonds'] > 3) & \ (df_all['TPSA'] > 60)) print(len(df_all)) # 3034076 print(df_all_both.sum()) # 563000 print((df_all_both.sum() / len(df_all)) * 100) # 18.55589642447981
元々の分子の総数 3034076、両基準を満たす分子の数563000（19%）となりました。
(最初の投稿では NumHA 10以下ではなく間違って5以下としていたため、このさらに半分の約30万個にまで減っていました。（01/28修正）)
両基準で絞り込んだ化合物数は約56万個で、 もともと300万個の化合物があったので、1/5程度にまで減らすことができたことになります。
かなり無茶苦茶な話をしているので怒られてしまいそうですが、とりあえずの結果としてこの基準をもとに進めていきたいと思います。
解釈や使用方法等に誤りがたくさんあると思うのでご指摘いただければ幸いです。
