magattacaのブログ

日付以外誤報

オープンソースで逆合成解析 ~AiZynthFinder~ part 3. CLI

オープンソース自動逆合成解析ツール、AiZinthFinder。インターフェースとしてGUIコマンドラインの2種類が用意されています。

Command Line Interface (CLI)の利点はバッチ処理ができることです。ポチッとやって多数の分子を処理できたら楽そうですね!

ツールは以下のGitHubに、

github.com

ドキュメントはこちらにあります。

molecularai.github.io

課題設定

逆合成解析のバッチ処理の利用として、AI-drivenな創薬のDMTA探索サイクルが想定されているのだと思います。(生成モデルでデザインされた化合物をそのまま逆合成解析、合成ロボットに流す、みたいな?) 

ですが、私はどうして良いのかわからないので、もう少し簡単な課題を設定したいと思います。

つまるところ、効率よく多種類の化合物を手に入れたいという課題です。そのためにはちょうど良い共通中間体を設定し、できるだけ少ない反応の数で異なる目的物が得られるようにすればOKです。

もう少し遡ると、化合物群に共通して必要となる試薬を確保する必要がある、ってなことになりそうです。

と、いうわけで、今回の目標は、「論文記載の化合物群をAiZynthFInderのCLIでまとめて逆合成解析し、それらの合成に必要な共通の試薬をみつけよう!」です。

データの準備

AiZynthiFinderのCLIで解析を行うには「SMILESが1行1化合物分格納されたテキストファイル」必要になります。早速準備していきましょう。

対象のデータとしては、py4chemoinformatics 第8章で用意してくださっているOrexin Receptorアンタゴニストを利用させていただきたいと思います。

Table 1から下記のMerckの文献3つを選び、ChEMBLから化合物をSDFでダウンロードしてきました。

Doc ID Journal
CHEMBL3098111 Bioorg. Med. Chem. Lett. (2013) 23:6620-6624
CHEMBL3867477 Bioorg Med Chem Lett (2016) 26:5809-5814
CHEMBL3352684 Bioorg. Med. Chem. Lett. (2014) 24:4884-4890

RDKitでSDFをテキストファイルに変換していきたいと思います。

まずはどんな化合物たちか見てみます。

CHEMBL3098111はこんな化合物群、、、

from rdkit import Chem
from rdkit.Chem import AllChem, PandasTools

df_8111 = PandasTools.LoadSDF('./data_oxr/CHEMBL3098111.sdf')
df_8111.head(2)

f:id:magattaca:20201018190637p:plain

CHEMBL3867477はこんな感じ、、、

df_7477 = PandasTools.LoadSDF('./data_oxr/CHEMBL3867477.sdf')
df_7477.head(2)

f:id:magattaca:20201018190746p:plain

中央に位置する環がピリジンだったり4,4-ジフルオロピペリジンだったりするみたいですね。

では、テキストファイルに変換する関数をつくってみます。

一部、塩酸塩が含まれているようでしたので、SaltRemoverを使って脱塩処理をします。*1
また、重複した化合物を除くためにset()でリストからユニークなSMILESのみを残します。

こんな感じに書いてみました。

# SaltRemoverで脱塩処理を行うための設定
from rdkit.Chem import SaltRemover
remover = SaltRemover.SaltRemover(defnData="[Cl]")

def desalt_smis(sdf_data):
    # 空のリストを用意
    desalted = []

    Suppl = Chem.ForwardSDMolSupplier(sdf_data)
    
    for mol in Suppl:
        # 塩酸塩の脱塩処理
        res = remover.StripMol(mol)
        # SMILESに変換
        smi = Chem.MolToSmiles(res)
        desalted.append(smi)
    
    # unqueなSMILESのみ保持する    
    unique_smi = set(desalted)
    
    return desalted

3つのSDFに適用してテキストファイルを作成します。「data_oxr」というディレクトリにSDFを保存したのでos.pathでパスを指定しました。

CHEMBL_IDs =  ['CHEMBL3098111', 'CHEMBL3352684', 'CHEMBL3867477']

import os

for chembl_id in CHEMBL_IDs:
    # SDFデータのパスを指定
    sdf_file = os.path.join("data_oxr", chembl_id + '.sdf')
    
    # 作成した関数でSMILESに変換
    chembl_smis = desalt_smis(sdf_file)
    
    # 出力テキストファイルを指定
    txt_path =  os.path.join("data_oxr", chembl_id + '.txt')
    
    with open(txt_path, mode='w') as f:
        for smi in chembl_smis:
            f.write(smi + '\n')

1行1化合物とするため、改行\nでつないでいます。

が、もっと簡単の方法があるらしい!

その名もSmilesWriter! @iwatobipen先生がこちらの記事で紹介してくださっていました。

