magattacaのブログ

日付以外誤報

MD計算の基本的な用語とパラメータをお勉強したメモ

前回の記事Google Colab上で分子動力学(MD)計算を行うMaking-it-rainで遊んでみました。計算したい系のパラメータをいれるだけなのでとても便利です。

・・・「でもそのパラメータが何かよくわからない!」「なんならMD計算が何かもよくわかってない!」

というわけで、パラメータの理解に必要そうなMD計算の基本的な用語をお勉強したメモです。せめてパラメータが何を言っているかわかりたい!

参考資料

参考にさせていただいた書籍はこちらです。

www.kyoritsu-pub.co.jp

いろいろとググった結果、本棚にあったこちらの書籍に戻ってきました。2009年出版の書籍ですが、知識ゼロから読めて「なんかこの単語見たことあるぞ!」って状態になれるので個人的にとてもオススメです。*1

また、以下のフリーでアクセスできる資料も参考にさせていただきました。*2

  1. 古明地 勇人, 上林 正巳, 長嶋 雲兵 生体分子の分子動力学シミュレーション(1)方法
    Journal of Chemical Software 2000(6)1-36
  2. 古明地 勇人, 田島 澄恵, 原口 誠, 高橋 伸幸, 上林 正巳, 長嶋 雲兵 生体分子の分子動力学シミュレーション (2) 応用
    Journal of Chemical Software 2001(7)1-28

1. MD計算の基本

1-1. MD計算って?

とりあえずWikipedia

分子動力学法(ぶんしどうりきがくほう、英: molecular dynamics、MD法)は、原子ならびに分子の物理的な動きのコンピューターシミュレーション手法である。(分子動力学法-Wikipedia

というわけで、分子の動的な性質がわかります。

タンパク質を対象としたMD計算は、1977年カープラス先生らの論文に端を発するそうです。 *3

www.nature.com

ウシ膵臓トリプシン阻害剤(Bovine Pancreatic Trypsin Inhibitor, BPTI)の9.2ps(ピコ秒)のシミュレーションで、これにより「タンパク質は動的な系である」という概念を提唱し、それまでの「タンパク質は固い」という概念を転換したそうです。新しい手法を適応してものの見方を変えてしまう研究、2013年のノーベル化学賞受賞につながっているのも納得ですね。*4

さて、この記事ではMD計算の用語について大体下図みたいに追っていきます。

f:id:magattaca:20220117204614p:plain

1-2. 動的な性質がなぜわかるの?

まず「どうやって動的な性質がわかるか?」という点ですが、「ある条件において、立体構造が時間に依存してどのように変化するか?」シミュレーションを行います。この構造変化の時系列のことをトラジェクトリ(軌跡, trajectory)とよびます。

f:id:magattaca:20220117204653p:plain

より具体的には、ある時点での「立体構造の状態(座標、速度)」をもとに、ニュートン運動方程式を解いて次の時点で立体構造の状態を導きます。この作業を繰り返していくことで各ステップの構造のスナップショット("静止画")が得られます。スナップショットを繋げば動画になるので、動的な状態を観察することができます。なお、ステップ間の時間間隔のことを時間刻み(time step)とよぶそうです。

1-3. 運動方程式で位置と速度を求める

MD計算ではニュートン運動方程式を解いて、次のステップの構造(座標、速度)を求めます。図のように「段階的に時間が進む」ことを時間発展(time evolution)とよび、「位置や速度を1ステップごとに求めて時間発展する」ことを時間積分(time integration)とよびます。

f:id:magattaca:20220117205038p:plain

上図のように、ある時刻の位置と速度は一つ前のステップの情報をもとに求まります。新しい位置と速度で構造を更新していくことで構造変化の時系列(トラジェクトリ)が得られます。

ただし、上記は2次の微少量(dvdt)を無視しているので、コンピュータによる数値計算では時間発展の間に誤差が蓄積して計算が破綻するそうです。誤差の蓄積を防ぐ時間積分法として「差分方程式から質点の位置と速度を求める方法」があり、Verlet法速度Verlet法蛙飛び(leap-frog)といった方法があります。

1-4. 分子力場

時間積分法で速度を求める際に各原子に働く力(F)を利用していますが、これは分子力場(Force field)から求めることができます。

力のベクトルはポテンシャルエネルギー(U)の勾配偏微分してマイナス(-∇U))で求まります。古典論に基づくMD計算では、原子の座標関数に代入することでポテンシャルエネルギーもとめます。分子のポテンシャルエネルギーを与える関数とその係数が分子力場です。

分子力場の例としてはAMBERCHARMMGROMOSといったものがあります。たとえばAMBER力場の関数形式は以下のような感じです(AMBER (分子動力学)-Wikipedia)。

f:id:magattaca:20220117205257p:plain

5つの項のうち、最初の3つは原子間の結合に関係する項(結合項)で、残り2つは結合の有無に関わらず関係する項(非結合項)です。

MD計算では、基本的に各原子が剛体球ばねで繋がれているとするモデルを前提としています。高校物理でよくみたやつですね!(格好良くいうと調和振動子?)

結合項のうち距離(r)や角度(θ)は調和ポテンシャルとなっていて、理想的な状態(平衡結合距離(r_0)、平衡結合角(θ_0))からズレると2乗に比例してエネルギーが大きくなります。また、二面角は回転ポテンシャルに基づいていて、「n : 二面角を作る4原子の対称性」「γ : 位相」です。

非結合項はいずれも距離に依存しますが、減衰の速いファンデルワールス・エネルギーと比較して、静電エネルギーは距離の逆1乗に比例するため減衰が遅く長距離の相互作用にも寄与します。

非結合項では(結合項と異なり)計算する原子のペア数が系の大きさに応じて膨大に増えるため、計算時間に大きく影響するという課題があります。そこで、計算を簡単にするため「一定の距離(カットオフ, cut-off)以内にある場合のみを考慮する」(カットオフ法)という方法が以前はよく用いられたそうです。ところが、上記の通り「静電エネルギーは長距離にも影響する」ので、一定の距離でバサッと切ってしまうと相互作用を無視してしまい、計算精度が低下してしまいます。現在では、これを改善するための効率的で精度の高い手法が色々と実装されているそうです。

1-5. 境界条件と(長距離)静電相互作用の計算

静電相互作用の計算の課題について述べたので、もう少し解決方法を眺めます。まず、前提知識として境界条件にふれておきます。

1-5-1. 境界条件

MD計算の計算系には孤立系周期系があります。孤立系では「計算系の外は真空」となっており、周期系ではずらっとくっついて並んだ「周期的な容器(セル)」を考えます。

孤立系で溶媒を明示的に扱う(explicit solvent)モデルでは、境界条件という「球状の容器の中にタンパク質やリガンド、水分子を閉じ込める方法」がよく用いられるそうです。この系では境界付近で水の挙動が不自然になってしまうという課題があります。

f:id:magattaca:20220117205701p:plain
上図は「タンパク質計算科学」第4章 図4-28および関連する記述をもとに作成

周期系では「周期的な箱」のなかに分子を詰めます。中心を基本セル、それ以外をイメージセルとよび、セルの形状は立方体直方体が一般的だそうです。

f:id:magattaca:20220117205732p:plain
上図は「タンパク質計算科学」第4章 図4-28および関連する記述をもとに作成

上図のようなイメージで、セルから外に飛び出た分子は反対側からセルに入ってきます*5。周期性を取りいれることで、「実験系よりも圧倒的に少ない原子数」を取り扱うだけで、「マクロな系の性質を求められる」というのが基本的な考え方のようです。周期的境界条件では、球境界条件の境界付近における不自然な挙動が解消できるという利点があります。

1-5-2. (長距離)静電相互作用の計算

先に見た通り、静電相互作用は考えるペアの組み合わせが多く、MD計算で計算時間がかかる部分で、なんとかして効率的に解きたい問題です。

基本的な解決方法の考え方は「近いところ」と「遠いところ」に領域を分けて、寄与の大きい前者は真面目に厳密計算し、少ない後者は近似計算することです。極端なものは「遠いところ=寄与ゼロ」とするもので、先のカットオフ法です。

もう少し真面目に計算しよう!ということで出てきたEwald法です。Ewald法は周期的境界条件周期性をうまく利用した高速計算方法です。うまいこと式を分割して収束を速くする感じです。具体的な式は書籍やWikipediaエバルトの方法-Wikipedia)をご参照ください。*6

f:id:magattaca:20220117205904p:plain
表式は「タンパク質計算科学」第4章 式(4.42)より引用

Ewald法をさらに高速化したものとして、Particle Mesh Ewald法(PME法)があります。こちらは波数空間の式(E_wave)を「電荷分布をグリッド上に内挿する形」に書き換え、高速フーリエ変換(FFT)することで高速化しています。

Ewald法は周期境界条件を利用したものでしたが、境界条件でも適用できる計算方法として多重極展開法があります。「遠いところ」にある電荷を持つ点群を多重極展開により近似する方法で、代表点との相互作用のみの計算に単純化することで、考慮する相互作用ペアを大幅に削減することができます。この手法を使った例に、セルマルチポール法(Cell Multipole Method, CMM)やファストマルチポール法(Fast Multipole Method, FMM)といったものがあるそうです。*7

以上が計算コストの高い静電相互作用についての工夫でした。

1-6. マクロな視点で(統計アンサンブル)

ここまでで、MD計算では運動方程式を解いて系の逐次変化をたどっていくこと、運動方程式の力はポテンシャルエネルギーから決まり、エネルギーを求めるために分子力場が使われることがわかりました。最後にもう少しマクロな視点でのお話にふれておきます。

多くの実験は定温定圧条件でおこなわれ、また、生体内で起こる現象も近似的には定温・定圧条件で起こります。これらと比較し生物学的な現象を解き明かす上で、MD計算結果をより自然(リアル)な系のシミュレーションとするために、さまざまな系の制御が行われます。

1-6-1. アンサンブル

「運動量や座標のある値の組」で指定される系の状態を微視的状態とよび、「多数の微視的状態の集まり」をアンサンブルとよびます。MD計算ではトラジェクトリがアンサンブルに対応します。

アンサンブルとしては以下があります。

f:id:magattaca:20220117211008p:plain
上図は「タンパク質計算科学」第4章 図4-22および関連する記述をもとに再構成して作成

先に見た通り、定温や定圧といった条件が意味のあるシミュレーションとして大事になってくるので、さまざまな温度制御圧力制御の工夫によりカノニカルアンサンブルNPTアンサンブルを得ることが目指されます。

1-6-2. 温度制御(カノニカルアンサンブル)

系の温度は「系を構成する原子の運動エネルギーの総和」で決まるので、原子の速度が制御できれば温度を制御できることになります。なかでも自然な速度分布に制御するには、統計力学的に妥当な条件を満たす必要があります。

いろいろな温度制御法が提案されていますが、共通して「仮想的な巨大熱浴系との相互作用を想定して系の温度を制御している」という特徴があるそうです。

f:id:magattaca:20220124232631p:plain
「タンパク質計算科学」4.4.2 の記述、Web公開の講義資料*8をもとに作成*9

温度制御法の例はこんな感じです。*10

f:id:magattaca:20220117211050p:plain

難しくてよくわからなかったので正しい情報は書籍や、自然科学研究機構 分子科学研究所の奥村先生が研究室HPで公開してくださっている資料(分子シミュレーション研究会会誌 アンサンブル 連載解説「分子動力学シミュレーションにおける温度・圧力制御」)などをご参照ください。*11

1-6-3. 圧力制御(NPTアンサンブル)

温度制御では速度が大事でしたが、圧力制御では系の体積を増減させることで系の圧力がコントロールされます。代表的な圧力制御の方法としてアンダーソン法があります。温度制御の能勢法(能勢ーフーバー法)と組み合わせた能勢-アンダーソン法では、NPTアンサンブルが実現できます。

f:id:magattaca:20220117211326p:plain
「タンパク質計算科学」4.4.5の記述を参考に作成

・・・詳細はよくわかりませんでした。

以上、アンサンブルについてでした。細かな理論については理解できませんでしたが、MD計算で得られるトラジェクトリ(アンサンブル)が意味のあるもの(熱力学的な平衡状態からのサンプリング)とするために、一定の条件を保つよう制御する工夫が行われている、ということがわかりました。

1-7. 古典MDと第一原理MD(と機械学習MD)

以上でMDの基本的な用語はおしまいです。ところで、ここまで見てきたものは古典分子動力学とよばれるもので、ニュートン運動方程式古典力学)に従い、経験的なポテンシャル(分子力場)を利用して原子間の相互作用を求めています。

これに対して、原子間相互作用を量子力学に基づいて計算する第一原理分子動力学というものもあるそうです。量子力学で電子の状態を求め、原子核、電子のすべての相互作用を考慮して原子どうしの相互作用を計算します。このような量子論に基づくMDは1985年に提案されたカー - パリネロ(Car - Parrinello, CP)に遡るそうです。

古典MDでも「スパコンを使って計算!」みたいな大変そうな感じなのに、時間ステップごとに第一原理計算をして力を求めるなんてめちゃくちゃ大変そうですね。

「計算が大変なところには機械学習がやってくる!」というわけで、この力の計算部分を機械学習に基づくポテンシャルニューラルネットワークポテンシャル)に置き換えて計算を高速化した機械学習分子動力学というのもあるそうです。*12

また、第一原理MDと厳密に等しい結果が得られる自己学習ハイブリッドモンテカルロ(SLMC)というのもあり、機械学習MDほど計算コストは下がりませんが、学習結果に依存せず高い計算精度が得られるという利点があるそうです。*13

2. MD計算の流れと設定

2-1. MD計算の流れ

MDの基本的な内容がだいたいわかったので、具体的な計算の流れをたどりたいと思います。だいたいこんな感じです。*14

f:id:magattaca:20220117212136p:plain

シミュレーションボックスの準備で「系の中和」と書いていますが、これは周期的境界条件でシミュレーションを行う際に扱う時に、電荷が中和されていないと静電相互作用の計算がうまくいかなくなることに関係しているそうです。

2-2. 力場の選択

分子力場の例としてAMBERCHARMMGROMOSといったものがあることをみました。*15

ただし、それぞれの力場の中でも、どのような種類の分子を対象に計算を行うかによって選択肢が変わってくるようです。例えば、Amberの Webサイトの説明では以下のようになっていました。

