静電ポテンシャルからタンパク質-リガンド複合体を眺めよう〜パート②：Tyk2阻害剤への適用

リガンド-タンパク質複合体を静電ポテンシャル図で眺める記事の後編です。前編で使い方をお勉強したESP-DNNを実際のより複雑な例に適用してみましょう！

計算対象としては②TYK2キナーゼ阻害剤 (ΔG_expデータセット)を用います。

1. 計算対象
- 1-1. Tyk2キナーゼ阻害剤データセット
- 1-2. データセットの取得元と変換
2. ESPを計算して活性差を考察しよう
3. おわりに

1. 計算対象

1-1. Tyk2キナーゼ阻害剤データセット

今回取り上げる例はTYK2キナーゼ阻害剤(ΔG_expデータセット)です。

自由エネルギーとリガンド結合親和性との関係に焦点をあてたデータセットで、もともとは以下の文献に由来するもののようです。

Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field
J. Am. Chem. Soc. 2015, 137, 7, 2695–2703

自由エネルギー計算の分野では有名なものらしく、Open Force Fieldのベンチマーク（protein-ligand-benchmark）にも含まれています。

github.com

また、機械学習を自由エネルギー計算に応用した下記プレプリントでも対象として取り上げられています。(bioRxiv, CC-BY 4.0)

www.biorxiv.org

日本語解説：Preferred Netowrks Blog Neural Network Potentialを用いたタンパク-薬剤結合自由エネルギー計算の改善

具体的には以下のようなキナーゼ阻害剤化合物とデータです。*1

共通の骨格を有し、Rグループが異なる化合物群です。それぞれについて結合親和性のエネルギー（ΔG_exp）が求められています。

置換基の構造で「-7.4 kcal/mol ~ -11.7 kcal/mol」と広い幅で親和性が変化していてとても興味深いですね。

この活性の差を静電ポテンシャルで考察できたら面白そう！！！

1-2. データセットの取得元と変換

さて、Tyk2キナーゼ阻害剤データセットの化合物構造ですが、今回はOpen Force FieldのGitHubに格納されているものを利用させていただくことにしました。*2

github.com

上記にタンパク質（pdb）およびリガンドの構造データ（sdf）が含まれています。

注意点として、ESP-DNNの対応入力ファイルはPDBです。そこでリガンドデータをRDKitでそれぞれPDBに変換して用いました。

from rdkit import Chem

m = Chem.SDMolSupplier("lig_ejm_54.sdf")[0]
mH = Chem.AddHs(m, addCoords=True)
Chem.MolToPDBFile(mH, "lig_ejm_54.pdb")

リガンドのsdfファイルは一つずつ別々のディレクトリに入っています。以下の図に構造とGitHubのディレクトリ名との対応を書いたのでご参考までに。。。

sdfは各化合物ディレクトリのcrdディレクトリの中です。

以上が今回の計算対象とデータの取得元です。

2. ESPを計算して活性差を考察しよう

2-1. 計算対象の絞り込みと傾向

「早速計算してみよう！」といきたいところですが、テーブルの化合物数がそこそこ多いです（16個）。

そこで、データセットの中でも構造が似通ったシクロプロピル基を持つもののみに注目してみます（下図、青枠）。

F置換（14）> 無置換（4）~ Cl置換（11） > CN置換（2）~ メチル置換（7）の順に活性が下がっています。

電子求引基の方が活性が良さそう（F vs. H, Cl > Me）ですが、大きくなると活性に悪そうです（H > CN, Me）。

2-2. PyMolで眺める

２次元の構造式で見た傾向が３次元ではどうなるか？PyMolで描画してみました。

複合体図の下半分はシクロプロピル基周辺の拡大です。近接するアミノ酸残基「Tyr、Arg、Leu」をSpheres表示にしています。

拡大してみると、シクロプロピル基上の置換基を側鎖が取り囲んでいます。「大きい置換基は好まれない」というのは複合体構造でも納得できそうです。

また、親水性のArgがあるので、疎水性のメチル基よりも「分極した電子求引基の方が好まれる」といった考察も可能かもしれません。

・・・このままではよくわからないので、上の構造を静電ポテンシャル図で描画してみましょう！

2-3. ESP-DNNの結果で親和性を考察しよう

タンパク質（Tyk2）と上記3つのリガンド構造をESP-DNNで計算してみました。Google Colab上で計算し、NGL viewerで描画しています。具体的な計算と描画の方法は前回の記事をご参照ください。

2-3-1. タンパク質

まずタンパク質の静電ポテンシャル図はこんな感じです。

ESP図の赤い部分の原子は負の(部分)電荷をもち、青い部分の原子は正の(部分)電荷をもちます。

真ん中左寄りのくぼみがリガンドの結合するポケットです（黄色丸枠）。ポケットを拡大（右下）すると右奥が青色、溶媒に面した左部分に青と赤色が混在しています。

リガンドのシクロプロピル基周辺は溶媒に面した向かって左に相当します。この付近のESPの色がどの残基に対応するか？、眺めてみましょう。

凡そですが、青い箇所は先に見たArgやTyrの側鎖に対応し、赤い箇所は主鎖（Pro）のカルボニル酸素に対応していそうです。*3

水素結合のドナーになりそうな水素原子近辺が青（+）、アクセプターになりそうな箇所が赤（-） と解釈することができるかもしれません。

2-3-2. リガンドの比較

つぎにリガンドのESP図を見てみましょう。こんな感じです。

リガンドの右側、ジクロロベンゼンをみるとCl基以外の水素原子が青色（+） になっています。ベンゼン環上で置換基導入が他の位置に及ぼす電子的な影響が視覚的でわかりやすいですね！

シクロプロピル基は左側、黄色の丸枠で囲んだ部分です。シクロプロピル基の置換基について以下の特徴が見えます。

1. フッ素置換（lig_jmc_23）ではFが赤色
1. 無置換（lig_ejm_46）では青みがかった白
1. メチル置換（lig_jmc_28）では青みがかった白の大きな構造

F以外ではシクロプロピル基全体がプラスっぽく、Fではマイナスっぽさが追加されている感じでしょうか？

置換基とは反対側の面からシクロプロピル基を眺めてみましょう。

こんな感じ。

面白いことに、F置換は、無置換やメチル置換と比べて青色が濃くなっています。強い電子求引基で引っ張られた逆側にプラス性が強く現れるというのはわかりやすいですね！

2-3-3. リガンドとタンパク質の比較

タンパク質とリガンドをそれぞれ確認したので、組み合わせるとどうなるか？重ねてみましょう。

まずは一番シンプルな無置換の構造です。

リガンド右側、ジクロロベンゼンの青い箇所がタンパク質側の赤い静電ポテンシャル面と接しているのが見えます。相補的な赤と青がタンパク-リガンドで接しており、相互作用に適していると考察できそうです。

また、リガンド中央の赤い箇所はシクロプロピル基に隣接したアミド結合のカルボニル酸素です。こちらもタンパク質側の青みがかった部分と接しており、相補的な組み合わせになっています。

データセットで「共通の骨格部分がタンパク質との相互作用に適した構造で一定の親和性を担保している」と考えると興味深いですね。

ついで活性の良いフッ素置換です。

共通骨格は無置換と同じなので、シクロプロピル基（左側）に着目します。

F置換により赤色の箇所が増えています。ここと接するタンパク質の表面は青色となっています（先に見たArg側鎖 etc.の領域に相当）。

つまり、F置換で相補的な静電ポテンシャル面が増加しており、これによって親和性が向上したと考察できそうです。

またリガンド上で、F置換は置換基と反対側のシクロプロピル基表面をより濃い青色に変化させていました。タンパク質ポケットを見ると該当箇所は赤色となっています（上図左の緑丸）。

重ねた図では見えなくなっていますが、「赤色（タンパク）と濃い青色（リガンド）が接している」ため、この箇所でも相補性が良くなっていると推測されます。

次に、活性の落ちるメチル基置換です。

導入されたメチル基は静電ポテンシャルがほぼニュートラル（白〜薄い青色）です。複合体の左側、緑色の丸で囲んだ箇所を見ると、タンパク質の青色表面とメチル基の表面が接しています。

「（部分）電荷を帯びている箇所にニュートラルなものが当たっている」という状態で、相補性が良いとは言い難くなっています。

つまり、「置換基の導入によりタンパクとぶつかる面が増えたが、静電ポテンシャルの相性でメリットがないため、立体障害により親和性が落ちた」と考察できそうです。

「立体障害が親和性にネガティブに働きそう」というのは（２次元構造式の）構造活性相関やPyMolの描画で見た内容とも合いそうですね。

以上をまとめると、小さな置換基でタンパク質と相補的な静電ポテンシャル表面をつくることができる電子求引性のフッ素導入は親和性の向上に最適だった、と言えそうです。

3. おわりに

以上、今回はESP-DNNで計算した静電ポテンシャル図をもとに、Tyk2キナーゼ阻害剤の親和性の差について考察してみました。

定性的な議論ではありますが、「タンパク質とリガンドの表面の色が相補的（赤 - 青）となることと、親和性が良い（ΔG_expが小さい）」ことに関係性がありそうなことが見えました。

（2D/3Dの）構造式や複合体を普通に描画して眺めるだけではわからない観点（ESP）をとりいれることで、一歩踏み込んで活性について議論できそうなことがわかってよかったです。

ESPをまともに計算しようとすると量子力学計算が必要でどうして良いかわかりませんが、ESP-DNNを使うことで簡単かつ高速で計算ができました。

論文のタイトル通り、まさしく実用的（Practical）なツールですね。

Practical High-Quality Electrostatic Potential Surfaces for Drug Discovery Using a Graph-Convolutional Deep Neural Network J. Med. Chem. 2020, 63, 16, 8778–8790

すばらしいツールを公開してくださっている研究者の皆様に感謝です！

さあ、ESP-DNNで観察眼を鍛えて、創薬人格付けチェックに正解できる能力者（エスパー）に君もなろう！

おしまい。

（おまけ）ダジャレ師匠の素晴らしい記事
「私のキャリアパス紹介　〜アカデミアから民間企業に転職して」小寺正明（Preferred Networks）

*1:bioRxiv より図を引用

*2:ライセンス：MIT, CC-BY-4.0

*3:描画の関係で帰属があいまいです。すみません。

2022-06-11

静電ポテンシャルからタンパク質-リガンド複合体を眺めよう〜パート①：ESP-DNNの使い方〜

もはや元旦恒例となっている創薬ちゃん（@souyakuchan）さん主催の創薬人格付けチェック！

元旦恒例 #創薬人格付けチェック！
今回はタンパク質側も実験構造と予測構造を混ぜてあり余計に訳分からなくした。画像(平行法ステレオ図) に示す４つの結合様式のうち、１つが概ね正しい複合体構造だ。
正しいと思うのは何番だろうか？次のツイートで投票を取る。#souyakuchan_quiz #格付けチェック pic.twitter.com/9SL7bd1sBN
— 叢雲くすり (創薬ちゃん) (@souyakuchan) 2022年1月1日

皆さんもちろん正解してますよね？・・・私は毎年外しています（解説を読んでもよくわからない）。

ぶっちゃけリガンドとタンパク質の複合体の図を見ても「当てはまりが良い！」とか、「ここが変！」とかさっぱりわかりません。

「観察眼を鍛えたい！もう一歩踏み込んだ視点が欲しい！」ということで、今回は静電ポテンシャル図（ESP）を描いて複合体を眺めてみたいと思います。

利用するツールは①ESP-DNNで、計算対象としては②TYK2キナーゼ阻害剤データセットを用います。

この記事ではまず簡単な分子（１置換ベンゼン）を例にESP-DNNの使い方を把握します（複合体の描画は次回）。*1

1. ESP-DNN
- 1-1. ESPの復習
- 1-2. ESP-DNNって？
2. ESP-DNNをGoogle Colabで
3. NGL viewerで可視化
4. 拡張子PQRと謎の原子
- 4-1. PQR？
- 4-2. 謎の原子？
5. おわりに

1. ESP-DNN

1-1. ESPの復習

静電ポテンシャル（Electrostatic Potential, ESP）は「ある分子の周りに＋１価の電荷をおいた際に、その電荷が感じるポテンシャルエネルギー」です。

これを分子について可視化するために分子表面として描いたのが静電ポテンシャル図（Electrostatic Potensial Surface）です。等電子密度面に静電ポテンシャル(電位)を色で表現したもので、分子の静電的な特徴を簡単に把握できます。

学部の有機化学で１置換ベンゼンの反応性（ortho-, meta-, para- 配向性）の説明で出てきたりするやつですね。

こういうの。

赤い部分の原子は負の(部分)電荷をもち、青い部分の原子は正の(部分)電荷をもちます。置換基による芳香環上の変化が目に見えるのが面白いですね。*2

ところで以前、どこかで「ニトロベンゼンとピリジンはESPを描くとめっちゃ似てる！」って力説してる先生がいらっしゃった気がします。確かに似てるかも・・・

ってなわけで（？）、ESPをつかって可視化すると構造式で単純に描く以上に分子の性質を考察することができます。生理活性物質への応用を考えると、「リガンドとタンパク質の静電的な相補性の観点から、親和性・活性について議論できる！」ってなりそうです。

つまり、ESP力を鍛えれば相互作用の観察力が上がって創薬人格付けチェックに正解できるはず！

・・・知らんけど

1-2. ESP-DNNって？

さて、とても魅力的なESPによる可視化ですが、きちんと計算しようとすると量子力学（QM）が必要です。例えば「化学の新しいカタチ」さんは以下の記事でPsi4を使った計算方法を解説してくださっています。

future-chem.com

・・・QM計算しんどいですね。

「面倒な計算はディープラーニングで置き換えよう！」ということで、QMのデータセットを学習して短時間の計算で置き換えられるようにしたのがESP-DNNです。

Practical High-Quality Electrostatic Potential Surfaces for Drug Discovery Using a Graph-Convolutional Deep Neural Network
Prakash Chandra Rathi, R. Frederick Ludlow, and Marcel L. Verdonk
J. Med. Chem. 2020, 63, 16, 8778–8790 (Open Access)

Astex Pharmaceuticalsの研究者によって開発されたプログラムで、GitHubで公開してくださっています。（Apache-2.0ライセンス）

github.com

また、ESP-DNNについては@iwatobipen先生が以下の記事で取り上げてくださっています。

Electrostatic Potential Surface(ESP) calculation with GCNN #RDKit #chemoinformaticsiwatobipen.wordpress.com

github.com

上記GitHubでは、Python 2系のコードで書かれていたESP-DNNをPython３系にアップデートしたものを公開してくださっています。

これは遊ぶしかないですね！

・・・でも残念ながら私のMacではうまく動きませんでした。

2. ESP-DNNをGoogle Colabで

OSが合ってないならGoogle Colabを使えばいいじゃない！ということでちょっとだけ書き換えました。

github.com

conda createで環境構築できるようにenvironment.ymlをPython3（@iwatobipen先生の記事で公開されている設定）に合わせただけです。

Google Colab上で以下を実行すればESP-DNNを使うためのconda環境（esp-dnn-env）ができます。

# Linux用のMinicondaの環境を作成
!wget -c https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
!chmod +x Miniconda3-latest-Linux-x86_64.sh
!bash ./Miniconda3-latest-Linux-x86_64.sh -b -f -p /usr/local

# GitHubからESP-DNNをもってくる
!git clone https://github.com/magattaca/ESP_DNN.git
    
# 作業ディレクトリを移動
%cd ESP_DNN

# condaの仮想環境をyamlから作成
!conda env create -f ESP_DNN/environment.yml

これで計算環境の準備はおしまいです（すこし時間かかります）。