折角なので使ってみます。SDFを入力としてテキストファイルを出力する関数を書き直してみます。

from rdkit.Chem.rdmolfiles import SmilesWriter

def desalt_smi_write(sdf_data, text_file):
    # 空のリストを用意
    desalted = []

    Suppl = Chem.ForwardSDMolSupplier(sdf_data)
    
    for mol in Suppl:
        # 塩酸塩の脱塩処理
        res = remover.StripMol(mol)
        desalted.append(res)
   
    # SDWriterで書き出し
    writer = SmilesWriter(text_file)
    
    for mol in desalted:
        writer.write(mol)
    writer.close
    
    return None

for chembl_id in CHEMBL_IDs:
    # SDFデータのパスを指定
    sdf_file = os.path.join("data_oxr", chembl_id + '.sdf')
 
    # 出力テキストファイルを指定
    txt_path =  os.path.join("data_oxr", chembl_id + '.txt')
    
    # SDFからテキストに変換する関数
    desalt_smi_write(sdf_file, txt_path)

できました!3つとも無事テキストファイルとして出力されました。

今回は短いのであまりご利益を感じられないかもしれませんが、SMILESに変換してリストに格納する操作がいらなかったり、プロパティをもたせて書き出したりできるようなので便利そうです。

但し、SmilesWriterで書き出したテキストファイルは、各行のSMILESのあとに「空白 + インデックス」も追加されていました。

処理されるうえでは無視されるものかもしれませんが、AiZynthFinderで問題が起きないかわからないので、今回は無難に先に出力したファイルを利用することとします。

AiZynthFinder コマンドライン

データが用意できたのでCLIバッチ処理を行います。

データを格納したディレクトリ「data_oxr」に以下の4つも合わせて入れて起きます。

  1. config.yml
  2. zinc_stock_17_04_20.hdf5
  3. full_uspto_03_05_19_rollout_policy.hdf5
  4. full_uspto_03_05_19_unique_templates.hdf5

これらはデータと同じディレクトリにある必要はありませんが、configのパスを直すのが面倒なので入れてしまいました。・・・無精ですみません。

コマンド自体は簡単!ターミナルで作業ディレクトリを移動し、以下を実行するだけ*2

export KMP_DUPLICATE_LIB_OK=TRUE

aizynthcli --config config.yml --smiles CHEMBL3098111.txt

最初の1行は「OMP: Error #15」を回避するためのものなので、他の環境では必要ないかもしれません*3

で、こんな感じでターミナルに出力されていきます。

Loading policy model from full_uspto_03_05_19_rollout_policy.hdf5 and templates from full_uspto_03_05_19_unique_templates.hdf5 to full_uspto
Loading stock from zinc_stock_17_04_20.hdf5 to zinc
Selected stocks: zinc
Compounds in stock: 17422831
Selected policy is full_uspto
Done with Cc1cc(C)cc(-c2cnc(-c3cccnc3)c(C(=O)NCc3ccc4c(c3)OCO4)c2)c1 in 70.4 s
Done with Cc1cc(C)cc(-c2cnc(-c3cccnc3)c(C(=O)NCc3ccc(F)c(C)c3)c2)c1 in 69.8 s
Done with COc1ccc(CNC(=O)c2cc(-c3cccc(Cl)c3)cnc2-c2cccnc2)cc1OC in 75.0 s

分子にもよりますが1化合物2分前後なので、10化合物で20~30分くらいをみておけば良いでしょうか? Macがウィンウィン唸るのを放っておくと、「output.hdf5」というhdf5ファイルが出力されます。

ファイル名が重なると上書きされてしまいそうなので、それぞれChEMBL IDを後ろにつけて別のフォルダに保存しておきました。

では、解析結果を見てみましょう!

バッチ処理結果解析

CHEMBL3098111

出力はhdf5なので、PandasのDataFrameに読み込むことができます。*4

import pandas as pd
data_8111 = pd.read_hdf('./data_oxr/output_files/output_CHEMBL3098111.hdf5', "table")

こんな感じ