f:id:magattaca:20220117212214p:plain

2-3. 明示的な水分子モデル

水分子の配置で明示的非明示的とかきましたが、これは水も分子として扱うか、それとも連続体モデルのように扱うか、という選択です。

明示的な水分子としてはSPC*16TIP系(TIP3P、TIP4P、TIP5P)、OPC*17といったモデルがあり、Wikipedia水モデル項にも詳しく説明されています。

f:id:magattaca:20220117212252p:plain

TIP系の数字(3、4、5)は、それぞれ「いくつの点をつかって水分子をモデル化するか?」と関係します。単純に考えると、酸素原子1つと水素原子2つ3点モデル(TIP3P、SPC)となりそうですが、4点モデルではH-O-H角垂直二等分線上の点Mもさらに考慮し(TIP4P)、5点モデルでは酸素原子ローンペア方向の2点Lを考慮する(TIP5P)モデルとなっています。

以上が、MDの基本でふれていなかった設定項目についてです。

3. Making-it-rainのパラメータ振り返り

MDでやりたいことと、計算の流れがだいたいつかめたので、前回の記事で遊んだMaking-it-rainチュートリアルのパラメータを振り返ってみましょう!

3-1. シミュレーションボックス(系)の準備

まずは「①シミュレーションボックス(系)の準備」に相当する段階で、Amberトポロジーファイルを作成している部分です。

f:id:magattaca:20220117212644p:plain

タンパク質対象なので力場(Force_field)として「ff19SB」を選択しているのは先にAmberのマニュアルで見た通りですね。また、水分子のモデル(Water_type)としては明示的な「OPC」を選択しています。

選択項目の中にありませんが、コードの中身を見るとMaking-it-rainは静電相互作用の計算として「PME(Particle Mesh Ewald)」を採用しているようです。ですので、境界条件周期的境界条件となっていて、「Size_box」でセルの大きさを設定することになります。またEwald法と関連するので系を中和するための設定(Ions, Concentration)も入っています。

3-2. エネルギー最小化と平衡化

続いて「②エネルギー最小化」と「③平衡化」に相当する段階です。パラメータ設定はEquilibration protocolでまとめられています。

f:id:magattaca:20220117212706p:plain

平衡化シミュレーションの時間を設定する項目で、「Time」の他に「Integration_timestep」があります。MD計算で運動方程式を解くときは、微少時間の時間積分を考えるという話だったので、微少時間に相当するパラメータとなりそうです。(上図では0.2nsのシミュレーション中、2fs間隔で積分計算を行う感じ?)

また、計算で得られるアンサンブルに関係するパラメータとして温度圧力があります。Making-it-rainはNPTアンサンブルのみとなっています。

MD計算で欲しいのは「一定の時間間隔ごとの構造のスナップショット(トラジェクトリ)」ですが、「Write_the_trajectory」がその時間間隔のパラメータになっています。運動方程式の計算はフェムト秒間隔ですが、すべての計算結果を保存したらデータ量が膨大になるので、ピコ秒オーダーで出力する感じでしょうか?

3-3. シミュレーション本番

最後に「④シミュレーション本番」に相当する部分です。ここのパラメータは平衡化の部分とほぼ一緒です(Making-it-rain特有の時間設定パラメータ(Stride)が増えているくらい)。

f:id:magattaca:20220117212732p:plain

ところで、Making-it-rainはNPTアンサンブルと書きましたが、「どのように温度・圧力制御を行うか?」の設定の選択肢はありませんでした。

コードの中身を見ると、MonteCarloBarostat(pressure, temperature)LangevinIntegrator(temperature, friction, dt)というそれっぽい箇所がみられます。Making-it-rainがMD計算に使用しているOpenMMのAPIを参照してみましょう。

MonteCarloBarostatは「モンテカルロアルゴリズムで周期的なセルのサイズを調整し、一定の圧力の効果をシミュレーションする」とあります。これが圧調整器(バロスタット)みたいですね。ただし「温度を維持する機能はないので、LangevinIntegratorAndersenThermostatといったものと合わせて使うように」とあります。

また、LangevinIntegratorは、「ランジュバン動力学を使って系のシミュレーションを行う積分器」とあります。「ランジュバン動力学 - Wikipedia」によると、「サーモスタットを用いた時のように温度を制御できるため、正準集団(カノニカルアンサンブル)を近似できる。」とのことですので、こちらが一定の温度に保つための温度調整器(サーモスタット)の役割を果たしているようです。

ところで、先に見た温度制御の例で、「サーモスタットはだいたい速度をコントロールしている」「ランジュバン方程式類似で摩擦力を使う」といった内容を見ました。LangevinIntegratorの引数に摩擦係数(friction coefficient)があるのを見ると、確かに摩擦力で速度制御して温度をコントロールしてそうだな、という印象です。*18

以上、Making-it-rainのパラメータ振り返りでした。MD計算の考え方や手順をぼんやり眺めた後だと、パラメータの意味が少しわかった気がします!

4. Making-it-rainトラジェクトリ解析の振り返り

Making-it-rainには基本的なトラジェクトリ解析もありましたが、これらもさっぱり意味がわかりませんでした。ついでに振り返っておきましょう。

一般的な解析手法について素晴らしい資料がありました!理化学研究所 計算科学研究センターのeラーニング「計算生命科学の基礎6 生体系分子シミュレーションの新展開」です!

横浜市立大学 池口満徳先生(研究室HP)の講義動画(YouTube)とスライド資料(SlideShare)を公開してくださっています。うれしい!*19

4-1. RMSDとRMSF

まずは、根平均二乗変(Root Mean Square Deviation, RMSD)と根平均二乗揺らぎ(Root Mean Square Fluctuation, RMSF)です。どっちもズレを評価している感じ。

RMSDについてMaking-it-rainの結果と説明を並べてみるとこんな感じ。

f:id:magattaca:20220117213521p:plain

上の例だと、PDBから取得した初期構造(参照構造)からのズレを可視化できていることになりそうです。また、トラジェクトリ前半でRMSDは増加傾向にあるものの、後半では一定の範囲に落ち着いているので、確かに収束具合がわかります。

RMSFについてはこんな感じ。

f:id:magattaca:20220117213809p:plain

RMSDとRMSFの用語が似すぎていて違いがよくわかっていませんでしたが、参照構造との比較か平均構造との比較かという違いだったんですね。横軸が違う(time, residue)意味もやっとわかりました。

4-2. 回転半径

回転半径はこんな感じ。

f:id:magattaca:20220117213918p:plain

実験で得られるデータと関係しているというのは知りませんでした。上記の例は時間も短く、フォールディングの検証でもないのであまりこの指標には合わないのかもしれません。

4-3. 主成分解析(PCA)

PCAはこんな感じ。

f:id:magattaca:20220117214102p:plain

タンパク質を相手にPCAを行う意味がよくわかっていなかったですが、揺らぎの解析だったんですね。平面プロット(PC1, PC2)で近くにあるCαは、一緒に動いている残基の可能性が高いということでしょうか??

4-4. 相関行列

Making-it-rainではPearson's Cross Correlation(CC)となっていましたが、出力結果を見た感じ相関行列(Correlation Matrix)や動的相関図(Dynamical Cross Correlation Map, DCCM)と呼ばれているものと同じでしょうか?

f:id:magattaca:20220117214124p:plain

この解析にはいくつか問題点も指摘されているそうです。式の通り2つの残基の平行な協同運動しか観測できず、また中心付近のコアドメイン内の協同運動を観測することもできません。また、揺らぎの相関なので収束が遅く、再現性が低いそうです。

詳しい分子内協同運動を調べるには、相関行列よりも先の主成分解析やEssential Dynamics(ED)といったものの方が良いそうです。

以上、トラジェクトリ解析の振り返りでした。解釈は難しいですが、とりあえず何を目的とした解析なのかを把握できたのでヨシ!

5. おわりに

以上、「今回はMD計算のパラメータがわからない!っていうか、そもそもMD計算て何をやっているんだ??」ということで、MD計算の基本的な内容とパラメータ等の用語のお勉強をしてみました。

ついでに、前回の記事で遊んだMaking-it-rainのチュートリアルを振り返って、パラメータや操作の意味・解析方法について復習してみました。

「計算がうまくいっているのか?」「妥当な結果なのか?」といったところまではわかりませんでしたが、とりあえず何をやろうとしているかぼんやりわかったのでヨシ!

書籍の参考部分が多くなってしまったので、最後にもう一度おすすめしときます。

タンパク質計算科学 ―基礎と創薬への応用―(CD-ROM付き) / 神谷 成敏 肥後 順一 福西 快文 中村 春木 著 | 共立出版

今回も色々と間違いが多そうなのでご指摘いただけると嬉しいです。(初心者向けの本もお薦めいただけると嬉しいです。)

ではでは!

*1:いろいろな資料を参考にさせていただいていますが、ベースは書籍を参考にしているので先に宣伝しときます。PDBファイルの読み方を含めてタンパク質をコンピューターで扱う方法がざっとさらえて素敵な書籍です。基本的な内容は変わっていないと思うので、まだまだ現役の書籍とは思いつつ、アップデートを期待している書籍でもあります。(もちろん私は全くの部外者で、アフィリエイトでもないので(安心して?)リンクをクリックしてしてください)。なお本記事中の図表の作成にあたって参考にさせていただいた書籍の対応箇所はキャプションに記載しています。

*2:こっちはさらに前の2000年、2001年の資料です。日本語の素晴らしい資料が20年経ってもフリーでアクセスできるのは最高ですね!

*3:McCammon, J. A., Gelin, B. R. & Karplus, M., Dynamics of folded proteins Nature 267, 585-590 (1977).

*4:なお、この記念すべきシミュレーションの動画はカープラス研究室HP(Karplus group HP)でみることができます。

*5:SF映画でドアを開けたら元の部屋につながっている感じ??

*6:フーリエ級数という言葉が出てきた瞬間に理解を諦めました。

*7:湯どうふさんのYouTube動画が簡潔でわかりやすいです。多重極展開法のYouTube動画はこちらで、セルマルチポール法のYoutTube動画はこちらです。

*8:東京大学大学院農学生命科学研究科 アグリバイオインフォマティクス教育研究ユニット 講義「分子モデリングと分子シミュレーション」 2019年5月23日の資料「分子動力学法とモンテカルロ法」の中の図表と記述を参考にしました

*9:イメージはいらすとやさんの図から

*10:下図の特徴に書いている部分がわかりにくいのは私が全く理解していないからです。ごめんなさい。

*11:J-STAGEでもフリーで読めます。

*12:たとえばAENET(The Atomic Energy Network)というオープンソースソフトウェアがあるそう(生体分子というよりもマテリアル寄り?)。

*13:原子力機構システム計算科学センター 永井佑紀先生が講演資料をGitHubで公開してくださっていて面白かったのでご紹介しておきます(スライドPDFのURL(GitHub)

*14:参考:@Ag_smith先生 Qiita記事 MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜 より

*15:Wikipediaの「力場(化学))」項に他にもいろいろな力場の例が載っています。

*16:単純点電荷, simple point-chargeの略?

*17:Oはoptimal?参考: J. Phys. Chem. Lett. 2014(5)3863–3871

*18:ランジュバン方程式・動力学が全くわかっていないので間違っていたらごめんなさい。

*19:以下の説明では池口先生の資料と、冒頭にあげた古明地先生らの資料の記述を参考に引用させていただいています。

Google CouldでMD計算をお試し!ーMaking it rainー

タンパク質がぐにゃぐにゃ動く分子動力学(MD)計算かわいいですよね。下はD.E. Shaw Research("DESRES")のツイートで、彼らの開発したスパコンアントン - wikipedia)で計算されたものらしいです*1

「お高いパソコンがないとMDで遊べないのかな?」と思ってましたが、Twitterで「低分子ならGoogle Colabでも結構できるよ!」「Making it rainが参考になるよ!」って教えていただきました。@ylogtさん、ありがとうございました!

というわけで、とりあえずMakig-it-rainをそのまま実行した結果のご報告です。*2

f:id:magattaca:20220108220841p:plain
CloudからMDの雨を降らせ!

1. Making it rain?

1-1. 概要

Making-it-rainは計算リソースが必要な分子動力学(MD)計算を、誰でもお手軽にできるようにしよう!というプロジェクトです。GoogleアカウントがあればOK!

Maiking-it-rainではOpenMMというオープンソースのMDシミュレーションライブラリをGoogle Colab上で実行することができます。ノートブックに従って順番に実行することで、①MD計算の準備、②実行、③結果の解析まで行えます。また、Google Driveとの連携されているので、計算結果やノートブックの保存も容易です。

論文は以下。

GitHubレポジトリは以下。

github.com

MD計算をハンズオンで学習するための教育や、資金の少ない研究室でもマイクロ秒スケールのシミュレーションができるようにするため開発されたそうです。ありがたいですね!

1-2. 論文の紹介(Zenodo公開版)

査読済みの論文はオープンアクセスではありませんが、プレプリントChemRxivZenodoで公開してくださっています*3。計算の中身やコードの詳細というよりも、プロジェクトのモチベーションや大枠の構成の説明と、実行例の紹介といった内容になっています。

Making it rainの基本的構成は、① 計算パワーの必要なMD計算をGoogle ColabGPUで実行し、② 入手力ファイルや設定ファイル、結果といったデータ保存・管理Google Driveで行うものとなっています。

Google Colab、Google DriveともにGoogleアカウントがあれば無料で利用できますが、有償版ではより長時間高容量の利用ができます。*4

f:id:magattaca:20220108221331p:plain

計算時間の目安として、リゾチーム(lysozyme)を題材に1μsのシミュレーションを行った例が論文に書かれています。より高性能のGPUが使える有償版(Colab Pro)では、平均して約230ns/dayのシミュレーションが可能で、1マイクロ秒には6日あればOKだったそうです。

Colabには時間制限(有償版でも1日)があるため、Making it rainは長時間の計算を分割して実行できるように作られています。論文の例では、1マイクロ秒全体のトラジェクトリ20n秒ごと、50個のストライドに分割しています。各ストライドの終わりに、その時点での系の状態を保存したファイル(.xml)をOpenMMが出力してくれるので、出力ファイルを元に再度計算を実行することで連続したシミュレーションができる、とのことです。例では5回実行することで計1μsの結果を得ています。

