化学構造生成器 Chemical VAE の使い方 - ご注文はリード化合物ですか？〜医薬化学録にわ〜

この記事は創薬 (dry) Advent Calendar 2020 の 12 日目の記事です。
一週間以上遅刻してすみません・・・（というか早く wet の方も書かねば）。
adventar.org

ケモインフォマティクスで盛んに研究されているものの一つに、構造生成器というものがあります。
構造生成器とは、任意の条件を入力すると、構造式を自動で出力するアルゴリズムのことです。
生物活性値などの情報を入力すると、その活性値を持つ医薬品の構造式が出てくるのが理想的ですが、実際には、scaffold を入力して、それを構造展開した構造が出力される、というのが現実的な使い方になると思います。
創薬化学者が中々デザインできない構造を自動で出力することを目標に、様々な構造生成器が（比喩ではなく）毎週のように報告されています。

構造生成器は、大まかに、building block タイプと生成モデルタイプに大別されます。
Building block タイプは、部分構造のパーツを予め用意しておき、レゴブロックのように組み合わせることで、新しい構造を生成します。化学構造としては妥当なものを生成しやすいことが長所として考えられます。
一方、生成モデルタイプは、変文オートエンコーダ（variational autoencoder: VAE）や敵対的生成ネットワーク（generative adversarial network: GAN）など、深層学習を用いて構造を生成します。学習データの構造にはない、新規な構造を提案することが期待されています。

生成モデルタイプの構造生成器の先駆けとなったのが、 ~~Twitter おじさん~~Aspuru-Guzik 教授らのチームが作成した、Chemical VAE です。
https://pubs.acs.org/doi/10.1021/acscentsci.7b00572

Chemical VAE は、SMILES を用いて作成した VAE で、潜在空間座標を decoder で復元することで、新たな化合物を提案することが可能となります。
Chemical VAE のコードは、Guzik 先生の GitHub で公開されているので、モデルをすぐに試すことができます。
github.com

実行環境も enviromment.yml が提供されているのですが、これを展開しても私のパソコンではコードが動かず・・・
仕方がないので、自分で環境を以下のように作成しました。

conda create -n chemvae python=3.7 anaconda
source activate chemvae
conda install -c rdkit rdkit
pip install tensorflow==1.13.1
pip install keras==2.0.5
git clone https://github.com/aspuru-guzik-group/chemical_vae.git

chemvae は作成した環境名です。Tensorflow のバージョンは GitHub に書かれているバージョンと違いますが、これで動きました（GitHub に書かれているバージョンが古過ぎる）。

git コマンドでコードをダウンロードすると、chemical_vae というフォルダが新しく出来ていると思います。
examples フォルダの中に、intro_to_chemvae.ipynb というファイルがあります。
このコードが、Chemical VAE のintroduction の役割を果たしています。
早速開いて実行したいのですが、そのまま実行すると、chemvae がない、と怒られてしまいます。
そこで、パスを通すため、最初のセルの一番上に、2 行書き込んであげます。

import sys
sys.path.append("../")

後のセルは特に追記することなく、そのまま動かすことができます。
ただし、最後の t-SNE は非常に遅いので、一定のマシンスペックがある方以外は実行しないことをオススメします（t-SNE 自体が計算コストの重い手法）。

このコードをそのまま動かしてもいいのですが、少し自分で構造を生成した感じがしない（分かりにくい）ので、 example フォルダ内に、自分で構造生成した感が出るコードを作成しました。

# 必要なライブラリを読み込む
import sys
sys.path.append("../")

# tensorflow backend
from os import environ
environ['KERAS_BACKEND'] = 'tensorflow'
# vae stuff
from chemvae.vae_utils import VAEUtils
from chemvae import mol_utils as mu
# import scientific py
import numpy as np
import pandas as pd
# rdkit stuff
from rdkit.Chem import AllChem as Chem
from rdkit.Chem import PandasTools
# plotting stuff
import matplotlib.pyplot as plt
import matplotlib as mpl
from IPython.display import SVG, display
%config InlineBackend.figure_format = 'retina'
%matplotlib inline

# VAE モデルを読み込む
vae = VAEUtils(directory='../models/zinc_properties')

# 初期条件の SMILES を入力する
input_smiles = "c1ccccc1"

# 入力構造を one-hot vector に変換する
input_tensor = vae.smiles_to_hot(input_smiles, canonize_smiles=True)

# 潜在空間座標に変換
latent_coordinate = vae.encode(input_tensor)

# 潜在空間座標をランダムにサンプリング
np.random.seed(0)
generated_coordinates = np.random.multivariate_normal(latent_coordinate[0], np.eye(196), 100)

# 生成した座標を SMILES に復元する
output_tensor = vae.decode(generated_coordinates)
generated_smiles = vae.hot_to_smiles(output_tensor, strip=True)

# 生成した文字列のうち、構造として成立しているものだけを選択
valid_smiles = []
for string in generated_smiles:
    if Chem.MolFromSmiles(string) != None:
        valid_smiles.append(string)

# Valid な SMILES を表示
print(valid_smiles)

# 生成した分子の構造を表示
Chem.MolFromSmiles(valid_smiles[0])

input_smiles を別の SMILES に書き換えることで、好きな化学構造を構造展開することができます。
また、発生させる潜在空間座標の数は、generated_coordinates の 3 つ目の引数を変えることで調整できます。
今回は 100 座標発生させましたが、SMILES として正しい文字列は 3 つしか出てきませんでした。しかもうち一つは、二重結合を二つ含む 4 員環を含んでおり、不安定だと考えられます。ただ、出てきた構造の一つはアニリンだったので、ちゃんと入力のベンゼンを構造展開した化合物が得られました。

Chemical VAE は正直性能は低く、論文の比較対象としてよく用いられていますが、変分モデルを構造生成器に持ち込み、コードを公開し、一大分野を築いたという点で、その業績が色褪せることはないでしょう。