ESP-DNNをつかって静電ポテンシャルを計算するコマンドは、計算対象が「リガンドか?タンパク質か？」によって少し変わります。

いずれにせよ入力はPDBファイルなので適当に準備してください（RDKitならChem.MolToPDBFileです）。

リガンド（低分子）の場合は、以下のコマンドで「①仮想環境をアクティベートして、②ligand_filesディレクトリにある複数のPDBファイルをまとめて計算」できます。

%%bash
source activate esp-dnn-env
python -m esp_dnn.predict -m ligand -i ligand_files

タンパク質の場合は引数-m（計算モード, mode）をproteinにすればOKです。

%%bash
source activate esp-dnn-env
python -m esp_dnn.predict -m protein -i protein_files

上記でprotein_filesディレクトリにあるPDBファイルが計算されます（ディレクトリ名は適宜修正してください）。

計算結果はpqr形式のファイルで、入力のPDBファイル名のあとに.pqrが追加されたファイルが同じディレクトリ内に作成されます。

これで計算したものが先の１置換ベンゼンです。簡単ですね！

3. NGL viewerで可視化

ESP-DNNの計算結果の描画はNGL viewerがオススメされています。

NGL viewer : http://nglviewer.org/ngl/

上記のURLにアクセスして「File → open」で目的の「.pqrファイル」を選びます。もしくはファイルをドラッグ&ドロップでもOKです。

例えばニトロベンゼンのpqrファイルだとこんな感じ。

PDBファイルを基にしているので結合次数の情報が入っていない（ベンゼンが全て単結合）であることにご注意ください。

ESP図を表示する手順は「① 表面(Surface)描画を追加して、② 描画方法を静電ポテンシャルに変更」です。

まず表面描画を追加します。

1. 右枠上の横三本線マークをクリックしてStructure オプションのメニューを開く
1. Representationにsurfaceを追加
1. 分子表面が描画され右枠にsurfaceメニューがあらわれる

ついで描画方法を静電ポテンシャルに変更します。

1. surfaceメニューの横三本線マークをクリックしてRepresentation オプションのメニューを開く
1. 各項目を変更
  - surface typeでavを選択
  - scale factorを4.0に変更
  - radius Typeをexplicitに変更
  - opacityをスライドバー調整（表面が透けて骨格が見える）
  - colorSchemeをelectrostaticに変更
  - colorScaleをRed-White-Blueに変更

以上の設定で以下の図ができます。

「目のアイコン」をクリックすると描画の表示をON/OFFできるので、複数の分子を１つのViewerで見るときに便利です。

4. 拡張子PQRと謎の原子

4-1. PQR？

さて、さらっと流しましたがPQRファイルとはなんでしょうか？

PQRはPDBをベースにしたフォーマットで、「PはPDB、Qは電荷（charge）、Rは半径（radius）」を意味するそうです。

この意味に合わせて、保持する情報が以下のように変更されています。*3

① PDB : 占有率（occupancy）→ PQR : 原子の電荷 (atomic charge, Q)
② PDB：温度因子（B factor） → PQR：原子半径（atomic radius, R）

但し、カラム（列）の位置はPQRとPDBでズレていることも多いそうです。

先のニトロベンゼンのファイルで比較してみましょう。

こんな感じ

各行の座標（xyz）情報は同じですが、それ以降の数字が変わっていますね。

4-2. 謎の原子？

ところで上図の４行とも炭素原子（C）ですが、面白いことにPQRファイルの原子半径は「C1, C2, C3 → 1.700」と「C4 → 1.780」で異なっています（黄色枠で示した箇所）　。

このC4はニトロ基置換のイプソ位になります。もう一度描画を確認すると、炭素原子のうちC4だけ白い球が追加されています。他の炭素との違いを反映して原子半径も異なっているのかもしれません（？？？）

この白い球の原子名はHpとなっていて、ニトロ基の窒素原子（N）、酸素原子（O）上にも同様にHpという名前の球が追加されています。また、酸素原子上には他にもHlpという原子が追加されています。

Hとついていますが水素原子ではなさそうです。追加の様子を見る限り非共有電子対や（芳香族以外の）二重結合の電子を表すダミーのようにも見えますが、PQRフォーマットやESP-DNNの説明をみてもそれらしい記述がありませんでした。

お分かりの方、ご教示いただけると嬉しいです。*4

5. おわりに

以上、今回はタンパク質-リガンド複合体を眺めるための観察方法「静電ポテンシャル図」を試す準備として、簡易計算ツール ESP-DNN の使い方と描画方法を確認しました。

ESPの計算には量子力学が必要で大変そうですが、ディープラーニングで置き換えたESP-DNNを使えばとても簡単に計算できました。素晴らしいツールを公開してくださっている皆様に感謝です。

結果の確からしさ（量子力学計算との比較）は検証できませんでしたが、そもそも私はQM計算できない・・・。できない人向けの簡易ツールということで許してください。

次回、実際にリガンドやタンパク質を計算してみたいと思います。

間違いが多いと思うのでご指摘いただけると嬉しいです。ではでは！

*1:ESP-DNNの図はGitHubのFigureより

*2:静電ポテンシャルの説明と見方についてははこちらの記事「静電ポテンシャルマップの誤解」を参考にさせていただきました。ずっと誤解してた。。。

*3:参考：Jmol/JSmol wiki File formats/Coordinates: PQRより

*4:なお、ESP-DNNのexamplesにあるPQRファイルにも同様にHpやHlpがあるので計算が失敗しているわけではなさそうです。

2022-05-27

RDKitで３次元構造に水素原子を付加するときの設定（addCoords）

RDKitの基本的な操作についてのメモです。

「Chem.AddHs()の引数のaddCoords=Trueを今更知った！」って話です。

よくやるミス

立体構造の座標が入ったMolファイル、SDFファイルから化合物を読み込んだ後で、「水素原子を付け足したいなー」って時があります。

まず、molファイルからmolオブジェクトを作ります。

from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole

m = Chem.MolFromMolFile('test.mol')

IPythonConsole.drawMol3D(m)

で、何も考えずにChem.AddHs()とやって、こんな構造がでてがっかりします。

mH = Chem.AddHs(m)

IPythonConsole.drawMol3D(mH)

遠近法みたいですね。

addCoords=Trueをつけると？

違う、そうじゃない！ということでaddCoords=Trueします。

true_mH = Chem.AddHs(m, addCoords=True)

IPythonConsole.drawMol3D(true_mH)

そう、これが欲しかった！

何が起きてる？

何が起きてるのか座標を確認してみましょう。

addCoords=Trueを設定しなかった場合のMOLブロックはこんな感じです。

print(Chem.MolToMolBlock(mH))

"""
    EFZ
         RDKit          3D
    
     30 32  0  0  1  0  0  0  0  0999 V2000
       -4.6650  -32.8110   24.1290 Cl  0  0  0  0  0  0  0  0  0  0  0  0
        0.7540  -37.4060   23.9780 F   0  0  0  0  0  0  0  0  0  0  0  0
       -1.2240  -36.9050   24.5760 F   0  0  0  0  0  0  0  0  0  0  0  0
       -0.8430  -37.4960   22.5660 F   0  0  0  0  0  0  0  0  0  0  0  0
        1.5560  -35.7230   20.0590 O   0  0  0  0  0  0  0  0  0  0  0  0
        0.9480  -35.2860   22.1560 O   0  0  0  0  0  0  0  0  0  0  0  0
       -0.6250  -35.1620   20.4090 N   0  0  0  0  0  0  0  0  0  0  0  0
       -1.6090  -34.5740   21.2770 C   0  0  0  0  0  0  0  0  0  0  0  0
       -2.8080  -34.0130   20.7730 C   0  0  0  0  0  0  0  0  0  0  0  0
       -3.7390  -33.4660   21.6780 C   0  0  0  0  0  0  0  0  0  0  0  0
       -3.4980  -33.4840   23.0470 C   0  0  0  0  0  0  0  0  0  0  0  0
       -2.3320  -34.0240   23.5560 C   0  0  0  0  0  0  0  0  0  0  0  0
       -1.4060  -34.6050   22.6750 C   0  0  0  0  0  0  0  0  0  0  0  0
       -0.0960  -35.2860   23.1970 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.4550  -34.6070   24.3910 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.8900  -34.0500   25.3200 C   0  0  0  0  0  0  0  0  0  0  0  0
        1.5190  -33.2180   26.6670 C   0  0  0  0  0  0  0  0  0  0  0  0
        2.2540  -32.3400   26.2030 C   0  0  0  0  0  0  0  0  0  0  0  0
        1.0630  -32.0720   26.5920 C   0  0  0  0  0  0  0  0  0  0  0  0
       -0.3570  -36.7960   23.5830 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.6810  -35.4120   20.8040 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.0000    0.0000    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
      1 11  1  0
      2 20  1  0
      3 20  1  0
      4 20  1  0
      5 21  2  0
      6 14  1  0
      6 21  1  0
      7  8  1  0
      7 21  1  0
      8  9  2  0
      8 13  1  0
      9 10  1  0
     10 11  2  0
     11 12  1  0
     12 13  2  0
     13 14  1  0
     14 15  1  0
     14 20  1  0
     15 16  3  0
     16 17  1  0
     17 18  1  0
     17 19  1  0
     18 19  1  0
      7 22  1  0
      9 23  1  0
     10 24  1  0
     12 25  1  0
     17 26  1  0
     18 27  1  0
     18 28  1  0
     19 29  1  0
     19 30  1  0
    M  END
"""

水素原子（H）の座標にご注目ください。複数ある水素原子のXYZ座標が全て「0.0000」になっています。

Chem.addHs()はデフォルトでは水素原子の座標を原点にするようです。付加した水素原子がすべて原点に集まっているので、３D描画では原点から飛び出る結合が描かれていたみたいです。

addCoords=TrueのMOLブロックはこうなっています。

print(Chem.MolToMolBlock(true_mH))

"""
    EFZ
         RDKit          3D
    
     30 32  0  0  1  0  0  0  0  0999 V2000
       -4.6650  -32.8110   24.1290 Cl  0  0  0  0  0  0  0  0  0  0  0  0
        0.7540  -37.4060   23.9780 F   0  0  0  0  0  0  0  0  0  0  0  0
       -1.2240  -36.9050   24.5760 F   0  0  0  0  0  0  0  0  0  0  0  0
       -0.8430  -37.4960   22.5660 F   0  0  0  0  0  0  0  0  0  0  0  0
        1.5560  -35.7230   20.0590 O   0  0  0  0  0  0  0  0  0  0  0  0
        0.9480  -35.2860   22.1560 O   0  0  0  0  0  0  0  0  0  0  0  0
       -0.6250  -35.1620   20.4090 N   0  0  0  0  0  0  0  0  0  0  0  0
       -1.6090  -34.5740   21.2770 C   0  0  0  0  0  0  0  0  0  0  0  0
       -2.8080  -34.0130   20.7730 C   0  0  0  0  0  0  0  0  0  0  0  0
       -3.7390  -33.4660   21.6780 C   0  0  0  0  0  0  0  0  0  0  0  0
       -3.4980  -33.4840   23.0470 C   0  0  0  0  0  0  0  0  0  0  0  0
       -2.3320  -34.0240   23.5560 C   0  0  0  0  0  0  0  0  0  0  0  0
       -1.4060  -34.6050   22.6750 C   0  0  0  0  0  0  0  0  0  0  0  0
       -0.0960  -35.2860   23.1970 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.4550  -34.6070   24.3910 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.8900  -34.0500   25.3200 C   0  0  0  0  0  0  0  0  0  0  0  0
        1.5190  -33.2180   26.6670 C   0  0  0  0  0  0  0  0  0  0  0  0
        2.2540  -32.3400   26.2030 C   0  0  0  0  0  0  0  0  0  0  0  0
        1.0630  -32.0720   26.5920 C   0  0  0  0  0  0  0  0  0  0  0  0
       -0.3570  -36.7960   23.5830 C   0  0  0  0  0  0  0  0  0  0  0  0
        0.6810  -35.4120   20.8040 C   0  0  0  0  0  0  0  0  0  0  0  0
       -0.9030  -35.4093   19.4485 H   0  0  0  0  0  0  0  0  0  0  0  0
       -3.0092  -34.0045   19.6916 H   0  0  0  0  0  0  0  0  0  0  0  0
       -4.6687  -33.0188   21.2963 H   0  0  0  0  0  0  0  0  0  0  0  0
       -2.1321  -33.9989   24.6374 H   0  0  0  0  0  0  0  0  0  0  0  0
        1.6124  -33.9674   27.4668 H   0  0  0  0  0  0  0  0  0  0  0  0
        3.1187  -31.9803   26.7800 H   0  0  0  0  0  0  0  0  0  0  0  0
        2.8248  -32.2051   25.2724 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.2807  -31.8348   25.8560 H   0  0  0  0  0  0  0  0  0  0  0  0
        0.6089  -31.3601   27.2970 H   0  0  0  0  0  0  0  0  0  0  0  0
      1 11  1  0
      2 20  1  0
      3 20  1  0
      4 20  1  0
      5 21  2  0
      6 14  1  0
      6 21  1  0
      7  8  1  0
      7 21  1  0
      8  9  2  0
      8 13  1  0
      9 10  1  0
     10 11  2  0
     11 12  1  0
     12 13  2  0
     13 14  1  0
     14 15  1  0
     14 20  1  0
     15 16  3  0
     16 17  1  0
     17 18  1  0
     17 19  1  0
     18 19  1  0
      7 22  1  0
      9 23  1  0
     10 24  1  0
     12 25  1  0
     17 26  1  0
     18 27  1  0
     18 28  1  0
     19 29  1  0
     19 30  1  0
    M  END
"""

水素原子のXYZ座標に数値がはいっていました。

以上です。

おわりに

３次元構造に水素原子を付け足すときに、以下の記事に紹介してくださっているような方法をとったりしていました。

qiita.com

「①一度SMILESに戻して、②molオブジェクトを作り直し、③水素原子を付加して、④３次元座標を生成」といった手順です。

これだと元のファイルの折角の３次元座標を捨ててしまっていて、もったいないなーと思っていました。簡単な設定で解決できることがわかってよかったです。

・・・むしろなぜ今まで気づかなかったのか。

ところで今回例につかった化合物はエファビレンツです。

この不斉をどう合成するか？？・・・答えはアートオブプロセスケミストリーで！

アートオブプロセスケミストリー: メルク社プロセス研究所での実例

作者:安田修祥
化学同人

Amazon

おしまい。

2022-05-05

OpenMMをステップバイステップで　〜Part 10 全体のまとめと感想、あと失敗例〜

OpenMMの使い方を順番にたどるシリーズのまとめです。

この記事はシリーズ全体の総括と各記事へのリンクまとめが主目的です。

また、備忘録としてうまくできなかったことものせておきます。ProLIFの踏み込んだ設定やParmEdのファイル変換機能にふれています。

1. 全体の総括
2. きっかけと感想
3. このシリーズでできなかったこと
4. （おまけ）試したけどうまくいかなかったもの
- 4-2. ProLIFの相互作用パラメータを変えた再解析
- 4-3. OpenMMからAmberのファイルに変換してMMPBSA
5. おわりに

1. 全体の総括

「OpenMMをステップバイステップで」というタイトルでMDシミュレーションの計算の設定や実行、結果の解析方法について複数の記事にまたがって眺めてきました。

前回で「タンパク質ー低分子複合体のシミュレーション解析」まで確認し、（Part 1で設定した）当初の目標「TeachOpenCADDのT019 · Molecular dynamics simulationとT020 · Analyzing molecular dynamics simulations の内容」を全てたどり終わりました。区切りが良いのでこのシリーズは一度終わりとして、全体の総括をしたいと思います。

