「ベイズ反応最適化」ツールで遊んだ話 ~EVML, auto-QChem, EDBO~
前回、化学合成のための反応条件最適化ツールの文献をご紹介しました。
Shields, B.J., Stevens, J., Li, J. et al. Bayesian reaction optimization as a tool for chemical synthesis. Nature 590, 89–96 (2021).
こちらの文献、データ・コードともに公開してくださっています。
- auto-QChem : Quantum Chemistry
- EDBO : Experimental Design via Basian Optimization
- EVML : Expert versus Machine Leaning
せっかくなので今回はこちらで遊んでみたいと思います。
EVMLのゲームを自分で遊んだ後、EDBOに従って遊んだらどうなるか検証してみます。
反応最適化ゲーム(EVML)
概要
まずは、EVML(Expert versus Machine Leaning)(GitHub)です。
これは、開発された反応条件最適化ツール(EDBO)の性能を比較するためのWebゲームアプリです。
化学者に以下のような条件最適化ゲームをしてもらい、その結果をソフトウェアと比較し、ベンチマークとしています。
Reaction optimization game
ボスから司令のあった化合物がパラジウム触媒による直接アリール化で合成できることがわかった。
パラメータ(base, Ligand, solvent, concentration, temperature)を変えて、反応収率を最適化したい。
ボスから与えられた期限は1ヶ月。1日にこなせる実験のキャパシティーは5実験(1 batch)。
つまり最大100実験のうちに最適化しなければならない(全組み合わせの約6%)。
使えるのは、自分の知識と、各実験から得た情報のみ(ネットや本をみたらダメ)。
さあ、やってみよう!*1
このゲーム、実際の実験データをつかっているのがすごいですよね。
正直、配位子(Ligand)は色々開発されすぎててよくわかりません。・・・トリフェニルホスフィンの安心感。
塩基(Base)も溶媒(Solvent)も研究室ではみなかったものが多いです。製薬企業さんではよく使われるものなんでしょうか???
やってみる
早速あそんでみましょう!
GitHubのファイル(aws_setup.TXT)には、アマゾンウェブサービス(AWS)でのインストール方法が書いてありますが、 よくわからないのでMacで試します。
RのShinyパッケージでつくられたWebアプリのようなので、Rで実行すれば良さそうです。
以下のRのパッケージが必要そうなのでインストールしました。*2
- shiny
- dplyr
- ggplot2
- reshape2
- rmarkdown
- DT
EVMLのGitHubからフォルダをダウンロードしてきたら、「arylation_appフォルダ」の「app.R」ファイルを実行すればOKです。*3
実行すると、アプリが起動し以下のようなページが開きました!
左側のパネルに、ゲーム利用者の情報をいれたり、反応条件をいれたりします。
プルダウンから条件を選んで「add」し、5つ条件を選んだら「Run Experiments」をクリック!
選んだ条件の収率(yield)は「Experimet Results」タブに表示されます。 収率を見ながら条件を変えて最適化を目指しましょう!
GitHubの「results」フォルダ「analysis.ipynb」をみると、50人の化学者による結果がまとめられています。
こんな感じです。下図、上が全体で、下3つは属性で分けたサブセットです。
全体では、5バッチ前後で最高収率(100%近く)に到達している方が多い印象でしょうか?
所属の違いをみるとfacultyやEngineerは比較的直線的な改善を見せているのに対し、 Med ChemistやProcess Chemistは横ばいとなっているところがあるのが面白いですね。 条件検討の進め方にも考え方の違いがありそうです。
ベイズ最適化ソフト(EBDO)で50回行った結果も紹介されています。以下の通りです。
こちらも5バッチ前後で最高収率に到達しているものが多そうです。ソフトは毎回Maxに到達するまで検討しつづけるのが、ヒトとの違いが見えて面白いですね。
詳細な分析と結果の考察は文献本文をご参照ください。
上から目線でコメントしましたが、私は100%近く(>99%)になるまで7バッチでした。完全にソフトに負けています。
文献で一度結果を読んでいたはずなのに。。。
合成反応の記述子DB(auto-QChem)とベイズ最適化ツール(EDBO)
続いて、ベイズ最適化を実施するためのパッケージauto-QChemとEDBOです。 前者は、関連する化学物質を記述子に変換・保管するデータベースで、後者はベイズ最適化パッケージです。
EDBOのGitHub、examplesやexperimentsディレクトリのipynb
ファイルには、
論文中のデータや解析がコードとともに掲載されています。ありがたいですね!
以下の記事はだいたいそこからのコピペ切り貼りです。
この記事の目標は、上記の反応最適化ゲームのデータを題材にして各ツールを使ってみることです。 論文の記述子は量子化学計算(Gaussian)によるものですが、Gaussianは使えないしよくわからないので、 より一般的なケモインフォの記述子(Mordred)を用いたいと思います。
auto-QChem
まずは、合成反応を入力に変換し保管するauto-QChemです。MongoDBタイプのデータベースで、ユーザーガイドやインターフェースについては以下で確認できます。
- auto-QChem GitHub
- auto-QChem GitHub documentation
- auto-QChem Database User Guide
- auto-QChem landing page *4
auto-QChemは名前の通り量子化学(Quantum Chemistry)記述子向けで、計算の実行を簡単にする(入力ファイルの作成と出力の解析)というものです。 Mordredを使うのであれば不要ですが、せっかくなのでさわってみます。
インストール
インストールの仕方はGitHubのinstall.mdに書いてあります。
condaで仮想環境を作って、関連するパッケージをインストールした後、GitHubから落としてきたauto-QChemのsetup.py
を使えばOKです。
# 新しい仮想環境を作ってactivate conda create --name autoqchem python=3.7 conda activate autoqchem # パッケージのインストール conda install jupyter pandas scipy matplotlib pymongo pyyaml fabric xlrd appdirs openpyxl conda install -c conda-forge openbabel=2.4.1 # v.3.0.0はよくないらしい conda install -c rdkit rdkit python -m pip install imolecule==0.1.13 #3dで描画するのに使う
ついで、git clone
などでauto-QChemをローカルにもってきたあと、保存したディレクトリに移動し、setup.py
を実行します。
python setup.py install
おしまい!*5
分子の取り扱い(moleculeクラス)
auto-QChemでは分子はmolecule
クラスで扱います。OpenBabelのOBMol
クラスを継承したクラスになっています。
デフォルトではmolecule
クラスの入力は分子のSMILESで、3次元構造(conformer)を構築します。
コンフォメーションはOpenBabelのOBConformerSearch()
で、遺伝的アルゴリズムを使って生成されているようです。
試しにリガンドを一つ読み込んでみましょう!
論文でフィーチャーされてたCgMe-PPhというリガンドです。PubChemではこちらの名前では引っかからずmeCgPPhがシノニムにありました(PubChem CID: 53427790)
from autoqchem.molecule import molecule # Cg-MePPhのSMILESから分子を構築 smiles_str = "CC12CC3(OC(O1)(CC(O2)(P3C4=CC=CC=C4)C)C)C" mol = molecule(smiles_str) print(type(mol)) # <class 'autoqchem.molecule.molecule'>
draw()
という関数が用意されており、jupyter notebook上でぐりぐり動かせる3D描画が見られます。
# mol.draw()
かわいい。PubChemと描画スタイル同じですね。PubChemもimolecule
使っているんでしょうか?
コンフォマーを複数発生させることもできます。また、logging
でログを表示させられるので、動作確認しておきたいときは表示させても良いかもしれません。
# ログの表示 import logging logging.basicConfig(level=logging.INFO) # 最大30のコンフォマー作成 mol_conf = molecule(smiles_str, max_num_conformers=30) # INFO:autoqchem.molecule:Initializing molecule with canonical smiles: CC12CC3(C)OC(P2c2ccccc2)(CC(O1)(O3)C)C # INFO:autoqchem.molecule:Creating initial geometry with option 'best'. # INFO:autoqchem.molecule:Initial geometry created successfully. # INFO:autoqchem.molecule:Conformer Search generated 6 conformations of CC12CC3(C)OC(P2c2ccccc2)(CC(O1)(O3)C)C molecule
6つのコンフォマーができたみたいです。comformer_num
で各コンフォマーにアクセスできます*6。
違いがわからない。。。
その他の使い方や、Gaussianのインプットファイルの作成方法等はFunctional documentationをご参照ください。 また、データベースのインターフェースはFlaskを使っているようですがよくわからないので飛ばします。
Mordredで記述子計算
ベイズ最適化パッケージに進む前に、入力の記述子を用意しておきたいと思います。必要な分子(Ligand、Base、Solvent)のSMILESは以下です。
Ligand_smiles = ["CC(C)C1=CC(C(C)C)=C(C(C(C)C)=C1)C2=C(P(C3CCCCC3)C4CCCCC4)C(OC)=CC=C2OC", "CC(C)(C)P(C1=CC=CC=C1)C(C)(C)C", "CN(C)C1=CC=CC(N(C)C)=C1C2=CC=CC=C2P(C(C)(C)C)C3=CC=CC=C3", "P(C1CCCCC1)(C2CCCCC2)C3CCCCC3", "P(C1=CC=CC=C1)(C2=CC=CC=C2)C3=CC=CC=C3", "CC(C1=C(C2=CC=CC=C2P(C3CCCCC3)C4CCCCC4)C(C(C)C)=CC(C(C)C)=C1)C", "P(C1=CC=CO1)(C2=CC=CO2)C3=CC=CO3", "CP(C1=CC=CC=C1)C2=CC=CC=C2", "CC(OC1=C(P(C2CCCCC2)C3CCCCC3)C(OC(C)C)=CC=C1)C", "FC(F)(F)C1=CC(P(C2=C(C3=C(C(C)C)C=C(C(C)C)C=C3C(C)C)C(OC)=CC=C2OC)C4=CC(C(F)(F)F)=CC(C(F)(F)F)=C4)=CC(C(F)(F)F)=C1", "C[C@]1(O2)O[C@](C[C@]2(C)P3C4=CC=CC=C4)(C)O[C@]3(C)C1", "CP(C)C1=CC=CC=C1"] Base_smiles = ["O=C([O-])C.[K+]", "O=C([O-])C(C)(C)C.[K+]", "O=C([O-])C.[Cs+]", "O=C([O-])C(C)(C)C.[Cs+]"] Solvent_smiles = ["CCCCOC(C)=O,CC1=CC=C(C)C=C1", "CCCC#N", "CC(N(C)C)=O"]
それぞれMordredで記述子を計算した後、エラーを取り除きcsvで書き出します。auto-QChemのnotebook Fast featurize with Mordred.ipynbを参考にして進めていきます。
まず、必要なライブラリをインポート。
import pandas as pd import numpy as np import mordred from rdkit import Chem from mordred import Calculator, descriptors
数が少ないSolventを例に進めます。PandasのDataFrameに記述子を格納したいので、SMILESのlistをDataFrameにしておきます。
ついでにrdkitのrdmol
オブジェクトに変換して列を加えておきます。
solvent_df = pd.DataFrame(Solvent_smiles, columns={'smiles'}) solvent_df['rdmol'] = solvent_df['smiles'].map(lambda x: Chem.MolFromSmiles(x)) print(solvent_df.shape) # (4, 2)
ModredのCalculator
で記述子を計算します。
calc=Calculator(descriptors, ignore_3D=True) md=calc.pandas(solvent_df['rdmol']) print(md.shape) # (4, 1613)
1613個の記述子が計算されました。エラーが出ている部分をNaN
に変換して、dropna
で列を削除します。
# エラーをNaNで置き換え md=md.applymap(lambda x: np.nan if type(x) in [mordred.error.Missing, mordred.error.Error] else x) # NaNを含む列の削除 md=md.dropna(axis=1) print(md.shape) # (4, 1275)
1275個残りました。あとはcsvとして保存するだけです。
SMILESのDataFrameと記述子のDataFrameをconcatenate
で結合します。rdmol
オブジェクトの列は不要なので削除しておきます。
また、他の要素と区別するため、列名にprefix
「solvent_」をつけておきます。
# 保存用のDataFrameを作成。rdmolは削除しておく。 df_to_save = pd.concat([solvent_df.drop('rdmol', axis=1), md], axis=1) # prefixをつける df_to_save = df_to_save.add_prefix('solvent_') # csvに保存 df_to_save.to_csv("solvent_desc.csv", index=False)
無事「solvent_desc.csv」ファイルが作成できました!
同様にして「ligand_desc.csv」「base_desc.csv」を作成しました。これで記述子の準備はおしまい。
ベイズ最適化パッケージ(EDBO)
ではいよいよ本命、EDBO(Experimental Design via Basian Optimization)です。
これはベイズ最適化を化学合成に適応するための実用的なパッケージです。
- GitHub : EDBO GitHub
- ドキュメント : edbo documentation
インストール
まずはGitHubのREADMEに従ってインストールします。
EDBO自体はpip
でインストールできます。必要なライブラリはrdkit、Mordred、PyTorchのようです。
指示に従ってconda
で新しい仮想環境(edbo
)を作成し、諸々インストールしました。
# https://github.com/b-shields/edbo # (0) Create anaconda environment conda create --name edbo python=3.7.5 # (1) Install rdkit, Mordred, and PyTorch conda activate edbo conda install -c rdkit rdkit conda install -c rdkit -c mordred-descriptor mordred conda install -c pytorch pytorch=1.3.1 # (2) Install EDBO pip install edbo # Running Notebooks conda install jupyterlab
おしまい。
計算の流れとモデルの中身
EDBOで反応を最適化する流れは以下の通りです。GitHubの例(Bayesian Optimization of a Mitsunobu Reaction)を参考にしています。
- 反応をエンコード、入力の記述子を取り込み:
edbo.utils
- 反応空間(reaction space)のパラメーター設定:
components``encodings``descriptor_matrices
- ベイズ最適化の設定(ebdoインスタンス化):
edbo.bro.BO_express
- 計算ループの実施
- 初期化(最初の実験セットを提案):
init_sample()
- 実験結果の追加 :
add_results()
- モデルの適合、獲得関数の最適化、次の実験の提案 :
run()
- 提案実験の出力 :
export_proposed()
- 提案実験の結果を追加し、収束するまで 2 -> 4を繰り返す
- 初期化(最初の実験セットを提案):
上記手順「3. ベイズ最適化の設定」でedbo.bro.BO_express
というオブジェクトを利用しています。
これは、パラメーターの自動的な読み込みや、全実験空間の組み合わせの計算、必要な前処理の実施といった作業を行ってれるもので、
edboをよりユーザーフレンドリーなプログラムにしてくれています。
BO_express
はベイズ最適化のモデルもデフォルトで設定してくれています。デフォルトの設計は以下となっています。(class edbo.bro.BO_express)
ガウス過程はGPyTorch
を使って実装されています。*7関連するやつを貼っておきます。*8*9
では順番に計算していきましょう!
EDBOでEVMLをプレイ
今回は「EDBOの指示に従ってEVMLをプレイする」という体で、計算を実行したいと思います。つまり、EDBOの計算で提案された実験をEVMLで実施し、結果をEDBOにフィードバックして次の実験を再度提案させ、収率の向上を目指します。
ではまず反応のエンコードです。
入力として使う分子の記述子を読み込みます。edbo.utils
を使って、Mordredの計算結果を含むcsvファイルを読み込みます。
計算実行
まず反応のエンコードです。
入力として使う分子の記述子を読み込みます。edbo.utils
を使って、Mordredの計算結果を含むcsvファイルを読み込みます。
# インポート import pandas as pd from edbo.utils import Data # 記述子計算のファイルを読み込み ligands_desc = Data(pd.read_csv('data/ligand_desc.csv')) solvents_desc = Data(pd.read_csv('data/solvent_desc.csv')) bases_desc = Data(pd.read_csv('data/base_desc.csv'))
読み込んだデータはdata
で確認でき、PandasのDataFrameとなっています。
print(type(ligands_desc)) print(type(ligands_desc.data)) # <class 'edbo.utils.Data'> # <class 'pandas.core.frame.DataFrame'>
次に、反応空間(reaction space)のパラメーター設定を行います。今回の合成反応条件で最適化するパラメータ(concentration、temperature)を含めて要素を指定する辞書を作成します。
# 反応空間のパラメーター components = {'Ligand':'mordred', 'Solvent':'mordred', 'Base':'mordred', 'Concentration' : [0.057, 0.100, 0.153], 'Temperature' : [90, 105, 120]} # Encodingの指定。指定しない場合、自動的にOne-Hot-Encoding (OHE)になる。 encoding = {'Concentration' : 'numeric', 'Temperature' : 'numeric'} # 記述子の設定 (descriptor_matrices) mordred = {'Ligand': ligands_desc.data, 'Solvent' : solvents_desc.data, 'Base' : bases_desc.data}
descriptor_matricex
の設定はEDBO外部で計算した記述子を使う場合に設定する項目で、Gaussianで計算したDFT記述子の場合などに使います。Mordredの計算はEDBOでも行えるので、実際にはここで指定しなくてもよかったのですが、今回は前もって計算したデータを使うことにしました。
次に、EDBOの BO_express
オブジェクトをインスタンス化します。上で指定した要素やモデルの設定など、もろもろを書き込みます。
from edbo.bro import BO_express # BO object bo = BO_express(reaction_components = components, encoding=encoding, descriptor_matrices=mordred, acquisition_function='EI', # 獲得関数はEI init_method='rand', # 初期化はランダム batch_size=5, # 1ラウンドあたりの実験数は5 target='yield' ) # yieldを最適化 # 事前分布の指定。文献の設定に合わせる from gpytorch.priors import GammaPrior bo.lengthscale_prior = [GammaPrior(2.0, 0.2), 5.0] bo.outputscale_prior = [GammaPrior(5.0, 0.5), 8.0] bo.noise_prior = [GammaPrior(1.5, 0.5), 1.0]
下準備が完了したので計算ループに進みます。
まずは初期化!最初にする実験をランダムに選択します。
選択された実験はbo.export_proposed()
でcsvファイルとして書き出すことができ、またnotebook上でもbo.get_experiments()
で確認できます。
bo.init_sample(seed=0) # 初期化 bo.export_proposed('data/init.csv') # CSVファイルに書き出し bo.get_experiments() # ノートブック上で確認
5つ実験が選ばれました!この条件で実験を行い、結果を追加して次のループに進むことになります。
実際の環境では、実験に時間がかかるのでPCを付けっ放しにしておくのはもったいないです。そのような場合に設定を保存しておけるよう、EDBOではpickel
化するメソッドedbo.bro.BO.save()
も用意されています。
bo.save() # optimizerの保存
作業ディレクトリに「BO.pkl」が作成されました。再度読み込みたいときはedbo.bro.BO.load()
を使えばOKです。
from edbo.bro import BO_express bo = BO_express() bo.load()
さて、話戻して、出力されたcsvファイルを見ると、以下のようにyield列があり「<Enter Response>」となっています。ここに収率を書き込めば、次のループで使用する入力ファイルとできます。
今回取り上げている反応の結果は、EBDOのGitHub 「experiments/data/direct_arylation/experiment_index.csv」に記載されています。もしくはEVMLで確認することもできます。
今回は「EDBOの指示に従ってEVMLをプレイする」ということにしているので、EVMLに条件を入力し、得られた結果をcsvに追記し、次のラウンドに回しました。
結果を記入したファイルはedbo.bro.BO_express.add_results()
メソッドで読み込みます。ついでrun()
を実行すれば、ガウス過程モデルのfit
と獲得関数(EI
)の最適化を行い、次にすべき実験を提案してくれます。
# 結果をboオブジェクトに追加 bo.add_results('data/init_results.csv') # 計算の実施 bo.run() # 次にする実験の取り出し bo.export_proposed('data/round0.csv') bo.get_experiments()
新しい提案の実験条件が手に入りました(round0.csv
)。これをEVMLに入れて結果を得て、再度EDBOで計算、次の条件を提案を繰り返します。
・・・結果・・・・
round6の提案実験で、収率100%の反応条件が2つ見つかりました!
最初のinitバッチを含めて計7バッチ(35実験)です。だめ押しでround7も実験しましたが、収率は下がったのでここで止めました。
EVML上で確認できる「条件探索の過程」は以下の図のようになりました。
右から左に向かって新しい実験です。赤色がEBDOの提案した実験の収率で、緑色はランダムに選んだ場合の実験収率です。
最初は10%にも満たなかった収率が徐々に上がってゆき、最終的に収率100%を達成している様が見て取れます。
結果の確認
EDBOは結果解析用のツールも提供しています。
例えば、最適化ループを回す過程でoptimizer
のconvergence
がどのように変化したかプロットすることができます。
bo.plot_convergence()
順調に収束している様子がわかりますね!
GitHubのnotebookで使われている解析も適用してみましょう!
実施した全40実験(8バッチ)の中で、最も収率の良かった上位5実験について反応条件を提示してみます。
全ての結果を読み込んで統合したDataFrameを作成します。カラムは反応条件(Ligand, Solvent, Base, Concentration, Temperature)と収率(yield)で、yieldで降順に並べ替えます。
results = pd.DataFrame(columns=bo.reaction.index_headers + ['yield']) for path in ['init', 'round0', 'round1', 'round2', 'round3', 'round4', 'round5', 'round6', 'round7']: results = pd.concat([results, pd.read_csv('data/' + path + '_results.csv', index_col=0)], sort=False) results = results.sort_values('yield', ascending=False) results.head()
できました!
上位5つのリガンド、塩基はどれも同じですね。温度はより高温、濃度は濃い方が良さそうですが、溶媒の影響もありそうです。
DMAc溶媒では薄め(0.057)でも収率96.6%なのが興味深いですね。
棒グラフも書いて見ます。
import matplotlib.pyplot as plt fig, ax = plt.subplots(1, len(results.columns.values[:-1]), figsize=(30,5)) for i, feature in enumerate(results.columns.values[:-1]): results[feature].iloc[:5].value_counts().plot(kind="bar", ax=ax[i]).set_title(feature) plt.show()
また、EDBOは構造式描画のためのChemDraw
というモジュールも用意されています。
上位5反応に出てきた分子の構造式を描画しましょう。
from edbo.chem_utils import ChemDraw for col in results.iloc[:,:3].columns.values: print('\nComponent:', col, '\n') cdx = ChemDraw(results[col].iloc[:5].drop_duplicates()) cdx.show()
Component: ligand_smiles_index
Component: solvent_smiles_index
Component: base_smiles_index
なるほど便利!!
やっぱり構造式でみるのがが一番しっくりくる!!
感想
以上、今回は前回の記事で取り上げたNatureの文献のツールで遊んで見ました。データとツールが公開されているというのはとてもありがたいですね!
特にベイズ最適化ツールEDBOは変更が必要な設定項目も少なく、コピペするだけで簡単に使えてとても便利な印象でした。
EDBOの実践として、記述子にMordredを使い、EVMLをプレイして見ました。最初は低い収率で類似の実験を検討しており、「大丈夫かなー、本当にこれでいいのかなー」と半信半疑でしたが、収率上昇の兆候が見え始めると一気に条件検討が進んでいき最終的に収率100%の条件を見つけることができました!
私は自分でEVMLをプレイした際、99%までしか行かなかったので、100%の結果を見てびっくりしました。・・・完敗です。
ところで、文献中では最終的に記述子として量子化学計算によるものが用いられていました。今回はMordredでも良い結果に到達できることが確認できたので良かったです。Gaussianなんて一生触れそうそうにないですからね!
溶媒の沸点を考慮して反応温度に拘束をかけるとか、追加の設定を加えたらより探索効率上がるかもなー、とか思いましたがその辺りどうなんでしょう???
いずれにせよ、面白いツールがどんどん出てきて公開されるのは見ていて楽しいですね!
ほとんどGitHubのコードそのままコピペした記事ですみませんでした。。。ではでは!
*1:READMEをそれっぽっくかえてみました
*2:参考「aws_setup.TXT」の「5. Install required packages」
*3:私はRStudioで実行しました。Rはほとんど触ったことがなかったのでRStudioの使い方に手間取りました・・・
*4:ランディングページ(LP)というのは、ユーザーがWebで最初にアクセスするページのことらしいです。知らなかった。
*5: 仮想環境でjupyter notebookの立ち上げにエラー(Bad config ~~)が出るときは、「pip install environment-kernels」をすれば動くようになるかもしれません。参考
*6:コンフォマーの数以上の数字を与えると、最後のコンフォマーが選ばれるそうです
*7:GPyTorchについてはHELLO CYBERNETICSさんの記事「GPyTorchでガウス過程を見てみよう」がとても参考になります。
*9:参考ガウス過程回帰の基礎 赤穂昭太郎 システム/制御/情報 2018(62)390