1-3. ノートブックの種類

Making it rainではGoogle Colabで実行可能なノートブックが複数用意されており、それぞれ計算内容計算方法が異なっているようです。論文では3つですが、現在はさらに3つアップデートされ、計6つのノートブックが利用可能です。

いずれもGitHubレポジトリのトップ(README)から利用可能で、こんな感じ。

f:id:magattaca:20220108221409p:plain

入力ファイルを作成するソフトウェアや利用する力場によってノートブックが分かれているようですが、MD計算を勉強したことがないので全く違いがわかりません。すみません。

よくわからないけどタンパク質モデリング(AlphaFold2+MD)やリガンドとの相互作用(Protein-ligand simulation)といったトピックもおさえてあるのが魅力的!!

以上、Making it rainの概要でした。

2. 使ってみた

2-1. AMBERノートブックの流れ

ざっとMaking it rainの目的が把握できたので早速遊んでみましょう!今回はノートブックのうち「AMBER」をそのままデフォルトで実行してみました。*5

AMBERノートブックの流れは以下のような感じでした。

f:id:magattaca:20220108221439p:plain

「ノートブックを順番に操作するとどうなるか?」ビデオ形式のデモまで著者のPablo Arantesさんが用意してくださっています。ノートブックを開くとすぐに見られる(2分ちょっと)ので、流れを把握するのに便利です。親切!!*6

2-2. やってみた 〜MD計算の実行〜

「AMBERノートブック」のデフォルトの題材は「ニワトリ卵白由来のリゾチームの水中における挙動」のシミュレーションでした。MD計算実行部分がどんな感じだったかご紹介します。

まず、初期構造の準備です。PDBからPDB ID: 1AKIを取得しています。手持ちの構造が使いたければGoogle Driveに入れておけば参照できます。

f:id:magattaca:20220108221602p:plain

次にシミュレーションの系の設定です。水で満たされたボックスタンパク質が一つ配置されます。さらにイオン(NaCl)を加えて系を中和しておきます。

f:id:magattaca:20220108221622p:plain

上がトポロジーファイルのパラメータ設定です。力場(force field。上ではff19SB)や水の種類、ボックスサイズ、イオンといった項目があります。これを元にMDシミュレーションのための系のパラメータがつくられます。

@Ag_smith先生の以下の記事に詳しいですが、AmberではtleapというモジュールでPDBファイルからトポロジーファイルを作成するそうです。実際にMaking it rainのコードを表示するとtleapを呼んでいます。

qiita.com

トポロジーファイル作成後、「シミュレーションの系がどのようにできたか?」、ノートブック上で3Dで可視化して確認できます。MD計算は時間が長いのできちんと準備できているか心配になりそうですが、ビジュアルで確認できると安心できていいですね!

f:id:magattaca:20220108221739p:plain

系が準備できたら平衡化(equilibration)を行います。MD計算ではシミュレーションの間、系の温度や圧力が一定の値を保つように計算を行うそうです。そのためにシミュレーションの本番前に短い計算で系をならしておくようです。*7

設定はこんな感じ。

f:id:magattaca:20220108221800p:plain

先の@Ag_smith先生の記事から関連する解説を引用させていただきます。*8

ここからは水溶液中に存在するタンパク質を物理演算で動かしていくことになります。しかし、これまでの操作で作り出した初期座標では水溶液の様子を再現しているとは言えません。・・・(中略)・・・そこで、まずは系のエネルギーをある程度小さくしてから(=エネルギー最小化)、各原子に(温度に対するBoltzmann-distributionに従う範囲で)ランダムに速度を与え、徐々に系全体を設定温度まで引き上げていく(=平衡化)という予備動作を行うことで水溶液の状態に近づけていくという操作を行います。

①エネルギー最小化、②平衡化の順に実施するそうです。確かに、設定の最初にminimizationの文字が見えますね。