各記事でとりあげた内容とツールを以下の表にまとめています。最初からタンパク質-低分子複合体のシミュレーションを行うのはハードルが高かったので、タンパク質単独の前処理方法から始めて、順番に計算の設定・実行、解析を眺めました。

記事リンク	おもな内容	ツール
Part 1	GUIツールを使ったタンパク質構造の前処理	OpenMM Setup
Part 2	Part 1の前処理操作をPythonで実行	PDBFixer
Part 3	GUIツールを使ったシミュレーション条件の設定	OpenMM Setup
Part 4	Google Colab上でMDシミュレーションを実行	OpenMM
Part 5	シミュレーション実行スクリプトの仕組みと書き方(Python)	OpenMM
Part 6	トラジェクトリ解析① (RMSD,RMSF,回転半径)	MDAnalysis, NGLView
Part 7	トラジェクトリ解析② (PCA, cross correlation)	MDAnalysis, pytraj
Part 8	低分子-タンパク質複合体の計算系の設定と計算	OpenMM, OpenFF, MDTraj
Part 9	低分子-タンパク質複合体のトラジェクトリ解析	ProLIF, MDAnalysis

メインに利用したツールはOpenMMでしたが、他にも様々なオープンソースのツールが必要でした。出てくるツール群の概要はだいたい以下の記事にまとめてあります。

magattaca.hatenablog.com

また、追加の情報として「OpenMMのよくある質問（FAQ @ Github）」を以下に翻訳しました。MD計算の設定や結果について書かれていて参考になりました。

magattaca.hatenablog.com

2. きっかけと感想

さて、今回MDシミュレーションを試してみようと思ったのは言うまでもなくTeachOpenCADDという素晴らしい教材が用意されていたからですが、それに加えてTwitterでMaking it rainというプロジェクトを教えていただいたことも大きかったです。

Making it rainについてはこちらの記事で取り上げました。

magattaca.hatenablog.com

MD計算にあこがれはあったものの、難しそうだし高性能なコンピュータがないとどうしようもないと思っていましたが、Google Colabを利用して計算できること、使いやすいオープンソースのプログラムが色々と開発されていると知って「これは遊べるぞ！！」となりました。

OpenMMにGUIツール（OpenMM Setup）が提供されていることを知ったのも大きかったです。「① GUIで使い方の流れと出力を把握 → ② Pythonのコードで同じことができるか確認」という手順を踏めたことで、私でも凡その流れを把握することができました。

実際に計算・解析してみて、「おおっ本当にタンパク動くんだな！」「X線だと綺麗にはまっているリガンドもMDだと割とパタパタ動くんだな！」となって楽しかったです。

特に、複合体のMDトラジェクトリで元の複合体構造にはなかった水素結合がリガンドとタンパクの間に形成されていたのは驚きでした。「相互作用も動的なもので、共結晶構造だけが唯一の正解の結合ポーズではないんだなー。」ってことが実感できたのが何よりでした。

さて、優秀な先人の皆様のおかげで学習する環境は整っているものの、私にとっては難しく、予想以上に長期のシリーズ記事になってしまいました。ここまでお付き合いいただきありがとうございました。

取り上げている計算内容はTeachOpenCADDとほとんど一緒なのに長さは約5倍！！・・・この無駄な長さの中に少しでも役立つ情報があればとても嬉しいです。

3. このシリーズでできなかったこと

MD計算は非常に範囲が広く、深い分野なのでこのシリーズで扱えなかった内容はとてもたくさんあります。

例えば、タンパク質の構造前処理一つとっても@Ag_smith先生のこちらの記事ほどは踏み込めていません。

qiita.com

ジスルフィド結合や金属イオンなど、難しいポイントはたくさんありそうです。

また、低分子の力場についてもOpen Forcefieldで開発されている新しい力場をうまく適用することができませんでした。

低分子-リガンド複合体のシミュレーションということで、「相互作用のエネルギーはどうなるのか？」「結合親和性の見積もりにはどのように活かせば良いのか？」といった点も非常に興味深いところではありますが、私の能力では手に負えず扱うことができませんでした。

AmberではMMPB-SA法（？）といったツールがあるようです。（以下の記事ではシミュレーション時間 500nsくらいは欲しい、といった記載がありました。Google ColabのGPUでも1nsに1時間くらいかかったのに・・・）

blog.livedoor.jp

OpenMMを使った結合自由エネルギー計算ツールとしてはChodera Labの開発したYANKといったものがあるようです。

まだまだ他にもトピックがたくさんありますね！もっといい方法があるよー、などなどご教示いただけると嬉しいです。

4. （おまけ）試したけどうまくいかなかったもの

総括は以上として、以下にはうまくいかなかった操作を備忘録として置いておきます。

4-2. ProLIFの相互作用パラメータを変えた再解析

4-2-1. 課題

記事（Part 9）では「低分子 - タンパク質複合体」のMDシミュレーションのトラジェクトリ解析で水素結合の解析を行いました。その中で、ProLIFとMDAnalysisに水素結合検出の差がありました。この点についてProLIFの相互作用の定義を変更することで解析結果が変わるか検討しました。

4-2-2. おさらい

まずおさらいですが計算対象の構造はPDB ID: 3POZでEGFRタンパクのキナーゼドメインとリガンド TAK-285のX線共結晶構造です。

図の通り共結晶構造からはEGFRのヒンジ領域「Met793の主鎖NH」と「リガンドのN原子」との間に水素結合ができても良さそうです。

実際、MDAnalysisによるトラジェクトリ解析では、デフォルトの定義「ドナー(D) - アクセプター(A) 距離 < 3.0Å」、「D - H - A角度 > 150°」を満たすフレームがいくつか見つかりました。

一方で、ProLIF を用いたInteraction Fingerprint解析（IFP）では、この領域に水素結合は検出されませんでした。

以上の原因を探るためProLIFの相互作用検出の定義を確認し、期待した相互作用を検出できるように定義を変更できるか検証します。

4-2-3. ProLIFの相互作用の定義

まず、ProLIFの相互作用の定義について、元論文の記述を確認しましょう。

論文：Bouysset, C., Fiorucci, S. ProLIF: a library to encode molecular interactions as fingerprints. J Cheminform 13, 72 (2021). （オープンアクセス、CC BY 4.0）

判定基準の距離と角度はTable 2に書かれています。

図中に枠で示した通り、水素結合（HBAcceptor、HBDonor）は「距離 ≦ 3.5 Å、角度 130° 〜 180°」です。MDAnalysisの定義よりも緩いですね。こちらは問題なさそうです。

それでは判定基準の原子パターンはどうでしょうか？こちらはSMARTSパターンで定義されているようです。

HBAcceptorは「（N or O or F or アニオン）and (カチオンではない)」です（多分）。「中性か負電荷を帯びているN原子、O原子、F原子」と考えると水素結合のアクセプターとして納得です。

