Python で量子化学計算
明けましておめでとうございます。
今年も、創薬化学、ケモインフォマティクス を中心にブログを更新していこうと思います。
量子化学計算といえば、Gaussian や GAMESS を利用している方も多いと思いますが、Python には量子化学計算をしてくれるライブラリがあります。
psicode.org
conda install psi4 psi4-rt python=3.7 -c psi4
で使えるようになります。
使い方ですが、基本的な使い方はこちらのブログにまとまっています（この記事もこちらを参考にして書いています）。
future-chem.com
また、psi4 を使いやすくした psikit なるものもあります（この記事もこちらを参考にして書いています）。
github.com
これだけではこの記事の新規性がないので、いくつかの分子をまとめて量子化学計算できるコードを作成したので、これを元に psi4 の基本的な使い方をまとめます。
import numpy as np import pandas as pd from rdkit import Chem from rdkit.Chem.rdDistGeom import ETKDGv3, EmbedMolecule from rdkit.Chem.rdForceFieldHelpers import MMFFHasAllMoleculeParams, MMFFOptimizeMolecule from rdkit.Chem.Draw import IPythonConsole import psi4 import datetime import time # ハードウェア側の設定 psi4.set_num_threads(nthread = 2) psi4.set_memory("1GB") # 入力する分子 smiles_list = ["C", "CC", "CCC", "CCCC", "CCCCC"] for smiles in smiles_list: t = datetime.datetime.fromtimestamp(time.time()) psi4.set_output_file("../result/psi4_log/{}{}{}_{}{}_{}.log".format(t.year, t.month, t.day, t.hour, t.minute, smiles)) # SMILES から三次元構造を発生させる mol = Chem.MolFromSmiles(smiles) mol = Chem.AddHs(mol) params = ETKDGv3() params.randomSeed = 1 EmbedMolecule(mol, params) # MMFF で構造最適化 MMFFOptimizeMolecule(mol) conf = mol.GetConformer() # Psi4 に入力可能な形式に変換 mol_input = "0 1" for atom in mol.GetAtoms(): mol_input += "\n " + atom.GetSymbol() + " " + str(conf.GetAtomPosition(atom.GetIdx()).x)\ + " " + str(conf.GetAtomPosition(atom.GetIdx()).y)\ + " " + str(conf.GetAtomPosition(atom.GetIdx()).z) geometry = psi4.geometry(mol_input) # 汎関数、基底関数を設定 level = "b3lyp/6-31G*" # 構造最適化計算(単位： au) energy, wave_function = psi4.optimize(level, molecule = geometry, return_wfn = True) LUMO_idx = wave_function.nalpha() HOMO_idx = LUMO_idx - 1 # HOMO を表示（単位： au） homo = wave_function.epsilon_a_subset("AO", "ALL").np[HOMO_idx] lumo = wave_function.epsilon_a_subset("AO", "ALL").np[LUMO_idx] with open("../result/psi4_summary.txt", "a") as writer: output_string = "{}\nEnergy: {}\nHOMO: {}\nLUMO: {}\n\n".format(smiles, str(energy), str(homo), str(lumo)) writer.write(output_string)
最初に、ハードウェアの設定をします。スレッドやメモリを大きくすれば、計算速度も上がると思いますが、手元のマシンのスペックと相談することになります。
出力するログファイルも設定できますが、当然、分子ごとに異なるので、for 文の中で名前を決めてあげます。
量子化学計算では当然三次元構造が必要なので、SMILES から立体配座を作成します。
水素原子を付加した後、ETKDG というアルゴリズムで構造を発生させます。イメージ的には、distance geometry 法に近いと思いますが、アルゴリズムを色々改良しているようです。
コンフォメーションを発生させてから、MMFF94 で構造最適化をしていますが、やった方がいいのか、厳密な検討はしていません（通常の DG 法は最適化する）。
余談ですが、MMFF94S もありますが、結晶構造の場合、こちらを使った方がいいらしいです（厳密な比較はしていない）。
立体配座を作成したら、いよいよ、psi4 にインプットするデータを作成します。データは、Gaussian や GAMESS へのインプットファイルをイメージすると分かりやすいです。最初の行に形式電荷などの情報、その下に元素記号と座標を書いていきます。String として作成したら、psi4.geometry で psi の Molecule オブジェクトに変換しましょう。今回は XYZ 形式（座標をそのまま使う）でインプットファイルを作成していますが、Z-matrix 形式にした方が計算が収束しやすい気がします（case by case ですが）。
計算の基本設定は簡単で、汎関数/基底関数を決めて、psi4.optimize に引数として入力すると計算をしてくれます。psi4.optimize はデフォルトだとエネルギーしか返してくれないので、return_wfn を True にして、色々な情報が入っている wave_function オブジェクトを取り出しましょう。
一点計算を行いたい時は、optimize を energy に変更すれば同じようにできます。
wave_function.nalpha でHOMO の軌道の番号が取れますが、1 始まりなので、wave_function.epsilon_a_subset('AO', 'ALL') で軌道のエネルギーを取り出すときは、番号が一つずれる（0 始まり）なので、気をつけてください。
epsilon_a_subset の引数を変えることで、他の情報も取り出せそうですが、psi4 のドキュメントが思ったより充実していないので、試行錯誤する必要がありそうです・・・
計算結果はログファイルにまとめられていますが、分量が多く、見にくいので、エネルギーを txt ファイルに出力しておくと、結果を確認しやすいです。
データの単位ですが、RDKit で構造を発生させた場合、座標の単位はオングストローム(10^-10 m) です。
psi4 へのインプットファイルも同様です。
今回は出力していませんが、psi4 の出力座標の単位はボーアになっているので、長さのスケールには注意してください。
また、エネルギーの単位は原子単位（ハートリー）になっています。
単位はちゃんとドキュメントに書いていて欲しい・・・
psi4 もまだまだ理解していないことも多いので、随時勉強したことを紹介していければと思います。
とりあえず原子レベルでの情報を取り出せるようになりたい。
触媒医薬品：新しい医薬品の形と落ちこぼれの夢
この記事は創薬 (dry) Advent Calendar 2020 の 13 日目の記事です。
クリスマスに書く advent calendar とは
クリスマスですし、飲みながら書いているので、偶にはポエミーな記事でも。
医薬品のモダリティの多様化が叫ばれて久しいですが、ものすごく大まかに分類すると、低分子、バイオ医薬（特に抗体医薬）、中分子に大別されると思います。 低分子はご存知の通り、昔から盛んに研究がなされていますが、ターゲット分子の枯渇や、rule of 5 の制約、動態、毒性の制御などの観点から、その重要性は下がると言われています（chemical space の広さを考えると、まだまだ可能性はいくらでもありそうですが）。 バイオ医薬（ワクチンも含む？）は、低分子医薬がターゲットにしないようなタンパク質も狙えることから、選択性も相まって、近年盛んに研究され、実際に認可されている薬も増えてきています。 最も、コストがとんでもなく、薬価の高騰にも寄与しているという・・・（バイオ医薬ってどうやって設計しているんでしょうか？詳しい方に御教授頂きたいです） 中分子は、 Ro5 に縛られない、官能基高密度型天然物やペプチド医薬などを扱います。これまで余り扱ってこなかったので、可能性はありますが、低分子医薬の発想の延長ではないか、という気もします（穿ち過ぎか？）。
これらの医薬品、モダリティは違えども、共通点があります。それは、医薬品とそのターゲットの化学構造は変化しない、ということです。 代謝などを除けば、そんなことは当たり前です。しかし、こう考えることもできるはずです。
ターゲット分子（特にタンパク質）と相互作用するだけが医薬品の形だろうか？
例えば、タンパク質をぶっ壊すことで、疾患を治療できないか？
高校時代（10 年前）、初めて化学で酵素を勉強したとき、基質特異性という性質に感動したことを今でも覚えています。そして、こう思ったことも。
酵素を自由自在に、天然アミノ酸以外で作れたら面白いだろうなぁ
大学では薬学部にいましたが、暗記ばかりの授業に嫌気がさしていたところ、何となくネットサーフィンをしていたら、一つのワードに出会いました。その衝撃は、もしかしたらもう二度と体験できないかもしれません。
触媒医療
触媒の機能を医薬品として応用する。化学反応を薬にする。私のやりたいことを全て詰め込んだようなワードでした。 この研究を絶対したい！これなら自分の人生をかけられる！若き日の私は触媒医薬品の研究をしたい、という思いでいっぱいでした。
お気づきの方も多いと思いますが、この研究は、東大薬学部の金井先生の研究室でなされているものです。
触媒医療というワードに出会って 10 年以上経っていますが、残念ながら、触媒医療の研究を私はできておりません。有機化学が好きなはずなのに、何故か今は実験すらせず、機械学習というよく分からない怪しげなものに手を出している始末です。実験をせずして何が有機化学者か！
しかし、私は触媒医療の研究をするチャンスはいずれは訪れるのではないか、と信じています。触媒医薬を設計する際には、普通の医薬品とは異なり、基質と生成物両方について毒性などを考慮する必要があります。当然、三次元的な相互作用や遷移状態なども考慮する必要があるでしょう。そういった場面で、計算化学の力は必要になってくるはずです。
自分が生きているうちにはなし得ないかもしれない、でも、その基盤を作るのに貢献したい。無能な研究者でも、夢を追いかけることはできるはずだから。
（余談ですが、金井先生はベンチャー企業を立ち上げられているようです。すげぇ。計算化学もしてるっぽい。）
化学構造生成器 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 は正直性能は低く、論文の比較対象としてよく用いられていますが、変分モデルを構造生成器に持ち込み、コードを公開し、一大分野を築いたという点で、その業績が色褪せることはないでしょう。
Deep Graph Library によるグラフ畳み込みネットワークの基本（追記予定）
最近、論文を書いてたり会議が増えたり人事関連で色々あったり投資を始めたりごちうさを観たりで忙しく、中々記事を書けていませんでしたが、久々にブログを更新しました。
AI ブーム、深層学習ブームは化学分野でも続いていますが、化学者としては、構造式を取り扱いたいことが多く、構造式を AI に入力をして、物性とか薬理作用とか出てこないかなぁ、と思うこともあります。
それを実現するのがグラフ畳み込みネットワーク（Graph convolutional neural network: GCNN）です。
GCNN を使った研究は時々あるものの、GCNN って何をしているのか、どうやって実装するのか、意外と日本語でまとまっていないので、以下のサンプルコードを元に解説していきます。
なお、オリジナルのコードは、pen さんのブログにあり、それを劣化させました。
Try GCN QSPR with pytorch based graph library #RDKit #Pytorch #dgliwatobipen.wordpress.com
import numpy as np import pandas as pd import torch import torch.nn as nn import torch.optim as optim import torch.nn.functional as F from torch.utils.data import DataLoader from torch.autograd import Variable from sklearn.metrics import r2_score from rdkit import Chem from rdkit.Chem.Draw import IPythonConsole import dgl import dgl.function as fn from dgl import DGLGraph # 原子の種類のリスト element_list = ["C", "N", "O", "S", "P", "F", "Cl", "Br", "I", "B", "Si", "Se", 'H', "Unknown"] atom_features_dim = len(element_list) + 6 + 5 + 1 MAX_ATOMNUM =60 BOND_FDIM = 5 MAX_NB = 10 # 特徴リストの one hot vector を算出 def atom_feature_encoding(feature, feature_list): # 特徴リストに含まれない場合、 "Unknown" 表記にする if feature not in feature_list: feature = feature_list[-1] # 返り値は True or False で構成されるリスト return [feature == f for f in feature_list] # 各原子の特徴量を出力 # 原子の種類（23 次元）+ 結合次数（6 次元）＋電荷（5 次元）＋ 芳香族性（1 次元） def atom_features(atom): return (atom_feature_encoding(atom.GetSymbol(), element_list) + atom_feature_encoding(atom.GetDegree(), [0,1,2,3,4,5]) + atom_feature_encoding(atom.GetFormalCharge(), [-1,-2,1,2,0]) + [atom.GetIsAromatic()]) # 結合の種類（単結合、二重結合、三重結合、共役系、環構造 それぞれ 1 or 0 (True of False) で表現） def bond_features(bond): bond_type = bond.GetBondType() return (torch.Tensor([bond_type == Chem.rdchem.BondType.SINGLE, bond_type == Chem.rdchem.BondType.DOUBLE, bond_type == Chem.rdchem.BondType.TRIPLE, bond_type == Chem.rdchem.BondType.AROMATIC, bond.IsInRing()])) # Mol class の分子データを GCNN に入力する DGL Graph 型に変換 def mol2dgl_single(mols): # 各分子の DGL Graph 型を入れるリスト cand_graphs = [] n_nodes = 0 n_edges = 0 # 各分子ごとに処理 for mol in mols: # DGL Graph の型を用意 g = DGLGraph() # 各原子（グラフのノード）の特徴量を格納するリスト atom_feature_list = [] # 原子の情報（Atom class）を atom_feature_list に突っ込む for atom in mol.GetAtoms(): atom_feature_list.append(atom_features(atom)) # 分子中の原子数を DGL Graph 型に入力 g.add_nodes(mol.GetNumAtoms()) # 各原子の特徴量ベクトルを DGL Graph 型に "h" という名前で入力 g.ndata["h"] = torch.Tensor(atom_feature_list) bond_begin_list = [] bond_end_list = [] # 結合毎に処理、結合の向きも考慮 for bond in mol.GetBonds(): # 結合の始点と終点の原子のインデックスを取得 begin_idx = bond.GetBeginAtom().GetIdx() end_idx = bond.GetEndAtom().GetIdx() # 結合の特徴量を算出 features = bond_features(bond) # インデックスに関する情報をリストとしてまとめる bond_begin_list.append(begin_idx) bond_end_list.append(end_idx) # （無向グラフなので）向きを変えて同じ操作を行う bond_begin_list.append(end_idx) bond_end_list.append(begin_idx) # 隣接するノードのインデックスを DGL Graph 型に入力 g.add_edges(bond_begin_list, bond_end_list) # 出力用リストに入力 cand_graphs.append(g) # DGL グラフ型として出力する return cand_graphs
少し長いですが、上が構造式を GCNN に入力できるよう変換するためのコードです。
最初に、いつものようにライブラリを読んでいきます。
今回は深層学習ライブラリとして PyTorch を用います。Tensorflow よりも GCNN 周りのライブラリが充実している気がします。
GCNN のライブラリは Deep Graph Library (DGL) を用います。
PyTorch と GCNN は terminal から pip で簡単に入れられます。
pip install torch torchvision pip install dgl
構造式中の各原子（グラフの用語ではノードとも呼ぶ）ごとに、特徴量を binary vector（0 or 1 で構成されたベクトル）として算出します。
今回は、原子の種類は何か、結合次数と形式電荷はいくつか、芳香属性はあるか、これらの特徴を取り出します。
atom_feature_encoding という関数は、第一引数の feature （文字でも数字でも良い）が第二引数に入力しているリストに含まれるどの要素に対応しているかを確認します。返り値は引数に用いたリストと同じ長さの binary vector で、対応する要素のみ 1 、他の要素は 0 となっています。
atom features と bond features という関数は atom_feature_encoding を実行するだけの関数で、引数が RDKit の atom クラス bond クラスを用いているかという違いだけで、役回りは一緒です（ただ、pen さんのオリジナルのコードだと、bond クラスで取得した情報がその後のモデルで使われていない気がします）。
mol2dgl_single が本丸で、RDKit の mol オブジェクトを DGL Graph という GCNN に入力可能な型に変換する関数となります。
引数の mols は mol オブジェクトのリストです（Chem.SDMolSupplier で直接読み込んだものは使えないので注意。[mol for mol in supplier] といった変換が必要）。
cand_graphs が出力用のリストです。このリストの中に、各分子の DGL Graph オブフェクトが入ります。
for 文で各分子の特徴量を算出と DGL Graph オブフェクトの用意をしていきます。
最初に出力する DGL Graph インスタンスを g と定義して用意します。atom_features で原子の特徴量を計算したら、g.ndata["h"] に Tensor 型として入力します。この操作で、DGL Graph オブジェクトに "h" という名称でノード（原子）の情報を定義します。入力する Tensor 型は原子数 * 原子の特徴量数の二次元ベクトルです。
（以下、追記予定）
生物活性物質データベース ChEMBL の使い方
QSAR モデルを構築するためには、多くのデータが必要です。 化合物とその生物活性についてまとめられたデータベースは色々ありますが、よく使われるものの一つに ChEMBL というデータベースがあります。
アルツハイマー病治療薬であるドネペジルのターゲットでもあるアセチルコリンエステラーぜ（AChE）を例に、ChEMBL から QSAR モデル構築に使うためのデータを取得してみましょう。
最初に、トップページから、ターゲットとなるタンパク質を検索します。
検索結果が出てきますが、左上の Targets を選択すると、タンパク質のアッセイデータセットが複数出てきます。
カラムの右側の方に、「Compounds」という項目があるので、右側の矢印を押して、データ数の多い順にソートしましょう。一回だと、データ数の少ない順にソートされるので、二回押す必要があります。 Organism のカラムには細胞の由来が記載されているので、目的に応じてデータを選択してください（今回は Homo sapiens を選びます）。Name にはタンパク質の名前があるので、確認しましょう（例えば EGFR で検索すると、EGFR1 と EGFR2 両方のデータが出てきます）。 欲しいデータセットを見つけたら、緑色で書かれている ChEMBL ID をクリックしてください（今回は ChEMBL220 を選択してみます）。
Target Report Card のページが出てくるので、少し下に行くと、Associated Bioactivities と書かれた円グラフが出てきます。円グラフ中の欲しいデータの領域をクリックすると、そのデータだけ得られます（今回は IC50 を選択します）。タンパク質によりますが、IC50 が最も多く、次に Ki が多数を占める場合が多いです。
Browse Activities に移動したら、中央上部に緑で小さく書かれた Select All を押して、左側の Table をクリック、最後に、右側の緑の CSV or TSV ボタンを押すと、化合物と生物活性のデータを落とすことができます。
原因は不明ですが、うまくデータがダウンロードできない場合があります。そのときは、ブラウザを変更するなどの対応をすると、落とせる場合があります（直近だと何故か Safari で落とせず、Google Chrome に変更したら落とせた）。 ダウンロードしたデータは、化合物の重複があったり、濃度が違ったりするデータも含まれるので、必要なデータの前処理をしていく必要があります。
タンパク質名以外にも、アッセイ方法などで検索をかけることもできます。 色々なデータを実際に触ってみて、ケモインフォマティクスで使われるデータの感覚に慣れてみてください。
新しい化合物の表記法 SELFIES
色々と書きたいことは多いのですがサボり癖が強く忙しかったり疲れてたり、少し考えて文章を書きたかったりで、少し久々の更新ですが、新しい知識を最近得たのでシェアも兼ねて記事を書きます。
化合物を扱うとき、我々化学者は構造式を使うことが一般的です。構造式は化合物の特徴を一目で理解することができますが、こういった「絵」をコンピューター上で扱うのは一般に困難です。化合物を文字列表記することで、コンピューター上でも扱うことが可能となります。
最も広く使われているのは SMILES という記法です。 SMILES のメリットとしては、広く使われているので、多くのデータベースに対応している、ということでしょうか。また、文字列表記でも意味を理解しやすいという点もメリットです。 例えば、エタノールの SMILES は CCO、ベンゼンの SMILES は c1ccccc1、S- イブプロフェンの SMILES は CC(C)Cc1ccc(C@@HC(=O)O)cc1 といった具合です。 デメリットとしては、一つの化合物でも、複数の表記法がある、という点です。エタノールでさえ、OCC や C(C)O といった表記も考えられます。Canonical SMILES という、一つの分子に対し一つの記法が対応する、というルールもありますが、詳細は DayLight 社が公表していないため、正確なルールやアルゴリズムはよく分かっていません。実務的には、RDKit の Chem.MolToSmiles(mol) を使えば、「RDKit で定義された」canonical SMILES が出てくるので、これを使うのが良いと思います。 また、当然ながら、文法が不正確だと、分子として成立していません。CC((O やN)12CO などは、カッコや数字がおかしく、分子として成立していません。
SMILES 以外にも分子の表記法はあり、有名どころでは SMARTS というものがあります。SMARTS は SMILES を拡張したもので、部分構造の検索などで用いられます。SMILES と違って、 A や * などのワイルドカードを利用することができるのが特徴です。
DeepSMILES という、ニューラルネットワークに利用しやすい（と主張している）ものもありますが、少し使ってみた限りではそんなにメリットはないような気がします・・・
2019 年に Twitter おじさん Chemical VAE を作ったことで有名な Aspuru-Guzik 教授が、新たな化合物表記方として、SELFIES というものを考案し、報告しています。
この SELFIES、どんな表記をしても必ず分子として意味をなすという、とんでもない性質を有しています。論文中の validity が 100% なのはもはやギャグかと・・・ 論文には SELFIES の文法について書かれていますが、正直複雑で、私も理解し切れていません。ただ、SELFIES のライブラリは GitHub に公開されていて、SMILES に簡単に変換できるので、実務上は理解していなくても問題ないと思います。
インストールも pip install selfies でおしまいなので、簡単です。ただ、one-hot encoding してくれる関数が書いてありますが、実際にはないような・・・ （最も、one-hot encoding する関数は自作した方が良いと個人的には考えています。この話も、いずれしていこうと思います）
既に SELFIES を使った VAE も報告されており、これから広く使われてくるのではないかと思います。
https://www.biorxiv.org/content/10.1101/2020.05.23.112201v2.full
ちなみに、SELFIES は「自撮り」という意味です。「笑顔」で「自撮り」するとは、Guzik も茶目っ気がある人のようです。
Jupiyter Notebook 上での構造式描写
ケモインフォマティクスと他の情報分野との最大の差異は、構造式を扱うことであると思います。
有機化学にとって構造式は切っても切れない関係にあり、それはプログラミング上においても例外ではありません。
構造式の情報を確認するのにも、RDKit は有用なライブラリとなります。
以下、構造式の描写方法を具体例を挙げて示しますが、先行のブログに詳細は書かれているので、詳しくはそちらをご覧頂ければと思います。
Draw molecules as SVG in horizontal layout #Drawing #RDKit #memoiwatobipen.wordpress.com
import numpy as np import pandas as pd from rdkit import Chem from rdkit.Chem.Draw import IPythonConsole mol = Chem.MolFromSmiles("CC(C)CCc1ccc([C@@H](C)C(=O)O)cc1") mol
上記のコードのキモは、IPythonConsole を読み込むことです。一度読み込んでしまえば、mol クラスを構造式として表示してくれるようになるので、データがどのような構造なのか、確認がしやすくなります。RDKit を読み込むときは、一緒に書いておくことをお勧めします。
構造式を表示させるだけでなく、画像データとして出力させたい場合、rdMolDraw2D というクラスを用いると便利です。
import numpy as np import pandas as pd from rdkit import Chem from rdkit.Chem.Draw import IPythonConsole, rdMolDraw2D from IPython.display import SVG # 化合物の mol データを用意 mol = Chem.MolFromSmiles("CC(C)CCc1ccc([C@@H](C)C(=O)O)cc1") # rdMolDraw2D クラスを用意して構造式を表示させる用意をする drawer = rdMolDraw2D.MolDraw2DSVG(300, 300) drawer.DrawMolecule(mol) drawer.FinishDrawing() svg = drawer.GetDrawingText() # 構造式を svg 形式で出力 with open("structure.svg", 'w') as writer: writer.write(svg) # 構造式を表示 SVG(svg)
RDKit のクラスの挙動は結構分かりにくいので、上記の書き方を覚えるなり、自分で色々試したりして慣れていくのが結果的に早く身に付くと思います。
ここで重要なのは、MolDraw2DSVG で Jupyter Notebook 上で表示させる枠を用意し、mol データを入力するということです。svg 形式のデータはテキストと形式が近いらしいので、GetDrawingText で出力用のデータ形式を用意することで、一般的なテキスト形式のデータ出力と同じ要領で構造式のデータを出力することができます。ちなみに、経験上、jpeg 形式や png 形式よりも、svg 形式の方が word や powerpoint 上で構造式をよりきれいに扱うことができます。SVG(svg) はセルの最後に書かないと、構造式を描写してくれないので、少し注意が必要です。
ここから少し応用編です。
構造式中の特定の原子をハイライトしたい、ということも RDKit ならできてしまいます。
import numpy as np import pandas as pd from rdkit import Chem from rdkit.Chem.Draw import IPythonConsole, rdMolDraw2D from IPython.display import SVG # 化合物の mol データを用意 mol = Chem.MolFromSmiles("CC(C)CCc1ccc([C@@H](C)C(=O)O)cc1") # rdMolDraw2D クラスを用意して構造式を表示させる用意をする drawer = rdMolDraw2D.MolDraw2DSVG(1000, 1000) drawer.drawOptions().addAtomIndices = True drawer.drawOptions().addStereoAnnotation = True drawer.DrawMolecule(mol, highlightAtoms = [0, 2, 9]) drawer.FinishDrawing() svg = drawer.GetDrawingText() # 構造式を svg 形式で出力 with open("structure_mark.svg", 'w') as writer: writer.write(svg) # 構造式を表示 SVG(svg)
drawOption で色々な設定ができます。addAtomIndice を True にすることで、mol クラスの原子のインデックスを図に表示することができます。このインデックスは、mol ファイル固有のものなので、化学的に等価な原子でも異なる整数値をとります。Python なので、インデックスは 0 始まりであることに注意して下さい。addStereoAnnotation で、不斉炭素原子に R, S のタグをつけることができます。DrawMolecule の highlightAtoms で、インデックスのリストを指定することで、対応するインデックスの原子に色をつけることができます。色付けしたい原子が一か所でも、リストとして指定して下さい。
Atom を Bond に書き換えれば、結合のインデックス表示や色付けを行うこともできます。