data_8111.head()
target search_time number_of_nodes max_transforms max_children top_score is_solved number_of_steps number_of_precursors number_of_precursors_in_stock precursors_in_stock precursors_not_in_stock top_scores trees
0 Cc1cc(C)cc(-c2cnc(-c3cccnc3)c(C(=O)NCc3ccc4c(c... 70.419213 545 7 19 0.986553 True 3 4 4 OB(O)c1cccnc1, Cc1cc(C)cc(B(O)O)c1, NCc1ccc2c(... 0.9866, 0.9866, 0.9866, 0.9866, 0.9866, 0.9866... [{'type': 'mol', 'smiles': 'Cc1cc(C)cc(-c2cnc(...
1 Cc1cc(C)cc(-c2cnc(-c3cccnc3)c(C(=O)NCc3ccc(F)c... 69.782567 592 7 25 0.986553 True 3 4 4 OB(O)c1cccnc1, Cc1cc(C)cc(B(O)O)c1, Cc1cc(CN)c... 0.9866, 0.9866, 0.9866, 0.9866, 0.9750 [{'type': 'mol', 'smiles': 'Cc1cc(C)cc(-c2cnc(...
2 COc1ccc(CNC(=O)c2cc(-c3cccc(Cl)c3)cnc2-c2cccnc... 74.974893 585 7 23 0.986553 True 3 4 4 OB(O)c1cccc(Cl)c1, OB(O)c1cccnc1, COc1ccc(CN)c... 0.9866, 0.9866, 0.9866, 0.9866, 0.9750 [{'type': 'mol', 'smiles': 'COc1ccc(CNC(=O)c2c...
3 Cc1cc(C)cc(-c2cnc(-c3cccnc3)c(C(=O)NCc3ccc(N(C... 71.364885 590 7 22 0.986553 True 3 4 4 OB(O)c1cccnc1, Cc1cc(C)cc(B(O)O)c1, CN(C)c1ccc... 0.9866, 0.9866, 0.9866, 0.9866, 0.9866 [{'type': 'mol', 'smiles': 'Cc1cc(C)cc(-c2cnc(...
4 COc1ccc(CNC(=O)c2cc(-c3cc(F)cc(F)c3)cnc2-c2ccc... 63.741923 587 7 22 0.986553 True 3 4 4 OB(O)c1cc(F)cc(F)c1, OB(O)c1cccnc1, COc1ccc(CN... 0.9866, 0.9866, 0.9866, 0.9750, 0.9634 [{'type': 'mol', 'smiles': 'COc1ccc(CNC(=O)c2c...

どのような合成ルートが提案されたか、画像ファイルとして出力することができます。

ドキュメントをコピペ・・・・

from aizynthfinder.analysis import ReactionTree

# 全化合物について全てのツリーのリストを取得
all_trees = data_8111.trees.values

# 最初の化合物についての情報を取り出す
trees_for_first_target = all_trees[0]

# 各ツリーを画像ファイル(png)として書き出していく
for itree, tree in enumerate(trees_for_first_target):
    imagefile = f"route{itree:03d}.png"
    ReactionTree.from_dict(tree).to_image().save(imagefile)

実行すると最初のターゲット化合物について11個の合成ルートがそれぞれpngファイルとしてフォルダに出力されました。

こんな感じです。

f:id:magattaca:20201018191630p:plain

鈴木-宮浦クロスカップリングが2回使われています。ClとBrがあって先にClを反応させるルートになっています。反応性としては逆になってしまいそうです。

では、最初にもどって課題「論文記載の化合物群の合成に必要な共通の試薬」を検証してみたいと思います。

逆合成解析を進めた終点、試薬はDataFrameのprecursors_in_stockというcolumnに格納されています。

一つ目の例を取り出してみます。

precursors_8111 = data_8111['precursors_in_stock']
print(type(precursors_8111[0]))
print(precursors_8111[0])

#  <class 'str'>
#  OB(O)c1cccnc1, Cc1cc(C)cc(B(O)O)c1, NCc1ccc2c(c1)OCO2, O=C(O)c1cc(Br)cnc1Cl

str型で一つの文字列中に複数のSMILESが「カンマ+空白」(, )で区切られて格納されています。

構造式に起こしてみましょう。文字列をsplitで区切って、SMILES一つずつに分解した後、Molオブジェクトに変換しています。

precursor_mols = []

l = precursors_8111[0].split(', ')
for smi in l:
    m = Chem.MolFromSmiles(smi) 
    precursor_mols.append(m)

Chem.Draw.MolsToGridImage(precursor_mols, molsPerRow=4)

f:id:magattaca:20201018191820p:plain

合成ルートと見比べてみてきちんと目的の原料が取り出せているようです。

では、全ての化合物についてそれぞれで提案されているルートの出発原料を取り出し、集計してみます。
集計結果で出現頻度が多くなった原料は、複数の化合物の合成に必要となる試薬なので、優先順位高く確保する必要がありそうです。

具体的な処理としては、原料のSMILESをリストとしたのち、Python標準ライブラリcollectionscounterを使って出現回数に基づく処理を行います。*5

# 全化合物の全ての原料のSMILESを格納するためのリストを用意
precursors_list_8111 = []

for i in range(len(precursors_8111)):
    # 原料のSMILES 一つずつを要素としたリスト
    l = precursors_8111[i].split(', ')
    
    # extendでリストの要素を追加
    precursors_list_8111.extend(l)

import collections

# Counterで各要素の出現頻度を数える
c_8111 = collections.Counter(precursors_list_8111)

# most_common()で出現頻度で並べ替え
mc_8111 = c_8111.most_common()

print(mc_8111)

#   [('OB(O)c1cccnc1', 17), ('O=C(O)c1cc(Cl)cnc1Cl', 10), ('Cc1cc(C)cc(B(O)O)c1', 9), ('COc1ccc(CN)cc1OC', 8), ('O=C(O)c1cc(Br)cnc1Cl', 6), ('OB(O)c1cc(Cl)cc(Cl)c1', 2), ('NCc1ccc2c(c1)OCO2', 1), ('Cc1cc(CN)ccc1F', 1), ('OB(O)c1cccc(Cl)c1', 1), ('CN(C)c1ccc(CN)cc1', 1), ('OB(O)c1cc(F)cc(F)c1', 1), ('CNCc1ccc2nccnc2c1', 1), ('CNCc1ccc2c(c1)OCO2', 1), ('COc1ccc(CN)cc1', 1), ('Cc1cc(C)cc(B2OC(C)(C)C(C)(C)O2)c1', 1), ('CN1CCc2cc(CN)ccc21', 1), ('COc1cccc(CN)c1', 1), ('NCc1ccc(C(F)(F)F)c(F)c1', 1), ('Cc1cc(F)cc(B(O)O)c1', 1), ('COC(=O)c1cc(Cl)cnc1Cl', 1), ('Nc1cc(C(F)(F)F)cc(C(F)(F)F)c1', 1), ('OB(O)c1cc(F)cc(Cl)c1', 1)]

(要素, 出現回数)という形のタプルが出現回数順に並んでいます。

では、マストバイアイテム3つを見てみましょう!

must_buy_reagents = []
occurence = []

for i in range(3):
    smi = mc_8111[i][0]
    m = Chem.MolFromSmiles(smi)
    must_buy_reagents.append(m)
    occurence.append(str(mc_8111[i][1]) + ' times')

Chem.Draw.MolsToGridImage(must_buy_reagents, legends=occurence)

f:id:magattaca:20201018192113p:plain

図の下、数字はこれらの試薬が出現した回数です。全17化合物ですので、3-ピリジンボロン酸はすべてに共通する試薬のようです。

とりあえずこの3つを購入してカップリングで骨格を作っておけば、様々なフラグメントとアミド結合をつくることで多種類の化合物を合成することができそうです。

カップリングの選択性は・・・頑張る!

CHEMBL3867477

CHEMBL3352684はCHEMBL3098111と構造が似ていたので割愛し、CHEMBL3867477の結果を見てみます。

まずは、提案ルート例から。コードは同じなので割愛します。

f:id:magattaca:20201018192251p:plain

で、出現頻度の高い上位3試薬はこちら、、、

f:id:magattaca:20201018192412p:plain

例に挙げた合成ルートでは、アルコールとフェノールの光延反応という感じですが、頻度の多い試薬を見るとフッ素原子を脱離基として求核置換のルートが多く提案されているようですね。

全49化合物ということを考えると、これらの試薬は半数前後で必要となる試薬です。これは買い、でしょうか???

まとめ

今回はAiZynthFinderのコマンドラインバッチ処理を行ってみました。複数の化合物の解析結果、共通する原料・試薬があれば、多種類の化合物を作る上で重要となりそう、ということで検証してみました。

論文というストーリーのある化合物群なので、互いに類似した化合物群があげられている、という前提はありますが、頻度の高い重要そうな試薬が結構いい感じでみつかったかな、という印象です。

さて、なぜ今回こんな解析をしてみたか? 

ここで思い出話を・・・

私は学生の頃、構造活性相関のようなテーマを行なっていました。
ダメな学生だったので指導教官から何度か急かされていたものの、ノラリクラリとかわしつつのんびり実験をしていました。

ある土曜の午後でした。TLCを焼いていると、指導教官が無言で足早に近づいてきます。

「早くテーブルの化合物作って!!!!」

温厚な指導教官も流石に堪忍袋の緒が切れたようでした。

やべぇやべぇ。至急、合成しなければ。。。

と、いうわけで多種類の化合物を急いで作らないといけないシチュエーションでした。
まずは共通の試薬を確保せねば!と、・・・必死でラボの試薬庫と冷蔵庫あさりました。

こんな時、今回みたいに重要な試薬がわかったら何を探せば良いか分かって楽だったかもしれませんね!

まあ、ちゃんと実験してればお叱り受けなかったですけどね!!

今回も色々と間違いが多いと思うので、ご指摘いただければ幸いです。

ではでは!