一方でHBDonorは原子番号を使って定義をされていて「窒素(#7)-H or 酸素(#8)-H or 硫黄(#16)-H」です。硫黄が入っているところに生体分子をターゲットとしたツールであることを感じますね。

ここで注目したいのはHBAcceptorの原子の種類の指定方法です。「N」のように大文字のアルファベットで書いた場合、SMARTSでは「アリファティックな窒素原子」を意味するはずです*1。ですので、芳香属性の窒素原子はこのアクセプターの定義では検出できない可能性があります。

4-2-4. RDKitでSMARTSパターンの認識をチェック

ProLIFの定義をふまえて「芳香族性の原子が水素結合アクセプターからはじかれている」可能性が出てきました。

「SMARTSパターンでどのような構造が認識されるか？」、RDKitで検証してみましょう。ProLIFも構造式の取り扱いはRDKitをベースとしているようなので、RDKitで検証しておけば挙動が予測できるはず・・・多分。

今回計算対象のリガンド（TAK-285）の構造を作って、HBAcceptorのSMARTSパターンと部分構造検索でマッチするか検証します。

from rdkit import Chem
from rdkit.Chem import Draw

# SMILESからTAK-285を構築
smi_tak285 = "CC(C)(CC(=O)NCCn1ccc2c1c(ncn2)Nc3ccc(c(c3)Cl)Oc4cccc(c4)C(F)(F)F)O"
rd_tak285 = Chem.MolFromSmiles(smi_tak285)
rd_tak285H = Chem.AddHs(rd_tak285)

Draw.MolToImage(rd_tak285H)

SMARTSからMolオブジェクトを作ってGetSubstructMatchesでマッチする部分構造を見つけ、ハイライトして描画します。

# SMARTSのMolオブジェクト化
HBA_query = Chem.MolFromSmarts('[N,O,F,-{1-};!+{1-}]')

# 部分構造検索と該当箇所の取り出し
HBA_tak285 = rd_tak285H.GetSubstructMatches(HBA_query)
matched_HBA = [a[0] for a in HBA_tak285]

# 描画
Draw.MolToImage(rd_tak285H, highlightAtoms = matched_HBA)

少し小さいですが赤い丸で示されているのがSMARTSに適合した部分です。構造式の左から時計回りにOH基、アミド、アニリンN、エーテルのO、CF3基のFと言ったものが認識されています。

ここで芳香族ヘテロ環（ピロロピリミジン環？）の窒素原子はクエリにかかっていないことに注目です。

MDAnalysisで解析した対象のリガンド窒素原子は、ProLIFではそもそも水素結合アクセプターを取りうる原子として認識されていなかった可能性が示唆されました。

ちなみに、HBDonorの定義のように原子番号でSMARTSを定義するとどうなるでしょうか？

# N → #7に変更して再検証
# SMARTSのMolオブジェクト化
HBA_query_2 = Chem.MolFromSmarts('[#7,O,F,-{1-};!+{1-}]')

# 部分構造検索と該当箇所の取り出し
HBA_tak285_2 = rd_tak285H.GetSubstructMatches(HBA_query_2)
matched_HBA_2 = [a[0] for a in HBA_tak285_2]

# 描画
Draw.MolToImage(rd_tak285H, highlightAtoms = matched_HBA_2)

今度はリガンドの全ての窒素原子がハイライトされました。水素結合ドナーになりそうにない窒素原子（ピロール様のN？）まで含んでしまっているところが問題ではありますが、、、適切なクエリを作るのは難しいですね。

4-2-5. ProLIFの相互作用パラメーターを変えよう

RDKitによる検証で水素結合ドナーの定義を変えてやれば、前回の宿題の水素結合も認識される可能性が出てきました。

ProLIFでは検出する相互作用のパラメーターを変えたり、新しい相互作用を定義することができます。

ドキュメント①：2.1 Changing the parameters for an interaction

ドキュメント②：2.2. Writing your own interaction

ProLIFのデフォルトの各相互作用のクラスはドキュメントのprolif.interactionsに記載があり、水素結合アクセプターは「prolif.interactions.HBAcceptor(donor='[#7,#8,#16][H]', acceptor='[N,O,F,-{1-};!+{1-}]', distance=3.5, angles=(130, 180))」となっています。

アクセプターのパラメーターを変えた新しいクラス（HBAcceptor_allN）とを作ってみましょう！

import prolif as plf

class HBAcceptor_allN(plf.interactions.HBAcceptor):
    def __init__(self):
        super().__init__(acceptor="[#7,O,F,-{1-};!+{1-}]")

元々のクラスを継承した新しいクラスを作成し、super()で親クラスの__init__メソッド実行時のパラメータを書き換えています*2

ここではアクセプターのSMARTSのNを#7に変更しました。

同様にして水素結合ドナー（prolif.interactions.HBDonor(donor='[#7,#8,#16][H]', acceptor='[N,O,F,-{1-};!+{1-}]', distance=3.5, angles=(130, 180))）についても新しいクラスを作っておきましょう。

class HBDonor_allN(plf.interactions.HBDonor):
    def __init__(self):
        super().__init__(acceptor="[#7,N,O,F,-{1-};!+{1-}]")

これらの新しい水素結合クラスを使って、ProLIFの水素結合認識範囲はひろがるでしょうか？？？

4-2-6. MDトラジェクトリの再解析

それではMDトラジェクトリをもう一度ProLIFで解析してみましょう。

まずはMDAnalysisで読み込んで解析の準備をします。

import MDAnalysis as mda
from MDAnalysis.analysis import align

# topology（.pdb）とtrajectory（.dcd）からUniverseオブジェクトを構築
md_universe = mda.Universe("complex_simulation_data/complex_prep.pdb", "complex_simulation_data/trajectory.dcd")

# Cα原子でアラインメントを取る
alignment = align.AlignTraj( mobile = md_universe, reference = md_universe, select="name CA", in_memory=True)
alignment.run()

リガンドとタンパクを選択します。

# MDAnalysisでリガンドとタンパクを選択
lig = md_universe.atoms.select_atoms("resname LIG")
prot = md_universe.atoms.select_atoms("protein")

新しくつくったProLIFの相互作用（['HBAcceptor_allN', 'HBDonor_allN']）を指定したFingerprintを作ります。

# 解析したい相互作用のFingerprint
fp = plf.Fingerprint(['HBAcceptor_allN', 'HBDonor_allN'])

# トラジェクトリを渡して解析実行
fp.run(md_universe.trajectory[::10], lig, prot, residues = "all")

df_fp = fp.to_dataframe()
df_fp

検出されたのは前回と同じくASP155との相互作用だけでした。解析対象をトラジェクトリの全フレームにしてもASP155以外は見つかりませんでした。

定義を変えてacceptorの窒素原子を全て（アリファティック以外も含む）に変更しても、目的のヒンジ領域との相互作用は検出できないようです。残念。

4-3. OpenMMからAmberのファイルに変換してMMPBSA

先にAmberToolsにはエネルギー計算（MMPB-SA.py）するためのツールがあると書きました。

Google検索しているとOpenMMで計算した結果をAmberの対応ファイルに変換してエネルギー計算しようとしている例などが出てきました。

そこでまずはファイル変換方法を試しました。

4-3-1. ParmEdのインストール

ファイルの変換ができるツールとしてParmEdがあります。（Making it rain でも使われていました）

ParmEd ドキュメント

ParmEd GitHub

インストールはcondaでOKです。

conda install -c conda-forge parmed

4-3-1. ParmEdのファイル変換でエラー

ParmEdでファイル変換を行うには、目的の構造をParmEdのStructureオブジェクトに変換したのち.save()で望みのファイル形式を選択します。

OpenMMからParmEdのStructureを作るには、load_fileで「①トポロジーと②系（SystemのXML）」を読み込みます。

import parmed as pmd

# トポロジーのためにPDB読み込み
test_pdb = pmd.load_file("complex_simulation_data/complex_prep.pdb")

# OpenMMで作ったsystem.xml情報の読み込み
test_system_xml = pmd.load_file("complex_simulation_data/system.xml")

# ParmEdのStructure作成
test_struct = pmd.openmm.load_topology(test_pdb.topology, test_system_xml)

ところがこれをAmberのファイル形式（prmtop or parm7）に変換しようとするとエラーが出ました。

test_struct.save("test_amb.prmtop")

"""エラ-メッセージ
〜省略〜
       1579 for bond in self.bonds:
    -> 1580     bond.type.used = True
       1581 self.bond_types.prune_unused()
       1582 data['BOND_FORCE_CONSTANT'] = [type.k for type in self.bond_types]

    AttributeError: 'NoneType' object has no attribute 'used'
"""

以上のように、「AttributeError: 'NoneType' object has no attribute 'used'」というエラーが出ます。bond_typeがNoneTypeになっているのがいけないようです。

4-3-2. エラーの原因と回避策

同様のエラーの報告例があり、回避策も２通りありました。

GitHub Issue: Bond-related error saving AMBER prmtop from ParmEd structure #930

Issues Antenna: Error writing out prmtop file from OpenMM system

上記の議論によると、OpenMMでforcefieldからcreateSystemで系を構築する際に、拘束条件（constraints）が設定されているとbondのパラメータを消してしまうそうです。そのためParmEdがbond typeをアサインできずNoneのエラーになるとのこと。

回避策①として、createSystemの引数にflexibleConstraints=Trueを与えて系を構築する方法が挙げられていました。

例えば、Part 8で使ったコードを修正すると以下のようにします。

# flexibleConstraints=TrueでSystemをつくる
system = FF.createSystem(complex_pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1.0*unit.nanometers, 
                         constraints = HBonds, rigidWater=True, ewaldErrorTolerance=0.0005, flexibleConstraints=True)

# SystemをXMLで出力する
with open("flexConst_system.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(system))

以上のようにして作成したSystemのXMLをParmEdで読み込んでみます。

# 作り直したSystemのxmlを読み込み
test_system_xml2 = pmd.load_file("flexConst_system.xml")

# ParmEdのStructure作成
test_struct2 = pmd.openmm.load_topology(test_pdb.topology, test_system_xml2)

# Amber形式で出力
test_struct2.save("test_amb.prmtop")

作り直したSystemの情報を使うと無事「test_amb.prmtop」ファイルが作成されました！

また、OpenMMのSystemを作り直すのが面倒な場合、回避策②としてダミーのbond typeを割り当てる方法も挙げられていました。

# ダミーのbond typeを割り当てる
bond_type = pmd.BondType(1.0, 1.0, list=test_struct.bond_types)
test_struct.bond_types.append(bond_type)
for bond in test_struct.bonds:
    if bond.type is None:
        bond.type = bond_type

# Amber形式で出力
test_struct.save('bond_type_assinged.prmtop')

こちらの方法でも無事「bond_type_assinged.prmtop」ファイルが作成されました！

4-3-3. MMPB-SAの失敗

Amber形式にファイル変換する方法がわかったので、MMPB-SA計算ができるかどうかを試します。

計算に必要なファイルは複合体、タンパク質のみ、リガンドのみのトポロジーのようです。

そこで先のParmEdのStructureからそれぞれ必要な部分構造を抜き出して以下のようにファイルを作成しました。

# タンパク質を残基番号 1-317 で選択して出力
pmd_protein = test_struct[':1-317']
pmd_protein.save('protein.prmtop')

# リガンドを残基名 LIG で選択して出力
pmd_ligand = test_struct[':LIG']
pmd_ligand.save('ligand.prmtop')

# 複合体を出力
pmd_complex = test_struct[':1-317, LIG']
pmd_complex.save('complex.prmtop')

print(pmd_protein)
print(pmd_ligand)
print(pmd_complex)

#  <Structure 5120 atoms; 317 residues; 5176 bonds; PBC (orthogonal); parametrized>
#  <Structure 63 atoms; 1 residues; 66 bonds; PBC (orthogonal); parametrized>
#  <Structure 5183 atoms; 318 residues; 5242 bonds; PBC (orthogonal); parametrized>

原子数を見る限り目的の構造を選択できていそうです。

これらのファイルを使って、Making it rainを参考にante-MMPBSA.py``MMPBSA.pyを使おうとしました。

ところが構造を作れないといったエラーがでてうまくいきませんでした。エラーの中身を見てもよくわかりませんでした。Amberの勉強をしないとわからなさそうです。。。

以上、試してうまくいかなかったこと２点でした。

5. おわりに

以上、今回は関連記事のまとめとMD計算を試した感想でした。また、失敗例は使い道はありませんが、検討中にわかったこと（ParmEdの使い方）もあったので備忘録として残しておきました。

Part 1から3ヶ月以上もかかってしまいましたが、素人のチュートリアルお試し記事にお付き合いいただきありがとうございました！

色々とご指摘・ご教示いただければ幸いです。

ではでは！

*1:参考：Daylight SMARTSの説明

*2:参考：ProLIF ドキュメント 2.1.2. Reparameterizing an interaction with another name

2022-05-04

OpenMMをステップバイステップで〜 Part 9 低分子-タンパク複合体のトラジェクトリ解析〜

OpenMMの使い方を順番にたどる記事のPart 9です。

前回の記事（Part 8）では「低分子 - タンパク質複合体」についてシミュレーション系の構築と計算の実行を行いました。この記事では結果の解析を行います。

一般的な解析手法についてはPart 6、Part 7で確認しました。この記事では「低分子リガンドとタンパク質の相互作用解析」に焦点を当てたいと思います。

コードは基本的にTeachOpenCADD T020：Analyzing molecular dynamics simulationsとMaking it rain : Protein-Ligand smulationsからお借りしています。

1. キナーゼと阻害剤
2. 低分子 - タンパク質複合体のトラジェクトリ解析
3. おわりに

1. キナーゼと阻害剤

具体的な解析を眺める前に、計算対象の構造の特徴について前提知識のおさらいをしておきます。

今回の構造はPDB ID: 3POZでEGFRタンパクのキナーゼドメインとリガンド TAK-285のX線共結晶構造です。

キナーゼは「ATPなどの高エネルギーリン酸結合を有する分子からリン酸基を基質あるいはターゲット分子に転移する（リン酸化する）酵素の総称」（Wikipedia - キナーゼ）です。なのでATP結合部位があることが多く、EGFRタンパクの場合以下のようなイメージになるでしょうか？*1

ATPのアデニンと水素結合を作っているのがヒンジ領域（Hinge region）と呼ばれる部分で、EGFRでは「Gln(Q) 791, Leu(L) 792, Met(M) 793」に相当します（あとで解析の際に残基番号を使います）。

キナーゼ阻害剤は結合様式によって分類され「ATP結合部位とその近縁のみ」に結合するものはType-Iと呼ばれるそうです。

木下誉富キナーゼを標的とした構造生物学および創薬の現状日本結晶学会誌 59，174-181(2017)

計算の対象としている複合体の構造では結合部位はこのようなイメージです（見やすさ重視で加工してます）。

アデニンに類似したピロロピリミジン環（？）がヒンジ領域のMet 793としっかり水素結合を作っていそうです。ただし、水素結合に関与しているピリミジン部位の窒素原子の位置は異なるようです。ATPとTAK-285の構造式を２次元でならべるとそっくりに見える部分構造でも、タンパク質との複合体では完全には重ならないというあたりに構造生物学の面白さを感じます。

なお、ピロロピリミジン環の先に伸びるフェニル環部分はEGFRポケットのゲートキーパー、疎水性ポケットへと伸びています。初期のキナーゼ阻害剤に対する獲得耐性変異としてT790Mなどが知られているそうです。確かに結合部位内で変異が起きてしまったらリガンドが結合しづらくなりそうですね。*2

ざっとですがキナーゼ阻害剤のおさらいは以上です。

2. 低分子 - タンパク質複合体のトラジェクトリ解析

2-1. トラジェクトリの動画

最初にPyMolで確認したトラジェクトリの動画を貼っておきます。（前回の記事の最後に貼ったもの）

リガンドが結構動いてますね。

2-2. MDAnalysisでRMSD解析

2-2-1. RMSDの推移

それでは具体的な解析を進めていきましょう！

まずはシミュレーションが崩壊していないかの確認も含めてトラジェクトリ全体のRMSDの推移を見てみます。

MDAnalysisを使うので、データを読み込んでUniverseオブジェクトを作成しアラインメントをとりましょう。

import MDAnalysis as mda
from MDAnalysis.analysis import align

# topology（.pdb）とtrajectory（.dcd）からUniverseオブジェクトを構築
md_universe = mda.Universe("complex_simulation_data/complex_prep.pdb", "complex_simulation_data/trajectory.dcd")

# Cα原子でアラインメントを取る
alignment = align.AlignTraj( mobile = md_universe, reference = md_universe, select="name CA", in_memory=True)
alignment.run()

トラジェクトリ全体のRMSDはMDAnalysis.analysis.rms.RMSD()で確認できます。

今回は主鎖（backbone）のRMSDを取ると同時に、タンパク質全体（protein）とリガンド（resname LIG）も追加で計算します。*3

追加の計算対象は「引数 groupselections」で指定します。selectで選んだメインのAtomGroupについて重ね合わせた状態で、追加のAtomGroupについて計算が行われます。

from MDAnalysis.analysis import rms

# フレームのリセット
md_universe.trajectory[0]

# RMSD解析（メイン "backbone", 追加 ["protein", "resname LIG"]）
RMSD_analysis = rms.RMSD(md_universe, reference = md_universe, select = "backbone", groupselections =  ["protein", "resname LIG"], ref_frame=0)
RMSD_analysis.run()

結果をプロットしましょう。

results.rmsdに「[[frame, time(ps), RMSD(A), RMSD(A), RMSD(A)], [...], [...]]」形式で入っているので、PandasのDataFrameにしてからmatplotlibでプロットします。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 結果のDataFrameへの変換
column_names = ['frame', 'time(ps)', 'backbone', 'protein', 'ligand']
rmsd_df = pd.DataFrame(np.round(RMSD_analysis.results.rmsd, 2), columns = column_names)

# 結果のプロット
plt.plot(rmsd_df['time(ps)'], rmsd_df['backbone'], label="backbone")
plt.plot(rmsd_df['time(ps)'], rmsd_df['protein'], label="protein")
plt.plot(rmsd_df['time(ps)'], rmsd_df['ligand'], label="ligand")
plt.legend(loc="upper left")
plt.xlabel("Time (ps)")
plt.ylabel("RMSD [$\AA$]")
plt.title("RMSD of protein and ligand")

plt.show()

主鎖が青色、タンパク全体がオレンジ色、リガンドが緑色です。タンパク質に比べるとリガンドの振れ幅が大きいですね。ただし、これはbackbone基準で合わせた場合なので、リガンドだけを考慮した場合とは異なることにご注意ください。

2-2-2. フレーム間のRMSD

先に見たRMSDは最初のフレーム（トラジェクトリ 0 ps）に対するものだったので、フレーム間のペアワイズのRMSDを見てみましょう。diffusionmap.Distancematrixを使います。

from MDAnalysis.analysis import diffusionmap

# proteinについてペアワイズRMSDを計算
pairwase_protein = diffusionmap.DistanceMatrix(md_universe, select = "protein")
pairwase_protein.run()
dist_matrix_protein = pairwase_protein.results.dist_matrix

# ligandについてペアワイズRMSDを計算
pairwase_ligand = diffusionmap.DistanceMatrix(md_universe, select = "resname LIG")
pairwase_ligand.run()
dist_matrix_ligand = pairwase_ligand.results.dist_matrix

結果は2Dのマトリックスです（ndarray）。ヒートマップで可視化しましょう！

proteinとligandでヒートマップの色のレンジを合わせるため、最初にフレーム間の距離の最大値をとっておいてvmaxに指定します。

# コードはTeachOpenCADD T20より
# フレーム間の距離の最大値を取得（heatmapのvmaxに使用）
max_dist = max(np.amax(dist_matrix_ligand), np.amax(dist_matrix_protein))

# プロットを作成
fig, ax = plt.subplots(1, 2)
fig.suptitle("RMSD between the frames")

# タンパクについてsubplot
img1 = ax[0].imshow(dist_matrix_protein, cmap="viridis", vmin=0, vmax=max_dist)
ax[0].title.set_text("protein")
ax[0].set_xlabel("frames")
ax[0].set_ylabel("frames")

# リガンドについてsubplot
img2 = ax[1].imshow(dist_matrix_ligand, cmap="viridis", vmin=0, vmax=max_dist)
ax[1].title.set_text("Ligand")
ax[1].set_xlabel("frames")

fig.colorbar(img1, ax=ax, orientation="horizontal", fraction=0.1, label="RMSD (Å)")

タンパク質のRMSDは時間推移で近いフレーム間では小さく、離れたフレーム間では大きくなっているようです。時間発展に伴って少しづつ全体の構造が変わっているようですね。

対して、リガンドではフレームが離れた箇所でもペアワイズのRMSDが小さくなる部分がありました。時間経過のプロットでも大きくRMSDが動いている箇所があったので、一定の範囲内で構造の揺り戻しが起きているのかもしれません。

2-3. Protein-Ligand Interaction Fingerprint（ProLIF）で相互作用解析

2-3-1. ProLIF？

つづいて具体的に「どの残基がリガンドと相互作用をしているか？」解析したいと思います。

ここではProLIF（Protein-Ligand Interaction Fingerprints）というツールを使います。タンパク質とリガンドとの相互作用をフィンガープリント形式にして解析してくれる素敵なツールです。

論文：Bouysset, C., Fiorucci, S. ProLIF: a library to encode molecular interactions as fingerprints. J Cheminform 13, 72 (2021).
ProLIF GitHub
ProLIF Documentation

ProLIFとは別のツールで恐縮ですが、動作のイメージをつかむためにPLIP（Protein-Ligand Interaction Profiler）の解析図を貼っておきます。*4

上図のようにPLIPではタンパク質とリガンド間にある（非共有結合の）相互作用を種類別に分けて解析してくれます。「ある残基との間にある種類の相互作用が有るか？無いか？」をフィンガープリントとして捉えるというイメージです。

PLIPの入力はPDBファイルで１点の静止画ですが、ProLIFを使うとMDトラジェクトリ（動画）を対象とした解析ができます。

ProLIFのインストールにはいくつか必要なライブラリがあります。ドキュメントのInstallationを参考にインストールしましょう。

推奨はcondaで専用の仮想環境を作った上で、githubからpipでインストールするみたいです。

# 仮想環境を作る
conda create -n prolif

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

# 主な依存ライブラリのインストール（これら以外は自動でインストールされる）
conda config --add channels conda-forge
conda install rdkit cython

# ProLIFのインストール
pip install git+https://github.com/chemosim-lab/ProLIF.git

※ 注意「pip install prolif」は使えないそうです。

私はMDお試し用の仮想環境に上記でインストールしました。

2-3-2. ProLIFでフィンガープリント解析

それではProLIFで相互作用フィンガープリントを解析してみましょう。

ProLIFはMDAnalysisで読み込んだトラジェクトリを対象として解析できます。リガンドおよびタンパク質の指定が必要なのでMDAnalysisでslect_atomsしてからProLIFのFingerprintを使いましょう。

import prolif as plf

# MDAnalysisでリガンドとタンパクを選択
lig = md_universe.atoms.select_atoms("resname LIG")
prot = md_universe.atoms.select_atoms("protein")

# デフォルトの相互作用解析
fp = plf.Fingerprint()

# トラジェクトリを渡して解析実行
fp.run(md_universe.trajectory[::10], lig, prot)

上記はトラジェクトリを10フレームごとに解析する場合です（スライスで指定）。今回は50フレームのトラジェクトリなので５セットの結果が出ます。

デフォルトではリガンドの周囲6.0Åの範囲にある残基を対象として解析を行います。

結果はPandasのDataFrameに出力することができます。

df_fp = fp.to_dataframe()
df_fp

ProLIFのデフォルトの解析で見つかった相互作用（True）は「0フレーム目でASP155.Aが水素結合ドナー（HBDonor）」として働いているというもののみでした。残りのフレーム、また他の残基との間に相互作用は検出できなかったようです。

今回のシミュレーション構造は「残基番号 700 〜 1017」なのでASP155.Aは「鎖AのD855」に相当します。ちょうどDFGループのDですね。

期待していたのはヒンジ領域（QLM）との相互作用でしたが出なかったものは仕方ない...残念。

2-3-3. ProLIFでLigand Interaction Network（LigPlot）

DataFrameのみだと味気ないのでもう少し凝った描画をしてみましょう。

ProLIFのLigNetworkを使うと、LigProtのようなお洒落な２Dダイアグラムが作成できます。

引数にはProLIFの「① 相互作用フィンガープリント解析結果（ifp）のデータフレーム」と「②リガンドをMoleculeオブジェクトにしたもの」を渡します。ただし、データフレーム（①）は「引数return_atoms=True」を設定している必要があるのでもう一度作り直しています。

from prolif.plotting.network import LigNetwork

# ① LigPlot用のifpのDataFrame
df_lpt = fp.to_dataframe(return_atoms = True)

# ② リガンドをMDAnalysisのオブジェクトからProLIFのMoleculeオブジェクトにしておく
lmol = plf.Molecule.from_mda(lig)

# ifpのDataFrameから2Dネットワークを作成
# フレーム 0から2D描画を作る場合
net = LigNetwork.from_ifp(df_lpt, lmol, kind='frame', frame=0)
net.display()

できました！アスパラギン酸（ASP155）の酸性側鎖とリガンド（TAK-285）の水酸基との間に水素結合ができているようです。シンプルでみやすいですね。

今回は最初のfram 0でのみ相互作用が検出（True）されていたので LigNetwork.from_ifpの「引数 kind='frame', frame=0」としました。

全てのフレームを合算した結果を描画したい時は「kind='aggregate'」とすればOKです。また、「threshold=0.3」とすれば合算した３割以上（閾値 0.3以上）の割合で見えている相互作用に限定して描画することもできます。

ちなみにPyMolでトラジェクトリの最初の方を確認するとこんな感じです。

ヒンジ領域も結構近いのに検出してくれないのかー。厳しい。。。*5

以上でProLIFを使った解析はおしまいです。

2-4. MDAnalysisで原子間距離の解析

つづいてMDAnalysisを使ってもう少し細かい解析を行いましょう。ProLIFで全体をざっとみただけでは拾いきれなかったものが見えるかもしれません。

まずは原子間距離を測ってみたいと思います。

シミュレーション時にリガンドの各原子につけた名前は以下のようになっています。

また、タンパク質は以下のように前処理をしているため、シミュレーション上の残基の番号は元の構造と700ずれます。

MDAnalysisでは各原子を「残基ID（resid） or 残基名（resname）」と「原子の名前（name）」の組み合わせで指定することができます（select_atoms）。

先の注意点を踏まえると「Met 793（→resid 93）の主鎖N原子」と「リガンド（resname LIG）のN4」は以下のように指定できます。

# EGFRのMet793のN
atomgroup_M793N = md_universe.select_atoms("resid 93 and name N")

# LigandのN4
atomgroup_LIG_N4 = md_universe.select_atoms("resname LIG and name N4")

この２原子の距離はdistances.distを使って求められます。

from MDAnalysis.analysis import  distances

dist_M793N_LigN4 = distances.dist(atomgroup_M793N, atomgroup_LIG_N4)
print(dist_M793N_LigN4)

#  [[93.        ]
#  [ 1.        ]
#  [ 3.30475404]]

戻り値はnumpy.ndarrayで「[[一つ目のAtomGroupの残基ID], [２つ目のAtomGroupの残基ID], [原子間距離(Å)]] 」となっています。

従ってMet残基とリガンドN4原子の距離は3.3Åとなります。トラジェクトリに沿った変化を眺めてみましょう。

# トラジェクトリの各フレームごとに原子間距離を求めてリストに格納
distlist_M793N_LigN4 = []

for _ in md_universe.trajectory:
    dist_t = distances.dist(atomgroup_M793N, atomgroup_LIG_N4)[2][0]
    distlist_M793N_LigN4.append(dist_t)

# 最小値、最大値
print("shortest length(A) : ", np.round(min(distlist_M793N_LigN4), 2))
print("longest length(A) : ", np.round(max(distlist_M793N_LigN4), 2))

# プロット
plt.plot(distlist_M793N_LigN4)
plt.ylabel("distance (Å)")
plt.xlabel("frame")
plt.ylim(0.0, 4.0)
plt.title("Distance between Met793's N and LIG's N4")
plt.show()

# shortest length(A) :  2.92
#  longest length(A) :  3.92

全体では「2.9Å〜 3.9Å」の範囲を行ったり来たりしているようです。

ちょっと長めですが最小値周りでは水素結合を形成する距離に入っていると言って良いかもしれません。

ところでこういう場合、結晶構造の分解能(PDB ID : 3POZ, Resolution: 1.50 Å）はどう考慮すれば良いのでしょうか？MD計算したから解消されていると考えるのか？計算後も誤差として分解能を考慮すべきか？・・・よくわかりません。ご教示いただけると嬉しいです。

2-5. MDAnalysisで水素結合の解析-1

さて、原子間距離の解析はできましたが水素結合を考慮するには角度も合わせて考える必要があります。

MDAnalysisには水素結合解析の機能も用意されています（hydrogenbonds.hbond_analysis））。

設定する引数が多いですが１つ１つはシンプルです。

① ドナー（D）原子（donors_sel(str)）
② 水素原子（H）本体（hydrogens_sel(str)）
③ アクセプター（A）原子（acceptors_sel(str)）
④ ドナーと水素原子（D-H）のカットオフ距離（d_h_cutoff(float)、デフォルト 1.2）
⑤ ドナーとアクセプター（D-A）のカットオフ距離（d_a_cutoff(float)、デフォルト 3.0）
⑥ ドナー - 水素 - アクセプター（D-H-A）の角度（d_h_a_angle_cutoff、デフォルト 150）

水素結合の定義そのものという感じですね。

それではドナーを「Met 793（→resid 93）の主鎖N原子」、アクセプターを「リガンド（resname LIG）のN4」として解析してみましょう。

from MDAnalysis.analysis.hydrogenbonds import hbond_analysis

# ① ドナー原子
donor_MetN = "resid 93 and name N"
# ② 水素原子
donor_MetH = "resid 93 and name H"
# ③ アクセプター原子
acceptor_LigN = "resname LIG and name N4"

hb_analyzer = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_MetN,
                                              hydrogens_sel = donor_MetH, 
                                              acceptors_sel = acceptor_LigN,
                                              d_h_cutoff = 1.2,
                                              d_a_cutoff = 3.0,
                                              d_h_a_angle_cutoff = 150)
hb_analyzer.run()

結果はresults.hbondsにnumpy.ndarrayとして格納されていて「[<frame>, <donor atom id>, <hydrogen atom id>, <acceptor atom id>, <distance>, <angle>] 」となっています。

DataFrameに変換しましょう。

column_names = ["frame", "donor_index", "hydrogen_index", "acceptor_index", "distance", "angle"]
hb_anal_df =pd.DataFrame(np.round(hb_analyzer.results.hbonds, 2), columns = column_names)
hb_anal_df.set_index("frame")

全50フレームのうち4フレームでのみ（デフォルトの）水素結合形成条件を満たしました。1割以下なのでProLIFで検出されなかったのも仕方ないかもしれませんね。

このままプロットしてもつまらないのでカットオフを緩めて再解析してみましょう（d_a_cutoff = 3.5）。解析結果のD-A距離とD-H-A角度は一つのプロットにまとめて描画します。

# ドナー - アクセプター カットオフ距離を3.5にして再解析
hb_analyzer_35 = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_MetN,
                                              hydrogens_sel = donor_MetH, 
                                              acceptors_sel = acceptor_LigN,
                                              d_h_cutoff = 1.2,
                                              d_a_cutoff = 3.5,
                                              d_h_a_angle_cutoff = 150)
hb_analyzer_35.run()

# DataFrame化
hb_anal_df_35 =pd.DataFrame(np.round(hb_analyzer_35.results.hbonds, 2), columns = column_names)
hb_anal_df_35.set_index("frame")

# プロットを作成
fig = plt.figure(figsize=(10, 6))

dist_plot = hb_anal_df_35["distance"].plot(style="o", x="frame")
dist_plot.set_ylabel("distance (Å)")
dist_plot.set_xlabel("frame")
dist_plot.set_ylim(0.0, 4.0)

angle_plot = hb_anal_df_35["angle"].plot(secondary_y=True, style="*", x="frame")
angle_plot .set_ylabel("angle (°)")
angle_plot .set_xlabel("frame")
angle_plot .set_ylim(0.0, 180)

fig.legend(loc="lower right", bbox_to_anchor=(1, 0), bbox_transform=dist_plot.transAxes)
plt.show()

D-A距離（distance）とD-H-A角度（angle）は一見連動して動いているように見えますが所々逆の動き（short length & wide angle）になっているのが面白いですね。

TeachOpenCADD T20では同様の解析を経て、距離のばらつきに対して角度は180度近くで安定化していることからリガンドの結合に重要な相互作用となっているのではないかと考察されていました。・・・なるほど、そう解釈するのか！

2-6. MDAnalysisで水素結合の解析-2

MDAnalysisによる水素結合解析の方法がわかったので、ProLIFで見つかった「ASP 855（→resid 155）の側鎖OD1/OD2原子」と「リガンド（resname LIG）のO3」との相互作用も解析してみましょう。

コードは同じことの繰り返しですが、アスパラギン酸側鎖のカルボキシ基には酸素原子が２つあるので両者で違いがでるか？興味があります。

水素結合のカットオフを緩めて解析してみました（d_a_cutoff = 4.0, d_h_a_angle_cutoff = 120）。

# ① ドナー原子
donor_LigO = "resname LIG and name O3"
# ② 水素原子
donor_LigH = "resname LIG and name H22"
# ③ アクセプター原子−１
acceptor_AspOD1 = "resid 155 and name OD1"
# ③ アクセプター原子−２
acceptor_AspOD2 = "resid 155 and name OD2"

# Asp OD1原子について解析
hb_analyzer_OD1 = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_LigO, hydrogens_sel = donor_LigH, 
                                              acceptors_sel = acceptor_AspOD1,
                                              d_h_cutoff = 1.2, d_a_cutoff = 4.0, d_h_a_angle_cutoff = 120)
