化学構造と生物学的活性における相関関係をStructure Activity Relationship(SAR)またはQuantative SAR(QSAR)と呼びます。一般的には似たような化合物は似たような生物学的活性を示すことが知られており、この相関関係を理解しドラッグデザインに活かすことが創薬研究において大変重要です。
また、このような問題には細胞の生死、毒性の有無といった化合物がどのクラスに入るのかを推定する分類問題と阻害率（%inhibition）といった連続値を推定する回帰問題の2つがあります。
ChEMBLからhERG阻害アッセイの73データを用いてIC50が1uM未満のものをhERG阻害あり、それ以外をhERG阻害なしとラベルします。
まずは必要なライブラリをインポートします。
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix, f1_score
from sklearn.ensemble import RandomForestClassifierChEMBLでダウンロードしたタブ区切りテキストの処理は8章とほぼ同じですが、今回は活性データが欲しいのでSTANDARD_VALUEという列を探して数値を取り出します。この値が1000nM未満であればPOSというラベルを、そうでなければNEGというラベルを振ります。最後にラベルをnumpy arrayにしておきます。
mols = []
labels = []
with open("ch09_compounds.txt") as f:
header = f.readline()
smiles_index = -1
for i, title in enumerate(header.split("\t")):
if title == "CANONICAL_SMILES":
smiles_index = i
elif title == "STANDARD_VALUE":
value_index = i
for l in f:
ls = l.split("\t")
mol = Chem.MolFromSmiles(ls[smiles_index])
mols.append(mol)
val = float(ls[value_index])
if val < 1000:
labels.append("POS")
else:
labels.append("NEG")
labels = np.array(labels)続いてmolオブジェクトをフィンガープリントに変換します。このフィンガープリントからhERG阻害の有無を予測するモデルを作成します。
fps = []
for mol in mols:
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
fps.append(arr)
fps = np.array(fps)データセットを訓練セットテストセットの2つに分けます。テストセットは作成した予測モデルの精度を評価するためにあとで使います。
x_train, x_test, y_train, y_test = train_test_split(fps, labels)予測モデルを作成するにはインスタンスを作成してfitメソッドで訓練させるだけです
rf = RandomForestClassifier()
rf.fit(x_train, y_train)先程分割しておいたテストセットを予測します。
y_pred = rf.predict(x_test)Confusion matrixを作成します。
Confusion matrixとはクラス分類の結果をまとめた表です。クラスを正しく分類できているかをわかりやすく視覚化できTP,TNが多くFP,FNが少ないほどよく分類できているということになります。
| Actual class | |||
|---|---|---|---|
Positive |
Negative |
||
Predicted class |
Positive |
True Positive(TP) |
False Positive(FP) |
Negative |
False Negative(FN) |
True Negative(TN) |
|
confusion_matrix(y_test, y_pred)
#array([[11, 1],[ 5, 2]])| 11 | 1 |
|---|---|
5 |
2 |
F1スコアを見てみましょう。
f1_score(y_test, y_pred, pos_label="POS")
#0.4あまりよくないですね。
|
Note
|
train_test_split関数がランダムに訓練セットとテストセット分割するので、Confution matrix,F1スコアの値は実行するたびに変わります。 |
-
正しいと予測されたもののうち本当に正しいものの割合を適合率と呼ぶ precision=TP/(TP+FP)
-
正しいものが正しいと予測された割合を再現率と呼ぶ recall=TP/(TP+FN)
F1スコアは適合率と再現率の調和平均となり
F1 = 2 * (precision * recall) / (precision + recall)
で計算されます。
回帰モデルは最初に説明したとおり、連続値を予測するモデルとなります。今回はRandomForestの回帰モデルを作成して、その精度をR2で評価します。データは分類問題で使ったhERGのアッセイデータを利用することにしましょう。最初に必要なライブラリをインポートします。
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from math import log10分類問題のときにはラベル化しましたが、今度は連続値を予測したいのでpIC50に変換します。(なぜpIC50にすると都合が良いのかはそのうち補足する)
pIC50s = []
with open("ch09_compounds.txt") as f:
header = f.readline()
for i, title in enumerate(header.split("\t")):
if title == "STANDARD_VALUE":
value_index = i
for l in f:
ls = l.split("\t")
val = float(ls[value_index])
pIC50 = 9 - log10(val)
pIC50s.append(pIC50)
pIC50s = np.array(pIC50s)データセットをトレーニングセットとテストセットの2つに分割します。フィンガープリントは分類モデルのときに作成したものを流用します。
x_train, x_test, y_train, y_test = train_test_split(fps, pIC50s)訓練します。Scikit-learnの場合はこの手順はどの手法でもほぼ同じメソッドでfitしてpredictです。
rf = RandomForestRegressor()
rf.fit(x_train, y_train)予測しましょう。
y_pred = rf.predict(x_test)予測精度をR2で出してみます。
r2_score(y_test, y_pred)
#0.52まずまずといったところでしょうか。
回帰の当てはまりの良さへの1つの評価指標としてよく利用されるもので決定係数とも呼ばれます。
今回紹介した手法は似たような化合物は似たような生物学的活性を示すという仮設に基づいて生成されるモデルです。もしトレーニングセットに似ている化合物が含まれなかった場合の予測精度はどうなるのでしょうか？
当然その場合は予測された値は信頼できませんよね。つまり、予測値にはその予測が確からしいか？という信頼度が常についてまわります。そのようなモデルが信頼できる、または適用できる範囲をapplicability domainと呼びます。これに関しては明治大学金子先生のモデルの適用範囲・モデルの適用領域が詳しいです。
Hugo Kubinyi先生は似ている化合物は果たして似た活性を示すのか？という疑問を、estradiolのOH基をMethoxy基に変換すると活性が消失する例を上げて説明されていました。その当時のジャーナルが見つからなかったのですが大体次のようなものでした。
Estrogen Receptorのリガンドであるestradiolの水酸基をMeO基に置換した右側の化合物は活性を示しません。両者の類似性は77%と非常に高いにも関わらずです。
この活性消失の理由は複合体結晶構造から綺麗に説明されます。水酸基はGLU353と水素結合を形成していますが、これをMeOに置換すると水素結合能が失われアゴニスト活性がなくなります。
以上が我々が定義した類似性の尺度が蛋白質の分子認識の尺度と一致しないという例です。上記のように必ずしも化合物の類似度から活性が予測できるわけではなく、極めて類似度が高いのに活性が消失してしまうことが多々あります。特にMMPの文脈で説明されたActivity Cliffはこのような事象にそれっぽい名前をつけたものですし、マジックメチルと呼ばれるメチル基導入の効果もコンフォメーションが大きく変化した結果立体的な類似性が減少したという事象に名をつけたような気がします。
applicability domainはそもそも、類似度（距離）という尺度が正しいと仮定をおいて、トレーニングセットの類似性からその予測が信頼できるかという確度を測る手法なのでその辺りは注意したほうがよいでしょう。リード最適化だと、似たような化合物を作るのが普通なのでapplicability domainの問題よりも上記の類似性の尺度が細かいところで狂ってくるほうが問題視されることが多いです。