設定が終わったら平衡化シミュレーションの計算を実行します。NPTアンサンブル物質の量(N)、圧力(P)、温度(T)が一定になるようにシミュレーションを行うそうです。「気温と大気圧に開放されているフラスコを用いた実験室条件に最も密接に対応している」らしい。。。(Wikipedia - 分子動力学法

f:id:magattaca:20220108221828p:plain

平衡化の長さは0.2nsとしていて、平均して約45ns/day(〜 0.03 ns/min)のスピードでシミュレーションできたようなので、計算時間は6分程度でした。"Time Remaining" で残り時間の目安がわかるのが親切ですね。

最後にシミュレーション本番(Production)を実行します。シミュレーションの長さはNumber_of_stridesStride_Timeで指定し、2つの掛け算が得られる長さとなります。

例)100n秒シミュレーション = (ストライド10) x(各ストライドの長さ 10n秒

f:id:magattaca:20220108221853p:plain

今回は全てデフォルトのまま2nsのシミュレーションを行いました。以下を実行し、計算時間は1時間ちょっとという感じでした。平衡化の10倍の長さなので、だいたい同じ感じのスピードで計算できた感じです。

f:id:magattaca:20220108221927p:plain

MD計算が終わったら、シミュレーションの最後にトラジェクトリ(trajectory)ファイルを一つにまとめる操作を行います。

系に含まれる各原子の3次元座標の時間変化のことをトラジェクトリー(trajectory、軌跡)と呼ぶそうです*9。パラメータ設定で「トラジェクトリファイルに書き出す頻度」を10psとしたので、2nsのシミュレーションでは200個のファイルができることになります。これを一つのファイルにまとめます。

f:id:magattaca:20220108221946p:plain

まとめると0.2GB弱でした。2nsで0.2GBなので1μsだと100GBくらいになりそうです。・・・でっかい。

2-3. やってみた 〜結果の解析〜

計算の実行とファイルの統合が終わったので、結果の解析部分を見てみましょう!

Making it rainのノートブックに用意されているので、引き続き実行結果をみていきます。

まずはトラジェクトリファイルを動画で見ます。「Load, view and check the trajectory」というセルを実行するだけ!

こんな感じ(以下はキャプチャのGIFです)

f:id:magattaca:20220108222020g:plain

ピクピクしてる!!かわいい!

「動画もいいけど、もう少し定量に結果を見たいよね!」というわけで、以下のような解析が用意されています。それぞれセルを実行するだけなので結果のプロットを貼っていきます。

まずは各アミノ酸残基のCα原子(中心炭素原子)についてこんな解析ができます。

f:id:magattaca:20220108222110p:plain

RMSD(根平均二乗変位, Root Mean Square Deviation)、②radius of gyration(慣性半径)、③RMSF(根平均二乗揺らぎ, Root Mean Square Fluctuation)です。①と②は時間に依存した変化分布、③は残基ごとの値についてそれぞれプロットされています。

まったくどう解釈していいかわからない!!

①RMSDが0.5nsくらいまで増加傾向にあるのをみると、production前の平衡化が足りていなかったのでしょうか?

③RMSFでは末端以外にも残基40~5070前後で少し大きくなっていますね。動画(GIF)で左側青色のループ部分に相当しそうなので、動きが大きそうな部位を数値化できてる感じでしょうか??。

残りの解析結果も貼っておきます。こんな感じ。

f:id:magattaca:20220108222131p:plain

2D RMSDPCAPearson's Cross Correlation (CC, ピアソン相関)。。。よくわかりませんでした。

結果の解析部分は以上です。

3. おわりに

以上、今回はGoogle Colab上でMD計算を実行できるMaking it rainと、実際に実行してみた結果のご紹介でした。

Making it rainではオープンソースのMD計算ソフトopenMMを利用しています。こちらはTeachOpenCADDでも使われているそうです。@iwatobipen先生が記事「Make mdtools for openmm」で紹介してくださっていました。*10

専門外の分野でも、オープンソースのツールがあったり使い方の解説があると「ちょっと遊んでみようかなー」となるので嬉しいですね!Making it rainはさらに低コストで高コストなMD計算を実施するフレームワークを提供してくれていて、とてもありがたいです。

残念ながら専門外すぎて、どういう設定で計算して、結果をどう評価していいのか全くわかりませんでした。不勉強ですみません。。。

ところで「Make it rain」という単語、直訳すると「雨を降らせる」ですが、Googleで検索すると強面のお兄さんたちと紙幣が舞う画像がでてきます。スラングで「大金をばら撒く、札束をまき散らす」といった意味があるとか。。。

クラウド(雲)コンピューティングから高価なデータの雨を降らせよう(たくさんお金がないとできないような計算をいっぱいやろう)なんて、意味が込められてたら面白いですね。

今回も色々と間違いが多そうです。ご指摘いただけると嬉しいです。

ではでは

追記:日本語化ノートブック

Making it rainのノートブックを日本語化したものをGitHubに入れておきました。ほとんどDeepL翻訳そのままで、動作確認もまだできていませんがご参考までにどうぞ。一応Colabでも動かせるはずです。

github.com

*1:引用したツイートと関連するMD計算の結果はDESRESのHPで公開されています(CC-4.0)。Molecular Dynamics Simulations Related to SARS-CoV-2

*2:下の図はGitHubのFigureといらすとやさんのイラストを利用させていただいています。

*3:ChemRxiv上のライセンスはCC BY NC ND 4.0です。ZenodoはCC BY 4.0になっていますがサイトの設定か論文本文かどちらかわからなかったです。この記事ではZenodo公開版を参照させていただいています

*4:無償版/有償版の費用や性能は変わる可能性があるので適宜検索してください。Google Driveの有償版としてはGoogle Oneの情報を記載しました。

*5:何もわからないので一番上にあるやつをクリックしただけ。。。

*6:デモ動画はYouTubeにありますが「限定公開」設定なのでここにはリンクを貼らないでおきます。ノートブックをひらけばすぐわかります。

*7:全然違う話題ですけど、平衡化を最初に行うのベイズ推定のマルコフ連鎖モンテカルロ法(MCMC)使ったパラメータ推定の計算みたいですね。

*8:ただし引用記事はこの部分でAmberではなくGROMACSというソフトに切り替わっています

*9:参考:古明地 勇人, 上林 正巳, 長嶋 雲兵 生体分子の分子動力学シミュレーション(1)方法 J. Chem. Software, 2000(6)1–36

*10:TeachOpenCADD、インストールで力尽きて新しいトークトリアルの中身をまだ見られていないです。すみません。

2Dも3Dも!Panel-Chemistryでまとめて描画しよう!

前回の記事Jupyter notebook構造式エディタを使う方法を調べました。

もっといい方法がありました!その名もPanel-Chemistry!!

2Dの構造式(JSME Editor)だけでなく、3Dの描画(NGL Viewer, Py3DMol)もカバーしてくれています。こりゃ便利!

1. Panel-Chemistry?

とりあえずデモが格好いいですよ!*1

f:id:magattaca:20220103233305g:plain

Panel-Chemistryは、化学分野の探索的データ解析・ビジュアリゼーションをサポートするためのプロジェクトです。Pythonでデータ可視化するためのオープンソースライブラリであるHoloViz - Panelに化学ドメインを組み込んで構築されているのでPanel-Chemistryと名付けられているみたいです。

Panel-Chemistryがサポートしている可視化ツールは3つです。

  • JSME Editor (化学構造式描画&編集、2D)
  • NGL Viewer (タンパク質、DNA/RNA etc.の描画、3D)
  • py3DMol (タンパク質、低分子化合物の描画、3D)

それぞれ対応可能なファイルフォーマットが異なります。2D、3Dの生物学・化学描画をカバーしてくれているのはありがたいですね!

ライセンスはApche-2.0GitHubは以下です。

github.com

pipcondaのどちらでもインストールできます。

# pipの場合
pip install panel-chemistry

# condaの場合
conda install -c conda-forge panel-chemistry

# もちろんmambaでもOK
mamba install -c conda-forge panel-chemistry

とりあえず使用感を知りたいと言う方は、GitHubレポジトリREADMEBinderへのリンクが用意されているのでそちらをご利用ください。

2. 使ってみよう!

GitHubの例(examples)ほぼそのままで恐縮ですが、どんな感じに使えるかご紹介します*2

サポートする可視化ツールそれぞれについて、①Jupyter notebook上で描画、あるいは②サーバーモードでアプリケーションとして利用する方法の2つの使い道があります。

2-1. JSMEの場合

2-1-1. Jupyter notebook上で使う

まずはJSMEをJupyter notebookで使う場合です。JSMEについては以下をご参照ください。

panelpanel_chemistryから使いたいツール(今回はJSMEEditor)をインポートしておきます。

import panel as pn
from panel_chemistry.widgets import JSMEEditor

次にpanelextension()jsmeを読み込みます*3。以下のpn.extension()実行前に、上記の通りpabel_chemistryを読み込んでおく必要があることに注意です。

pn.extension("jsme", sizing_mode="stretch_width")

panel可視化スペースのレイアウトはcolumn()を使って設定できます。JSMEEditorを以下のようにレイアウトに追加するだけでもう使えます。

editor = JSMEEditor(height=500)
pn.Column(editor)

f:id:magattaca:20220103234037p:plain

JSMEが起動しました!簡単!

さらに便利な機能は、JSME起動時に同時に指定した構造式が描画された状態で起動できることです。例えばSMILESで構造を指定するなら以下のようにすればOKです。

smiles="N[C@@H](CCC(=O)N[C@@H](CS)C(=O)NCC(=O)O)C(=O)O"
editor = JSMEEditor(value=smiles, height=500, format="smiles")

pn.Column(editor, editor.param.value)

f:id:magattaca:20220103234055p:plain

上図の通りグルタチオンが描かれた状態でJSMEが起動しました。またpn.Column()editorに加えてeditor.param.valueもレイアウトに追加しているので、上図下部のようにvalueに格納されたSMILESも表示されています。

valueは他にもmolsdf形式も使えます。例えば以下のようにRDKitで作成したMol形式の座標を読み込むこともできます。

from rdkit import Chem

toluene = Chem.MolFromSmiles("Cc1ccccc1")
toluene_mb = Chem.MolToMolBlock(toluene)
print(toluene_mb)
"""
    
         RDKit          2D
    
      7  7  0  0  0  0  0  0  0  0999 V2000
        3.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
       -0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
       -1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
       -0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
      1  2  1  0
      2  3  2  0
      3  4  1  0
      4  5  2  0
      5  6  1  0
      6  7  2  0
      7  2  1  0
    M  END
"""  
# mol形式の座標を読み込んで起動する
editor = JSMEEditor(value=toluene_mb, height=500, format="mol")
pn.Column(editor)

f:id:magattaca:20220103234203p:plain

トルエンのように単純な構造だと有り難みは少ないですが、より複雑な分子を描画する時、SMILESから構築するとツールに依存した描画スタイルになってしまいます。複雑な分子を好みの向き・スタイルで描きやすいと言う点で、座標情報のあるMol形式が使えるのは良いですね!

2-1-2. アプリケーションとして使う

つづいてアプリケーションとして使う場合です。ビューのレイアウトといった諸々の設定をJupyter notebook上に書き込んだあと、コマンドラインpanel serveコマンドを使うことでアプリケーションとして起動することができます。

この方法ではPanelFastLstTemplateテンプレートを利用してレイアウトを作成します。最後に.servabl()アノテーションをつけておくことでコマンドラインから起動できるようになります。

Panelのアプリケーションサーバーの仕組みはドキュメントのServer Deploymentをご参照ください。デフォルトのポートは5006なので、起動後http://localhost:5006にブラウザでアクセスすればアプリを利用できます。

設定例はこんな感じ。

# アプリケーションの見た目の色調設定
accent="#00A170"

# メイン描画部位の設定
values = pn.Param(editor, parameters=["value", "jme", "smiles", "mol", "mol3000", "sdf"], widgets={
    "value": {"type": pn.widgets.TextAreaInput, "height": 200},
    "jme": {"type": pn.widgets.TextAreaInput, "height": 200},
    "smile": {"type": pn.widgets.TextAreaInput, "height": 200},
    "mol": {"type": pn.widgets.TextAreaInput, "height": 200},
    "mol3000": {"type": pn.widgets.TextAreaInput, "height": 200},
    "sdf": {"type": pn.widgets.TextAreaInput, "height": 200},
    })

# サイドバーの設定
settings = pn.Param(editor,
    parameters=["height", "width", "sizing_mode", "subscriptions", "format", "options", "guicolor"],
    widgets={
        "options": {"height": 300},
    })

# テンプレートを利用したレイアウトの作成 + servable()
pn.template.FastListTemplate(
    site="Panel Chemistry", 
    title="JSME Editor", 
    sidebar=[settings],
    main=[editor, values],
    accent_base_color=accent, header_background=accent,
).servable();

上記設定が完了したら、コンソール(私はmacなのでターミナル)でpanel serve ノートブック名を打ち込めばサーバーアプリケーションとして起動します*4。「PanelChemistry_JSME.ipynb」というノートブックならこんな感じ。

panel serve PanelChemistry_JSME.ipynb

もちろんJupyter notebook上で最初に「!」をつけてコマンドを実行してもOKです*5

実行後「サーバーを起動したぜ!」「アプリケーションはここだぜ!( Bokeh app running at: http://localhost:5006/PanelChemistry_JSME )」みたいなことを言ってくるのでこのURLをブラウザにコピペすればアプリが使えます。

こんな感じ!

f:id:magattaca:20220103234249p:plain

緑を基調にしてるのが素敵ですね!

JSMEを使う場合は以上です。

2-2. NGL viewerの場合

2-2-1. Jupyter notebook上で使う

つぎにNGL viewerの場合です。NGL Viewerについては以下をご参照ください。

  • AS Rose, AR Bradley, Y Valasatava, JM Duarte, A Prlić and PW Rose. Web-based molecular graphics for large complexes. ACM Proceedings of the 21st International Conference on Web3D Technology (Web3D '16): 185-186, 2016. doi:10.1145/2945292.2945324.

  • AS Rose and PW Hildebrand. NGL Viewer: a web application for molecular visualization. Nucl Acids Res (1 July 2015) 43 (W1): W576-W579 first published online April 29, 2015. doi:10.1093/nar/gkv402.

やることは大体一緒です。まずはライブラリのインポートとngl_viewerのpanelへの読み込みです。

import panel as pn 
from panel_chemistry.pane import NGLViewer
from panel_chemistry.pane.ngl_viewer import EXTENSIONS
pn.extension("ngl_viewer", sizing_mode="stretch_width")

次いで、ビューワーの設定です。

描画する目的の構造はobjectパラメーターで設定し、①url like、②non-url like、③blobの3種類から選べます。①と②はPDBに対応するもので、①なら「'rcsb://1NKT'」、②なら「'1NKT'」のようにします。②はPDB IDだけで、自動的に「rcsb://」を保管してURLにして構造を持ってきてくれます。

blobの時は合わせて拡張子(extension)パラメーターを設定する必要があります。pdbcif等色々と使えるのでGitHubなどをご参照ください。

今回はGitHubのexample通りやってもRCSBからPDBファイルを持ってきてくれなかったので、別途ダウンロードした「1nkt.pdb」を使っています*6

# ビューワーと読み込む構造の設定
viewer = NGLViewer(object="1nkt.pdb",background="#F7F7F7", min_height=700, sizing_mode="stretch_both")

# パラメーター設定
settings = pn.Param(
    viewer,
    parameters=["object","extension","representation","color_scheme","custom_color_scheme","effect",],
    name="⚙️ Settings"
)

# 入力ファイルの設定
file_input = pn.widgets.FileInput(accept=','.join('.' + s for s in EXTENSIONS[1:]))

def filename_callback(target, event):
    target.extension = event.new.split('.')[1]

def value_callback(target, event):
    target.object = event.new.decode('utf-8')

file_input.link(viewer, callbacks={'value': value_callback, 'filename': filename_callback})

# レイアウトの設定
header = pn.widgets.StaticText(value='<b>{0}</b>'.format("&#128190; File Input"))
file_column = pn.layout.Column(header, file_input)

layout = pn.Param(
    viewer,
    parameters=["sizing_mode", "width", "height", "background"],
    name="&#128208; Layout"
)
pn.Row(
    viewer,
    pn.WidgetBox(settings, layout, width=300, sizing_mode="fixed",),
)

f:id:magattaca:20220103234330p:plain

できました!

2-2-2. アプリケーションとして使う

つづいてアプリケーションとして使う場合です。コマンドラインで起動する際にノートブック単位で指定するので、JSMEを起動したノートブックとは別のもの(PanelChemistry_NGL.ipynb)を作っています。

アプリケーションにおける設定は上のJupyter notebook上で使う場合と同様です。viewersettingsflie_columnlayoutといった変数の部分なのでそのまま引き継ぎます。

こんな感じで起動します。

# アプリケーションの見た目の色調設定
accent="#D2386C"

# テンプレートを利用したレイアウトの作成 + servable()
pn.template.FastListTemplate(
    site="Panel Chemistry", 
    title="NGLViewer", 
    sidebar=[file_column, settings, layout],
    main=[viewer],
    accent_base_color=accent, header_background=accent,
).servable();

あとはコンソールでコマンドを打ち込むだけです。

panel serve PanelChemistry_NGL.ipynb

ポートもデフォルト「5006」で変わらずなので「http://localhost:5006/PanelChemistry_NGL 」にブラウザでアクセスすればアプリを利用できます。

f:id:magattaca:20220103234359p:plain

できました!起動時にファイルを指定しているので構造を読み込んでいますが、左上のファイルを選択から新しく構造を読み込むこともできます。

2-3. Py3DMolの場合

2-2-1. Jupyter notebook上で使う

最後にpy3Dmolの場合です。py3Dmolについては以下をご参照ください。

また、化学の新しいカタチさんの記事「py3Dmolを使って化学構造をJupyter上で美しく表示する」で詳細に解説してくださっています。

今までと同様にまずはライブラリを読み込みます。

import py3Dmol
import panel as pn

from panel_chemistry.pane import Py3DMol
pn.extension()

基本的な使い方は以下です。1行目py3DMol.view()queryを設定するとデータベースから構造をダウンロードして描画してくれます。下記ではmmtf(macromolecular transmission format)形式でPDB ID 1nktの構造をとってきています。

p = py3Dmol.view(query="mmtf:1nkt")
p.setStyle({"cartoon": {"color": "spectrum"}})
pviewer=Py3DMol(p, height=400, sizing_mode="stretch_width", name="Basic")
pviewer

f:id:magattaca:20220103234436p:plain

できました!これ以外にもグリッド表示インタラクティブなプロット(低分子の動くやつ)の例がGitHub上にありますが割愛します。

2-3-2. アプリケーションとして使う

つづいてアプリケーションとして使う場合です。ここでも新しいノートブックを(PanelChemistry_py3Dmol.ipynb)を作っています。

今までと同様、設定を書き込んだあと起動するだけです。

# 設定
background = pn.widgets.ColorPicker(value="#e6f6ff", name="Background")
style=pn.widgets.RadioButtonGroup(value="sphere", options=["line", "cross", "stick", "sphere"], name="Style", button_type="success")

# アプリケーションの見た目の色調設定
accent = "#0072B5"

# テンプレートを利用したレイアウトの作成 + servable()
pn.template.FastListTemplate(
    site="Panel Chemistry", 
    title="Py3DMol Pane", 
    sidebar=[background, style],
    main=[pviewer],
    header_background=accent, accent_base_color=accent
).servable();

あとはコンソールでコマンドを打ち込むだけです。

panel serve PanelChemistry_py3Dmol.ipynb

ポートもデフォルト「5006」で変わらずなので「http://localhost:5006/PanelChemistry_py3Dmol 」にブラウザでアクセスすればアプリを利用できます。

f:id:magattaca:20220103234519p:plain

できました!

2-4. 3つのアプリを同時に

最後に、3つのアプリを同時に使う方法です。

それぞれのノートブックでservable()にした状態で、コマンドラインで以下の通り3つ並べて起動すればOKです。

panel serve PanelChemistry_JSME.ipynb PanelChemistry_NGL.ipynb PanelChemistry_py3Dmol.ipynb

この状態で「http://localhost:5006/ 」にアクセスすると以下のようになります。

f:id:magattaca:20220103234536p:plain

真ん中に並んでいるブロックがそれぞれJSMENGL viewerpy3Dmolそれぞれのアプリに対応するので、クリックすれば開けます。もちろん別々のタブで同時に開くことも可能です。

3. おわりに

以上、今回は描画ツールPanel-Chemistryのご紹介でした。

化学構造式を描画・編集できるJSMEと、3Dの描画が得意なNGL Viewerpy3Dmolを並べて扱えるので便利ですね!一度コピー&ペーストしてノートブックを作っておけば、気軽にアプリとして起動できそうです。

まあ私はサーバーとかアプリとか全くわかってないんですけどね!

備忘録がわりのメモをはさんでいるのでGitHubのexampleよりわかりづらくなっているかもしれません。間違いがあったらご指摘いただけると嬉しいです。

ではでは!

*1:GitHubのteaser用gifです。大きいのでサイズを小さくさせていただきました。

*2:あれこれ余計なコメントをつけていますが、自分のための備忘録ですのでご容赦ください。

*3:Bokehとか言うのと関係ありそうですがよくわからなかったです。

*4:蛇足ですが、作業ディレクトリとcondaで入れているなら仮想環境のアクティベートも忘れずに!

*5:jupyter notebook上でサーバーを起動するコマンドを実行した場合、アプリケーション利用中そのセルが実行中になって、次の他のセルの実行ができなくなることにご注意ください

*6:RCSB PDB : 1NKT

Jupyter Notebookでも構造式を編集したい! 〜JSMEとipywidgetsでMolを扱う話〜

RDKitとJupyter notebookの組み合わせとても便利ですよね!「プログラミングなんてよくわからないぜ!」っていうへっぽこにとって結果が逐一見られるのは最高です!

ですが、遊んでいて一つ困ることが・・・「RDKitの構造式を書き変えるの面倒!!ちょっと変えたいだけなのに。。。」

というわけで今回はJupyter Notebook上で構造式を編集する方法を試してみたいと思います。具体的にはRDKitMolオブジェクトからMolファイルで座標を取り出して、構造式を編集後、新しい座標から再度Molオブジェクトを作り直します。

f:id:magattaca:20220101142250p:plain

1. JSMEとipywidgetを使おう

1-1. 先例 〜SMILESを取り出す方法〜

調べてみると同様のニーズはあるみたいで、優秀な方々が検討結果を共有してくださっていました。今回は、lithium0003さんのGitHubレポジトリ(JSME_ipywidget)を参考にさせていただくことにしました。*1

こちらの方法によると構造式エディタとしてJSMEを用い、Jupyter notebook上で実効するためipywidgetsを使うとのことでした。

1-2. JSMEって?

JSMEはフリーで利用可能な構造式エディタで、2013年に下記の論文で報告されています。*2

jcheminf.biomedcentral.com

元々Javaで書かれていたJMEというエディタを、JavaScriptにアップデートしたもののようでBSDライセンスで配布されています。

Webページはこちら(→ JSME Molecule Editor)です。「Download the latest JSME distribution as zip file」をクリックすると最新版をダウロードできます。今回は少し古いバージョンを使います。

f:id:magattaca:20220101135442p:plain

今回使うバージョンは2017-02-26です。最新版ではうまく動かなかったので参考例に近いバージョンにしました。

f:id:magattaca:20220101135505p:plain

ダウンロードすると「JSME_2017-02-16ディレクトリ」が手に入りますが、このうち「jsmeディレクトリ」が今回必要となる部分です。これだけ作業ディレクトリに移しておきます。

f:id:magattaca:20220101135528p:plain

なお、JSMEの見た目はこんな感じです。以前こちらのブログでも紹介したZINC20データベースでの利用例をもってきました。

f:id:magattaca:20220101135548p:plain

ChemDrawMarvinといった商用ソフトほどの機能はなさそうですが、十分な機能があり、何よりフリーで公開してくださっていることに感謝です。

1-3. 先例でコードのお勉強

lithium003さんの公開してくださっている方法はGitHubレポジトリのノートブック(smiles_draw.ipynb)の通りです。このノートブックでは「JSMEエディタをjupyter notebook上で使用し、描いた構造式のSMILESを取り出す」ということが行われています。

こんな感じ(READMEより)

f:id:magattaca:20220101135654p:plain

でもさっぱりコードの中身がわからない!!Pythonも分からないのにJavaScriptも必要なの??

ということでGoogleで検索しながら順番に見ていきました。

全部で4つのセルからなります。

1-3-1. セル1

まず最初にipywidgetsの準備です。ipywidgetsはJupyter notebook上でインタラクティブなUIを手軽に作れるライブラリだそうです*3。ドキュメントはこちら(→ ipywidgetsドキュメンテーション

とりあえずセル1のコードを引用します。ライブラリをインポートした後、SmilesEditorというクラスを作成しています。

from traitlets import Unicode, Bool, validate, TraitError
from ipywidgets import DOMWidget, register

@register
class SmilesEditor(DOMWidget):
    _view_name = Unicode('SmilesView').tag(sync=True)
    _view_module = Unicode('smiles_widget').tag(sync=True)
    _view_module_version = Unicode('0.1.0').tag(sync=True)

    # Attributes
    value = Unicode('', help="SMILES value").tag(sync=True)

ここで行われているのはipywidgetsドキュメンテーションの「Low Level Widget TutorialWidget skeleton項」 と「Building a Custom Widget - Email widgetMaking the widget stateful項」を参考にすると様子がわかりそうです。

1行目のtraitletsはJupyter notebookで機能しているライブラリだそうで、Pythonの(動的に決まる)「クラスの属性の型をきちんと定めて、さらに細かいチェック機能を簡単に呼び出せる」ようにできるそうです*4

2行目でipywidgetsライブラリからインポートしています。DOMWidgetDOMDocument Object Modelの略だと思います。DOMはJavaScriptの「表示されたWebサイトを動的に書き換えることができる」という特徴を支えている仕組みで、「HTMLをJavaScriptで操作することが出来る」そうです*5

またregisterはビューが動的にアップデートするのに関わる関数のようです。*6

このセルの目的はJSMEを使えるようにするための箱(SmilesEditor)を準備することのようです。JavaScriptベースのJSMEを使うためのDOMWidgetのクラスSmilesEditorをつくっているんだと思います。たぶん。。。

1-3-2. セル2

セル2ではJSMEを使うための設定が書き込まれています。マジックコマンド%%javascriptでJupyter notebook上でJavaScriptが実行できるようにした上で処理が書かれています。

セル2を引用します。

%%javascript
require.undef('smiles_widget');
require(['jsme/jsme.nocache.js'])
define('smiles_widget', ["@jupyter-widgets/base"], function(widgets) {
    var SmilesView = widgets.DOMWidgetView.extend({
        // Render the view.
        render: function() {
            this.smiles_input = document.createElement('div');
            this.smiles_input.id = "jsme_container";
            this.smiles_display = document.createElement('div');
            this.smiles_display.textContent="SMILES : ";
            this.el.appendChild(this.smiles_display);
            this.el.appendChild(this.smiles_input);

            function myFunc(callback){
                let jsmeApplet = new JSApplet.JSME("jsme_container", "480px", "480px");
                jsmeApplet.setCallBack("AfterStructureModified", callback);
            };
            setTimeout(myFunc, 500, this.smilesChanged.bind(this));
        },

        smilesChanged: function(jsmeEvent) {
            console.log(this, jsmeEvent)
            let jsme = jsmeEvent.src;
            let smiles = jsme.smiles();
            this.smiles_display.textContent="SMILES : "+smiles;
            this.model.set('value', smiles);
            this.model.save_changes();
        },
    });
    return {
        SmilesView: SmilesView
    };
});

冒頭2行のrequireは「一般的にモジュール化されたJavaScriptファイルを読み込む」ために用いられるものだそうです*7

1行目のrequire.undef()はモジュールの定義を外すためのもので、内部の状態をリセットするための関数のようです*8。一度リセットしてから、2行目でJSME本体を呼びにいっている感じでしょうか??

2行目でJSMEを読み込むパスの通り、JSMEはノートブックと同じディレクトリに置かれた「jsmeディレクトリ」内の「jsme.nocache.jsファイル」を利用しています。requireはモジュール化されたJavaScriptファイルをよぶという説明通りですね。

つづいて3行目以降では具体的なSmilesViewの中身が書かれています。

thisというのは「JavaScriptに最初から用意されている特別な変数のこと」で、「呼び出した場所や方法によってその中身が変化する」という特徴があるそうです*9。よくわかりませんが、変数thisに色々格納して、JSME構造式描画ビューワーsmiles_input)やSMILES表示smiles_display)の機能を設定している感じです。変数を設定した後appendChildメソッドを使って要素を実際に追加していく感じ*10

準備ができたら、真ん中のfunctionから始まるブロックでJSEMのアプレットを読み込んでいます。AfterStructureModifiedというコールバックを使っているので、構造が編集されると変更後の構造を読み込んで反映する仕組みがつくられている感じです。・・・たぶん。

一番最後のブロックでは、構造式エディタに書かれた構造式からSMILESをとり出す仕組みが書かれているようです。letでは始まる2行はJSMEからSMILESの情報を取り出して、変数に格納している感じです*11jsmeEvent.src.smiles()はJSME特有の箇所でJSME API ドキュメントなどに記載があります。

とりだしたSMILES情報をもとに、「this.smiles_display.textContent="SMILES : "+smiles;」ではインタラクティブな編集に応答してSMILESを表示する仕組みが、「this.model.set('value', smiles);」ではmodelvalueというプロパティにSMILESを格納しています。

以上、セル2でJavaScriptのアプリ、JSMEを利用するための設定と書かれた構造からSMILESを取り出す仕組みの設定までが終わったようです。

1-3-3. セル3

セル3は、セル1とセル2で設定が完了したSmilesEditorを起動しています。

コードはシンプルです。

smiles = SmilesEditor()
smiles

f:id:magattaca:20220101135853p:plain

参考例のREADMEと同じものが立ち上がりました!

構造式を書き込むとEditor上部の「SMILES:」欄に対応するSMILES表記が反映されます。

1-3-4. セル4

セル4は構造式に対応するSMILESをJupyter notebookのセルに取り出すだけです。

コードはこちら

smiles.value
#  'OCc1ccccc1'

セル2でみたように、valueプロパティにSMILESを格納していたので.vlaueとすることでSMILESが取り出せました!

以上がlithium003さんが公開してくださっている方法です。

2. この記事で試す方針

コードを辿るだけで長くなってしまいましたが、今回私が行いたいのは「編集した構造式の新しい座標をベースにMolオブジェクトを作り直す」ことです。このためには「JSMEからSMILESではなくMolファイル形式の情報として編集結果を取得」すれば良さそうです。

具体的には以下を行います。

  1. RDKitのMolオブジェクトをMolFileで出力する
  2. JSMEをnotebook上で起動する
  3. MolFileを読み込んで構造式を編集
  4. 編集後の座標(Molファイル形式)を取得
  5. RDKit Molオブジェクトを作り直す(MolFromMolBlock

一時的にMolFileで出力するのが格好悪いですが、書き直すよりは楽かな???ということで。

3. RDKitで構造式を編集する正攻法

お試しの前に、RDKitで構造式を編集する真っ当な方法をご紹介しておきます。以下の日本語解説記事がとても参考になります。

sishida21.github.io

RWMolオブジェクトに変更して、②編集した後、③再度Molオブジェクトに戻すそうです。

例えばトルエンベンジルアルコールにしようとするとこういう感じです。

from rdkit import Chem
from rdkit.Chem import Draw

# トルエンのMolオブジェクト
toluene = Chem.MolFromSmiles("Cc1ccccc1")

# 編集前後でindexを保つためにAtom Mapping
for i, atom in enumerate(toluene.GetAtoms(), start=1):
        atom.SetAtomMapNum(i)

Draw.MolToImage(toluene)

f:id:magattaca:20220101140045p:plain

# RWMolオブジェクトに変換
rw_toluene = Chem.RWMol(toluene)

# 酸素原子(Chem.Atom(8))を炭素原子(C:1)に追加して結合をつくる(AddBond)
from_idx = rw_toluene.AddAtom(Chem.Atom(8))
to_idx = [atom.GetIdx() for atom in rw_toluene.GetAtoms() if atom.GetAtomMapNum() == 1][0]
rw_toluene.GetAtomWithIdx(from_idx).SetAtomMapNum(8)
rw_toluene.AddBond(from_idx, to_idx, Chem.BondType.SINGLE)

# 編集後にMolオブジェクトに戻す
benzylalcohol = rw_toluene.GetMol()
Chem.SanitizeMol(benzylalcohol)

Draw.MolToImage(benzylalcohol)

f:id:magattaca:20220101140101p:plain

できました!けど面倒。。。分子ごとにatom indexを確認しないといけないのも一手間です。

4. JSMEとくみあわせてRDKit Molを編集しよう

では、JSMEと組み合わせて楽になるか???試してみましょう。

4-1. RDKitのMolオブジェクトから座標抽出

まずは、RDKitのMolオブジェクトから構造式の座標Molファイル形式temp.mol)で出力しておきます。

from rdkit import Chem

# トルエンを作ってMolFileで出力する
toluene = Chem.MolFromSmiles("Cc1ccccc1")
Chem.MolToMolFile(toluene, 'temp.mol')

4-2. JSMEの設定(MolFile抽出)

JSMEを利用する箇所はほとんどlithium003さんのコードを利用させていただきます。ただし、ここで欲しいのは編集後のMol形式の座標なので、SMILESではなくMolFileを取り出す様に変更します。

目的にあわせて名前もSmilesEditorからMolEditorにしてみました。

準備のセル1はクラス名称等をかえていますが基本そのままです。

from traitlets import Unicode, Bool, validate, TraitError
from ipywidgets import DOMWidget, register

@register
class MolEditor(DOMWidget):
    _view_name = Unicode('MolView').tag(sync=True)
    _view_module = Unicode('mol_widget').tag(sync=True)
    _view_module_version = Unicode('0.1.0').tag(sync=True)

    # Attributes
    value = Unicode('', help="Mol value").tag(sync=True)

JavaScriptのセル2では、よりシンプルにするためSMILES確認用ビューワーを削除しています。最後のブロックでMolFileを取り出すように変更しています。

%%javascript
require.undef('mol_widget');
require(['jsme/jsme.nocache.js'])
define('mol_widget', ["@jupyter-widgets/base"], function(widgets) {
    var MolView = widgets.DOMWidgetView.extend({
        // Render the view.
        render: function() {
            this.mol_input = document.createElement('div');
            this.mol_input.id = "jsme_container";
            this.el.appendChild(this.mol_input);

            function myFunc(callback){
                let jsmeApplet = new JSApplet.JSME("jsme_container", "480px", "480px");
                jsmeApplet.setCallBack("AfterStructureModified", callback);
            };
            setTimeout(myFunc, 500, this.molChanged.bind(this));
        },

        molChanged: function(jsmeEvent) {
            console.log(this, jsmeEvent)
            let jsme = jsmeEvent.src;
            let mol_data = jsme.molFile();
            this.model.set('value', mol_data);
            this.model.save_changes();
        },
    });
    return {
        MolView: MolView
    };
});

準備完了!

4-3. JSMEの実行と構造編集

JSMEは立ち上がるでしょうか?

MolEditorクラスをインスタンス化して実行します。

mol_editor = MolEditor()
mol_editor

f:id:magattaca:20220101140253p:plain

立ち上がりました!

右上の上下三角アイコンをクリックすると、構造式を種々の形式で読み書きできるプルダウンが開きます。

ここから「Paste Mol or SDF or SMILES」を選択すると新しい画面が開きます。「ファイルを選択」という箇所を選択するとファイルの読み込みができるので、先に出力しておいたMolファイル(temp.mol)を選択して取り込みます。「Accept」で描画スペースに構造式が反映されるので、あとは好きに編集しましょう!

f:id:magattaca:20220101140356p:plain

適当に書き加えてみました。編集後の構造式の座標はvalueプロパティから無事取り出せるでしょうか???

print(mol_editor.value)

"""
    Cc1ccc(C(C)O)nc1
    JME 2017-02-26 Sat Jan 01 12:59:19 GMT+900 2022
    
     10 10  0  0  0  0  0  0  0  0999 V2000
        5.6001    1.2124    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        4.2000    1.2124    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        3.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        2.1000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        1.4000    1.2124    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        2.1000    2.4249    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        3.5000    2.4249    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
        6.3001    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
        6.3001    2.4249    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    1.2124    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
      1  2  1  0  0  0  0
      2  3  2  0  0  0  0
      3  4  1  0  0  0  0
      4  5  2  0  0  0  0
      5  6  1  0  0  0  0
      6  7  2  0  0  0  0
      7  2  1  0  0  0  0
      1  8  1  0  0  0  0
      1  9  1  0  0  0  0
      5 10  1  0  0  0  0
    M  END
"""    

できました!

最後の仕上げにこの情報からRDKit Molオブジェクトを再度作り直します。

# 編集後の構造式からMolオブジェクトを作る
modified_mol = Chem.MolFromMolBlock(mol_editor.value)

Draw.MolToImage(modified_mol)

f:id:magattaca:20220101140618p:plain

再構築できました!

JSMEを使った構造式編集は、Atom indexとか考えなくても直感的に操作できるのでやりやすいですね!

5. 試したけど上手くいかなかったこと

やりたかったことは以上ですが、他にも試して上手くいかなかったことがいくつかあります。改善方法をご存知の方がいらっしゃったらご教示いただければ幸いです。

5-1. 失敗例1:JSMEバージョン依存

JSMEは継続的に開発されており、最新版は2021-07-13です。こちらを使って上記コードを試してみましたが上手くいきませんでした。

具体的にはJupyter notebok上にJSMEエディタ画面が立ち上がりませんでした。なので、今回は古いバージョンの2017-02-26を利用しています。

5-2. 失敗例2:Molファイルのクリップボードからのペースト

記事ではRDKitのMolオブジェクトを一度Molファイルとして出力し、JSMEにファイル読み込みという操作をしています。

「余計なファイルを増やすのは嫌だ」ということで、ファイル出力せずにMol情報だけコピー&ペーストしようとしました。

具体的にはpyperclipというライブラリを試しました。コード内でクリップボードコピーしたりペーストしたりできるそうです*12

pipでもcondaでもインストールできるそうなのでmambaしました*13

mamba install -c conda-forge pyperclip

RDKitのMolオブジェクトからクリップボードへのコピーまでは以下の通りできました。

import pyperclip

m = Chem.MolFromSmiles('CC')
mb = Chem.MolToMolBlock(m)

pyperclip.copy(mb)

このあと「Cntrol + V」すると、以下の通りペーストできます。

     RDKit          2D

  2  1  0  0  0  0  0  0  0  0999 V2000
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.2990    0.7500    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0
M  END

問題はJSMEでペーストできなかったことです。Jupyter notebook上で起動したJSMEの「Paste Mol or SDF or SMILES」で枠内で「Cntrol + V」しましたが、全く反映されず書き込めませんでした。JSMEのテストページでは問題なく貼り付けと構造式の描画ができたので、私のコード(or Jupyter notebook上での挙動)に問題があるようです。

5-3. 失敗例3:JSME起動時に構造反映

最後に、全く方法がわからなかったのが「JSMEを起動する際に、編集したい元の構造式が書かれた状態で起動する」という方法です。

「構造式の読み込み作業が面倒だなー」と思ったので、RDKitのMolオブジェクトからMolBlockで取り出した情報を、JSME起動時に書き込みたかったのですが、ipywidgetもJavaScriptもよくわからないので無理でした。

6. おわりに

以上、今回は「Jupyter notebook上でも構造式を編集したい!」ということで、先行例を参考にJSMEとの連携を試してみました。

途中でふれた通りRDKitの構造編集はノンプログラマーの私にはハードルが高いです。「直感的に構造式を編集したい!」という点でエディターを使えるととても便利な気がします。

商用ソフトはお高そうなので個人では手を出しづらかったのですが、フリーでつかえるJSMEと組み合わせられることがわかってよかったです。開発・公開してくださっている方々に感謝です!

結局、コードの中身をあまり理解しないままコピー&ペーストしてしまいました。今回も間違いがたくさんありそうです。ご指摘、また改善方法等ご教示いただけると嬉しいです。

ではでは!!今年もよろしくお願いいたします!!

Img2Molでピクトグラムを構造式に!

今年も残りわずかですね。2021年も色々なことがありました。

素人がお馬鹿をさらしているだけのブログですが、コメントいただけたり、twitterでリプライをいただいたりと皆様に暖かく接していただけて本当に感謝しています。

というわけで、今年最後のお遊びです!Img2Molをつかったピクトグラムの構造式化!!!

Img2Mol

以前こちらのブログでImg2Molというソフトウェアをご紹介しました。化学構造式の画像からOCRでSMILESを取り出すというものです。

magattaca.hatenablog.com

ナノプシャン

ところで有機化学の分野には、化学構造式でお絵描きしてしまおう!面白い構造式の化合物を作ろう!という素敵な研究があります。

有名なところでは分子の小人、ナノプシャンがありますね!

f:id:magattaca:20211231234833p:plain
Wikipediaより

ピクトグラム

2021年いろいろなことがありましたが東京オリンピックが印象的でしたね。賛否両論あったものの、この大変な状況で素晴らしいパフォーマンスを発揮されたスポーツ選手の皆様に勇気をもらった方も多かったのではないかと思います。

さて、オリンピックのでは開会式のパフォーマンス、ピクトグラムを再現!というのが話題になりましたね。

ピクトグラムを化学構造式に?

お馬鹿な私はこう思いました。

ピクトグラムをImg2Molに投げ込んでめば、ピクトグラムを化学構造式で表現する方法がわかるのでは???」

Img2Molでオリジナルなナノプシャンをつくってやろうではありませんか!!

というわけで雑に作ったピクトグラムを投げ込んでみました!

はしれ

まずは走る人です。こんなピクトグラム作ってみました

f:id:magattaca:20211231234906p:plain

Img2Molによる解釈はこちら

f:id:magattaca:20211231234920p:plain

足が表現されてる!宇宙人感!

ソーシャルディスタンス

距離を保つのが大事ですよね

f:id:magattaca:20211231234938p:plain

構造式では?

f:id:magattaca:20211231234950p:plain

Brはどこから??

レディオ体操

開催直前に色々ありましたが、小林賢太郎さんの所属するユニット ラーメンズのコント「バニー部」は独特の魅力がありますよね。

レディオ体操第一!

f:id:magattaca:20211231235015p:plain

構造式では?

f:id:magattaca:20211231235045p:plain

はじめて窒素原子がはいった!

おしまい

というわけでなんの中身もない記事でした。

こんな感じのくだらないブログですが、来年もこりずにお付き合いいただければ幸いです。

おしまい

TeachOpenCADD 2021をインストールした話

大幅にアップデートされたTeachOpenCADD 2021がリリースされました!

とりあえずインストールした!!

ちょっとつまずいた!!

のでインストール方法のメモです。

1. TeachOpenCADDって?

TeachOpenCADDは「計算機支援ドラッグデザイン(CADD, computer-aided drug design)の基本的なタスクをどのように実行すれば良いか?」学習・教育するためのプラットフォームで、ケモインフォマティクス構造バイオインフォマティクスについて学ぶ事ができる教材です。

Volkamer研究室が主体となって開発されており、オープンソースCC BY 4.0)で公開されています。

2019年に公開・論文発表されており、以前こちらのブログでも取り上げさせていただきました。

jcheminf.biomedcentral.com

magattaca.hatenablog.com

今回、大幅にアップデートされたTeachOpenCADD 2021が公開されましたので改めてご紹介します。

2. TeachOpenCADD 2021

TeachOpenCADDはCADDのトピックをそれぞれ取り上げたトークトリアルtalk + tutorial = talktorial)と呼ばれる教材で構成されています。