hb_analyzer_OD1.run()

# Asp OD1原子の解析結果データフレーム
hb_anal_OD1_df =pd.DataFrame(np.round(hb_analyzer_OD1.results.hbonds, 2), columns = column_names)

# Asp OD2原子について解析
hb_analyzer_OD2 = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_LigO, hydrogens_sel = donor_LigH, 
                                              acceptors_sel = acceptor_AspOD2,
                                              d_h_cutoff = 1.2, d_a_cutoff = 4.0, d_h_a_angle_cutoff = 120)
hb_analyzer_OD2.run()

# Asp OD2原子の解析結果データフレーム
hb_anal_OD2_df =pd.DataFrame(np.round(hb_analyzer_OD2.results.hbonds, 2), columns = column_names)

まずはD-A距離についてAspの2つの酸素原子（OD1/OD2）の結果を一緒にプロットしてみましょう。

# D-A距離についてOD1/OD2を同時にプロット
plt.scatter(hb_anal_OD1_df["frame"], hb_anal_OD1_df["distance"], color = 'red', marker="o", label="with OD1")
plt.scatter(hb_anal_OD2_df["frame"], hb_anal_OD2_df["distance"], color = 'blue', marker="*", label="with OD2")

plt.xlabel("frame")
plt.ylabel("distance (Å)")
plt.ylim(0.0, 4.0)

plt.legend()
plt.show()

Asp855側鎖のカルボキシ基のうち、OD2原子（red circle）はリガンド水酸基と4.0Åの距離に入ったのは3/50フレームのみでいずれもギリギリ4.0以下でした。

一方で、OD1原子（blue star）は26/50フレームで4.0Å以下の距離に入り、3.0Å以下の点も多くみられます。

同様にD-H-A角度についてもプロットしてみましょう。

# D-H-A角度についてOD1/OD2を同時にプロット
plt.scatter(hb_anal_OD1_df["frame"], hb_anal_OD1_df["angle"], color = 'red', marker="o", label="with OD1")
plt.scatter(hb_anal_OD2_df["frame"], hb_anal_OD2_df["angle"], color = 'blue', marker="*", label="with OD2")

plt.xlabel("frame")
plt.ylabel("angle (°)")
plt.ylim(0, 180)

plt.legend()
plt.show()

角度についても、OD2原子（red circle）は3フレームとも150°よりも小さく、水素結合の基準を満たしそうにありませんでした。

また、OD1原子（blue star）も半数近くは150°以下となっており、距離、角度ともに基準を満たすものは10フレーム前後（2割程度）に落ち着きそう、という結果になりました。

以上の比較から、今回の複合体トラジェクトリではAsp855側鎖カルボキシ基の一方の酸素原子(OD2)とは水素結合を形成するものの、もう一方の酸素原子（OD1）に結合が切り替わることは観察されませんでした。

もう少しつっこむと、30フレーム（600 ps）以降では距離、角度ともに（デフォルトの）水素結合要件を満たすものはありませんでした。つまり、1nsの短いシミュレーションの初期でたまたま見られた偶然の結合の可能性もあります。

ご参考までに計算対象の元のPDB（3POZ）の構造を貼っておきます。

PyMolのデフォルトではリガンド（TAK-285）の水酸基はArg841の主鎖カルボニル酸素およびAsn842の側鎖アミドカルボニル酸素と相互作用しています。一方、タンパク（EGFR）のAsp855はK745と水分子を介したネットワークを形成しているようです。

この共結晶構造解析だけからリガンド OH基とタンパクのAsp855が水素結合を形成する可能性を考慮するのは難しそうです。今回のMDシミュレーションで見られた相互作用をランダムに与えられた初速の産んだ偶然と捉えるか？、X線共結晶では捉えられなかった相互作用と捉えるか？、繰り返しのMDシミュレーションやより長時間のシミュレーションを行えばわかるのかもしれませんが、私の手には負えません。。。。

3. おわりに

以上、今回は低分子ータンパク複合体のMDシミュレーションについてトラジェクトリ解析を試してみました。

特に（イメージしやすい古典的な）水素結合形成の解析に焦点を当てて、ProLIFやMDAnalysisを使った解析を追いかけてみました。・・・といってもほとんどTeachOpenCADD T20とMaking it rainのコードのCopy&Pasteですが。。。

個人的にはProLIFの解析で当初想定していたリガンド TAK-285とタンパク Hinge領域との相互作用が検出できなかったものの、もともとの共結晶構造（PDB ID：3POZ）を眺めるだけでは想像し難いリガンド OH基とタンパク D855 側鎖カルボキシ基との間に水素結合を形成する可能性が示唆された点が非常に興味深かったです。

・・・といっても、MD計算で初期に与えるランダムな速度分布の偶然の結果かもしれません。

より突っ込んだ解析、考察は門外漢の私には手に余るので専門家の皆様にお任せします（丸投げ）。

「MD計算の解釈で元の構造の分解能をどう考慮すれば良いか？」など、色々とわからないことが多かったです。今回も間違いが多いと思うのでご指摘をいただければ幸いです。

ではでは！

*1:PDB ID: 2TIXのリガンドは正確にはATPそのものではなくアナログの「Phosphoaminophosphonic acid-adenylate ester（ANP）」と呼ばれるものです。あくまでイメージとしてご容赦ください。

*2:Structure Based Drug DesignでTAK-285をベースにした構造から展開し、T790M/L858R変異に対しても活性を有する化合物を得る文献もありました。ACS Med Chem Lett. 2012(4)201。PubMedで読めます。

*3:リガンドの名前（LIG）は前回シミュレーション系構築の際につけました

*4:PLIPの文献はこちら Melissa F Adasme, Katja L Linnemann, Sarah Naomi Bolz, Florian Kaiser, Sebastian Salentin, V Joachim Haupt, Michael Schroeder, PLIP 2021: expanding the scope of the protein–ligand interaction profiler to DNA and RNA, Nucleic Acids Research, Volume 49, Issue W1, 2 July 2021, Pages W530–W534,

