この記事では化合物をベクトルに変換する方法について説明し、このベクトルを利用して (1)化合物間の類似度を計算する方法と (2)散布図に化合物をプロットする方法を紹介します。なお、この記事は創薬 Advent Calendar 2017 (#souyakuAC2017) の5日目の記事です。

化合物空間とフィンガープリント

ケモインフォマティクスの論文を呼んでいると、しばしば化合物空間 (chemical space) という言葉が出てきます。これは、世の中に存在しうるすべての化合物が配置された仮想的な空間のことです。概念的には原子のつなぎ方次第で無限の化合物が存在し得ますから、化合物空間は途方もない広さをもっていると言えます。創薬においては、この空間の中から医薬品になる化合物を見つけ出すことが求められます¹。

ここまでは非常に概念的な話でしたが、実際に作業をするにあたっては化合物空間を数学的な空間に落としこむ必要があります。より具体的に言えば、各化合物を一定のルールで固定長のベクトルに変換する必要があります。このベクトルはフィンガープリント (fingerprint)と呼ばれ、さまざまなフィンガープリントの生成方法が提案されています。

最後に紹介するコードで実際に使用しているのはMACCS Keysと呼ばれるフィンガープリント(生成手法)です。MACCS Keysは、166種類の部分構造からなる辞書を用い、各構造が計算対象の化合物に含まれているか否かを0/1で表しただけのシンプルなフィンガープリントです。従って、このフィンガープリントは166次元のバイナリのベクトルとなり、各次元が特定の部分構造と1対1対応することになります。

化合物間の類似度を計算する

化合物の類似度を評価するには、それぞれの化合物についてフィンガープリントを計算し、そのフィンガープリント間の類似度を計算すればいいことになります。化合物の類似度の指標としてよく使われるのはTanimoto係数 (Tanimoto coefficient) です。フィンガープリントAとBの間のTanimoto係数は下記のように計算されます。:

T(A, B) = \frac {\sum_{i}(A_i \wedge B_i)} {\sum_{i}(A_i \vee B_i)}

早い話が、AとBで共通に立っているビット (ANDをとったときの立っているビット) の数をAかBいずれかあるいは両方で立っているビット (ORをとったときの立っているビット)の数で割っただけの話です。従って、Tanimoto係数は0から1の値をとり、前者が共通する部分構造が1件もない場合なく全く似ていない場合で、後者がフィンガープリントが完全一致して最も類似度が高い場合に対応します。

実際にRDKitを使って、化合物間の類似度を計算してみましょう。

compare_compounds.py

from rdkit import Chem, DataStructs
from rdkit.Chem import MACCSkeys

smilesDict = {
        "salicylic acid": "c1ccc(c(c1)C(=O)O)O",
        "acetylsalicylic acid": "O=C(Oc1ccccc1C(=O)O)C",
        "acetaminophen": "CC(=O)NC1=CC=C(C=C1)O"}

molDict = {name: Chem.MolFromSmiles(smiles) for name, smiles in smilesDict.items
()}
fpDict = {name: MACCSkeys.GenMACCSKeys(mol) for name, mol in molDict.items()}

for nameA, fpA in fpDict.items():
    for nameB, fpB in fpDict.items():
        tanimoto = DataStructs.FingerprintSimilarity(fpA, fpB)
        print("%s - %s: %0.2f" % (nameA, nameB, tanimoto))

compare_compounds.pyの実行結果

salicylic acid - salicylic acid: 1.00
salicylic acid - acetaminophen: 0.45
salicylic acid - acetylsalicylic acid: 0.74
acetaminophen - salicylic acid: 0.45
acetaminophen - acetaminophen: 1.00
acetaminophen - acetylsalicylic acid: 0.42
acetylsalicylic acid - salicylic acid: 0.74
acetylsalicylic acid - acetaminophen: 0.42
acetylsalicylic acid - acetylsalicylic acid: 1.00

FingerprintSimilarityty関数が、化合物間のTanimoto係数を計算している部分です。引数で指定すれば別の種類の類似度を計算することもできます。実行結果においては、サリチル酸とアセチルサリチル酸、アセトフェノンの3化合物の間でのTanimoto係数が出力されています。同一の化合物間では1.00、サリチル酸とその誘導体であるアセチルサリチル酸の間では比較的高い0.74という値になっています。一方、人間の目で構造式を見ても明らかに似ていないアセトアミノフェンとサリチル酸、アセトアミノフェンとアセチルサリチル酸の間では、0.45と0.42という比較的低い類似度が出力されており、直感に即した計算結果になったことが確認できました。

PCAによる化合物空間の可視化

PCAはprincipal component analysisの略で、日本語では主成分分析と呼ばれています。PCA自体は、ケモインフォマティクスに限らず、数理解析の分野で古くから広く使われている手法です。ここでPCAを行う目的は、166次元という高次元のフィンガープリントを無理やり2次元に落とし込んで可視化することにあります。下記がRDKitで生成したMACCS Keysを元に、scikit-learnの機能を使ってPCAを実行し、matplotlibで散布図を描くプログラムです。なお、解析対象の化合物として3日目の記事でPubChemから取得したyes1_inhibition.sdfを使っています。

pca.py

import numpy as np
from matplotlib import pyplot as plt


def main():
    from rdkit import Chem
    from rdkit.Chem import MACCSkeys
    spl = Chem.SDMolSupplier('yes1_inhibition.sdf')

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

    fps = [MACCSkeys.GenMACCSKeys(mol) for mol in mols]
    fpMtx = np.array([fp2arr(fp) for fp in fps])

    plotPCA(fpMtx)

# convert rdkit fingerprint to numpy array
def fp2arr(fp):
    from rdkit import DataStructs
    arr = np.zeros((1,))
    DataStructs.ConvertToNumpyArray(fp, arr)
    return arr


# plot each compound using PCA
def plotPCA(fpMtx):
    from sklearn.decomposition import PCA

    pca = PCA(n_components=2)
    res = pca.fit_transform(fpMtx)
    # extract each component
    pc = res.T

    cum_cr = sum(pca.explained_variance_ratio_)
    print("cumulative contribution ratio=%.2f" % cum_cr)

    plt.figure()
    plt.scatter(pc[0], pc[1], marker=".")
    plt.xlabel("PC1")
    plt.ylabel("PC2")
    plt.legend()
    plt.savefig("pca_plot.png")

if __name__ == "__main__":
    main()

pca.pyの実行結果

cumulative contribution ratio=0.20

標準出力に表示されているのが、累積寄与率と呼ばれる数値です。これは、簡単に言うと、次元を圧縮した結果の第1主成分と第2主成分の2軸で、元々のデータをどの程度説明できるかを示した値です。このプロットで2割の情報しか説明できていないということになりますから、飽くまで参考程度のものということになります。まだまだ工夫する余地はありそうですね。とはいえ、化合物をフィンガープリントというベクトルに落とし込み初歩的な解析を行うことで、概念的なものに過ぎなかった化合物空間の"一端"を可視化にすることができました。

もっとも概念的に存在するすべての構造が合成可能というわけではありませんし、医薬品を想定した場合には分子量等の制限もありますから、文字通りに全化合物空間を探索するということはありません。 ↩