2019年版ではT1 ~ T10の10回のトークトリアルでしたが、2021年版はT22までと大幅に内容がアップデートされています。

ChemRxivに論文(preprint)が公開されています。

  • Sydow D, Rodríguez-Guerra J, Kimber TB, Schaller D, Taylor CJ, Chen Y, et al. TeachOpenCADD 2021: Open Source and FAIR Python Pipelines to Assist in Structural Bioinformatics and Cheminformatics Research. ChemRxiv

Figure 1. をご覧いただけば扱われている内容がつかめると思います。

f:id:magattaca:20211231000639p:plain

すごい!!

もう少し概要を知りたい方はVolkamer Labのブログ(TeachOpenCADD 2021 release is out!)や、RDKit UGM2021でのご発表などもご参照いただくと良いかも。。。

  • UGM2021のライトニングトーク(関連発表は9分50秒くらいからスタート(10分弱))

www.youtube.com

3. アップデートはコンテンツだけじゃない!Webサイトも充実!

盛りだくさんにコンテンツがアップデートされたTeachOpenCADDですが、2021のアップデートはさらに豪華です!

格好いいWebサイトまでできました!!

f:id:magattaca:20211231000729p:plain

素敵だと思ったのは「Talktorials by collection」。関連するトークトリアルをグループ分けしたものが確認できます。