*5:なお、ProLIFはトラジェクトリ以外のデータも解析できます。Fingerprintオブジェクトを作成後、run()ではなくgenerate()を使えば１対１、run_from_iterable()を使えば複数を対象にした解析ができます。ドッキングの結果の複合体に対して適用したりすることができるようです。参考：Working with docking poses instead of MD simulations

2022-05-02

OpenMMをステップバイステップで　〜Part 8 低分子-タンパク複合体シミュレーション〜

OpenMMの使い方を順番にたどる記事のPart 8です。

Part 7まででタンパク質単体のMDシミュレーションについて、「①準備（前処理）、②実行、③結果の解析」の流れがなんとなくわかりました。第１ステージはクリア！…ということにして、第２ステージ「低分子 - タンパク質複合体」のシミュレーションを試したいと思います。

この記事ではTeachOpenCADD T019：Molecular dynamics simulationで行われている「複合体シミュレーション系の構築と計算」の理解を目指します。題材にする構造はPDB ID: 3POZでEGFRタンパクとリガンドの複合体です。

低分子リガンドのための力場としてopen forcefieldを利用し、OpenMMのシミュレーション系を構築します。計算はGoogle Colabで行いました。*1

1. 低分子を含むMDシミュレーションに向けた準備
2. 計算の対象構造について
3. タンパク質-低分子複合体の系をつくろう！
4. 計算条件の設定と実行
5. おわりに
6. おまけ

1. 低分子を含むMDシミュレーションに向けた準備

1-1. タンパク質単独と何が違う？

最初にタンパク質単独と低分子-タンパク質複合体のシミュレーション準備で違う点です。

一番大きな違いは力場です。タンパク質の構造は基本的に（限られた種類の）アミノ酸からなりますが、低分子リガンドはより多様な原子タイプ・結合パターンを含んでいます。なのでより広い領域をカバーした力場パラメータが計算に求められます。

OpenMMにデフォルトで含まれる力場パラメータは生体分子（タンパク質、DNA、RNA ...）や水モデルには対応していますが、低分子には対応していません。そこで別に低分子の力場を作成して取り込んでやる必要があります。

ここではOpenMM Forcefieldsを経由してopen forcefieldの力場を利用します。

1-2. OpenMM Forcefields?

OpenMM Forcefields はOpenMMの追加の力場を提供しているツールです（GitHub）。

サポートしている力場は① AMBER、② CHARMM、③ Open Force Field Initiative force fields です。

低分子についてはAMBERのGeneral AMBER Force Field（GAFF）やOpen Force Fieldのopenff、smirnoffといった力場が使えるようになります。

インストールはcondaでOKです。

conda install -c conda-forge openmmforcefields

1-3. Open Force Field Initiative？

似たような名前が多くて混乱しますが、Open Force Field InitiativeはOpenMMとは別の分子力場に焦点を当てたプロジェクトです。以前こちらのブログ記事でも取り上げました。

継続的な分子力場の改善を目指したイニシアティブで、基になるデータ管理から力場の作成までをオープンかつ再現性よく行うため、自動化した仕組みづくりを目指しているそうです。

現在（2022/05）、以下の３系統の力場が提供されています。

左から右に新しくなります。 SMIRNOFF（スミノフ？）はSMIRKS-native Open Force Fieldの略で、「OpenFF 1.y.x」系統のコードネームがParsley（パセリ）、「OpenFF 2.y.x」系統のコードネームがSage（セージ）です。*2

GitHubへのリンクはそれぞれ以下です。

smirnoff99Frosst

Parsely (OpenFF 1.3.1)

Sage (OpenFF 2.0.0)

利用するにはOpen Force Field ToolkitをインストールすればOKです。condaがあれば簡単です。

conda install -c conda-forge openff-toolkit

使える力場パラメータのバージョンは以下で確認できます。

from openff.toolkit.typing.engines.smirnoff import forcefield
print(forcefield.get_available_force_fields())

#  ['smirnoff99Frosst-1.0.2.offxml', 'smirnoff99Frosst-1.0.0.offxml', 'smirnoff99Frosst-1.1.0.offxml', 'smirnoff99Frosst-1.0.4.offxml', 'smirnoff99Frosst-1.0.8.offxml', 'smirnoff99Frosst-1.0.6.offxml', 'smirnoff99Frosst-1.0.3.offxml', 'smirnoff99Frosst-1.0.1.offxml', 'smirnoff99Frosst-1.0.5.offxml', 'smirnoff99Frosst-1.0.9.offxml', 'smirnoff99Frosst-1.0.7.offxml', 'openff-1.0.1.offxml', 'openff-1.1.1.offxml', 'openff-1.0.0-RC1.offxml', 'openff-1.2.0.offxml', 'openff-1.3.0.offxml', 'openff_unconstrained-2.0.0-rc.1.offxml', 'openff_unconstrained-1.3.1.offxml', 'openff_unconstrained-1.2.1.offxml', 'openff-2.0.0-rc.2.offxml', 'openff_unconstrained-1.0.0-RC2.offxml', 'openff_unconstrained-1.1.0.offxml', 'openff_unconstrained-1.0.0.offxml', 'openff-2.0.0.offxml', 'openff_unconstrained-2.0.0.offxml', 'openff_unconstrained-2.0.0-rc.2.offxml', 'openff-1.1.0.offxml', 'openff-1.0.0.offxml', 'openff-1.0.0-RC2.offxml', 'openff-1.3.1.offxml', 'openff-1.2.1.offxml', 'openff-1.3.1-alpha.1.offxml', 'openff_unconstrained-1.0.0-RC1.offxml', 'openff_unconstrained-1.2.0.offxml', 'openff_unconstrained-1.3.0.offxml', 'openff-2.0.0-rc.1.offxml', 'openff_unconstrained-1.0.1.offxml', 'openff_unconstrained-1.1.1.offxml', 'openff_unconstrained-1.3.1-alpha.1.offxml']

ファイル形式はoffxmlのようです。OpenFFのxmlという感じでしょうか？

OpenFF toolkitがあるとRDKitのMolオブジェクトをOpenFFで使えるように変えたりすることもできます。

1-4. MDTrajのインストール

準備の最後にMDTrajをインストールしておきます。

こちらはTrajとついている通りトラジェクトリ解析用のツールです。前回の記事でご紹介した埼玉大学　松永先生のコロキウム（生体分子系の分子動力学シミュレーションデータの解析入門）でも、ファイルの拡張子変換機能（mdconvert）が使いやすいと紹介されていました。

今回の記事では、別々に前処理したタンパク質と低分子リガンドを一つの複合体にまとめる（マージする）ツールとして利用します。

インストールはcondaでOKです。

conda install -c conda-forge mdtraj

これでツールの準備はおしまいです。

（省略しましたがRDKitも使います。*3）

2. 計算の対象構造について

具体的なシミュレーションの準備を行う前に、計算の対象の構造を確認しておきます。

今回計算の対象とする低分子-タンパク質複合体の構造はTeachOpenCADD T019：Molecular dynamics simulationを踏襲してPDB ID: 3POZです。

EGFRタンパクのキナーゼドメインとリガンド TAK-285のX線共結晶構造です。*4

文献はこちら

Aertgeerts K, Skene R, Yano J, et.al., Structural analysis of the mechanism of inhibition and allosteric activation of the kinase domain of HER2 protein. J Biol Chem. 2011(286)18756-65. doi:10.1074/jbc.M110.20619351133-2/fulltext)

複合体構造の中でのTAK-285のLigand IDは03Pです。あとでPDBファイルから構造を取り出す際に使います。*5

3. タンパク質-低分子複合体の系をつくろう！

それでは具体的な処理を進めていきましょう！

この記事では操作の理解を重視するため、TeachOpenCADDでは関数にまとめて定義してある部分を解体して眺めていきます。

3-1. 処理全体の流れ

最初に処理の流れを確認します。以前の記事（part 5 シミュレーションプログラムの中身を理解しよう）で取り上げましたが、OpenMMのシミュレーションプログラムはざっくり以下のようなクラス構成です。

今回、低分子リガンドを加えるにあたって処理を増やすのは左半分です。

追加する処理を加えてみました（水色）。

タンパク質-リガンド複合体のPDBファイルから「①リガンド部分を抜き出し」て、「② Open Forcefieldで力場を設定」し、「③ タンパクとマージ」して複合体に戻します。

より具体的に、処理と使うツールは以下のようになります。

PDBファイルには低分子リガンドの結合次数や水素原子の情報が含まれていないので、これらを補完するためにRDKitを使います。Open Forcefieldを使って低分子の力場を補充し、タンパク-低分子の複合体の系を構築し、最後にシミュレーションの箱（水、イオン追加）を用意して設定完了です。

3-2. Step 1：タンパク質の前処理

Step ①はタンパク質部分の前処理です。ここで行う作業は以前の記事（「Part 2 PDBFixerでタンパク質の前処理」）で取り上げた内容とほぼ同じです。

PDBFixerを使って以下のような作業をします。

水やイオンの追加はここではまだ行わないことにご注意ください*6。

こんな感じのコードです。*7

# 処理したいPDBファイルの指定
pdb_file = 'Data/3poz.pdb'

# 必要なライブラリのインポート
import pdbfixer

# PDBFixerでPDBファイルを読み込む
fixer = pdbfixer.PDBFixer(pdb_file)

# 不要な構造の削除
fixer.removeHeterogens()

# 欠けている残基のチェック（欠損原子の確認のためにも必要）
fixer.findMissingResidues()

# タンパク質末端の欠けている残基を取り除く処理
chains = list(fixer.topology.chains())
keys = fixer.missingResidues.keys()
for key in list(keys):
    chain = chains[key[0]]
    if key[1] == 0 or key[1] == len(list(chain.residues())):
        del fixer.missingResidues[key]

# 非標準な残基が含まれているか確認、あれば標準的なものに置き換える
fixer.findNonstandardResidues()
fixer.replaceNonstandardResidues()

# 欠けている原子の確認、あれば追加する
fixer.findMissingAtoms()
fixer.addMissingAtoms()

# 水素原子の付与（pHを設定する）
ph = 7.0
fixer.addMissingHydrogens(ph)

print(type(fixer))
#  <class 'pdbfixer.pdbfixer.PDBFixer'>

得られたPDBFixerのオブジェクトには処理後のタンパク質の原子の座標（positions）とトポロジー(topology)が含まれています。

3-3. Step 2：リガンド構造の修正

続いてリガンド構造の準備です。

PDBから構造や座標を取り出しますが、結合次数や水素原子の情報がPDBファイルには含まれていません。そこでこれらの情報をRDKitを利用して追加します。

まず、RDKitでPDBファイルを読み込んで残基ごとにスプリット（SplitMolByPDBResidues）し、リガンド（ID：03P）を抜き出します。

from rdkit import Chem
from rdkit.Chem import AllChem

# PDBファイルを読み込んで残基ごとにスプリット
rdkit_pdb = Chem.MolFromPDBFile('Data/3poz.pdb')
rdkit_pdb_split = Chem.SplitMolByPDBResidues(rdkit_pdb)

# 残基 03Pの抜きだし
pdb_tak285 = rdkit_pdb_split["03P"]

Chem.Draw.MolToImage(pdb_tak285)

単結合だけだと寂しい感じ・・・。

SMILES情報からTAK-285のリファレンス構造を作成し、リファレンスをもとに結合次数を追加しましょう（AssignBondOrdersFromTemplate]）。

SMILESはRCSB PDBのリガンドページ（03P）から持ってきました。

# 中途半端に残ってる水素原子を消しておく
pdb_tak285_noH = Chem.RemoveHs(pdb_tak285)

# SMILESからリファレンス構造を作成
smi_tak285 = "CC(C)(CC(=O)NCCn1ccc2c1c(ncn2)Nc3ccc(c(c3)Cl)Oc4cccc(c4)C(F)(F)F)O"
ref_tak285 = Chem.MolFromSmiles(smi_tak285)

# リファレンスをもとに結合次数をアサイン
prep_tak285 = AllChem.AssignBondOrdersFromTemplate(ref_tak285, pdb_tak285)
prep_tak285.AddConformer(pdb_tak285.GetConformer(0))

# 水素原子の付加
prep_tak285 = Chem.AddHs(prep_tak285, addCoords=True)

Chem.Draw.MolToImage(prep_tak285)

無事、結合次数のアサインと水素原子の付加ができました！*8

3-4. Step 3：リガンド用の力場の準備

リガンド構造が用意できたので力場を準備します。

最初にOpenMMでタンパク質（amber14-all.xml）と水モデル（amber14/tip3pfb.xml）の力場を設定しておきます。

from openmm import *
from openmm.app import *

FF = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')

このForceFieldに低分子の力場を追加します。

まず、リガンドを「RDKitのMolオブジェクト」から「OpenFFのMoleculeオブジェクト」に変更します。openff toolkitのMolecule.from_rdkitを使えばOKです。*9

from openff.toolkit.topology import Molecule

off_tak285 = Molecule.from_rdkit(prep_tak285, allow_undefined_stereo=True, hydrogens_are_explicit=True)

print(type(off_tak285))
#  <class 'openff.toolkit.topology.molecule.Molecule'>

OpenMMのForceFieldに低分子の力場を追加するには、低分子の構造から① Templateを作ってそれを②ForceFieldにregisterします。

今回はOpen Forcefieldに用意されているGAFFを力場に使います。この場合は「① GAFFTemplateGenerator」と「② registerTemplateGenerator」を使うことになります。

from openmmforcefields.generators import GAFFTemplateGenerator

gaff_temp = GAFFTemplateGenerator(off_tak285)
FF.registerTemplateGenerator(gaff_temp.generator)

なお、GAFFはAMBER力場の仲間で、低分子医薬品デザインのために作られたものだそうです。タンパク質側をAMBER力場で設定しているときに互換性の良い低分子力場となる感じでしょうか？

今回使用したGAFFのバージョンや、Open forcefieldにインストールされているGAFFのバージョン一覧は以下のように確認できます。

print("version info : ", gaff_temp.gaff_version)
print("installed versions : ", GAFFTemplateGenerator.INSTALLED_FORCEFIELDS)

# version info :  2.11
# installed versions :  ['gaff-1.4', 'gaff-1.8', 'gaff-1.81', 'gaff-2.1', 'gaff-2.11']

使われたバージョン 2.11のようです。

これで力場の準備が終わりました。

3-5. Step 4：リガンドのオブジェクト変換

つづいてリガンドのオブジェクトを変換します。タンパク質とマージするための準備です。

力場を作る際にRDKitからOpenFFのMoleculeオブジェクトにしましたが、これをさらにOpenMMのオブジェクトにします。

まず、OpenFF Moleculeオブジェクトに、分子および原子の名前を追加します。原子の名前は「元素記号 + 数字」形式です。

# 分子の名前（"LIG"）を追加
off_tak285.name = "LIG"

# 原子の名前を追加
# 辞書を用意しておいて、元素記号＋数字形式にする
element_counter_dict = {}
for off_atom, rdkit_atom in zip(off_tak285.atoms, prep_tak285.GetAtoms()):
    element = rdkit_atom.GetSymbol()
    if element in element_counter_dict.keys():
        element_counter_dict[element] += 1
    else:
        element_counter_dict[element] = 1
    off_atom.name = element + str(element_counter_dict[element])

RDKitのGetSymbolで原子の元素記号を一つずつ取り出して辞書に格納しています。同じ元素記号があったら数字を+１するので、炭素なら「C1、C2・・・」と順番に原子に名前がつけられていきます。

名前をつけたらOpenFFからOpenMM（topology、positions）に変換します。ひとつ注意点として、OpenMMの単位はnm基準なので、PDBファイルのオングストロームから変更する必要があります。

# OpenMMの単位系導入
from openmm import unit

# トポロジー（topology）の取り出しと変換
off_top = off_tak285.to_topology()
omm_top = off_top.to_openmm()

# 座標（positions）の取り出し
omm_pos = off_tak285.conformers[0]

# 単位をÅからnmへ
for atom in omm_pos:
    coords = atom / atom.unit
    atom = (coords / 10.0) * unit.nanometers

# TopologyとPositionsをまとめてOpenMMのオブジェクトに変換
omm_tak285 = app.Modeller(omm_top, omm_pos)

最後にtopologyとpositionsをまとめてOpenMMのModellerオブジェクトにしています。

これでリガンド側での準備はもろもろ完了です。タンパク質と合体しましょう！

3-6. Step 5：リガンドとタンパク質のマージ

ここまででタンパク質のPDBFixerオブジェクトとリガンドのModellerオブジェクトが用意できました。これらをまとめて一つにします。

topologyとpositionsを別々に処理します。

まずtopologyですが、MDTrajを利用します。タンパク質とリガンドをそれぞれMDTrajのオブジェクトに変換（Topology.from_openmm）したのちjoinで結合し、OpenMMオブジェクトに戻します（to_openmm）。

import mdtraj as md

# タンパク質とリガンドのトポロジーをそれぞれMdTrajに変換
md_protein_top = md.Topology.from_openmm(fixer.topology) 
md_ligand_top = md.Topology.from_openmm(omm_tak285.topology) 

# joinでひとつのトポロジーにまとめる
md_complex_top = md_protein_top.join(md_ligand_top)

# OpenMMオブジェクトにする
omm_complex_top = md_complex_top.to_openmm()

ついで、positionsです。

各原子の位置の値と単位（nm）を含むタプルの配列を作成します。OpenMMのQunatityクラスとなります。

import numpy as np

# タンパクとリガンドを合わせた総数
total_atoms = len(fixer.positions) + len(omm_tak285.positions)

# positionsを格納するQuantity配列の準備
complex_positions = unit.Quantity(np.zeros([total_atoms, 3]), unit=unit.nanometers)

# タンパク質のpositionsを追加
complex_positions[: len(fixer.positions)] = fixer.positions
# 後に続けてリガンドのpositionsを追加
complex_positions[len(fixer.positions) :] = omm_tak285.positions

topologyとpositionsがそれぞれ処理できたので、これらをまとめたOpenMMのオブジェクトにしましょう。

comp_model = app.Modeller(omm_complex_top, complex_positions)

できました！

3-7. Step 6：水、イオンを追加し系を構築

前処理の最後として系に水とイオンを追加します。

水のpaddingの厚みは「1.0 nm」、イオンは「NaCl 0.15 M」とします。水モデルはTIP3Pで、forcefiled（amber14/tip3pfb.xml）で指定されています。

comp_model.addSolvent(FF, padding=1.0 * unit.nanometers, ionicStrength=0.15 * unit.molar)

処理後の状態をPDBファイルで保存しておきましょう。

# 処理後の状態（トポロジー、原子の位置）をPDBファイルで出力
top = comp_model.getTopology()
pos = comp_model.getPositions()
app.PDBFile.writeFile(top, pos, open('Data/gaff_complex_processed.pdb', 'w'))

こんな構造になりました。

4. 計算条件の設定と実行

4-1. 計算条件の設定

計算系の準備が完了したので、シミュレーションを実行するための設定を行いましょう。

タンパク質単独でシミュレーションした過去記事（「Part 3」）と同様の設定にしたいと思います。

こんな感じです。

コードにするとこんな感じです。

from openmm import *
from openmm.app import *
from openmm.unit import *

# シミュレーションの準備  
## ForceFiledからSystemを構築する
system = FF.createSystem(comp_model.topology, nonbondedMethod=PME, nonbondedCutoff=1.0*nanometers, constraints = HBonds, rigidWater=True, ewaldErrorTolerance=0.0005)

## 圧力制御の機能をaddForceする
system.addForce(MonteCarloBarostat(1.0*atmospheres, 300*kelvin, 25))

## 積分計算の設定
integrator = LangevinMiddleIntegrator(300*kelvin, 1.0/picosecond, 0.002*picoseconds)
integrator.setConstraintTolerance(0.000001)

## 計算環境をPlatformで設定
platform = Platform.getPlatformByName('CUDA')
platformProperties = {'Precision': 'single'}

## simulationの構築
simulation = Simulation(comp_model.topology, system, integrator, platform, platformProperties)
simulation.context.setPositions(comp_model.positions)