f:id:magattaca:20211231000803p:plain

構造バイオインフォマティクスってインスタ映えしそうですね。

4. まだまだ増えてるコンテンツ

TeachOpenCADD 2021はトークトリアル T22までのアップデートと書きましたが、実はコンテンツは継続的に開発されていてWebサイト上では最新コンテンツも確認できます。

Volkamer Labのブログ(TeachOpenCADD Kinase edition is out!, 2021.12.14)によると、キナーゼに焦点を当てたスペシャルエディションができたそう。

Webサイトでチェックだ!

f:id:magattaca:20211231000857p:plain

新しいトークトリアル(T23 ~ T28)がアップデート済みです!ますます充実していますね!

5. オンライン上で遊ぼう!(インストール不要)

TeachOpenCADDで遊ぶにはオンライン上で実行する方法と、ローカル環境にインストールする方法の2通りがあります。

オンラインはBinderで実行する形式で、セットアップ不要で手間要らずです。TeachOpenCADDのGitHub上にある「lauch / binder」ボタンをクリックするだけ!

BinderのURLは「 https://mybinder.org/v2/gh/volkamerlab/TeachOpenCADD/master 」なので、毎回GitHubに行くのが面倒な方はブックマークしておくと便利かもしれません。

こんな感じでBinderが起動して環境を準備してくれます。

f:id:magattaca:20211231000931p:plain

準備にちょっと時間がかかります(5~10分くらい?)が、ワンクリックで使えるお手軽さは最高です。

準備できるとJupyter Labが立ち上がります。こんなの。

f:id:magattaca:20211231000959p:plain

上図のように「teachopencaddディレクトリ → talktorialsディレクトリ → 見たいトークトリアルのディレクトリ(例:T001_query_chembl)」とたどって目的のノートブック(talktorkal.ipynb)を見つけます。

開くとWebサイト上で閲覧できるものと同じものが見られます。こんな感じ。

f:id:magattaca:20211231001017p:plain

Webサイトとの違いは実際にノートブックのセルを実行できる事です。また編集できるので気になるコードを修正して結果の違いを見比べたり、試したい解析を追加したりと、いろんな遊び方ができそうです。

6. ローカル環境で動かしたい!〜インストール①〜

毎回Binderを待つのが嫌!オフラインでも試したい!という方はconda環境があればローカルにインストールすることもできます。

WebサイトのRun locallyページにインストール方法が書かれています。

Condaパッケージからのインストールは下記で、mambaを使った方法が書かれています。

f:id:magattaca:20211231001104p:plain

mambaは「C++で実装された高速なcondaパッケージマネージャー」だそうです*1。conda環境を使われているなら、mambaを追加インストールしておけば諸々のパッケージインストールが速くなるみたいです。(mambaオフィシャルドキュメント)

mambaは以下でインストールできます。私はMacなのでターミナルで以下を実行しました。

conda install mamba -c conda-forge

インストールできたら今まで「conda install ~~~」としてきたのを「mamba install ~~~」とするだけでOKです。

なお、TeachOpenCADDのインストール方法では新しい仮想環境をつくる工程が省略されているようです。通常通りconda create -n <仮想環境名>とすればOKですが、仮想環境の作成はmambaでもできるそうです。

やってみた!

# 仮想環境の作成
mamba create -n teachopencadd

できた!

f:id:magattaca:20211231001240p:plain

上図の通りできた仮想環境(teachopencadd)のアクティベート、ディアクティベートの方法はcondaと変わらず同じコマンドです。

なお、mambaで作成した仮想環境もconda info -eで確認できる環境一覧に出てきますし、削除したくなったら以下でOKです。

conda remove -n teachopencadd --all

仮想環境ができたら、アクティベートしてTeachOpen CADDをインストールします。(Webサイトの順番だとできない気がするので逆にしてます)

# 仮想環境のアクティベート
conda activate teachopencadd

# TeachOpenCADDのインストール 
mamba install teachopencadd

インストールできたら以下のコマンドを実行しましょう!

teachopencadd start .

カレントディレクト.)に新しいワークスペースが作成され、ターミナルにこんなのが表示されます。

f:id:magattaca:20211231001318p:plain

ロゴがかわいいですね。使えるトークトリアルの一覧とJupyterLabの起動コマンドが書かれています。

jupyter lab teachopencadd-talktorials

指示通り上記コマンドを実行すればBinderで確認したのと同じページがローカルで立ち上がります。

f:id:magattaca:20211231001349p:plain

無事インストールできました!

なお、ワークスペースを作成する場所を指定したい時はteachopencadd start .コマンドのピリオドの部分にお好みのパスを入れてあげてください。

7. もっとローカル環境で動かしたい! 〜インストール②〜

「リリースを待ちきれない!」という方や、「自分で新しいトークトリアルを作るぜ!」という方は最新の開発版をインストールされるのも良いかもしれません。

一応ご紹介しておきます。

f:id:magattaca:20211231001432p:plain

以下のコマンドでGitHubのTeachOpenCADDレポジトリに用意されているテスト環境用のtest_env.ymlを利用して仮想環境teachopencaddが作成されます。

mamba env create -f https://raw.githubusercontent.com/volkamerlab/TeachOpenCADD/master/devtools/test_env.yml

次に、以下でGitHubのレポジトリが作業ディレクトリ内にZIPファイルとしてダウンロードされます。

wget https://github.com/volkamerlab/teachopencadd/archive/master.zip -O teachopencadd.zip

unzipで解凍すればOKですが、以下では新しくDocumentsディレクトリを作ってその中に解凍しています。

# 作業ディレクトリ(test)の下にDocumentsディレクトリを作成
mkdir -p ./Documents

# Documentsディレクトリ内にファイルを解凍
unzip teachopencadd.zip -d ./Documents

これでteachopencadd-masterディレクトリができました。この中のtalktorialsディレクトリに各トークトリアルのディレクトリが入っています。

f:id:magattaca:20211231001509p:plain

あとはtalktorialsディレクトリに移動してJupyterLabを立ち上げればおしまいです!仮想環境(teachopencadd)のアクティベートも忘れずに!

# talktorialsディレクトリに移動
cd ./Documents/teachopencadd-master/teachopencadd/talktorials

# 仮想環境のアクティベート
conda activate teachopencadd

# JupyterLabの立ち上げ
jupyter lab

mamba installしたときと同じJupyterLab画面が立ち上がりました!

先に見たインストール方法①ではワークスペース作成コマンドの実行があり、talktorialsフォルダができました。開発環境版を使う場合はワークスペース作成の代わりに、GitHubレポジトリを丸ごと持ってきてtalktorialsディレクトリをそのまま利用する、という違いがあるようです。

8. おしまい

以上、新しくパワーアップしたTeachOpenCADDについて簡単な紹介とインストール方法のご紹介でした。

コンテンツの大幅なアップデートだけでなく、Webサイトもとても見やすくなっていて素晴らしいですね!

Webサイトに書かれているインストール方法そのままだと、どうも私の環境(mac)ではうまくいかなかったので、ちょっとだけ手順をかえたものを記事にしてみました。色々と理解しないまま行っているので間違っているかもしれません。ご指摘いただけると嬉しいです。

なお、2021年版ではアップデートがなかったようですが、2019年版はKNIME版もありました。

こちらナイメストさんがnote上で使い方を解説してくださっています。

note.com

KNIME各ノードについて初歩から丁寧に説明してくださっていて、ほとんどKNIMEを使ったことのない私でも「ぽちぽちクリックしてたら動いたけど、背景ではこんなことが行われてるのねー」ととても勉強になるのでおすすめです。