# XML結果出力ための設定
with open("system.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(system))
with open("integrator.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(integrator))

# エネルギー最小化の実行
print('Performing energy minimization...')
simulation.minimizeEnergy()

# 平衡化のステップ
print('Equilibrating...')
simulation.context.setVelocitiesToTemperature(300*kelvin)
simulation.step(50000)

# シミュレーション本番の実行の前にreporterオブジェクトを追加する
print('Simulating...')

dcdReporter = DCDReporter('trajectory.dcd', 10000)
dataReporter = StateDataReporter('log.txt', 1000, totalSteps=500000, step=True, speed=True, progress=True, potentialEnergy=True, temperature=True, separator='\t')
checkpointReporter = CheckpointReporter('checkpoint.chk', 10000)

simulation.reporters.append(dcdReporter)
simulation.reporters.append(dataReporter)
simulation.reporters.append(checkpointReporter)
simulation.currentStep = 0

# シミュレーション本番の実行
simulation.step(500000)

# シミュレーション最後の状態をXMLで書き出す
simulation.saveState("final_state.xml")

# 最後の状態をPDBx/mmcifで出力する設定
state = simulation.context.getState(getPositions=True, enforcePeriodicBox=system.usesPeriodicBoundaryConditions())
with open("final_state.pdbx", mode="w") as file:
    PDBxFile.writeFile(simulation.topology, state.getPositions(), file)

項目が多いので長いですが、前回同様にメインは「①系の構築　→ ② エネルギー最小化 → ③ 平衡化のステップ → ④ トラジェクトリのプロダクションステップ」を行っているだけで、間に結果出力（Reporter）の設定が入っています。

それぞれの操作の意味については以前の記事（「Part 5：シミュレーションプログラムの中身を理解しよう」）をご参照いただければと思います。

4-2. TeachOpenCADDの設定との違い

以上の計算条件の設定は基本的にはTeachOpenCADD T19と同じですが以下の点が違います。

① アンサンプル：TeachOpenCADDは圧力制御を加えていない
② 非結合相互作用の拘束条件（constraints、rigidWater）
③ シミュレーション時間（平衡化の有無、プロダクションの時間）
④ 出力ファイル形式の設定

また、TeachOpenCADDでは変数としてforcefieldやmodellerなどを使っています。これらは同じ名前のメソッドがOpenMMに含まれている場合があり、最初のimportで読み込むmoduleを変えると誤認識される恐れがあります。重複のない変数名に変えた方が無難だともいます。

4-3. GoogleColabで実行する際の注意点

以上の計算条件の設定を実行すればMDシミュレーションが行えるはずですが、1nsの計算をするのは大変です。前回（「Part 4：ColabでOpenMMシミュレーションを実行」）と同様に、今回も計算はGoogle Colabで実行しました。

ひとつ注意点として、前回は「①minicondaをインストール」し、そこに「② condaで色々なソフトウェアをインストール」する手順を取りました。

ですが、今回この方法ではRDKitがうまく動きませんでした*10。

回避策としてTeachOpenCADDを参考にcondacolabを利用しました。

こんな感じ。

!pip install condacolab
import condacolab
condacolab.install()

!conda install -q -y -c conda-forge mdtraj openmm openmmforcefields openff-toolkit pdbfixer rdkit

condacolab.install()を実行すると一度「kernel restart」しますが、気にせず次の「!conda install ~~~」に進んでOKです。

より良い方法をご存知の方教えていただけると嬉しいです。

4-4. シミュレーション動画

結果はこんな動画になりました！

めっちゃ動くやん！

解析は次回！

5. おわりに

以上、今回はMDシミュレーションお勉強の第２ステージということで「低分子 - タンパク質複合体」計算系の設定と実行を行ってみました。

…といってもTeachOpenCADD T19の内容をひとつひとつ確かめただけですが、、、余計にわかりにくくなっていたらごめんなさい。

タンパク質単独のシミュレーションと違って、複合体では低分子リガンドの構造情報の取り扱いや専用の力場の設定が必要になりました。そのため複数のツールを組み合わせて行ったり来たりしなければならず難しかったです。

また、低分子力場を扱うためにOpen Forcefieldを使いました。本当はGAFFではなくSMIRNOFFやOpenFFを使ってみたかったのですが、私の能力では無理でした。

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

ではでは！

6. おまけ

低分子の力場としてGAFFではなくOpenFFの利用を試みるとこうなりました。

# 力場の設定・追加
from openmmforcefields.generators import SMIRNOFFTemplateGenerator

smirnoff_temp = SMIRNOFFTemplateGenerator(off_tak285)
forcefield.registerTemplateGenerator(smirnoff_temp.generator)
print("forcefield version : ", smirnoff_temp.forcefield)
# forcefield version :  openff-2.0.0

最新のOpenFF 2.00（Sage）が使われました。

これを使って進めるとOpenMMで系に「水、イオンを追加（addSolvent）」する際に以下のような出力が出ました。

# 〜〜〜略〜〜〜
<Improper class1="[H][O][C]([C]([H])([H])[H])([C]([H])([H])[H])[C]([H])([H])[C](=[O])[N]([H])[C]([H])([H])[C]([H])([H])[N]1[C]([H])=[C]([H])[c]2[n][c]([H])[n][c]([N]([H])[c]3[c]([H])[c]([H])[c]([O][c]4[c]([H])[c]([H])[c]([H])[c]([C]([F])([F])[F])[c]4[H])[c]([Cl])[c]3[H])[c]21$C22#29" class2="[H][O][C]([C]([H])([H])[H])([C]([H])([H])[H])[C]([H])([H])[C](=[O])[N]([H])[C]([H])([H])[C]([H])([H])[N]1[C]([H])=[C]([H])[c]2[n][c]([H])[n][c]([N]([H])[c]3[c]([H])[c]([H])[c]([O][c]4[c]([H])[c]([H])[c]([H])[c]([C]([F])([F])[F])[c]4[H])[c]([Cl])[c]3[H])[c]21$N5#32" class3="[H][O][C]([C]([H])([H])[H])([C]([H])([H])[H])[C]([H])([H])[C](=[O])[N]([H])[C]([H])([H])[C]([H])([H])[N]1[C]([H])=[C]([H])[c]2[n][c]([H])[n][c]([N]([H])[c]3[c]([H])[c]([H])[c]([O][c]4[c]([H])[c]([H])[c]([H])[c]([C]([F])([F])[F])[c]4[H])[c]([Cl])[c]3[H])[c]21$C20#27" class4="[H][O][C]([C]([H])([H])[H])([C]([H])([H])[H])[C]([H])([H])[C](=[O])[N]([H])[C]([H])([H])[C]([H])([H])[N]1[C]([H])=[C]([H])[c]2[n][c]([H])[n][c]([N]([H])[c]3[c]([H])[c]([H])[c]([O][c]4[c]([H])[c]([H])[c]([H])[c]([C]([F])([F])[F])[c]4[H])[c]([Cl])[c]3[H])[c]21$C23#30" periodicity1="2" phase1="3.141592653589793" k1="1.5341333333333336"/>

# 〜〜〜略〜〜〜

<Element Improper at 0x13f32bbc0>

# 〜〜〜略〜〜〜

エラーではないようですが、improperとなっているのが気になったので今回は採用しませんでした。

推測に過ぎませんが、OpenFFのSMIRNOFFはSMIRKSベースで２面角の力場を設定しているようなので、適切なクエリに該当しない場合improper判定されるのかもしれません。ご存知の方がいらっしゃったらご教示いただければ幸いです。

*1:ロゴは書くソフトウェアのWebサイトより。イラストはいらすとやさんより。

*2:セージで検索すると薬用サルビアが出てきました。ウォッカからハーブ系統に路線変更したんでしょうか？

*3:このブログを読んでくださっている方は既にインストール済みと思います。

*4:3DビューはRCSB PDB提供のNGL Viewer (AS Rose et al. (2018) NGL viewer: web-based molecular graphics for large complexes. Bioinformatics doi:10.1093/bioinformatics/bty419)

*5:PDBで「Ligand ID：03P」をで検索すると他に3RCD、3W2Oがヒットします。見比べても面白かもしれません。

*6:リガンドを入れる前に水を入れると後でぶつかるから？でしょうか？？

*7:詳しい中身は過去記事をご参照いただければと思います。

*8:TeachOpenCADDのコードに従って、結合次数を調整した後、再度AddConformerで元の3D座標を追加していますが、これが必要かはよくわかりません。この操作をしなくてもGetConformer().GetPositions()で確認すると座標は変わらなかったです。

*9:TeachOpenCADDでは引数にhydrogens_are_explicit=Trueが設定されていませんが、デフォルトではFalseで自動的にCham.AddHsを実行することになります。あらかじめ水素原子を付加しているのでTrueにしておきました。

*10:インストール後にimport Chemしようとしたら「必要なGLIBCXX（？）が見つからない」といったエラーでimport出来ませんでした

2022-04-30

OpenMMをステップバイステップで　〜Part 7 MDAnalysisでトラジェクトリ解析（②）〜

オープンソースのプログラムでMD計算をお勉強する記事のPart 7です。

今回は前回（part 6）に引き続き、MDAnalysisを使ったシミュレーション結果のトラジェクトリ解析（パート②）を行います。また、解析ツールとしてPyTraj も利用します。

使用するデータはEGFRタンパク質（PDB ID: 3POZ）についてOpenMMで1ns計算したものです。

1. この記事で試す解析
2. 主成分分析（PCA）
3. 相互相関解析（cross correlation）
4. おわりに

1. この記事で試す解析

前回の記事では①RMSD、②RMSF、③回転半径の計算を行いました。後半戦では④主成分分析と⑤相互相関解析を行います。

ざっくりと「タンパク質のどのあたりが一緒にうごいているか？」、おおまかな動きの特徴をつかめる解析のようです。

2. 主成分分析（PCA）

2-1. MD計算の主成分分析って？

まずは主成分分析です。統計っぽい話でよくみますがMD計算においてはどのように使われるのでしょうか？

前回の記事でも取り上げさせていただいた埼玉大学 松永康佑先生のコロキウム動画わかりやすく解説してくださっています（動画の33分あたりからです）。*1

www.youtube.com

こういう感じのお話？

そのままではごちゃごちゃ動いていてよくわからないのが、主成分分析することでノイズを減らして全体としての動きの傾向をつかみやすくなるようです。

上記動画の35分過ぎからは粗視化モデルをつかった主成分分析例が紹介されています。非常にわかりやすいのでおすすめです。

2-2. MDAnalysisで解析の準備（おさらい）

それではMDAnalysisを使った解析を試してみましょう！

前回のおさらいになりますがMDAnalysisはUniverseオブジェクトを中心にデータ解析を行います。topologyとtrajectoryを読み込んで構築します。

import MDAnalysis as mda

# topologyとtrajectoryからUniverseオブジェクトを構築
md_universe = mda.Universe("simulation_data/3poz-processed.pdb", "simulation_data/trajectory.dcd")
print(type(md_universe))

また、解析に先立ってトラジェクトリのタンパク質のアラインメントをとっておきます（念のため・・）。

from MDAnalysis.analysis import align

# トラジェクトリの位置を最初のフレームにリセットする
md_universe.trajectory[0]

# align.AlignTrajでアラインメントを取る
alignment = align.AlignTraj( mobile = md_universe, reference = md_universe, select="protein", in_memory=True)
alignment.run()

準備完了！

2-3. MDAnalysisで主成分分析

2-3-1. PCA計算の流れ

主成分分析のやり方については「MDAnalysis User Guide：Principal component analysis of a trajectory」に詳しいチュートリアルがあります。

MDAnalysis.analysis.pca.PCA()を使えば求めることができ、以下のような流れで計算するようです。

それではチュートリアルに従って順番に確認してみましょう。

2-3-2. PCAの実行と累積寄与率の確認

MDAnalysisでPCAを行うにはAnalysisモジュールのMDanalysis.analysis.pca.PCA()を使います。引数にUniverseオブジェクトと目的のAtomGroupの選択基準を渡します。

Cα原子を選択してPCAを行うとこんな感じ。

import MDAnalysis.analysis.pca as pca

# Cα原子のPCA
Ca_pca = pca.PCA(md_universe, select ='name CA')
Ca_pca.run()

PCAの結果の主成分はp_componentsに格納されていて、「原子数 x 3」の配列が成分数分入った配列です（デフォルトで実行した場合は成分数も「原子数 x 3」です。）。

# 主成分を含む配列のサイズ
components = Ca_pca.results.p_components
print(components.shape)

#  (951, 951)

シミュレーションの対象としているタンパク質の構造は317残基（Cα原子：317個）なのでちょうど３倍の951になりました。

分散（variance）および累積寄与率（cumulated_variance）を確認してみます。どちらもnumpy.ndarrayなのでインデックスを指定して取り出せます。

# 寄与率
print("1st varinace : ", Ca_pca.results.variance[0])
print("2nd varinace : ", Ca_pca.results.variance[1])

# 累積寄与率
print("1st cumulative varinace : ", Ca_pca.results.cumulated_variance[0])
print("2nd cumulative varinace : ", Ca_pca.results.cumulated_variance[1])

# 1st varinace :  63.373490214896336
# 2nd varinace :  14.112964500184088
# 1st cumulative varinace :  0.39144314860346985
# 2nd cumulative varinace :  0.4786156120633174

第１主成分と第２主成分を合わせた累積寄与率は48%でした。

2-3-3. 累積寄与率のプロット

累積寄与率が95% に達するには何番目まで考慮すれば良いでしょうか？

cumulated_varianceはndarrayなのでnp.where()で95より大きくなる要素のインデックスを取得し、その一番最初を取り出します。

import numpy as np
import matplotlib.pyplot as plt

# 累積寄与率が95%を超える成分
pc_95 = np.where(Ca_pca.results.cumulated_variance > 0.95)[0][0]
print(">95% component number : ", pc_95)

# 累積寄与率のプロット
plt.plot(Ca_pca.results.cumulated_variance[:pc_95])
plt.xlabel("Principal component")
plt.ylabel("Cumulative variance")
plt.show()

#  >95% component number :  35

第35主成分まで含めると95%を超えることがわかりました。少しずつ増えていく感じです。

2-3-4. トラジェクトリの次元削減とプロット

PCAを使ってトラジェクトリを次元削減してみましょう。transform()を使うと、指定した数（n_components）の主成分に変換されます。

戻り値はサイズ「（トラジェクトリのフレーム数, 主成分の数）」のndarrayです。

主成分の数を3としてCα原子について変換するとこんな感じです。

# CαについてPCAで次元削減
transformed = Ca_pca.transform(md_universe.select_atoms("name CA"), n_components=3)
print("shape of the array : ", transformed.shape)

# shape of the array :  (50, 3)

２次元にプロットして確認してみましょう。準備としてPandasのデータフレームに変換します。

import pandas as pd

# トラジェクトリの１フレームの時間(ps)
frame_time = md_universe.trajectory.dt

# データフレームに格納する
df = pd.DataFrame(transformed, columns=['PC1', 'PC2', 'PC3'])
df['Time (ps)'] = df.index * frame_time
df.head()

20ps 毎の各フレームで、主成分PC1〜PC3の重みに変換されました。

プロットしてみます。時間推移（ Time (ps)）で色づけしています。

import matplotlib.pyplot as plt
import seaborn as sns

pair_grid_plot = sns.PairGrid(df, hue='Time (ps)')
pair_grid_plot.map(plt.scatter, marker='.')

PC1とPC2の２次元プロットでは近い色の点が集まった意味ありげなパターンになりました。

・・・何の意味かはわかりません。

2-3-5. 主成分の動きを動画で

PCAを使ったトラジェクトリの基本的な動きの確認には、プロット以外にも動画で確認する方法もあるそうです。

元々のトラジェクトリを各主成分に投影します。

重み（w）は先のtransform()で得られていて、主成分の固有ベクトルはresults.p_components、平均座標はmeanでアクセスできます。上記式を実行しましょう。

# 第１主成分
pc_1 = Ca_pca.results.p_components[:, 0]

# 第１主成分についての重み
weight_1 = transformed[:, 0]

# 平均座標
mean_coord = Ca_pca.mean

# 投影
projected_1 = np.outer(weight_1, pc_1) + mean_coord.flatten()

# arrayの形を整える
coordinates = projected_1.reshape(len(weight_1), -1, 3)
print(coordinates.shape)

#  (50, 317, 3)

最後のreshapeでは、投影を計算後のarrayを「フレーム数、原子数（N）、xyz座標（3）」の次元となるように整えています。

各原子の座標の時間（フレーム）変化のデータが手に入ったので、これをもとにUniverseオブジェクトを作成します。

proj_univ = mda.Merge(md_universe.select_atoms("name CA"))
proj_univ.load_new(coordinates, order="fac")

#  <Universe with 317 atoms>

Merge()で既存のAtomGroupから新しいUniverseオブジェクトを作り、load_newで座標の情報を更新しています。

NGLviewを使ってトラジェクトリの動画を確認しましょう。

import nglview as nv

view = nv.show_mdanalysis(proj_univ.atoms)
view

全体が上下方向にクシャッととつぶれるような動きをピクピクしています。

ちなみにPCAで次元削減するまえの元々の動画はこんな感じです。

view = nv.show_mdanalysis(md_universe)
view

そのまま見ているよりもPCAの方が動きがシンプルでわかりやすいですね。

PCAのプロットで第１主成分とパターンをつっくていた第２主成分ではどうなるでしょうか？（コードは同じなので省略します。）

今度は左右方向に揺れている感じです。第１主成分と第２主成分で異なる動きの向きが取り出せていそうです。

楽しい！！・・・意味はわからないですが。

主成分分析は以上です。

3. 相互相関解析（cross correlation）

3-1. 相互相関？

続いて相互相関解析を行います。こういう感じの話です。

こちらも動きの解析ですが、PCAでは基本的な運動のモードを取り出す感じだったのに対し、分子内の運動の相関に焦点を当てているようです。

3-2. PyTrajのインストール

MDAnalysisにそれらしい項目がなかったので、Making it rainで解析に使われていた PyTraj を利用したいと思います。

PyTrajはAMBERのAmber Toolsに組み込まれているcpptrajというプログラムをPythonから扱えるようにしたフロントエンドパッケージだそうです

GitHub：PyTraj

ドキュメント：pytraj 2.0.2.dev0

インストールはcondaでできます（チャンネルはambermdです）*2。

conda install -c ambermd pytraj

ただし上記はPython バージョン3.9以上では使えません。Python 3.8以下の仮想環境を作るか、AmberToolsを丸ごとインストールすればPyTrajを使えます*3。

AmberToolsのインストールもcondaでOKです。

conda install -c conda-forge ambertools=21 compilers

こっちはPython 3.9の環境にも入りました。

3-3. PyTrajでトラジェクトリファイルの読み込み

PyTrajで解析を行うためにトラジェクトリを読み込みます。

load()にトラジェクトリファイル（ex. DCD）とトポロジーファイル（ex. PDB）を渡します*4。

import pytraj as pt

traj = pt.load("simulation_data/trajectory.dcd", "simulation_data/3poz-processed.pdb")
print(traj)

#  pytraj.Trajectory, 50 frames: 
#  Size: 0.064427 (GB)
#  <Topology: 57648 atoms, 17895 residues, 17881 mols, PBC with box type = cubic>

なお、PyTrajには他の解析ソフトとの相互利用を意識した機能もあり、MDAnalysisのUniverseオブジェクトからトラジェクトリを読み込むこともできます（PyTraj tutorial：Interface with mdanalysis）。

# MDAnalysisのUniverseオブジェクトから座標をnumpyのarrayに取り出す。
all_atoms = md_universe.select_atoms("all")
coords_from_mda = np.asarray([all_atoms.positions for frame in md_universe.trajectory], dtype='f8')

# 座標とPDBファイルからPyTrajに読み込む
traj_from_mda = pt.Trajectory(xyz=coords_from_mda, top= "simulation_data/3poz-processed.pdb")
print(traj_from_mda)

# pytraj.Trajectory, 50 frames: 
# Size: 0.064427 (GB)
# <Topology: 57648 atoms, 17895 residues, 17881 mols, non-PBC>

こちらはなぜか周期境界条件が認識されず、non-PBCとなってしまいました。原子数やトラジェクトリのフレーム数など合っていそうです。

3-4. PyTrajでアラインメント

ファイルから読み込んだトラジェクトリでアラインメントをとっておきます。

align()を使います。デフォルトでは全ての原子について最初のフレームを参照構造（ref）にしてアラインメントを取ります。

PyTraj（というよりAmber?）では、原子の選択はmaskと呼ばれるようです（Amber mask selection）。特定の原子団を選択したければ「mask=選択基準」といった形式で引数に渡すことになります。

例えば、Cα原子についてアラインメントを取るには引数mask='@CA'を渡します。

# Cα原子について最初のフレームを参照にアラインメントを取る
traj_align = pt.align(traj, mask='@CA', ref=0)

無事、アラインメントが取れました。

3-5. PyTrajで相互相関解析

PyTrajでの解析準備が終わったので相互相関解析を行います。

Matrixのcorrel()を使うことで相互相関の行列が得られます。Cα原子を選択して計算しましょう。

from pytraj import matrix

mat_cc = matrix.correl(traj_align, mask ='@CA')

print(type(mat_cc))
print(mat_cc.shape)

# <class 'numpy.ndarray'>
# (317, 317)

計算対象のタンパク質構造の残基数が317なので、「アミノ酸残基 x アミノ酸残基」の2次元行列が得られました。

あとはこれをプロットするだけです。seaborbでヒートマップを可視化します。

ax = sns.heatmap(mat_cc, vmin = -1, vmax = 1, cmap = 'PiYG_r', square=True, xticklabels=30, yticklabels=30)

ax.set_xlabel('Residues')
ax.set_ylabel('Residues')

① 60 ~ 90（残基番号 760 ~ 790）、② 120 ~ 160（残基番号 820 ~ 860）、③ 180 ~ 240（残基番号 880 ~ 940）あたりに相関が大きくなっている群がありそうです。

PyMolで少し色分けしてみました。

ヘリックスだったりシート構造だったりする部分のようです。「協同運動している箇所」と言われると「そうかも・・・？？？」という感じですが雑過ぎですね。。。ごめんなさい。

4. おわりに

以上、今回の記事ではMDトラジェクトリの解析（後半戦）を行いました。取り上げた解析は主成分分析と相互相関解析で、MDAnalysisとPyTrajを利用しました。

どちらも揺らぎの大きな特徴をつかむ解析のようですが、主成分分析では運動モードの種類を取り出す感じで、相互相関解析では分子内のアミノ酸残基の協同運動の傾向がでる感じでした。

個人的には第１主成分と第２主成分でそれぞれに投影した動画の運動の向きが変わっていたのが面白かったです。

短いシミュレーションでデフォルトの解析を実行しただけなので生物学的に意味のあることを言えるような記事ではありませんが、とりあえずの使い方がわかったのでヨシ！

以上、今回も間違いが多そうです。ご指摘いただければ幸いです。

ではでは！

*1:イベントの特設サイトはこちら（スパコンコロキウム）で、スライドPDFはこちら（直リンク）です。

*2:参考：pytraj - Installation

*3:Amberはライセンスが必要だそうですが、AmberToolsはフリーで使えるようです（GPL? LGPL?）。

*4:PyTrajのトラジェクトリファイル読み込みは、MDAnalysisでUniverseオブジェクトを作る時と引数の順番が逆であるのに注意です。

1. 計算対象

1-1. Tyk2キナーゼ阻害剤データセット

1-2. データセットの取得元と変換

2. ESPを計算して活性差を考察しよう

2-1. 計算対象の絞り込みと傾向

2-2. PyMolで眺める

2-3. ESP-DNNの結果で親和性を考察しよう

2-3-1. タンパク質

2-3-2. リガンドの比較

2-3-3. リガンドとタンパク質の比較

3. おわりに

1. ESP-DNN

1-1. ESPの復習

1-2. ESP-DNNって？

2. ESP-DNNをGoogle Colabで

3. NGL viewerで可視化

4. 拡張子PQRと謎の原子

4-1. PQR？

4-2. 謎の原子？

5. おわりに

よくやるミス

addCoords=Trueをつけると？

何が起きてる？

おわりに

1. 全体の総括

2. きっかけと感想

3. このシリーズでできなかったこと

4. （おまけ）試したけどうまくいかなかったもの

4-2. ProLIFの相互作用パラメータを変えた再解析

4-2-1. 課題

4-2-2. おさらい

4-2-3. ProLIFの相互作用の定義

4-2-4. RDKitでSMARTSパターンの認識をチェック

4-2-5. ProLIFの相互作用パラメーターを変えよう

4-2-6. MDトラジェクトリの再解析

4-3. OpenMMからAmberのファイルに変換してMMPBSA

4-3-1. ParmEdのインストール

4-3-1. ParmEdのファイル変換でエラー

4-3-2. エラーの原因と回避策

4-3-3. MMPB-SAの失敗

5. おわりに

1. キナーゼと阻害剤

2. 低分子 - タンパク質複合体のトラジェクトリ解析

2-1. トラジェクトリの動画

2-2. MDAnalysisでRMSD解析

2-2-1. RMSDの推移

2-2-2. フレーム間のRMSD

2-3. Protein-Ligand Interaction Fingerprint（ProLIF）で相互作用解析

2-3-1. ProLIF？

2-3-2. ProLIFでフィンガープリント解析

2-3-3. ProLIFでLigand Interaction Network（LigPlot）

2-4. MDAnalysisで原子間距離の解析

2-5. MDAnalysisで水素結合の解析-1

2-6. MDAnalysisで水素結合の解析-2

3. おわりに

1. 低分子を含むMDシミュレーションに向けた準備

1-1. タンパク質単独と何が違う？

1-2. OpenMM Forcefields?

1-3. Open Force Field Initiative？

1-4. MDTrajのインストール

2. 計算の対象構造について

3. タンパク質-低分子複合体の系をつくろう！

3-1. 処理全体の流れ

3-2. Step 1：タンパク質の前処理

3-3. Step 2：リガンド構造の修正

3-4. Step 3：リガンド用の力場の準備

3-5. Step 4：リガンドのオブジェクト変換

3-6. Step 5： リガンドとタンパク質のマージ

3-7. Step 6：水、イオンを追加し系を構築

4. 計算条件の設定と実行

4-1. 計算条件の設定

4-2. TeachOpenCADDの設定との違い

4-3. GoogleColabで実行する際の注意点

4-4. シミュレーション動画

5. おわりに

6. おまけ

1. この記事で試す解析

2. 主成分分析（PCA）

2-1. MD計算の主成分分析って？

2-2. MDAnalysisで解析の準備（おさらい）

3-6. Step 5：リガンドとタンパク質のマージ