正直、コードをみてもよくわからないのでKNIME版のアップデートが待ち遠しいです!!

なおなお、ついでのついでに皆様TalktorialsフォルダのトップがいずれもT000_templateとなっていることにお気づきになられましたか???

f:id:magattaca:20211231002040p:plain

明日のTeachOpenCADDを作るのは君だ!!!

おしまい。

医薬品と回転異性体 〜実例とちょっと計算の検証〜

こちらは創薬アドベントカレンダーに間に合わなかった記事です。毎度、初歩の初歩の記事で恐縮ですが折角なのでここに供養します。

さて、クイズ!以下の4つの医薬品!構造も作用機序も全然違いますが、ある共通点があります!一体なんでしょうか??

f:id:magattaca:20211226180044p:plain
これらの医薬品の共通点は?

答え・・・「回転異性体アトロプアイソマー)をもつ医薬品」。タイトルでネタバレしてましたけどね!

正直、私は構造式を見ただけではどこに回転異性が潜んでいるか分かりません。ということで、回転異性体判別力(?)をあげるべくエネルギー変化をシミュレーションして遊んでみたいと思います。

1. 回転異性と医薬品

まずは基本的な背景知識のおさらいです。

1-1. 回転異性の復習

低分子医薬品において不斉キラリティー)は重要な要素です。2021年ノーベル化学賞不斉有機触媒」の解説でも、不斉合成の重要性の例として医薬品が引き合いに出されていました。

この不斉の一つ、軸不斉は「分子が不斉中心を持たないが、キラリティーを有するキラリティーの特別な形式の一つ」です(wikipedia)。

軸不斉の例として回転異性アトロプ異性)があり、回転異性を持つ化合物の例としては野依先生が合成方法を確立されたBINAPがあります*1

f:id:magattaca:20211226180408p:plain

要するに、立体障害が大きくて自由回転できないから鏡像異性体ができちゃうぜってやつです。*2

1-2. 冒頭の例ではどこ?

「4つ異なる置換基がついた炭素」みたいな分かりやすい不斉中心がないので、回転異性体の有無はパッと見では分かりにくい気がします。

冒頭に挙げたそれぞれの化合物の不斉に関する軸はここです。

f:id:magattaca:20211226180442p:plain
回転異性体が生じる軸

いずれも(ヘテロ)アリールに結合する箇所です。BINAPと見比べると雰囲気がわかってきました。

ちなみにWikipediaによるとTelenzepine(+)-isomer(-)-isomerムスカリン受容体に対する活性に500倍以上差があるそうです(ラット大脳皮質)。鏡像異性体ってやっぱり生理活性に重要なんですね。

1-3. 回転異性体の安定性と分類

「回転異性体で活性差が大きいなら、良い方のだけ医薬品にすれば良いんじゃない?」となりますが、そんなに簡単ではないそうです。なぜなら、回転異性体の中には安定に分離できるものとできないものがあるからです。

回転異性体は軸周りの回転障壁の高さに由来します。障壁が十分に高ければ、異性体間での変換が十分に遅い(実質、変換しないとみなせる)ので、分離してそれぞれ別の化合物として扱う事ができます。一方で、障壁が中くらいであれば異性体を分離したところで、しばらくするとまた混ざってしまうので分離する意味がありません。

厳しい条件をクリアして医薬品として開発・臨床応用するためには、その化合物の均一性の保証も大事な要素です。生理活性も異なる2つの異性体が作るたびにバラバラに混ざってしまったら困ってしまいます。

以下の様に、回転障壁の高さに応じて回転異性体の安定性はざっと3つに分類できるそうです。(障壁のエネルギーは目安)*3*4

f:id:magattaca:20211226180753p:plain
回転異性体の大まかな分類

エネルギーの低いClass 1であれば、そもそも区別できるほど安定な回転異性体が生じないので一つの化合物として取りあつかえます。また、障壁が高いClass 3であれば年単位で安定に取り扱えます。面倒なのが中間のClass 2で、場合により、単体あるいは混合物として取り扱うようです。

1-4. 回転異性体への対処例

実際のところ、研究過程で回転異性体がでてきたらどうするのでしょうか?今年FDAに承認されたAmgen社のKras G12C阻害薬 Sotorasibの発見に至る経緯が面白いので参照しておきます。

オープンで読めます!*5

研究の途中で見つかった化合物 18は有望な性質を持っていたものの、回転異性体がみつかりました。障壁は26 kcal/molで、25℃での半減期8日。先の分類でいうとClass 2にあたりそうです。*6

f:id:magattaca:20211226185623p:plain

回転異性体に対する戦略として著者らがあげているのは以下の3つです。

  • ① もっと回りにくくして安定な異性体を分離できる様にする(Class 3にする)
  • ② 逆に回りやすくして回転異性体が生じない様にする (Class 1にする)
  • ③ 分子を対称形にして、回転しなくても鏡像異性体にならないようにする

検討の結果、戦略①がうまくいったそうです。軸不斉近傍にメチル基を導入することで回転障壁が30 kcal/molより大きくとなった化合物 24を見出しています。ここからさらに発展させて最終的にSotorasibを見出しています。

f:id:magattaca:20211226181151p:plain

メチル基一つで回転異性という現象を制御してしまうの格好いいですね!承認された医薬品の構造式だけを見ても分かりませんが、そこに至る経緯を辿ると科学者の戦いが垣間見えて楽しいです!

2. 2面角とエネルギー変化

前置きが長くなりました。

構造式だけを眺めるだけではイマイチよくわからない回転異性体ですが、軸周りに回転させたエネルギー変化をたどれば感触がつかめるかもしれません。オープンソースのツールを使ってシミュレーションしてみましょう!

試す手法は以下の2つです。

  • ① RDKitのコンフォマー生成による解析(Molecular Mechanics)
  • ② xTBの2面角スキャンによる解析(Semiempirical Extended Tight-Binding)

以下の2つの記事で解説してくださっている内容を参考にさせていただいています。というよりほぼコピペです。ごめんなさい。

まず試すのは以下のモデル化合物です。先のJ .Med. Chem. の文献を参考にしています。

f:id:magattaca:20211226183132p:plain

それでは2面角に依存してエネルギーがどのように変わるか? 眺めてみましょう!

2-1. モデル化合物 1の描画とatom indexの確認

まずはモデル化合物 1を描画して必要な情報を確認します。Molオブジェクトを作成し、考慮する二面角のatom indexをチェックします。

from rdkit import Chem
from rdkit.Chem import AllChem, Draw

# モデル構造作成
model_1 = Chem.MolFromSmiles("CC(C)c1ccccc1n3c(=O)ncc2cccnc23")

# 水素原子付加
model_1H = AllChem.AddHs(model_1)

Draw.MolToImage(model_1H)

f:id:magattaca:20211226181429p:plain

RDKitのブログ記事より、drawMolDraw2DaddAtomIndices=Trueとすることでatom indexを確認できます。*7

from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

d2d = rdMolDraw2D.MolDraw2DSVG(350,300)
d2d.drawOptions().addAtomIndices=True
d2d.DrawMolecule(model_1H)
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20211226181548p:plain

イソプロピル基の根元から軸周りのatom indexを今回の標的とします。小さくて見づらいですが3, 8, 9, 10が目的の2面角を定義する原子のatom indexです。

ハイライトして描画してみます(highlightAtoms)。

cp = Chem.Mol(model_1H)
d2d = rdMolDraw2D.MolDraw2DSVG(350,300)
d2d.drawOptions().setHighlightColour((0.8,0.8,0.8))
d2d.DrawMolecule(cp,highlightAtoms=[3, 8, 9, 10])
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20211226181759p:plain

結合(bond index)についても同じことはできます(addBondIndices=True)。見やすさのため水素原子を付加せずに確認します。

d2d = rdMolDraw2D.MolDraw2DSVG(350,300)
d2d.drawOptions().addBondIndices=True
d2d.DrawMolecule(model_1)
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20211226181850p:plain

着目部位のbond indexは19, 8, 9でした。

結合にコメントをつけて描画する事もできます(bondNote)。こんな感じ。

cp = Chem.Mol(model_1)
cp.GetBondWithIdx(8).SetProp("bondNote","Here!!")
d2d = rdMolDraw2D.MolDraw2DSVG(350,300)
d2d.drawOptions().setHighlightColour((0.8,0.8,0.8))
d2d.DrawMolecule(cp,highlightAtoms=[],highlightBonds=[19, 8, 9])
d2d.FinishDrawing()
SVG(d2d.GetDrawingText())

f:id:magattaca:20211226181928p:plain

atom indexはSetAtomMapNumを使っておくと途中でズレたりせず安心だそうです*8

また、これを使うと追加の描画の設定をしなくてもatom indexが表示されます。

for i, a in enumerate(model_1H.GetAtoms()):
        a.SetAtomMapNum(i)

Draw.MolToImage(model_1H)

f:id:magattaca:20211226182004p:plain

脱線しましたがatom indexがわかりました。

2-2. RDKitの制約付きコンフォマー生成

モデル化合物のMolオブジェクトとatom indexがわかったので、コンフォマーを生成してエネルギー変化を確認します。

目的の2面角に、180°10°ずつ回転した制約を加えたコンフォマーを作成します。各コンフォマーはMMFF力場で構造最適化します(分子力学法)。

# コード参考元 1; https://future-chem.com/rdkit-constrained-conformer/
# コード参考元 2; https://pschmidtke.github.io/blog/torsion/dihedral/oss/opensource/rdkit/xtb/energy/2021/02/16/torsion-angle-scans-xtb.html

# ベースになる3次元構造を生成(再現性のためrandomSeedを設定)
AllChem.EmbedMolecule(model_1H,randomSeed=10)

# ベースのコンフォマーを取得
conformer=model_1H.GetConformer(0)

# ベースが更新されない様にdeepcopyしておく
import copy
m_1 = copy.deepcopy(model_1H)

# MMFF力場の設定
m_1_p = AllChem.MMFFGetMoleculeProperties(m_1)

# エネルギーを格納するリスト
energy=[]

# 回転する角度の範囲
angles=range(0,180,10)

# コンフォマーのID = 0からスタート
confid=0

# 角度範囲でforループを回す
for angle in angles:

    # MMFF力場
    ff = AllChem.MMFFGetMoleculeForceField(m_1, m_1_p)

    # 目的の2面角に制約を設定(atom indiex : 3, 8, 9, 10)
    ff.MMFFAddTorsionConstraint(3,8,9,10, False, angle - .1, angle + .1, 10000.0)

    # 構造の最適化(minimize)
    ff.Minimize()

    # エネルギーを取得してリストに追加
    energy.append(ff.CalcEnergy())

    # 次の新しいコンフォマーを準備
    ## 新しいコンフォマーID
    confid+=1

    ## 原子数を引数で指定して新コンフォマーを作成
    new_conf = Chem.Conformer(model_1H.GetNumAtoms())

    ## 新コンフォマーに原子の位置を設定
    for i in range(model_1H.GetNumAtoms()):
        new_conf.SetAtomPosition(i, (m_1.GetConformer(-1).GetAtomPosition(i)))

    ## 新コンフォマーにコンフォマーIDを設定
    new_conf.SetId(confid)

    ## 新コンフォマーをベースに追加
    model_1H.AddConformer(new_conf)

コンフォマーができました。コンフォマーの数を確認します。

print(model_1H.GetNumConformers())
#  19

ベースとなるコンフォマーと、0° 〜 170° の10° 刻み18個の合わせて19個が含まれています。アザキナゾリノン骨格をコア構造にしたアラインメントを取って3次元で眺めてみましょう。

# 3次元を眺めるためにpy3Dmolを使う
from rdkit.Chem.Draw import IPythonConsole
import py3Dmol

# 重ねたいコア構造をatomIdで指定してコンフォマーのアラインメントをとる
AllChem.AlignMolConformers(model_1H, atomIds=[9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19])

# 3Dビューワーに表示(ベース以外の18個)
v = py3Dmol.view(width=600, height=600)
for cid in range(1, 19):
    IPythonConsole.addMolToView(model_1H, confId=cid, view=v)
v.setBackgroundColor('0xeeeeee')
v.zoomTo()
v.show()

f:id:magattaca:20211226182200p:plain

ぐちゃぐちゃですが、とりあえず見たい2面角周りの回転はできていそうです。アザキナゾリノン骨格の平面性がところどころ崩れていることに不安はありますが。。。。

角度とエネルギーの関係をプロットしてみます。

# matplotlibの設定
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.style.use('seaborn')

# 角度とエネルギーの関係
plt.plot(angles, energy, 'o')
plt.xlabel('dihedral angle')
plt.ylabel('energy')

f:id:magattaca:20211226182223p:plain

なるほど。よくわからん。。。

2-3. RDKitの制約付きコンフォマー生成~part 2~

先程3Dを眺めた際に骨格の平面性が崩れている感じがありました。どうもこれはRDKitのバージョンによって挙動が違うようです。

参考記事2によるとMMFFout of plane項の設定に問題があるそうです。記事のおすすめに従ってSetMMFFOopTerm(False)を使ってみましょう。

オプションの設定以外は同じコードです。

# molオブジェクト作り直し
model_1 = Chem.MolFromSmiles("CC(C)c1ccccc1n3c(=O)ncc2cccnc23")
model_1H = AllChem.AddHs(model_1)

AllChem.EmbedMolecule(model_1H,randomSeed=10)
conformer=model_1H.GetConformer(0)
m_1 = copy.deepcopy(model_1H)

# MMFF力場の設定
# さっきとここが違う!Out of plane項の設定を変える。
m_1_p = AllChem.MMFFGetMoleculeProperties(m_1)
m_1_p.SetMMFFOopTerm(False)

energy=[]
angles=range(0,180,10)
confid=0

for angle in angles:
    ff = AllChem.MMFFGetMoleculeForceField(m_1, m_1_p)
    ff.MMFFAddTorsionConstraint(3,8,9,10, False, angle - .1, angle + .1, 10000.0)
    ff.Minimize()
    energy.append(ff.CalcEnergy())
    confid+=1
    new_conf = Chem.Conformer(model_1H.GetNumAtoms())
    for i in range(model_1H.GetNumAtoms()):
        new_conf.SetAtomPosition(i, (m_1.GetConformer(-1).GetAtomPosition(i)))
    new_conf.SetId(confid)
    model_1H.AddConformer(new_conf)

新しく作ったコンフォマーを重ね合わせて眺めてみます。(コードは省略)

f:id:magattaca:20211226182432p:plain

悪化しました。平面性が全然保たれていなさそうです。

一応、こっちでも角度とエネルギーの関係をプロットしてみました。

f:id:magattaca:20211226182525p:plain

構造が崩れているので参考にはなりませんが、SetMMFFOopTerm()の設定で結果が大きく変わることは分かりました。

2-4. xTBによる2面角スキャン

RDKitのコンフォマー生成を用いた手法がうまくいかなかったので、今度はxTBを使った2面角のスキャンを試したいと思います。こちらも参考記事2で紹介されていたものです。

xTBオープンソースSemiempirical Extended Tight-Binding(半経験的拡張強結合近似(?))プログラミングパッケージです。Wikipedia-強結合近似によると、 「系の波動関数を各原子の場所に位置する孤立原子に対する波動関数の重ね合わせにより近似する手法」で「量子化学で用いられるLCAO法と密接な関係がある」そうです。*9

conda-forgeにあるのでcondaでインストールできます。*10

参考記事2で行われている方法は「RDKitで作成したコンフォマーをベースにしてxTBでエネルギー最適化を行おう!」というものです。

ここでは先に発生させたRDKitのコンフォマーを使います(SetMMFFOopTerm(False)していないもの)。インプットが悪すぎるかな?とは思いますが、今回はコードと使い方の学習重視ということでお許しください。

xTBで計算を行うには以下の点に注意すれば良さそうです。

  • 入力ファイルの準備: SDFでOK
  • 制約をかけたい2面角の原子のインデックスを設定:xTBは1始まりでRDKitの0始まりとズレるのに注意
  • 2面角の制約は.inpファイルに記載する
  • 出力ファイルはxtbopt.sdf

各角度の計算結果の構造を保持したいので、参考記事と少しコードを変えています。

# コード参考元 2; https://pschmidtke.github.io/blog/torsion/dihedral/oss/opensource/rdkit/xtb/energy/2021/02/16/torsion-angle-scans-xtb.html
# jupyter notebookから実行するためosを利用する
import os

# 検証したい2面角の角度
angles=range(0, 180, 10)

# エネルギーを格納するリスト
xtbenergy=[]

# 結果の構造を格納するリスト
output_structures = []

# rdkitで生成したコンフォマーを使ってループを回して計算
# xTBの入力に使うsdfを作成
for idx,deg in enumerate(angles):
    w = Chem.SDWriter('mol_1H.sdf')
    w.write(model_1H,confId=idx+1)
    w.close()

    # 2面角のatom indexを設定。1ずれるので(3, 8, 9, 10)ではなく(4, 9, 10, 11)
    atoms = '4,9,10,11'

    # 2面角制約の設定に使うxtb input fileを書く
    fh = open("dihedral_constraint.inp","w")
    fh.write("""$constrain
    force constant=1.00
    dihedral: {},{}
    $end""".format(atoms,float(deg)))
    fh.close()

    # xTBの実行
    # OMP_STACKSIZEは並列計算のための設定(計算環境に合わせて変更)
    # OMP_NUM_THREADSは計算に使えるスレッドの数(計算環境に合わせて変更)
    os.system("export OMP_STACKSIZE=4G && export OMP_NUM_THREADS=4,1 && xtb mol_1H.sdf --opt vtight --charge 0 --input dihedral_constraint.inp")

    # 出力のxtbopt.sdfに含まれる情報からエネルギーを取り出す。
    sdr=Chem.SDMolSupplier("xtbopt.sdf")
    for xtbmol in sdr:
        xtbenergy.append(xtbmol.GetProp("total energy / Eh"))

        # 構造を保持する
        output_structures.append(xtbmol)

結果のエネルギーをプロットしてみます。xtbenergyリストに格納されているのはstr型なのでfloatに変換しておきます。RDKitの計算結果とは単位が異なる(ハートリーエネルギー, Eh)ことにご留意ください。

xtb_energy = [float(i) for i in xtbenergy]

plt.plot(angles, xtb_energy, 'o')
plt.xlabel('angles')
plt.ylabel('xtb_energy')

f:id:magattaca:20211226182829p:plain

滑らかなプロットになりました。xTBの入力に使ったRDKitの計算結果(一つ目)と見比べると大きな違いです。

構造を眺めてみましょう。

f:id:magattaca:20211226182925p:plain

xTBで最適化した構造ではほとんどのコンフォマーで平面性が再現できているようです!!すごい!!計算手法でこんなに変わるんですね!

念のため、最適化に伴って2面角の制約が崩れていないか確認しておきます。

from rdkit.Chem import rdMolTransforms

angle_check = []

for m in output_structures:
    conf = m.GetConformer(0)
    ang = rdMolTransforms.GetDihedralDeg(conf, 3,8,9,10)
    angle_check.append(int(ang))
print(angle_check)

# [-3, 10, 22, 31, 41, 50, 60, 70, 79, 89, 99, 109, 119, 129, 138, 148, 158, 167]

少々ズレはありますが、10°刻みで0° ~ 170°の18ポイントとる事ができていました。

3. 別のモデル化合物2で検証

RDKitのコンフォマー生成とxTBを組み合わせる事でそれらしい結果が得られる事が分かりました。別のモデル化合物で同じ計算を行うとどうなるか、試してみましょう。

f:id:magattaca:20211226183309p:plain

# モデル構造作成
model_2 = Chem.MolFromSmiles("CC(C)c1cccc(C)c1n3c(=O)ncc2cccnc23")
# 水素原子付加
model_2H = AllChem.AddHs(model_2)

# Atom Indexの追加
for i, a in enumerate(model_2H.GetAtoms()):
        a.SetAtomMapNum(i)

Draw.MolToImage(model_2H)

f:id:magattaca:20211226183321p:plain

モデル化合物2で目的とする2面角のatom indexは3, 9, 10, 11でした。

xTBに入力するベースのコンフォマーをRDKitで生成します。

# モデル2でコンフォマーを生成
# 初期立体構造の作成
AllChem.EmbedMolecule(model_2H,randomSeed=10)

conformer=model_2H.GetConformer(0)
m_2 = copy.deepcopy(model_2H)

# MMFF力場の設定
m_2_p = AllChem.MMFFGetMoleculeProperties(m_2)

energy_2=[]
angles=range(0,180,10)
confid_2=0

for angle in angles:
    ff2 = AllChem.MMFFGetMoleculeForceField(m_2, m_2_p)

    # モデル2に合わせた2面角のatom indexを指定
    ff2.MMFFAddTorsionConstraint(3, 9, 10, 11, False, angle - .1, angle + .1, 10000.0)
    ff2.Minimize()
    energy_2.append(ff2.CalcEnergy())
    confid_2 += 1
    new_conf = Chem.Conformer(model_2H.GetNumAtoms())
    for i in range(model_2H.GetNumAtoms()):
        new_conf.SetAtomPosition(i, (m_2.GetConformer(-1).GetAtomPosition(i)))
    new_conf.SetId(confid_2)
    model_2H.AddConformer(new_conf)

立体構造とエネルギーはこんな感じでした。

f:id:magattaca:20211226183557p:plain

xTBの最適化をおこないます(indexの設定以外は同じです)。

# 検証したい2面角の角度
angles=range(0, 180, 10)

# エネルギーを格納するリスト
xtbenergy_2 =[]

# 結果の構造を格納するリスト
output_structures_2 = []

# rdkitで生成したコンフォマーを使ってループを回して計算
# xTBの入力に使うsdfを作成
for idx,deg in enumerate(angles):
    w = Chem.SDWriter('mol_2H.sdf')
    w.write(model_2H, confId=idx+1)
    w.close()

    # 2面角のatom indexを設定。1ずれるので(3, 9, 10, 11)ではなく(4, 10, 11, 12)
    atoms = '4, 10,11,12'

    # 2面角制約の設定に使うxtb input fileを書く
    fh = open("dihedral_constraint.inp","w")
    fh.write("""$constrain
    force constant=1.00
    dihedral: {},{}
    $end""".format(atoms,float(deg)))
    fh.close()

    # xTBの実行
    # OMP_STACKSIZEは並列計算のための設定(計算環境に合わせて変更してください)
    # OMP_NUM_THREADSは計算に使えるスレッドの数(計算環境に合わせて変更してください)
    os.system("export OMP_STACKSIZE=4G && export OMP_NUM_THREADS=4,1 && xtb mol_2H.sdf --opt vtight --charge 0 --input dihedral_constraint.inp")

    # 出力のxtbopt.sdfに含まれる情報からエネルギーを取り出す。
    sdr=Chem.SDMolSupplier("xtbopt.sdf")
    for xtbmol in sdr:
        xtbenergy_2.append(xtbmol.GetProp("total energy / Eh"))

        # 構造を保持する
        output_structures_2.append(xtbmol)

エネルギーはこんな感じ。

f:id:magattaca:20211226183803p:plain

モデル2ではモデル1と比べて置換基を増やしたためか、プロットに谷が2つできました。

立体の確認です。別々のmolオブジェクトとしてリストに格納しているのでAlignMolConformersではなくAlignMolを使っています。

v = py3Dmol.view(width=600, height=600)
for mol in output_structures_2:
    AllChem.AlignMol(prbMol=mol, refMol=model_1H, atomMap=[(i, i) for i in [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]])
    IPythonConsole.addMolToView(mol, view=v)
v.setBackgroundColor('0xeeeeee')
v.zoomTo()
v.show()

f:id:magattaca:20211226184110p:plain

今回もxTBの構造最適化により平面性の崩壊はかなりマシになっているようです。

4. 2つのモデル構造の比較

2つのモデル化合物について2面角とエネルギー変化の計算ができました。合わせてプロットしてみたいと思います。

絶対値での比較はズレが大きいので、それぞれのコンフォマーの最小値との相対的なエネルギー差をとってプロットします。

# それぞれのモデルの最小値と比べた相対的なエネルギーのリスト
rel_xtb_energy = [i - min(xtb_energy) for i in xtb_energy]
rel_xtb_energy_2 = [i - min(xtb_energy_2) for i in xtb_energy_2]

# プロット
plt.plot(angles, rel_xtb_energy, label="model_1")
plt.plot(angles, rel_xtb_energy_2, label="model_2")
plt.xlabel('angles')
plt.ylabel('relative_xtb_energy')

plt.legend()
plt.show

f:id:magattaca:20211226184128p:plain

モデル化合物1が青色、モデル化合物2が緑色です。

重ねてみるとどちらも最安定は軸に沿って直交した90°付近にあり、置換基の増えたモデル化合物2は0°、180°付近に近づくに従ってモデル化合物1よりもエネルギーが高く(回転障壁が高く)なっているようです。

それぞれの最大の相対的エネルギー値の差を比較してみましょう。

# 最大の相対的エネルギー値の差
diff_max = max(rel_xtb_energy_2) - max(rel_xtb_energy)
print(diff_max, "Eh")
#  0.007345450856000468 Eh

# 単位を変換(1 hartree = 627.51 kcal/mol)
diff_max_kcal = diff_max * 627.51
print(diff_max_kcal, "kcal/mol")
#  4.609343866648854 kcal/mol

PC Chem Basics.comさんの記事(1Hartreeは何kcal/mol?(単位換算と物理定数))を基に換算したところ、モデル化合物1と2の最大のエネルギー障壁の差4.6 kcal/molとなりました。

ところでこれらのモデル化合物はSotorasibの文献中の化合物(18と24)を簡略化したものです。比較するとこんな感じ。

f:id:magattaca:20211226184312p:plain

化合物18から24への変換でエネルギー障壁は4 kcal/mol以上上昇したとのことですが、今回のモデル化合物の検証結果は4.6 kcal/molでした。

ひょっとして結構良い結果でてる????

おわりに

以上、今回は医薬品と回転異性体をテーマに、その実例の紹介と回転障壁を題材とした計算を試してみました。

回転異性体は構造式をパッとみてわかる不斉中心が無いので分かりづらいように思います(私だけ??)。こんな時、気軽にちょっと計算してみてエネルギー差がわかったりしたら「あーはいはいそんな可能性もあるのねー」と納得しやすく便利な気がします。

計算の題材としては、近年の話題の医薬品Kras G12C阻害薬Sotorasib(AMG510)文献上の化合物をベースにしてみました。

インターネットの賢い人々の検討結果をコピー&ペーストしたところ、意外にも(?)それらしい結果が得られて正直びっくりでした。素晴らしい情報を公開してくださっている皆様に感謝です。

実際には色々な問題点に目をつぶった、無理矢理な議論をしてしまっています。xTBでマシになったとはいえ平面性が崩れたりしていました。。。

もっと良い方法があるよ!お前のやり方はあまりにもひどい!頭悪すぎ!等々、ご意見いただけると嬉しいです。

なお、本記事はtorsion driveって用語格好いい!と思って書きました。*11

Calculation drives new chemistry!!! 

知らんけど

*1:言わずもがな、こちらもノーベル化学賞(2001年)ですね

*2:アレン(allene)みたいなそもそも回転できないものもあります。念のため。

*3:参考: FutureMed.Chem.(2018)10(4),409–422
PuMedでFull textよめます。

*4:もう一つACSのwebinar資料を参考にしています。こちらは何故かインターネット上でPDFが閲覧できるものの出どころ不明なので引用は避けます。

*5:Creative Commonsでは無さそうなので図表の引用は控えます

*6:図中の構造式は論文のSupporting Informationで提供されているSMILESをもとにRDKitで描画しました。

*7:他にはこちらの記事のようにatom indexを確認する方法もあるそうです:RDKitで原子に番号を表示できなくて困っていた話

*8:参考:【RDKit】RWMol で分子(Mol instance) を編集する

*9:誰か説明してください

*10:mambaが使えるならそっちの方が速いかも…

*11:創薬アドベントカレンダーの@iwatobipen先生の記事 Compare torsion drive results between C and O linker of specific structure を読んで、回転って面白いなと思ったのがきっかけです。アドカレに間に合わなくてすみません。