「分子の表現についてのレビューと実践ガイド」のメモ ~part 4 マクロ分子の表現~

レビューをDeepLに投げ込みながら読んでいるのでそのメモです。

jcheminf.biomedcentral.com

ここまで３回で見て来たものは低分子に関するものでした。最後はマクロ分子の表現です。ケモインフォマティクスとバイオインフォマティクス の間の領域って感じのお話。

なお、順番は一部前後していますが、目次は論文そのままにしてます。また、 Creative Commons のもとで図表・内容を利用させていただいています。*1

Representations for macromolecules

取り上げられているのはバイオポリマーやオリゴマー、合成ポリマーといったポリマー構造を持つものです。

マクロ分子の例と表現例

ポリマーには以下の２種類があり、表現が難しくなっています。

1. monodisperse：構成するモノマーが同じ鎖長
2. polydisperse：確率的な (stochastic)性質があり鎖長が定義されない

またマクロ分子の表現はバイオインフォマティクスの要素が強くなります。
ケモインフォマティクスでは原子レベルで低分子を表現していたのに対し、ヌクレオチドやアミノ酸配列といった配列情報で表現するといった傾向があります。

Amino acid-based stuructures

ペプチド、タンパク質のビルディングブロックであるアミノ酸(AA)は1文字あるいは3文字の略記で表されます。

１文字表記の限界として、遺伝暗号に含まれる20種を表すのには十分ですが、 それ以外にも自然に存在するアミノ酸を表すには数が足りないということが挙げられています。

Peptides

ここでは大きく分けて2つの表記（CHUCKLES、HLEM)が紹介されています。

CHUCKLESはポリマーの配列からそのSMILESを導けるようにしたものです（逆方向の翻訳も可能）。 説明を読むとこんな感じ、、、

CHUCKES

CHUCKLESのオリゴマーにおける応用例としてBIOPEP-UWMというデータベース、また拡張版としてCHORTLESが挙げられています。

さて、より広範なマクロ分子を表現できる表記として以下の２つが挙げられています。

1. HELM (Hierarchical Editing Language for macromoleluces)：SMILESベース
2. SCSR (Self-Contained Sequence Representation)：v3000 Molfile フォーマットベース

以下ではPfizer社により開発されたというHELMの詳細についてもう少し見てみます。イメージを掴むために先に図を引用・・・

HELM

構成要素として複数のタイプの構造を持つ組み合わせを表現するシステムを作ることを目的として開発されたそうです。

HELMの表現と階層

HELM2ではさらに発展してポリマー混合物や自由形式のアノテーションも可能となっており表現できる幅が増えているそうです。

HELMは多くの製薬会社で実装されており、様々なデータベースやパッケージでもサポートされているそうです。 ChEMBLやRDKitにもあるんですって！知らなかった。。。

で、HELMについて字数を割いて詳しく説明されているの何でかな？と思ったのですが、生物と化学の融合地点としてBiocheminformatics表現が広がることで、 ペプチド医薬品といった低分子と生物学的製剤の間にあるような分野の表現が発展することを重視しているようです。

つまり、アトムベースのSMILESを配列ベースの表現に加えることで、天然のL-アミノ酸と非天然のD-アミノ酸の混ざったものといった表現の拡張ができる。非天然アミノ酸の組み入れはペプチドのバイオアベイラビリティの改善で重要なので、記述の幅が広がることでペプチド医薬の発展も記述できる、取り扱えるようになるというわけです。

なるほど！・・・それを踏まえた上で言わせてほしい！ Fig. 9の図、フェニルアラニンに修飾（Cl)入ってるからFで表現するのは違うと思うで！！

Proteins

タンパク質（50以上のアミノ酸残基）については、PDBデータベースとそこで使われている表現が挙げられています。

PDBにおける原子は以下の情報で表されています。
1. 連続する番号
2. 特定の原子名
3. 対応する残基の名前と番号
4. 鎖(chain)を特定する１文字のコード
5. 空間座標 (x, y, z)
6. 占有率 (occupancy)
7. 温度因子 (temperature factor)

2008年にはProtein Line Notation(PLN)というものがBiochemfusionにより開発されPubChemにも実装されているそうです。

PLNでは残基の構造を擬原子(pseudo-atom)で表すことで、化学表現と配列フォーマット間の変換をロスなくできるようになっているそうです。

Key macromolecules

グリカンや合成ポリマーといったマクロ分子についても見ていきます。

Glycans

グリカン(olygosaccharide、polysaccharide)にもデータベースがありますが、monosaccharideベースの表記となっています。

ドッキングといった相互作用解析ではアトムベースの表記が必要となるので、 グリカンを医薬品探索におけるリガンドとして扱うにはmonosaccharideベースの表記では不十分です。 そこで、アトムベースの表記へと変換するツールの開発が行われています。

また、Web3 Unique Representation of Carbohydrate Structures (WURCS)という表現が開発されています。

WURCS

WURCSは多くのデータベースに実装されていますが、ケモインフォマティクス のソフトウェアではほとんどサポートされていないそうです。

他のアプローチとして、ファーマコフォアに基づくものや言語モデルに基づくものといった研究もなされているようです。

Polymeric drugs

ポリマーの表現として最近、BigSMILES構文というものが開発されました。

以下のようなポリマーをエンコードすることができます。

BigSMILES

まだカノニカルなものはなく、応用例もないようですがさらなる研究が行われているそうです。

Graphical representations for molecules and macromolecules

先のセクションで見たマクロ分子の表現は、データの格納とケモインフォマティクス的解析を意図したものでした。 ここでは、マクロ分子とその物理化学的特性の可視化表現をみていきます。

2Dと3Dをそれぞれ扱いますが、何はともあれ見た方が早いので図を引用しておきます。

Fig.10より

2D depictions

分子の2D描写として、骨格構造のラスタ図、ベクトル図がよく使われます (Fig. 10a)。

以下のような点が問題になることがあり、2008年にはIUPACがオススメの標準描画方法を出したりしています。

各ソフトウェアそれぞれより良い描画方法を目指してアルゴリズムの改善等試みられているそうです。

構造そのものの描画以外にも、反応や相互作用の研究といった点でも２Dの描画が使われます (Fig. 10e)。 ここでは分子の置かれた環境と振る舞いの表現が重視されています。

Fig. 10には描かれていませんが、以下の３つも特徴のある例として取り上げられています。

3D depictions

3D描画のソフトウェアの例としてAvogadro、PyMOL、VMDといったものがあります。

表現方法の例として以下が挙げられています。

vdW球による可視化の応用例として、分子表面の描画による相互作用解析が挙げられています。
3D描画はドッキングや機構の研究で、2D描画は構造活性相関研究でよく使われているそうです。

以上がマクロ分子の表現でした。このあたりがAIにどのようにつながっていくのか、最後に応用例を見ていきたいと思います。

AI applications within drug discovery using molecular representations

Representations for macromolecules

マクロ分子の表現において人気のある応用はタンパク質構造の予測です。　　

また、グリカンの分野における研究例として、

1. 擬受容体(pseudo-receptor)モデルを使ったバーチャルスクリーニング
2. 機械学習によるグリカン-タンパク質相互作用の解析
3. グリコシル化部位予測

が取り上げられています。

Graphical representations for molecules and macromolecules

可視化については、さらなる高速化の研究が行われており、最近ではバーチャルリアリティや3Dプリンティングといった技術と紐づいています。

また、増え続けるデータを処理するための構造マイニング技術として、Optical Character Recognition (OCR)が挙げられています。 機械学習や確率的パターン認識技術により、2D描画の標準的な化学表現への変換が行われていますが依然として課題はあります。

変換が難しくなるものとして、以下のような項目の処理があげられています。

1. 画像の解像度
2. 化学的略記のコンピュータによる解釈
3. テキストに埋め込められた画像
4. 複数の構造を含む図に埋め込まれた画像
5. 反応パスウェイ中に埋め込まれた画像
6. 骨格の式あるいはMarukush構造で表された画像

マクロ分子の表現については以上です。

さて、４回にわたって分子の表現方法について眺めて来ましたが、いよいよまとめの段階に入ってきました。残りはDiscussionとConclusionです。

Discussion

ここまで様々な表現方法をみてきましたが、著者らは医薬品探索における問題を解決するためには、複数の表現を同時に活用することが必要と指摘しています。

タンパク質の構造予測を例として以下のような流れをあげています。

1. タンパク質の配列から始まり
2. 粗い3Dモデルの作成
3. 分子動力学法によるフォールディングメカニズムの理解
4. 最終的な配置と構造予測
5. ドッキング計算への応用

また、表現の選択に影響する要因として、2点あげられています。

1. 表現を生成する手法の複雑さ
2. オープンに利用可能なものか否か

・・・確かに。オープンソースツールで簡単に使えるものがあると遊んでみよう、ってなりますもんね！！

また、分子表現の変遷の歴史をたどってきたことを踏まえて、使われ続ける表現とそうでないものとの違い、要因を議論しています。

まず、変化の要因の一つとしてコンピューターテクノロジーの進化が挙げられています。

1. ストレージの増強
2. プロセッサーの質
3. パラレルプログラミング

かつての線形表記、IUPAC-Dyson、WLNは当時の状況としては妥当だったものの、主に人間が操作するもので、計算機による取り扱いが難しかったので、計算機がより使われるようになると廃れていきました。

現在ではよりコンピューターが単純に扱える表現が使われる傾向にあります(molecular string representation)。また、より計算が必要となる詳細な表現も現在では使えるようになりました(e.g. hashed fingerprints)。

別の要因としては、ケモインフォマティクスコミュニティにおける受容があげられています。より人間が読める(human readable)表現の方が好まれ、結果として長く使われるようになったかも、、、らしいです。

最後に、分野の違いの及ぼす影響があげられています。
異なる分野(chemoinfomatics, bioinfomatics, AI)ではグループ内の歴史的、継続性の理由から表記の選択、使われ続けるか否かが変わってくるのではないか、ということです。

継続は力なりですね！（違うか）

Concolusions

分子は複雑な構造なので、その表現においては様々な特性とともに、低分子とマクロ分子のそれぞれの異なる性質といった点も表現しなければなりません。ケモインフォマティクス 、バイオインフォマティクスの発展とともに、医薬品探索のプロセスが加速し、分子の振る舞いについての理解も進んで来ました。このレビューで取り上げられているような様々な表記方法を学び、AI駆動の医薬品探索の裏側を知ることでもっと有効活用できるようになると良いですね！

みたいな感じのことが書いてありました！

まとめにかえて

以上、懲りもせずに適当なメモを垂れ流してみました。ひたすらCreative Commonsにのっかっていくスタイル！

線形表現の歴史や、派生系、マクロ分子の表現等々、知らないことがたくさんあって個人的にはとても勉強になる文献でした。 途中謎の自作の表や無意味な箇条書きが出て来たのは、私が文章では単語の関係性が理解できない読解力低いマンだからです。 human-readableですら私には遠いですわー。

さあ、これでAI創薬に飛び込む準備ができた！・・・のか？

知らんけど。。。

色々と間違いがありそうなのでご指摘いただけると幸いです。

*1:https://creativecommons.org/licenses/by/4.0/

「分子の表現についてのレビューと実践ガイド」のメモ ~part 3 化学反応の表現~

レビューをDeepLに投げ込みながら読んでいるのでそのメモです。

jcheminf.biomedcentral.com

先の２回の記事では分子のグラフ表現と線形表現でした。単独の分子の表現方法を踏まえ、次のステップとして分子の変化、化学反応の表現方法へと話題がうつっています。知らない事が多かったので雑訳の羅列になってしまっていますがまあいつもの事ですね。

なお、順番は一部前後していますが、目次は論文そのままにしてます。また、 Creative Commons のもとで図表・内容を利用させていただいています。*1

Representations for chemical reactions

Harnessing reaction data for drug discovery

さて、化学反応は次のように理解できます。

1. 一群の分子のセットを
2. 関係するもう一つの分子のセットへと
3. 特定の条件のセットの下で
4. 相互変換すること

あらためて言葉で記述されると、なるほどそうなるのかーという感じです。

で、以下のような図式でよく書かれるわけですが、これが機械可読性が低い(machine-readableでない)のでどうにかしたい！、というわけです。

こういう感じのやつ

機械で処理できよるようにするには、前回、前々回見てきた分子の表現を応用する事ができます。（但し、各表現のもつ限界もそのまま引き継いでしまいます！予め注意喚起！）

また、フォーマットについてはRXNとRDファイルがよく使われます。グラフ表現のCTfileのところで出てきたやつですね。

以下、様々な化学反応の表現表現を見ていきます。下図のようなものが取り上げられています (Fig. 6より引用)。

化学反応の表現

Reaction SMILES and SMIRKS

Reaction SMILESはSMILESで化学反応が扱えるように拡張されたものです。

記載の仕方は以下のような感じです。

Reaction SMILES

尚、Atom mappingは反応剤には使えません。

また、Reaction SMILESそのものにはテキストによる情報（反応中心や反応条件 etc.）を追記することはできませんが、ファイルレベルで、RXNやRDファイルフォーマットが取り扱えるメタデータへの記載という形で情報を追加することができます。

次にとりあげるSMIRKSは「SMILESとSMARTSのハイブリッド」とも言えるものです。

SMIRKSは反応における分子の変化のパターンを一般化して記述します。

（※ SMARTSでは分子のパターンと部分構造が一般化されて記述されていました。）

記述の特徴はこんな感じ…

SMIRKS

また、図のようにSMIRKSはグラフ(reaction graph)に変換することもできます。

RInChI

RInChIはInChIを反応の表現に拡張したものです。

反応の「単一で、かつ順序の影響を受けない識別子 (unique & order invariant)」として導入されました。
Reaction SMILESは分子の順序が任意でしたが、RInChIでは同一の反応の表現が同じになることを意図されているので、反応の重複を見分ける事ができるという利点があります。

一方、文法がより複雑という課題があります。

カノニカルな表現を意図して導入されたInChIがSMILESと比べて人間にわかりづらいものだったこととそのまま対応する感じですかね？

RIhChIの記述はこんな感じ…

RInChI、RInChI key

「5. 反応の進行方向」については、例えば同一の反応がわずかな条件の差で異なる結果になった時に、それぞれの条件での進行方向が記載できることで、細かな条件間の比較を行うことが容易になる、という利点があります。

更に、RInChIの拡張としてProcAuxInfoというものが提案されており、メタデータ(収率、温度、濃度 etc.)を格納できるようになっています。

また、ハッシュ化されたRInChI keyを使うと反応データをインデックス化して効率的に検索する事ができます。

難点としては、文法の複雑さ以外には、SMARTSやSMIRKSに相当するものがなく、 部分構造検索や一般化した化学反応の変化の記述が難しいということがあげられています。

Condensed graph of reaction (CGR)

CGRはVernekらにより開発された表現です。こんな感じ…

CGR

図では少し見づらいですが、エチル基はカルボキシ基の酸素原子と緑色の線で繋がれています。「エステル化の仮想的な遷移状態」っていうことでしょうか？？

Bond electron matrices (BE-matrix)

Bond-electron matrixは反応を行列で表現しようとするものです。

BE-matrix

原料のBE-matrixに反応のR-matrixを足す事で生成物のBE-matrixが得られます。

行列表現で表せる情報を増やそうとしている点が、他の表現とは異なるアプローチのようです。

Hierarchical organization of reaction of reactions through attribute and education (HORACE)

HORACEは機械学習のアルゴリズムを利用して化学反応の分類を行うもので、化学反応の階層的な表現ができます。

HORACE

階層化された表現となっていることで、SMILESのような構造だけを表す反応表現よりも豊かな表現が可能となっています。

InfoChem CLASSIFY

InfoChem CLASSIFYはルールベースの合成計画法で使われるアプローチに発想を得た反応の表現です。

InfoChem CLASSIFY

各原子について生成されたハッシュコードは全ての出発物質と一つの生成物について合計され、反応中心のユニークな表現として使われます。

また、ハッシュコードは反応の分類や、様々な合成計画ツールで反応中心を抽出するアプローチとしても使われています。

文献中では反応中心を見つけ出すアプローチとしてよく用いられるものとして、「出発物質と生成物の間の最大共通部分構造(MCS)をみつける事」があげられています。
また、MCSの問題点として、NP完全で非決定的な多項式時間がかかるという事指摘されてます。
Inform CLASSIFYはMCSに代わるアプローチになるかも、という事が言いたいのでしょうか？ちょっとよくわからなかったです。

ちなみに取り出した反応中心は必要とする特異度(specificity)のレベルに応じて、隣接する原子にも範囲を拡張できるそうです。

範囲と特異度の関係は以下のようになっています。

Reaction fingerprints

Reaction fingerprintは反応のベクトル表現です。反応中心で起こる構造の変化を表します。

こうして作成したベクトル表現は例えば、機械学習でモデルの構築に用いることができます。

reaction fingerprintは反応中心を特定する別のアプローチとして期待されていますが、「反応グラフ(reaction graph)への変換が難しい」という問題点があります。

また、立体化学の取り扱いはまだできません（研究は盛んに行われているそうです）。

以上が化学反応の表現でした。つづいてマクロ分子の表現になりますが、ここで一旦区切ります。

AI applications within drug discovery using molecular representations

今回も途中をすっ飛ばして化学反応表現の応用について触れた箇所に飛んでおきます。

Representations for chemical reactions

化学反応表現の応用の主なものは逆合成と反応予測です。

これらはAI駆動創薬において「ループを閉じる」ので重要となります。(Design-Make-Test-Analysisのループですかね？)

提案された構造が合成可能なものかを判定し、適切な逆合成により合成ルートを提案することが大事なんですって！

D! M! T! A! P! D! C! A!

この辺りの話題についてより知りたい場合は以下の文献がオススメされていました。

• Mater AC, Coote ML (2019) Deep learning in chemistry. J Chem Inf Model 59:2545–2559
• Jensen KF, Coley CW, Eyke NS (2019) Autonomous discovery in the chemical sciences part I: Progress. Angew Chemie Int Ed

ACIEの文献はarxivにもあるっぽいです。

arXiv:2003.13754v1

70ページある・・・長い・・・

まとめ

以上、化学反応の表現でした。私はreaction SMILESくらいしか聞いたことがなかったので、行列表現やベクトル表現に相当するものがあるというのはなるほど！という感じでした。

分子の表現と反応の表現の対応関係はこんな感じですかね？

・・・HORACEとInfoChem CLASSIFYは対応するものがよくわからなかった。

また、「反応中心をどう見つけるか？」という点が議論に上がっていたのが興味深かったです。 確かに、複数の分子が記載されてる場合に、「どの分子を（目的物質に直接変化する）出発原料とするか？」は重要だし難しい問題だな、と思いました。

便利な試薬がたくさん開発されている分、複雑・高価な試薬が増えて、時と場合によっては原料よりも試薬の方が分子量が大きくて複雑な構造をしている、ということもあると思います。 そうすると、各々の分子だけからどれが原料かを判断するのは難しく、反応前後の差分、複数分子間の比較を取らないと分からない、というのも納得できます。

また、原料をどれに設定するかは、反応の収率の計算に直接影響するので、原料の設定次第で同じ反応に紐づいた数値データが変わってしまうことになります。 合成スキームとして書いてあったらすぐに気づくことができそうですが、機械可読な表現・数値に変換してしまうとこのあたりの微妙な差が埋もれて、 数字がひとり歩きしてしまいそうだな、と思ってしまいました。

（収率をb.r.s.m.で記載してたら良い反応に見えるけど、実際には・・・？？？みたいな問題とかもありますしねー・・・関係ないか。）

途中無意味な箇条書きを乱発しています。わかった気になるから！！わかんないけど・・・。

次回、最後！マクロ分子に続く！

色々と間違いがあると思います。ご指摘・ご教示いただければ幸いです。

*1:https://creativecommons.org/licenses/by/4.0/

「分子の表現についてのレビューと実践ガイド」のメモ ~part 2 線形表現~

レビューをDeepLに投げ込みながら読んでいるのでそのメモです。

jcheminf.biomedcentral.com

先の記事では分子のグラフ表現の箇所を取り上げました。続いて線形表記についてです。

なお、順番は一部前後していますが、目次は論文そのままにしてます。また、 Creative Commons のもとで図表・内容を利用させていただいています。*1

Linear notations for small molecules

線形表記(linear notation)は Ctab (connection table) をエンコードした文字列で、体系的なルールで解釈することができる分子の表現です。

線形表記には、行列と比べて以下のような利点があります。

1. サイズが小さい(ディスク容量的にも、サイズ的にも)
2. 操作が簡単 (Excelでコピペとか余裕)

Table 1 にアラニンを例に線形表記の例がのっています。

線形表記の例

The IUPAC quest for a universal notation

線形表記はルールで解釈できると書きましたが、そのルールをどう設定するか？については紆余曲折があったようです。（現在も進行中？）

IUPAC(International Union of Pure and Applied Chemistry)の戦いの歴史が書いてありました。ざっとこんな感じ・・・

歴史

現在ではDyson、WLN両者ともに下火になっているそうです。次節、より現代使われている表記法へと話が続きます・・・

が、その前に、IUPACが望ましいとした特性や表記法の分類がなるほど感があったので載せときます。

好ましい表記とは？

また、表記法の形式による分類の考え方はこんな感じ…

形式

この辺りを意識しながら以降の表記法の紹介を読んでいくと面白いかもしれません。

The advent of contemporary notations

Simplified Molecular Input Line Entry System

WLNよりも直感的な表記としてSMILES(Simplified Molecular Input Line Entry System)が1988年Weingingerらにより開発されました。
Daylight社のソフトウェアに組み込まれたことをきっかけに、現在でも同社が維持運営しているそうです。

SMILESには表記を生成するルールの違いや表される構造情報により複数の流派(?)があります。

生成のルールは以下のような感じです。次に引用したFig. 4と比較するとわかりやすいかもしれません。

SMILESの生成

色々とあるSMILESのうち、まずは生成ルールが異なる２種類（canonical smiles、randomized smiles）です。こんな感じ・・・

canonical or random?

細かい変換ルール(Cとcの違いなど)は割愛します。

randomized SMILESは表記が毎回変わってしまい困ってしまいそうですが、逆にデータ拡張(data augmentation)に使えるという利点があるそうです。

他にも派生したSMILESが色々とあります。

SMILESいろいろ

一口にSMILESっていってもいろいろあるんですね。

続いてもう少し別の特徴がある線形表記SMARTSとInChIへと続きます。

SMILES Arbitrary Target Specification (SMARTS)

SMARTSは部分構造検索のために開発されたSMILESの拡張です。

1. SMILESよりも利用できるシンボルが多く、より一般的な形で分子グラフを表現できる (→ コンピュータ科学における正規表現に類似)
2. 原子が一つ異なる、あるいは、結合の位置が異なる分子の集合を表現できる
3. 論理演算子 ("OR"、"NOT")も使える
4. 同位体(isotope)や結合のタイプ(aromatic, aliphatic)を明示できる
5. Recursive SMARTSでは原子の環境についての情報も含められる(ex. ortho, meta, para)

といった特徴があります。

SMILESは全てSMARTSとしても有効ですが、SMARTSがSMILESとして有効とは限りません。

International Chemistry Identifier

InChIは2006年 NISTにより標準的に自由に使える形式表現として導入されたものです（オープンソースのカノニカルな表記）。

様々な情報を含む多層構造で形成されており、下の図のような感じです。

InChI

InChIはSMILESと異なり元の分子グラフに戻すことができることは保証されていません。また、SMILESの方が人間には読みやすいという利点があります。

また、InChIのハッシュ化されたバージョンとしてInChIKeyがあります。

元のInChIの一義的な表現となるようになっていますが、InChIKeyが複数のInChIにマッピングされることがあります（InChIKey collision)。

Using chemical descriptors to represent molecules

ここまで見てきたのは原子に基づく(atom-based)ものです。

これに対して分子の特性をもとに表現するアプローチとして分子記述子(molecular descriptor)があります。

単純化するとこんな感じ…

アトムベース or 記述子

すごいたくさん種類があるらしい（Dragonソフトウェアでは4885種類！）ですが、以下では大きく２つの分類(Structural Kyes, Hashed fingerprints)が紹介されています。

Structural Keys

構造キー(structural key)は特定の化学構造が有る(1)か、無い(0)かにしたがってエンコードします。

２つの例が挙げられています。

構造キーの例

MACCSキーの注意点として、フラグメント166(or 960)種のセットに基づくといっても、ソフトごとに実装が違う場合があるので、異なるソフトでは同じ構造が異なるビット列に割り当てられる可能性がある、という点が指摘されています。

Hashed fingerprints

Chemical fingerprintはベクトルで、

1. 各要素は物理化学的特性あるいは構造の特徴をエンコード
2. それら要素がインデックス化（順序をつけられて）されて並べられて

います。

先の構造キーでは表す特徴のパターンがあらかじめ定義されていましたが、ハッシュ化されたフィンガープリントでは、各特徴が分子そのものから生成されるという点が異なります。

ビット列の長さを決めた上で生成し、分子の持つパターンをハッシュ関数によってビットに割り当てます。

２種類例が挙げられています。

フィンガープリント

Chemoinfomaticsでフィンガープリントがよく用いられるのは、 「分子のグラフ構造を直接的に素早くベクトル表現へと変換できるので、数理モデル化の入力として使える」ということにあるようです。

また、柔軟性のある表現なので、物理化学的特性を整数(integer, ex. 水素結合の数)や浮動小数点(float, ex. 分子量)にエンコードできるというのも魅力のようです。

線形表現については以上でお終いです。

つづいて話は変わって化学反応をどう表現するか？という話題が取り上げれています。ちょうど良いのでここで区切ります。

途中に出てきたIUPACによる分類指針を参考にここまでの分子の表現をざっくりまとめると次のようになりますでしょうか？

グラフ表現と線形表現の雰囲気

SMILESの派生系などを全て無視して単純化しすぎな分類ですが・・・

AI applications within drug discovery using molecular representations

今回も途中をすっ飛ばして線形表現の応用について触れた箇所に飛んでおきます。

Linear notations for small molecules

線形表現(SMILES、molecular fingerprint)は分子の特性予測やQSARでよく利用されます。

SMARTSは部分構造を定義し、関連する化合物を選択したり取り除いたりする目的で使われます。

ディープラーニングと関係する内容としては、

1. 自然言語処理(NLP)のツールを使ったde novo分子デザインの分野でSMILESが活躍
2. randomized SMILESを使ったデータ拡張
3. オートエンコーダーを使った潜在空間表現の学習による分子特性の予測

といったものが挙げられます。

ニューラルネットワークでは分子のベクトル表現を学習したうえで予測の利用しますが、 これは古典的な機械学習のアプローチでハッシュ化されたフィンガープリントを入力として使うことに類似しているので、 「学習されたフィンガープリント(lerned fingerprint)」とも呼ばれるそうです。

以上が線形表現の応用でした。

次回、化学反応の表現に続く！

途中浅い理解で雑な表にまとめたりしているので、色々と間違いがあると思います。ご指摘・ご教示いただければ幸いです

*1:https://creativecommons.org/licenses/by/4.0/

「分子の表現についてのレビューと実践ガイド」のメモ ~part 1 グラフ表現~

HTSなどの技術の発展に伴い実験データが増加し、計算機による解析が盛んになるにつれて、コンピューターが処理可能かつ、様々な分野の科学者が理解できるような分子の表現（Molecular representation） が必要とされるようになってきました。

たくさん開発されてきた表現方法の中で、特に医薬品探索で使われる電子的な分子の表現 (electronic molecular representation)と巨大分子の表現(macromolecular representation)にフォーカスを当てて解説しちゃうよ！

・・・って感じのイントロのレビューを読みました。オープンアクセス万歳！！(CC by 4.0)

David, L., Thakkar, A., Mercado, R. et al. Molecular representations in AI-driven drug discovery: a review and practical guide. J Cheminform 12, 56 (2020).

jcheminf.biomedcentral.com

グラフィカルアブストラクトからして素敵な感じ・・・

TOC

知らないことがたくさんあったのと、「グラフ表現 → 線型表現」という説明の順が面白かった（1D → 2D の説明が多い気がする）ので備忘録を残します。

と言っても英語だと覚えられない自分のための雑メモですが・・・

長いので４回にわけて見ていきます。１回目はグラフ表現！

なお、内容は一部前後しますが、目次は論文そのままにしてます。

Graph representations for small molecules

まずは分子グラフ（Molecular Graphs）をきちんと理解しましょう。

molecular graph

分子グラフは２次元のオブジェクトですが、3Dの情報(ex.原子座標、結合角、不斉点)を表すこともできます。 また、簡単に可視化できるソフトも多くあります。

注意点として、ノード間の関係性は属性(attribute)としてノード and/or エッジにエンコードしておく必要があります。

理由は、グラフのノードは数学的なオブジェクトであって（正式には）空間の位置を持っておらず、対ごとの関係性(pairwise)しか保持していないからです。

Mathmatical definition of graph

グラフはタプルで定義されます。

タプルによる定義

上は抽象的な数学の概念なので、コンピューターで扱えるように線形データ構造(linear data structure)にマッピングします。

線形データ構造

一般的な表現方法は隣接行列(adjacency matrixあるいはconnectivity matrix)です。

隣接行列の例

原子や結合の種類は特徴行列(feature matrix)の形式で表すことができます。 下の図ではone-hotエンコーディングで表していますが、要素の値を整数(integer)として、よりサイズの小さい行列とすることもできます。

特徴行列

Graph traversal algorithms

非線形データ構造(non-linear data structure)のグラフを線形の行列表現にするにあたって、ノードに番号をつける必要が生じますが、この番号の付け方によって表現が変わってしまいます。

一つの分子を一貫して同じ行列表現で表すために、番号をつけるアルゴリズム（graph traversal algorithm）が重要となります。

実際、ソフトウェアごとのグラフ検索における分岐の切り方の違いの影響で、得られる行列表現が変わってしまうことがあります。

ただし、表現の一貫性を重視しないような場合、例えばディープラーニングへの応用で、よりノイズの多いデータが欲しいといった場合には、ランダムサーチを用いることもできます。

Fig. 2ではよく使われるグラフ検索アルゴリズムが図示されています。

グラフ検索アルゴリズム

Advantages of graph representations

分子をグラフで表現することの利点は、

1. 3Dの情報を自然に埋め込むことができること
2. サブグラフが解釈可能であること

です。

利点①は特徴行列にあわらされている通りです。　　

利点②はサブグラフを部分構造として解釈できるということを意味するのでしょうか？？？比較として、文字列表現であるSMILESの一部(substirngs)を取り出しても、妥当性のあるグラフとなるとは限らない、ということをもとに理解すれば良いようです。

Limitations of molecular graph representations

Breakdown of graph model

反対に欠点として「グラフでは表せない分子」の例が挙げられています。

1. 非局在化した結合を持つ分子 (ex. 配位化合物、イオン結合・金属結合を持つ分子)
2. 3D空間で原子の配置が変わる分子 (ex. 互変異性体)

欠点①は、原子価結合法(valence bond theory)で説明できないもので、原子間のペアワイズな関係性だけ表すのが難しい、というのがその理由です。

①の解決策の一例としてhypergraphが紹介されています。
Hypergraphでは、エッジを原子のタプルではなく、少なくとも２つの原子からなる一群のセット(hyperedge)として表すので、多価結合(multi-valent bond)を表すことができます。

欠点②は、結合が切断されたりで生じたりする分子は、静的なグラフ表現１つでは表し切れない、表現してもあまり意味がない、といったことのようです。

Challenges of working directly with the graph representation

またその他の欠点として、「グラフ表現が（メモリ的にも、文字通りサイズ的にも）コンパクトではない」ということが指摘されています。

行列にしろ図にしろ、グラフは(1D)線形表現と比べて

1. 検索が困難
2. 必要なメモリが大きい

ため、多数の分子を扱う時に問題となるようです。

Connection tables and the MDL file formats

続いて分子のグラフ表現と関係の深いフォーマットとして、connection table*1とMDLファイルフォーマット*2があげられています。

## connection tables

connection table(Ctab)は以下の図のようなもので、６つの部分に分かれます。

connection table

connection tableは化学構造情報の取り扱いにおける標準的なフォーマットの一つで、Molファイルフォーマットの下地になっています。

atom block部分の注意点として、より実用的には水素原子を暗に (implicit) 取り扱うということがあります。
具体的には、

1. 水素原子をheavy atomのプロパティとして取り扱う
2. 利点はatom blockとbond blockのサイズを小さくすることができること
3. 水素原子は必要に応じて原子価モデル(valence model)のルールに基づいて再計算できる
4. 再計算のため、価数の情報を明示しておく必要がある

という感じです。

The Molfile format

MDLにより開発された、Molファイルフォーマットを含む一連のファイルをCTファイル(CTfiles、chemital table files)と呼びます。

CTファイルでは分子の構造を表現するのにconnection tableを使用します。つまり、connection table自体はファイルフォーマットではない、ということに注意してください。

Molファイルは拡張が可能で、SDファイル(structure/data file)では複数の分子の構造情報とその他のプロパティのデータを含むことができるようになっています。

CTファイル

以上が分子のグラフによる表現とそれを格納するためのファイル形式でした。

論文では続いて線形表記が紹介されていますが、長くなったので一度表現の紹介は中断します。

AI applications within drug discovery using molecular representations

論文の順序とは前後してしまいますが、この論文の主題は「医薬品探索におけるAIの利用と分子の表現」ですので、 忘れないうちにグラフ表現の応用可能性について触れられた箇所をご紹介しておきたいと思います。

Graph representations for small molecules

この文献で取り上げられた分子の表現は、それ自体がディープラーニングの表現学習(representation learning)に使われています。

表現学習の考え方は

1. 与えられたオブジェクト（e.g. 分子)の内部表現(e.g. ベクトル)を学習し
2. その表現を予想タスクに用いる

というものです。

表現学習においては、何を入力とするかが大事で

1. 入力に適した分子の表現をまずみつけること
2. その表現が問題を解くのに必要とされる情報を十分に有している事

が最初のポイントとなります。

古典的な機械学習(RFやSVM etc.)とは異なり、ディープニューラルネットワーク(DNN)を使ったアプローチでは、基本的に表現学習を行っています。

最近の医薬品探索の分野ではグラフニューラルネットワークの発達とともに、分子のグラフ表現が利用される流れが来ています。

以下のような応用例があります。

1. 分子の特性の予測 (property prediction)
2. 分子グラフの生成 (molecular graph generation, de novo design)
3. 合成ルートの提案 (synthesis prediction)

やり方の流れとしては大きく以下の２段階となっています。

1. グラフ表現学習の実施
   • グラフ表現全体から
   • グラフネットワークを使って
   • グラフの埋め込み(graph embedding)を取得
2. 学習した表現を入力として予測モデルを構築
   • RFやDNNなど
   • ここは古典的なフィンガープリントを使うときと一緒　　

以前はよりコンパクトな線形表現が好まれていたそうですが、 最近はちょっとずつ変わってきていて、（メモリを食う）グラフ表現も使われるようなってきているそうです。

分子グラフの他の利用としては、原子レベルのシミュレーションの入力（e.g. 分子動力学）が挙げられていました。 力場を使ったエネルギー計算に原子配置や結合の情報が必要となるため応用されているとのことです。

以上、グラフ表現のAIへの応用でした。

関連する内容として、文献では以下のレビューがオススメされていました。

• 化学と医薬品探索におけるディープラーニングの応用
  • Mater AC, Coote ML (2019) Deep learning in chemistry. J Chem Inf Model 59:2545–2559
  • Chen H, Engkvist O, Wang Y, Olivecrona M, Blaschke T (2018) The rise of deep learning in drug discovery. Drug Discov Today. 23(6):1241–1250
• 生成モデル
  • Sanchez-Lengeling B, Aspuru-Guzik A (2018) Inverse molecular design using machine learning: generative models for matter engineering. Science. 361(6400):360–365

では、次回線形表現に続く！

途中浅い理解で雑な表にまとめたりしているので、色々と間違いがあると思います。 ご指摘・ご教示いただければ幸いです。特にハイパーグラフはなんのことやら分からなかった・・・

*1:接続表？適切な訳語がわからない

*2:MDLは現在はBIOVIA

TeachOpenCADD トピック１１（パートC）〜オンラインAPI/サービスを使った構造に基づくCADD〜

こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル１１をもととしておりCC BY 4.0のライセンスに従っています。

トークトリアル１１は論文発表時のTeachOpenCADDには含まれておらず新しく追加されたものです。
４つのノートブックに分かれており、ここではその4つ目を訳出します。

パートCではパートBで行ったドッキング結果を解析しますが、パートBで使われていたオンラインサービスはすでに利用できなくなっていたため、ローカル環境ドッキングを行いました（実施方法はこちらの記事に記載いたしました）。
AutoDockVinaの２つの出力ファイルlog.txt、out.pdbqtをそれぞれファイル名を変更してlocal_vina.out、local_result.pdbqtとしてdataフォルダに格納し解析を行っています。

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

以下、訳。

トークトリアル 11 （パートC）

オンラインAPI/サービスを使った構造に基づくCADD

Developed at AG Volkamer, Charité

Dr. Jaime Rodríguez-Guerra, Dominique Sydow

このトークトリアルの目的

これは「オンラインWebサービス」についてのトークトリアルのパートCです。:

• 11a. キナーゼ阻害剤の候補をKLIFとPubChemで検索
• 11b. 11aで取得した候補化合物をターゲットタンパク質に対してドッキング
• 11c. 結果を評価し既知のデータと比較

入力構造を取得し、ドッキングが終わったので、結果が役に立つものかどうか評価します。

学習の目標

理論

• タンパク質ーリガンド相互作用
• ドッキングの偽陽性

実践

• 結果の可視化
• 自動化した解析の実行

レファレンス

• タンパク質ーリガンド相互作用に関するレビュー(Int. J. Mol. Sci. (2016), 17, 144)
• PDBデータベースにおける非共有性結合相互作用の系統的な解析 (M. Med. Chem. Commun. (2017), 8, 1970-1981)
• タンパク質ーリガンド相互作用が医薬品の作用にとって鍵となることについての次の書籍の中の章（Klebe G. (eds) Drug Design. Springer, Berlin, Heidelberg. DOI: 10.1007/978-3-642-17907-5_4）
• Notebookのためのインタラクティブな可視化ソフト、NGLView (Bioinformatics (2018), 34, 1241–124)
• 自動的にタンパク質ーリガンド相互作用を見つけるためのPLIP(Nucl. Acids Res. (2015), 43, W1, W443-W447)

---

理論

タンパク質ーリガンド相互作用

リガンドの結合は、リガンドとタンパク質ポケットの表面との間、あるいはタンパク質ータンパク質境界面における非共有結合性相互作用によって主に支配されています。このプロセスは、静電的相補性と形状相補性、誘導適合（inducued fitting）、脱溶媒和過程、その他に応じて変化するものです。

文献からいくつか引用します。

José L. Medina-Franco, Oscar Méndez-Lucio, Karina Martinez-Mayorgaから改変して引用します。

タンパク質ーリガンド相互作用（protein-ligand interactions, PLIs）とタンパク質ータンパク質相互作用（protein-protein interactions、PPIs）を理解することは、分子の認識における核となる部分で、多くの科学分野において重要な役割を果たします。PLIsとPPIsは医薬品探索において幅広い分野で実践的な応用があります。含まれるものをいくつかあげれば、ドッキング、構造に基づくデザイン、フラグメント化合物と低分子化合物、そして他の種類の化合物のバーチャルスクリーニング、複合体のクラスタリング、そしてアクティビティクリフの構造的な解釈といったものがあり、他にもあります。

もちろん、これらの相互作用はいくつかの方法によって合理的に説明することが可能で、ドッキング結果の系統的な解析への道を開きます。例えばMed. Chem. Commun., 2017,8, 1970-1981改変して引用します。

私たちはPDBから、タンパク質と複合体を形成している低分子化合物で、分解能2.5Å以下の全てのX線構造を抽出し、その結果11,016個の複合体のコレクションを得ました。リガンドとして考慮するためには、化合物は低分子である、医薬品化学の応用にとって興味深いものである、といったいくつかの判断基準を満たしていなければなりません（緩衝液や結晶作成用の混合物の一部は除外します）。このコレクションには750,873のリガンドータンパク質の原子のペアが含まれており、ここでは原子のペアは4Å以内の距離にある２つの原子として定義されています。最も頻度の高いトップ100のリガンドータンパク質原子ペアは7つの相互作用のタイプに分類することができました（下の図を参照してください）。最も頻度高くみられた相互作用の中には、疎水性相互作用（hydophobic contact）、水素結合そしてπ-スタッキングといった、リガンドデザインにおいてよく知られており、広く使われている相互作用がありました。それらに続いて、弱い水素結合、塩橋、アミドスタッキング、そしてカチオンーπ相互作用といったものがありました。

タンパク質ーリガンド相互作用を自動的に評価するプログラムはいくつかあります。最も人気のあるプログラムの一つはPLIPで、Webサーバーが公に利用可能で無料で使えるPythonライブラリがあるというのが人気の理由です。文献に付随するsupporting informationはタンパク質ーリガンド相互作用を簡潔に解説しており、非常に簡単に理解できるようになっています。導入としては、次のパラグラフで十分でしょう。

PLIPはタンパク質残基とリガンドの非共有結合性相互作用をみつけるためのルールベースのシステムを使っています。タンパク質とリガンドの接触している原子間の非共有結合性相互作用の特徴を見つけるために、特定の相互作用に寄与することができる化学官能基（例、水素結合のドナーに必要とされるもの）と相互作用の幾何学的特徴（例、距離と角度の閾値）についての文献に基づく情報を使います。各結合サイトについて、このアルゴリズムは最初に、特定の相互作用のパートナーとなりうるタンパク質とリガンドの原子、あるいは原子団を探します。ついで、相互作用を形成するタンパク質とリガンドの条件に適合したグループに対して幾何学的なルールが適用されます。

PLIPのSIからのFigure。ここでは芳香族スタッキングが描かれています。

より詳細な情報が必要なら、PDFドキュメントをチェックしてください！

# "ファイル > ライブPDFプレビューを使うためにこのNotebookを信頼する"をクリックする必要があるかもしれません
import requests
from IPython.display import HTML
pdf = "https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4489249/bin/supp_gkv315_nar-00254-web-b-2015-File003.pdf"
display(HTML(f"""
<iframe src="https://docs.google.com/viewer?url={pdf}&embedded=true"
     frameborder="0" webkitallowfullscreen mozallowfullscreen 
     allowfullscreen width="900" height="600">
</iframe>
"""))

訳注(2020/05)
ここにはPLIPの文献PDFが表示されていますが、サイズの問題で反映させていません。PubMedのこちらのページでSUPPLEMENTARY DATAをご確認ください。
訳注ここまで

偽陽性

ドッキングソフトウェアは最終的な答えをあたえるものではありません！全ての結果は厳密に検証しなければなりません。例えば、このレビューL.Yang, J.Zhang, X.Che, Y.Q.Gao, in Methods in Enzymology, 2016 では次の指摘がなされています。

タンパク質ーリガンド相互作用の研究においては、リガンドの自由度も重要な役割を果たす可能性があります。しかしながら、リガンド/タンパク質相互作用の研究においては、通常非常に限られた数のコンフォメーションしか考慮されません。一方で、リガンド/基質ータンパク質の相互作用を調べ、予測する上で、リガンドのコンフォメーションの集合を徹底的に調べることが極めて重要であることを示す例が増えています。さらに、構造ベースの医薬品デザインにおいて、テストされたリガンドのコンフォメーションが続いて行う計算の基礎となります。特に、リガンドは作用するタンパク質に結合する前と結合した後で、非常に異なるコンフォメーションをとることがあり、多くのリガンドが局所最小エネルギーのコンフォメーションではタンパク質に結合していないことを示す証拠があります。そのような場合、自由なリガンド単体のコンフォメーションの不適切な解析は、誤った結合エネルギー予測につながり、さらに酵素反応が含まれている場合には問題となる可能性があります。したがって、自由なリガンド単体の溶液中の構造を徹底的にサンプリングすることが必要で、それは適切な強化されたサンプリング手法と組み合わせたMDシミュレーションによって費用対効果良く実現することができます。

これはありうる失敗する点のうちの一つです。

---

実践

結果の可視化

nglviewを使いましょう！Webベースの分子ビューワーでJupyter Notebook上で実行することができます！また、ボックスから取り出したPDBQTファイルとも互換性があります（が、最初のモデルだけしか読み込みません。。。これをどうするかは後ほどとりあげます）。

nglviewをインストールするには次を実行してください。

!pip install nglview ipywidgets>=7.5  
!nglview install  
!nglview enable  

ノートブックでインタラクティブなGUIを作成するためにipywidgetsを使います。これにより、異なるリガンドをクリックし、クリックに応じてビューワーが更新されるようになります。特に、我々の小さなGUIには以下のことをさせたいと思います。

• Vinaの出力で報告されているように、ポーズと親和性のリストを表示する
• タンパク質構造をリボン表示、リガンドを Ball&Stick、そして周りの残基をlicorice(Stickのみ)で表示する
• 3Dの可視化がユーザーによるリストの異なるポーズの選択に応答するようにする

つまり、私たちは次のことをする必要があります。

1. NGLビューワーが適切な描画をするよう呼び出す
2. 結果のインタラクティブなテーブルを作成する（ヒント：ipywidgets.Selectを使ってください）
3. NGLビューワーとコミュニケーションできるイベントハンドラーを書く。つまりユーザーが新しいエントリをクリックすると、ディスプレイのリガンドと周囲の残基を更新するようにします（リボン表示は更新する必要がない）。

import pandas as pd
import time
import nglview as nv
# AppLayoutをつかうにはipywidgetsがv7.5+でなければなりません
from ipywidgets import AppLayout, Layout, Select, Button

Vinaの作成したPDBQTファイルはいくつかモデルを含んでいますが、nglviewは最初のものしか解析しません。この問題を回避する方法は単純です。ENDMDL行が見つかる度にスプリットすることでファイルを各モデルに分割します。

def split_pdbqt(path):
    """
    複数のモデルのPDBQTを別々のファイルに分けます
    """
    files = []
    with open(path) as f:
        lines = []
        i = 0
        for line in f:
            lines.append(line)
            if line.strip() == 'ENDMDL':
                fn = f'data/results.{i}.pdbqt'
                with open(fn, 'w') as o:
                    o.write(''.join(lines))
                files.append(fn)
                i += 1
                lines = []
    return files

Vinaの出力は結果のテーブルを含む単純なテキストファイルです。テーブルの解析は比較的簡単です。必要であれば、単純な可視化のためにpandas.DataFrameを返します。

def parse_output(out):
    """
    Vina出力ファイルからDataFrameを作成します
    """
    with open(out) as f:
        data = []
        for line in f:
            if line.startswith('-----+'):
                line = next(f)
                while line.split()[0].isdigit():
                    index, *floats = line.split()
                    data.append([int(index)] + list(map(float, floats)))
                    line = next(f)
    return pd.DataFrame.from_records(data, 
                                     columns=['Mode', 'Affinity (kcal/mol)', 'RMSD (l.b.)', 'RMSD (u.b.)'], 
                                     exclude=['Mode'])

Vinaの出力を解析

作成した関数をつかって次のテキストファイルを…

# !cat data/vina.out 
# T11bでオンラインの取得がうまくいかなかったためローカル環境で
# AutoDock Vinaを実施した結果を用いています。
!cat data/local_vina.out

#################################################################
# If you used AutoDock Vina in your work, please cite:          #
#                                                               #
# O. Trott, A. J. Olson,                                        #
# AutoDock Vina: improving the speed and accuracy of docking    #
# with a new scoring function, efficient optimization and       #
# multithreading, Journal of Computational Chemistry 31 (2010)  #
# 455-461                                                       #
#                                                               #
# DOI 10.1002/jcc.21334                                         #
#                                                               #
# Please see http://vina.scripps.edu for more information.      #
#################################################################

Detected 4 CPUs
Reading input ... done.
Setting up the scoring function ... done.
Analyzing the binding site ... done.
Using random seed: 267553816
Performing search ... done.
Refining results ... done.

mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1         -6.8      0.000      0.000
   2         -6.8      5.072      7.461
   3         -6.7      5.243      7.469
   4         -6.6      4.035      6.317
   5         -6.6      6.833      8.653
   6         -6.6      8.876     10.836
   7         -6.5      7.611      9.316
   8         -6.5      1.094      2.109
   9         -6.5      7.680      9.301
Writing output ... done.

綺麗にフォーマットされたpandas.DataFrameに変換することができます。

# parse_output("data/vina.out")
parse_output("data/local_vina.out")

	Affinity (kcal/mol)	RMSD (l.b.)	RMSD (u.b.)
0	-6.8	0.000	0.000
1	-6.8	5.072	7.461
2	-6.7	5.243	7.469
3	-6.6	4.035	6.317
4	-6.6	6.833	8.653
5	-6.6	8.876	10.836
6	-6.5	7.611	9.316
7	-6.5	1.094	2.109
8	-6.5	7.680	9.301

Vinaの出力には３つの列があります。

1. 予測された結合親和性
2. 対称性修正アルゴリズム（symmetry-corrected algorithm）を使った最良の解（親和性最小値）に対するRMSD
3. 対称性修正（symmetry correction）を用いていない最良の解（親和性最小値）に対するRMSD

私たちは親和性にだけ興味があるので、次のセルではpandas.DataFrameのために列を一つだけ取得しています。

---

それでは、NGLビューワーのインスタンスを作成します。各タンパク質ーポーズのペアについて新しいインスタンスを作成する代わりに、同じ描画面を通して再利用し、必要なリガンドを隠すか表示します。最初に全てを読み込み、リガンドをそれぞれの親和性でラベルづけします。ビューワーに読み込まれた各分子を「component」とよびます。最初にタンパク質を読み込むので、component_0となります。リガンドはcomponent_1から始まり続いていきます。

import time
def create_viewer(protein, ligands, affinities):
    """
    タンパク質と親和性のラベルをつけた全てのリガンドをもつnglview widgetを作成します
    """
    viewer = nv.NGLWidget()
    viewer.add_component(protein, ext="mol2")
    # 親和性のラベルを有する分子(@0)の最初の原子を選択します
    label_kwargs = dict(labelType="text", sele="@0", showBackground=True, backgroundColor="black")
    for ligand, affinity in zip(ligands, affinities):
        ngl_ligand = viewer.add_component(ligand, ext="pdbqt")
        ngl_ligand.add_label(labelText=[str(affinity)], **label_kwargs)
    return viewer

そして最後に、次の下のセルで実際のGUIを作成します！

ipywidgets.AppLayoutレイアウトを使って水平に配置された２つのwidgetsからなります。

• セレクター（ipywidgets.Select）
• NGLビューワーそのもの

ユーザーがセレクターで新しいエントリーをクリックすると、_on_selection_changeが呼ばれ以下のことを行います。

1. 新しい値が前のものと異なるかを確認します。異なる場合、次に…
2. 全てのリガンドを隠しビューをリセットします
3. 選択したものを表示し、500msの格好いいアニメーションとともにカメラの中心に合わせます
4. 新しいポーズの重心の5Å内にある残基のリストを更新するためにNGLビューワーにJavaScriptを実行します。

# Python widgetからは隠されているので、
# リガンドの周りの残基を更新するにはJavaScriptのコードが必要です。
# http://nglviewer.org/ngl/api/manual/snippets.html に基づきます
_RESIDUES_AROUND = """
var protein = this.stage.compList[0];
var ligand_center = this.stage.compList[{index}].structure.atomCenter();
var around = protein.structure.getAtomSetWithinPoint(ligand_center, {radius});
var around_complete = protein.structure.getAtomSetWithinGroup(around);
var last_repr = protein.reprList[protein.reprList.length-1];
protein.removeRepresentation(last_repr);
protein.addRepresentation("licorice", {{sele: around_complete.toSeleString()}});
"""

def show_docking(protein, ligands, vina_output):
    # 複数のPDBQTリガンドファイルを別々のファイルに分割する
    ligands_files = split_pdbqt(ligands)
    # 親和性（'Affinity'）の取得（DataFrameのこの列だけが必要です）
    affinities = parse_output(vina_output)['Affinity (kcal/mol)']
                                
    # viewer widgetの作成
    viewer = create_viewer(protein, ligands_files, affinities)
    
    # selection widgetの作成
    #   オプションは(text, value)のタプルのリストです。select上でクリックすると、
    #   `.observe(...)`に登録されている呼び出し可能なもの(callable)にvalueが渡されます。
    selector = Select(options=[(f"#{i} {aff} kcal/mol", i) for (i, aff) in enumerate(affinities, 1)],
                      description="",  rows=len(ligands_files), layout=Layout(width="auto"))
                 
    # GUI要素の配置
    # selectionボックスを左、viewerをウィンドウの残りの部分を占めるようにします。
    display(AppLayout(left_sidebar=selector, center=viewer, pane_widths=[1, 6, 1]))
    
    # これがイベントハンドラーです。ユーザーがselectionボックスでクリックした時に動作を取り出します。
    # viewerの変数を「見る」ことができるようにここで定義する必要があります。
    def _on_selection_change(change):
        # ユーザーが異なるエントリーでクリックした時だけ更新する
        if change['name'] == 'value' and (change['new'] != change['old']):
            viewer.hide(list(range(1,len(ligands_files) + 1)))  # 全てのリガンドを隠す（components 1-n）
            component = getattr(viewer, f"component_{change['new']}")
            component.show()  # 選択した一つを表示
            component.center(500)  # viewをズームする
            # リガンド周りの残基を表示するためにJSコードを呼ぶ
            viewer._execute_js_code(_RESIDUES_AROUND.format(index=change['new'], radius=5))
    
    # イベントハンドラーの登録
    selector.observe(_on_selection_change)
    # 最初の解に焦点をあてるために手動でイベントを起こす
    _on_selection_change({'name': 'value', 'new': 1, 'old': None})

    return viewer

上の関数ではPythonの応用をいくつか行なっています。もし興味があれば、下のヘッダーをクリックしてより詳細について見ることができます。

Python応用の説明

最初に、JavaScriptとPythonのコードを混ぜていることに気づいたかもしれません！これはNGLViewer.Viewer._execute_js_codeで提供されているように、ipywidgetsの基盤のおかげで可能になっています。ここで（JavaScriptのコードを含む）文字列を渡しておくと、後でNGL widgetのスコープで実行されます（thisはインタラクティブな描画面を差します）。パラメーター化するために、_on_selection_changeを呼ぶ度にフォーマット化されるテンプレート・プレースホルダーを追加しました。

この関数_on_selection_changeがユーザーのインタラクション（selectionボックス上でのクリック）をPythonの世界と繋げる糊の役割を果たします。selectionボックス上でクリックする度ごとに呼ばれ、引数changeを一つもちます。イベントのタイプと新旧両方の値を保持する辞書です（したがって比較し、それを取りあつかうことができます）。しかしながら、vewerを関数に含める必要もあります…そして、追加の値を渡すことはできません！

関数でviewerを使えるようにする一つの方法はshow_docking関数にその定義をネストすることです。そうすることで、Notebookの(global)scopeの中を探す必要なく外側の領域にアクセすることができるようになります。 このポストでscopeとclosureについて情報をもっと見ることができます。

ドッキング結果を見る

# viewer = show_docking("data/protein.mol2", "data/results.pdbqt", "data/vina.out")
viewer = show_docking("data/protein.mol2", "data/local_result.pdbqt", "data/local_vina.out")

訳注(2020/05)
以下は画像ですが、ノートブック上ではインタラクティブなビューワーが表示されます。
訳注ここまで

viewer.render_image(),

(Image(value=b'', width='99%'),)

# 静的な出力
viewer._display_image()

import time
time.sleep(10)

タンパク質ーリガンド相互作用の解析

PLIPは自動化された解析のためのWebサーバー を提供していますが、残念ながらAPIはありません。標準のWeb UIを使っているかのようにHTML形式を使ってみることもできますが、ライブラリ自体がPython-3で使えるように準備されており、pipで簡単にインストールすることができるので、簡単のためローカルで使うことができます。

!conda install -y -c openbabel openbabel lxml    
!pip install https://github.com/ssalentin/plip/archive/stable.zip --no-deps  

from plip.modules.preparation import PDBComplex
from plip.modules.report import BindingSiteReport
from glob import glob
import os
import pandas as pd

PILIPはPDBファイルだけを受け入れますが、PDBQTファイルはPDBの拡張された記法なので、相互の変換は列とレコードのうち正しい部分を選択するだけの問題です。この場合、ATOM（タンパク質原子）とHETATM（リガンド原子）レコードにだけ興味があるので、追加のフィールド（66列以降）を廃棄します。ついで、単に連結を行うだけでそれらを合わせてまとめることができます。

def pdbqt_to_pdbblock(pdbqt):
    lines = []
    with open(pdbqt) as f:
        for line in f:
            if line[:6] in ('ATOM  ', 'HETATM'):
                lines.append(line[:67].strip())
    return '\n'.join(lines)

def create_protein_ligand_pdbs(protein_pdbqt, ligands_pdbqt):
    ligands_pdbqt = sorted(list(glob(ligands_pdbqt)), key=lambda s: int(s.split('.')[-2]))
    protein_pdbblock = pdbqt_to_pdbblock(protein_pdbqt)
    for i, ligand_pdbqt in enumerate(ligands_pdbqt, 1):
        ligand_pdbblock = pdbqt_to_pdbblock(ligand_pdbqt)
        with open(f"data/docked_protein_ligand.{i}.pdb", 'w') as f:
            f.write(protein_pdbblock)
            f.write("COMPND    UNNAMED\n")
            f.write(ligand_pdbblock)

ジャジャーン!できました！

create_protein_ligand_pdbs("data/protein.mol2.pdbqt", "data/results.*.pdbqt")

これでPDBファイルをPLIPに渡し、その魔法が使えるようになりました。BindingSiteReportクラスはPDBComplexの検出された各結合サイトを処理し、興味のあるフィールドをもつオブジェクトを作成します。

• hydrophobic
• hbond
• waterbridge
• saltbridge
• pistacking
• pication
• halogen
• metal

これらのフィールドは（列名を保持する）<field>_featuresと（実際のレコードを保持する）<field>_infoに分けられます。オブジェクトを順に処理してgetattr()で正しい属性（attribute）の名前を取得すれば、概要をよくみるためにpandas.DataFrameに渡す辞書を構成することができます。以下をチェックしてください。

def analyze_interactions(pdbfile):
    protlig = PDBComplex()
    protlig.load_pdb(pdbfile)  # pdbの読み込み
    for ligand in protlig.ligands:
        protlig.characterize_complex(ligand)  # リガンドを見つけ相互作用を解析する
    sites = {}
    for key, site in sorted(protlig.interaction_sets.items()):
        binding_site = BindingSiteReport(site)  # 相互作用に関するデータをあつめる
        # *_featuresと*_infoのタプルをpandas dfに変換します
        keys = "hydrophobic", "hbond", "waterbridge", "saltbridge", "pistacking", "pication", "halogen", "metal"
        interactions = {k: [getattr(binding_site, k+"_features")] + getattr(binding_site, k+"_info") for k in keys}
        sites[key] = interactions
    return sites

interactions_by_site = analyze_interactions("data/docked_protein_ligand.1.pdb")
interactions_by_site

PLIP解析結果の辞書をみるにはここをクリックしてください

{'A:A:1': {'hydrophobic': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'DIST',
    'LIGCARBONIDX',
    'PROTCARBONIDX',
    'LIGCOO',
    'PROTCOO')],
  'hbond': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'SIDECHAIN',
    'DIST_H-A',
    'DIST_D-A',
    'DON_ANGLE',
    'PROTISDON',
    'DONORIDX',
    'DONORTYPE',
    'ACCEPTORIDX',
    'ACCEPTORTYPE',
    'LIGCOO',
    'PROTCOO'),
   (45,
    'LYS',
    'A',
    1,
    'A',
    'A',
    True,
    '2.64',
    '3.42',
    '133.60',
    True,
    431,
    'N3+',
    2953,
    'O3',
    (3.232, 17.45, 38.708),
    (4.782, 15.554, 36.326)),
   (150,
    'ASP',
    'A',
    1,
    'A',
    'A',
    False,
    '3.07',
    '3.43',
    '102.37',
    True,
    1470,
    'Nam',
    2955,
    'O3',
    (1.743, 18.722, 37.254),
    (0.523, 19.864, 34.258)),
   (45,
    'LYS',
    'A',
    1,
    'A',
    'A',
    False,
    '2.66',
    '3.10',
    '106.87',
    True,
    422,
    'Nam',
    2951,
    'O3',
    (5.392, 21.015, 40.948),
    (7.994, 19.469, 41.639)),
   (72,
    'LEU',
    'A',
    1,
    'A',
    'A',
    False,
    '3.65',
    '4.01',
    '103.83',
    True,
    694,
    'Nam',
    2945,
    'Npl',
    (3.257, 24.961, 34.128),
    (2.841, 27.919, 36.804)),
   (61,
    'MET',
    'A',
    1,
    'A',
    'A',
    False,
    '2.75',
    '3.76',
    '161.43',
    False,
    2945,
    'Npl',
    595,
    'O2',
    (3.257, 24.961, 34.128),
    (2.892, 28.611, 33.325)),
   (149,
    'THR',
    'A',
    1,
    'A',
    'A',
    True,
    '2.02',
    '2.89',
    '147.03',
    False,
    2955,
    'O3',
    1467,
    'O3',
    (1.743, 18.722, 37.254),
    (0.131, 21.066, 36.729)),
   (150,
    'ASP',
    'A',
    1,
    'A',
    'A',
    True,
    '2.70',
    '3.64',
    '162.21',
    False,
    2953,
    'O3',
    1478,
    'O3',
    (3.232, 17.45, 38.708),
    (1.632, 16.073, 35.738))],
  'waterbridge': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'DIST_A-W',
    'DIST_D-W',
    'DON_ANGLE',
    'WATER_ANGLE',
    'PROTISDON',
    'DONOR_IDX',
    'DONORTYPE',
    'ACCEPTOR_IDX',
    'ACCEPTORTYPE',
    'WATER_IDX',
    'LIGCOO',
    'PROTCOO',
    'WATERCOO')],
  'saltbridge': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'DIST',
    'PROTISPOS',
    'LIG_GROUP',
    'LIG_IDX_LIST',
    'LIGCOO',
    'PROTCOO')],
  'pistacking': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'CENTDIST',
    'ANGLE',
    'OFFSET',
    'TYPE',
    'LIG_IDX_LIST',
    'LIGCOO',
    'PROTCOO')],
  'pication': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'DIST',
    'OFFSET',
    'PROTCHARGED',
    'LIG_GROUP',
    'LIG_IDX_LIST',
    'LIGCOO',
    'PROTCOO')],
  'halogen': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'SIDECHAIN',
    'DIST',
    'DON_ANGLE',
    'ACC_ANGLE',
    'DON_IDX',
    'DONORTYPE',
    'ACC_IDX',
    'ACCEPTORTYPE',
    'LIGCOO',
    'PROTCOO')],
  'metal': [('RESNR',
    'RESTYPE',
    'RESCHAIN',
    'RESNR_LIG',
    'RESTYPE_LIG',
    'RESCHAIN_LIG',
    'METAL_IDX',
    'METAL_TYPE',
    'TARGET_IDX',
    'TARGET_TYPE',
    'COORDINATION',
    'DIST',
    'LOCATION',
    'RMS',
    'GEOMETRY',
    'COMPLEXNUM',
    'METALCOO',
    'TARGETCOO')]}}

この辞書は２つのレベルからなります。最初のレベルは見つかった結合サイトです。一つのリガンドだけをドッキングしているので、一つだけ見つかることが予想されます。各結合サイトについて、８個のリストを含むもう一つの辞書があり、リストは各特定の相互作用のためのものです。各リストは最初の列に列名を持ち、つづいて（該当するものがあれば）データが格納されています。これを踏まえて、次のようにPandas.DataFrameを構築することができます。

# オリジナルのノートブックと異なる解析結果を利用しているため、結合サイトの名前も異なっています。
# 結果に合わせてオリジナルのコード ["UNL:Z:1"] を['A:A:1']に変更しました。
# df = pd.DataFrame.from_records(interactions_by_site["UNL:Z:1"]["hbond"][1:], # データは列名の後に格納されています。
#                                columns=interactions_by_site["UNL:Z:1"]["hbond"][0]) # 列名は常に最初の要素です。

df = pd.DataFrame.from_records(interactions_by_site["A:A:1"]["hbond"][1:], 
                               columns=interactions_by_site["A:A:1"]["hbond"][0]) 

df

	RESNR	RESTYPE	RESCHAIN	RESNR_LIG	RESTYPE_LIG	RESCHAIN_LIG	SIDECHAIN	DIST_H-A	DIST_D-A	DON_ANGLE	PROTISDON	DONORIDX	DONORTYPE	ACCEPTORIDX	ACCEPTORTYPE	LIGCOO	PROTCOO
0	45	LYS	A	1	A	A	True	2.64	3.42	133.60	True	431	N3+	2953	O3	(3.232, 17.45, 38.708)	(4.782, 15.554, 36.326)
1	150	ASP	A	1	A	A	False	3.07	3.43	102.37	True	1470	Nam	2955	O3	(1.743, 18.722, 37.254)	(0.523, 19.864, 34.258)
2	45	LYS	A	1	A	A	False	2.66	3.10	106.87	True	422	Nam	2951	O3	(5.392, 21.015, 40.948)	(7.994, 19.469, 41.639)
3	72	LEU	A	1	A	A	False	3.65	4.01	103.83	True	694	Nam	2945	Npl	(3.257, 24.961, 34.128)	(2.841, 27.919, 36.804)
4	61	MET	A	1	A	A	False	2.75	3.76	161.43	False	2945	Npl	595	O2	(3.257, 24.961, 34.128)	(2.892, 28.611, 33.325)
5	149	THR	A	1	A	A	True	2.02	2.89	147.03	False	2955	O3	1467	O3	(1.743, 18.722, 37.254)	(0.131, 21.066, 36.729)
6	150	ASP	A	1	A	A	True	2.70	3.64	162.21	False	2953	O3	1478	O3	(3.232, 17.45, 38.708)	(1.632, 16.073, 35.738)

さて、これらの相互作用をNGLビューワーで表しましょう。相互作用ポイント（pandas.DataFrameのLIGCOOとPROTCOO）の間に円柱（cylinder）を描き、color_mapに表されているように色分けします（RGBタプルを使います）。

color_map = {
    'hydrophobic': [0.90, 0.10, 0.29],
    'hbond': [0.26, 0.83, 0.96],
    'waterbridge': [1.00, 0.88, 0.10],
    'saltbridge': [0.67, 1.00, 0.76],
    'pistacking': [0.75, 0.94, 0.27],
    'pication': [0.27, 0.60, 0.56],
    'halogen': [0.94, 0.20, 0.90],
    'metal': [0.90, 0.75, 1.00],
}

def show_docking_and_interactions(protein, ligands, vina_output):
    # 複数のPDBQTリガンドファイルを別々のファイルに分割します
    ligands_files = split_pdbqt(ligands)
    # 親和性（'Affinity'）を取得します（ DataFrameのこの列だけが必要です）
    affinities = parse_output(vina_output)['Affinity (kcal/mol)']
                                
    # selection widgetの作成
    #   オプションは(text, value)のタプルのリストです。select上でクリックすると、
    #   `.observe(...)`に登録されている呼び出し可能なもの(callable)にvalueが渡されます。
    selector = Select(options=[(f"#{i} {aff} kcal/mol", i-1) for (i, aff) in enumerate(affinities, 1)],
                      description="",  rows=len(ligands_files), layout=Layout(width="auto"))
    
                 
    # GUI要素の配置
    # selectionボックスを左、viewerをウィンドウの残りの部分を占めるようにします。
    app = AppLayout(left_sidebar=selector, center=None, pane_widths=[1, 6, 1], height="600px")

    # これがイベントハンドラーです。ユーザーがselectionボックスでクリックした時に動作を取り出します。
    # viewerの変数を「見る」ことができるようにここで定義する必要があります。
    def _on_selection_change(change):
        # ユーザーが異なるエントリーでクリックした時だけ更新する
        if change['name'] == 'value' and (change['new'] != change['old']):
            if app.center is not None:
                app.center.close()

            # NGLビューワー
            app.center = viewer = nv.NGLWidget(height="600px", default=True, gui=True)
            prot_component = viewer.add_component(protein, ext="pdbqt") # add protein
            time.sleep(1)
            lig_component = viewer.add_component(ligands_files[change['new']], ext="pdbqt") # add selected ligand
            time.sleep(1)
            lig_component.center(duration=500)

            # 相互作用の追加
            interactions_by_site = analyze_interactions(f"data/docked_protein_ligand.{change['new']+1}.pdb")
            interacting_residues = []                                    
            for site, interactions in interactions_by_site.items():
                for interaction_type, interaction_list in interactions.items():
                    color = color_map[interaction_type]
                    if len(interaction_list) == 1:
                        continue
                    df_interactions = pd.DataFrame.from_records(interaction_list[1:], columns=interaction_list[0])
                    for _, interaction in df_interactions.iterrows():
                        name = interaction_type
                        viewer.shape.add_cylinder(interaction['LIGCOO'], interaction['PROTCOO'], color, [0.1], name) 
                        interacting_residues.append(interaction['RESNR'])
            # 相互作用する残基の表示
            res_sele = " or ".join([f"({r} and not _H)" for r in interacting_residues])
            res_sele_nc = " or ".join([f"({r} and ((_O) or (_N) or (_S)))" for r in interacting_residues])
            prot_component.add_ball_and_stick(sele=res_sele, colorScheme="chainindex", aspectRatio=1.5)
            prot_component.add_ball_and_stick(sele=res_sele_nc, colorScheme="element", aspectRatio=1.5)
                                                        
                                                        
    # イベントハンドラーの登録
    selector.observe(_on_selection_change)
    # 最初の解に焦点をあてるために手動でイベントを起こす
    _on_selection_change({'name': 'value', 'new': 1, 'old': None})

    return app

show_docking_and_interactions はほとんどshow_docking関数とと同じですが、各選択イベントについてPLIP相互作用も計算します。この場合、イベントハンドラーはもう少しだけ関与しています。相互作用を示す円筒を手動で消すことができないので、各選択イベントについてNGLビューワーを作り直す必要があります。ですが、今回はJavaScriptは含まれていません。純粋なPythonだけです。ですが、ビューワーのwidgetを初めから作り直し、分子を読み込み、相互作用を計算するにはより時間が必要となり、体感としてはそれほどダイナミックなものにはなりません。

ドッキングと相互作用の表示

# app = show_docking_and_interactions("data/protein.mol2.pdbqt", "data/results.pdbqt", "data/vina.out")
app = show_docking_and_interactions("data/protein.mol2.pdbqt", "data/local_result.pdbqt", "data/local_vina.out")
app

訳注(2020/05)
以下は画像ですが、ノートブック上ではインタラクティブなビューワーが表示されます。
訳注ここまで

app.center.render_image(),

(Image(value=b'', width='99%'),)

# 静的なプレビュー
app.center._display_image()

ディスカッション

このトークトリアルの最後のパートでは、Vinaによるドッキングの結果を可視化し解析する方法を学びました。他のパートと異なり、このNotebookはローカル環境で実行され、Webサーバーは呼びませんでした。しかしながら、とてもたくさんのWebテクノロジーを取り扱っており、主にNGLビューワーとipywidgetsレイアウトに基づいています。

はじめに、選択したタンパク質ーリガンドポーズとその周囲でインタラクティブにビューワーを更新するために(observe()で登録されているように）イベントハンドラーを使いました。次に、ビューワーで相互作用を描くためにPLIPの自動化された解析を使って可視化を拡張しました。

クイズ

• 観察したドッキング解析において最も重要だったのはどのような相互作用でしたか？文献との一致は見られましたか？またそれはなぜでしょうか？
• 疎水性相互作用と水素結合の間の主な違いは何でしょうか？またどの点で類似しているでしょうか？
• パートBで学んだ内容を応用してPLIPのWebサーバーを使う関数を描いてください。メインのformエントリーとactionのエンドポイント、リクエストを投げるのに必要なinputタグの場所を見つけるためにHTMLソースコードをチェックする必要があります。

訳以上

追記

以上です。今回も色々と誤りがあると思います。ご指摘いただければ幸いです。

以下、今回のトークトリアルでよく分からなかった箇所とエラーについて調べた内容を記載します。

• 追記 1) AutoDockVinaの結果、RMSDのsymmetry correction に関して

RMSDについてsymmetry correctionを行ったか否かの２種類の記述があります。以下は、私の理解です。
リガンドの各原子を区別してラベルしている場合、分子に対称性があると、化学的には同じ構造でも原子のラベルとしては異なるという可能性があります（ex. ベンゼンの６つの炭素原子にラベルをつける場合）。ドッキングの結果の比較として、同じリガンドの異なる結合ポーズをRMSDで評価しようとする場合、各原子のラベルの対応を元に計算できます。ですが、対称性がある分子では、実際には重なりのよい（RMSDの小さい）ポーズでも、ラベルの付け方が複数通り可能なために原子同士が離れている（RMSDが大きくなる）と計算してしまいます。この問題を考慮するか否か、というのがsymmetry correctionを行うか否か、という違いだと思います。

• 追記2) PLIP実行時のエラーについて

PLIPでPDBを読み込む際に下記のエラーが出ました。　　

~/.pyenv/versions/anaconda3-5.3.1/envs/teachopencadd/lib/python3.6/site-packages/plip/modules/supplemental.py in read_pdb(pdbfname, as_string) in read_pdb(pdbfname, as_string)
    396     if os.name != 'nt':  # Resource module not available for Windows
    397         maxsize = resource.getrlimit(resource.RLIMIT_STACK)[-1]
--> 398         resource.setrlimit(resource.RLIMIT_STACK, (min(2 ** 28, maxsize), maxsize))
    399     sys.setrecursionlimit(10 ** 5)  # increase Python recoursion limit
    400     return readmol(pdbfname, as_string=as_string)
    
ValueError: current limit exceeds maximum limit   

使用するシステムのリソースを設定する箇所でのエラーです。こちらはどうもPLIPの問題ではなく、Mac OSとPython 3.6の組み合わせで生じる問題のようです。Pythonにissue（　課題34602:python3 resource.setrlimit strange behaviour under macOS ）　が出ていました。PLIPの該当箇所をコメントアウトしGitHubの修正履歴を参考に以下のように書き換えました。

def read_pdb(pdbfname, as_string=False):
    """Reads a given PDB file and returns a Pybel Molecule."""
    pybel.ob.obErrorLog.StopLogging()  # Suppress all OpenBabel warnings

    """エラー箇所
    if os.name != 'nt':  # Resource module not available for Windows
       maxsize = resource.getrlimit(resource.RLIMIT_STACK)[-1]
       resource.setrlimit(resource.RLIMIT_STACK, (min(2 ** 28, maxsize), maxsize))
    """

    # 変更ここから
    if sys.platform == 'darwin':
        try:
            import resource
        except ImportError:
            pass
        else:
            soft, hard = resource.getrlimit(resource.RLIMIT_STACK)
            newsoft = min(hard, max(soft, 1024*2048))
            resource.setrlimit(resource.RLIMIT_STACK, (newsoft, hard))
    # 変更ここまで
    
    sys.setrecursionlimit(10 ** 5)  # increase Python recoursion limit
    return readmol(pdbfname, as_string=as_string)

これでエラーが出なくなりました。

• 追記3) ビューワー上のリガンドに関して

ビューワーを画像にしてしまっているためわかりにくいですが、ドッキング結果のリガンドは結合におかしな点がみられます。ドッキング結果を出力したPDBファイルには原子の座標のみで、結合関係の情報が入っていないようでしたので、その辺りが影響しているのかもしれません。

終わりに

これでトークトリアル１１含め全てのトークトリアル終了となります。長々とお付き合いいただきありがとうございました！

私は、CADDに関連する項目についてはぼんやりと用語を聞いたことがある程度だったので、各トピックごとに理論と実践の両方が扱われていて、まとめて勉強する良い機会となりました。一つの標的EGFRに絞って、色々なトピックでアプローチすることで、様々な視点から情報を捉えるというのは面白いやり方ですね。恐らく内容としては初歩の初歩で、CADDのほんの入り口に過ぎないのだと思いますが、全体の流れと道しるべ的な専門用語を拾えたので、何が話題になっているか程度には理解できるようになっているといいなと思います。ではでは！

オフラインでドッキングする話（T11パートBの備忘録）

T11パートBはオンラインWebサービスを用いたターゲットタンパク質に対するドッキングでした。残念ながら用いられていたSwissDock上のドッキング、OPAL共に現在では利用ができなくなっていました。そこで、趣旨とは異なってしまいますがオフラインでのドッキングを行って見たいと思います。

目標はこんな感じです。

用いるデータ	概要	詳細
参照構造	EGFR-リガンド共結晶構造	PDB ID:3W33、但しKLIFSに登録されているデータ
ドッキングする化合物	リガンドと構造類似の化合物群	リガンド（ADP）をテンプレートとしてPubChemから取得

やること	T11パートB	この記事
結合サイトの位置	DoGSiteScorer	同左
タンパク質・リガンドの準備	AutoDockTools	ODDT
ドッキング	OPAL上のAutoDockVina	PCにインストールしたAutoDockVina

では、素人による適当解説やっていきますよ！

DogSiteScorerで良い感じのポケットを探す

Webサイトで遊ぶ

まずはドッキングを行うポケットの位置を決めます。ここはT11パートBでも実行できていたのでやり方を踏襲します。・・・が、やっぱりコードだけをみてもよくわからないのでWebサイトではどのような結果になるか見に行きたいと思います。

DoGSiteはProteinsPlusというサイトで提供されているソフトです*1。ProteinsPlusはハンブルク大学Computational Molecular Design Group（AMD、Matthias Rarey教授）で開発されたツールを利用できるように公開されているWebサーバーで、ポケットを探索するDoGSite以外にもタンパク質の構造を取り扱うためのツールが多数組み込まれています*2。

トップページ格好良い・・・

ProteinPlusトップページ

PDB IDでの検索も可能なようですが、今回はT11に従いKLIFSから取得したデータを使うため、Upload Proten(PDB format)にPDBファイルを投げ込みGo!。すぐに結果がでてこんなページが出てきます。

素敵・・・。

右側のパネルから利用できるサービスを選べるようになっています。DoGSiteScorerに行きましょう。色々とオプションがありますが今回はデフォルトのまま計算して見ます。

DoGSiteScorerの設定

1−２分で結果が出ます。今回は10個のポケットが見つかりました。

DoGSiteScorerによるポケット探索

ポケットのサイズに関する情報（体積、表面積）に加えてスコア（Drug Score、Simple Score）も出ています。スコアが高ければ医薬品の標的とするポケットとして良いのだと思います（論文を読んでいません。すみません）。確かにトップのポケット（P_0）は共結晶構造中のリガンドの位置と一致しています。

ポケットの情報はテーブルをスクロールした下部にあるDownloadから取得できます（zipファイル）。ありがたく頂戴して中身を見て見ましょう。

ダウンロードした情報の確認

ダウンロードしたディレクトリには以下のファイルが含まれています。

1. 「parameter.txt」・・・計算の設定　　
2. 「ジョブ名_desc.txtファイル」・・・ポケットの情報（Webサイト上に表示されているテーブル）。　　
3. 「PocXlsDescriptors.txt」・・・ポケット情報の省略語の説明
4. 「pockets」・・・各ポケットの座標等の情報のファイルを含むディレクトリ（.ccp4.gz、ポケットごとに別のファイル（今回は10個））
5. 「residues」・・・各ポケットの残基等の情報のファイルを含むディレクトリ（.pdb、ポケットごとに別のファイル（今回は10個））

Webサイト上で画像表示されていたもののベースになる情報がそっくり取得できるようです。

ポケット情報のテキストファイル（2）はタブ区切りのテーブルのようでした。Pandasの出番だすな！

import pandas as pd
df = pd.read_table('./DogSiteResults/Job_desc.txt')  # ファイルの名前は長いので変えています
print(df.shape)
df.head()
#  (10, 50)

	name	lig_cov	poc_cov	lig_name	volume	enclosure	surface	depth	surf/vol	lid/hull	...	MET	PHE	PRO	SER	THR	TRP	TYR	VAL	simpleScore	drugScore
0	P_0	0.0	0.0	NaN	1268.67	0.13	1422.18	26.50	1.121001	-	...	2	3	2	1	2	1	0	4	0.64	0.817267
1	P_1	0.0	0.0	NaN	371.84	0.17	681.06	13.45	1.831594	-	...	2	1	3	0	1	0	1	2	0.23	0.642846
2	P_2	0.0	0.0	NaN	338.56	0.07	524.19	13.65	1.548293	-	...	1	1	1	3	1	2	1	2	0.18	0.625701
3	P_3	0.0	0.0	NaN	234.94	0.10	473.53	10.68	2.015536	-	...	0	2	2	1	0	0	1	1	0.11	0.421320
4	P_4	0.0	0.0	NaN	207.36	0.12	353.27	11.60	1.703655	-	...	1	1	0	0	0	0	1	2	0.04	0.444417

5 rows × 50 columns

ポケットにどんな残基がいくつ含まれているかも含め、たくさんの情報が入っているようです。

具体的なポケットをPyMolで確認してみます。

ポケットの可視化

ccp4ファイルはObjectMapCCp4として認識され、立方体が読み込まれました。ポケットを単位胞のように認識しているのでしょうか？・・・不勉強なのでわかりません。すみません。

pdbファイルのHEADERにはポケットの中心座標の情報も記載されています。例えば今回のポケットP_０では・・・

• HEADER Geometric pocket center at 0.63 18.01 37.64 with max radius 22.90

でした。この情報があればドッキング計算でポケットの位置を指定できそうです。

APIで実行

同じ作業をAPIで行ってみます。といってもT11パートBのコードを切りはりするだけですが・・・。

PDBファイルをアップロードするの複雑なのでそのままPDB IDを使ってやってみます（PDB ID: 3W33）。ProteinPlusはREST APIのドキュメンテーションもあり、DogSiteScorerはこちらです。POSTでジョブを作成し、GETで結果を取得します(json形式)。

import requests

r = requests.post("https://proteins.plus/api/dogsite_rest",
                  json={
                      "dogsite": {
                          "pdbCode": "3W33",
                          "analysisDetail": "1",
                          "bindingSitePredictionGranularity": "1",
                          "ligand": "",
                          "chain": ""
                      }
                  },
                  headers= {'Content-type': 'application/json', 'Accept': 'application/json'}
                 )

r.raise_for_status() # エラーに備える
print(r.status_code)
# 200

HTTP ステータス：200なのでうまくいっているようです。ドキュメンテーションには各ステータスコードの場合に得られる結果も記載されています。

POSTが成功した場合の結果

実際にtext属性を確認してみると同じ結果となっています。

print(r.text)
#  {"status_code":200,"location":"https://proteins.plus/api/dogsite_rest/ジョブのID","message":"Job already exists"}

ジョブがうまく投げられたので計算結果を取得します。locationがジョブのIDなのでこちらを使います。

job_location = r.json()['location']
result = requests.get(job_location)
print(result.status_code)
# 200

こちらもうまくいったようです。結果は以下となっているそうです。

GETの結果

このうちresiduesの一番最初(インデックス0)が最もスコアの良いポケットだそうです。Webで結果をダウンロードした時と同じですね。
residuesはPDBファイル形式となっています。今欲しい情報は結合ポケットの中心の座標と半径で、こちらはHEADERのGeometric pocket center ~という行に含まれていました。こちらを取得します。

# PDBの取り出し
best_pocket_pdb = requests.get(result.json()['residues'][0]).text

# PDBから目的の情報を含む行を見つける
# スペースで区切ったインデックス5,6,7（[5:8]）が中心の座標（x,y,z）
# 最後（インデックス[-1]）が半径(radius)
for line in best_pocket_pdb.splitlines():
    line = line.strip()
    if line.startswith('HEADER') and 'Geometric pocket center at' in line:
        fields = line.split()
        center = [float(x) for x in fields[5:8]]
        radius = float(fields[-1])
        
print("center: ", center)
print("radius: ", radius)
# center:  [20.37, 29.52, 10.77]
# radius:  19.95

先ほどWebインターフェースで実行した場合と座標、半径が異なっています。KLIFSから取得したファイルとPDB IDを指定して行った違いのためでしょうか？同じ構造でも情報の取得減に応じて座標が異なる可能性があるというのは要注意そうです。

以降ではT11でKLIFSから取得したファイルを使うので、先にWebで得たポケットの情報を使うことにしたいと思います。

ドッキングの実行

ドッキングソフトの取得

結合ポケットの情報がわかったので実際のドッキングに進みたいと思います。ドッキングはAutoDockVinaで行いますが、AutoDockVinaを使うためにはタンパク質とリガンドをそれぞれ専用の形式（PDBQT）で準備する必要があります*3。T11パートAではAutoDockToolsを修正したものが使われていましたが、今回は別のツールODDT（Open Drug Dicovery Toolkit）を利用したいと思います。

ODDTは様々なオープンソースのツールの入出力方法を統合してもっと便利に使えるようにしましょう、という素敵なツールです。ツール多すぎて無理・・・って感じの私にはとてもありがたいです。詳細は以下をご参照ください。

• Wójcikowski, M., Zielenkiewicz, P., and Siedlecki, P. Open Drug Discovery Toolkit (ODDT): a new open-source player in the drug discovery field Journal of Cheminformatics 2015(7)26
• 文献URL
• ドキュメント
• GitHub

インストールはpipかcondaでできます。

# pipでインストールする場合  
pip install oddt  

# condaでインストールする場合  
conda install -c oddt oddt  

AutoDockVina自体はこちらのページのDownloadから取得できます。計算化学.COMさんのこちらのページ にインストールと使用方法が詳しく書かれいるのでご参照ください*4。

タンパク質の準備

ではまずタンパク質だけを含むファイル（mol2）をODDTで読み込んでみます。ODDTではタンパク質の変数名をproteinとするようにとドキュメントに書かれていたので合わせておきます。*5

import oddt
from oddt import toolkit

protein = next(oddt.toolkit.readfile('mol2', "data/protein.mol2"))
protein.protein = True
type(protein)
#  oddt.toolkits.ob.Molecule

タンパク質の情報を読み込みオブジェクトの生成ができたので、こちらをpdbqtファイルに変換します。ODDTのdockingモジュールからAutodockVina.write_vina_pdbqt(molオブジェクト, 出力先のディレクトリ)を実行します。タンパク質の場合は引数flexible=Falseとしておく必要があります。

from oddt import docking

oddt.docking.AutodockVina.write_vina_pdbqt(protein, "data", flexible=False)
# 'data/3w33_A.pdbqt'

タンパク質の情報（PDB ID:3w33のChain A）を反映したファイル名で出力されました。変数proteinの中身をこちらに切り替えておきます。

protein = next(oddt.toolkit.readfile('pdbqt', "data/3w33_A.pdbqt"))
protein.protein = True
print(type(protein))
# oddt.toolkits.ob.Molecule

実際にはODDTはAutoDockVinaのドッキングを行う際にpdbqtへの変換も同時に行うので予め変換しておく必要はありませんが、ここでは使い方の練習ということでご容赦ください。

リガンドの準備

リガンドはT11パートAで取得したPubChem化合物を使います。SMILESで保存しているので、RDKitで３次元構造におこします。AllChem.EmbedMolecule()で３次元にしますが、性能が良いというETKDGv2()を指定して構造を生成します。*6

# SMILESファイルの読み込み
SMILES_FILE = "data/similar_smiles.txt"
with open(SMILES_FILE) as f:
    smiles = [line.strip() for line in f]

# 必要なものをインポート
from rdkit import Chem
from rdkit.Chem import AllChem

# ３次元構造を生成してリストに格納
mols = []

for s in smiles:
    mh = Chem.AddHs(Chem.MolFromSmiles(s)) 
    AllChem.EmbedMolecule(mh, AllChem.ETKDGv2())
    mols.append(mh)

# 描画で確認
from rdkit.Chem import Draw

Draw.MolsToGridImage(mols, molsPerRow=5)

立体感でてますね！SDFに書き出しておきます。

writer = Chem.SDWriter("data/similar_ligands.sdf")
for mol in mols:
    writer.write(mol)
writer.close()

ODDTでリガンドを読み込みPDBQTファイルに変換します。

# SDFの読み込み
ligands = list(oddt.toolkit.readfile('sdf', 'data/similar_ligands.sdf'))

# PDBQTに一つずつ変換
ligand_id = 0

for ligand in ligands:
    file_name = "ligand_" + str(ligand_id)
    oddt.docking.AutodockVina.write_vina_pdbqt(ligand, "data", name_id=file_name)
    ligand_id += 1

dataディレクトリの中に10個のリガンドpdbqtファイルができました。次のドッキングでは一番最初のリガンドを使って試そうと思うので読み込んでおきます。

ligand = next(oddt.toolkit.readfile('pdbqt', 'data/ligand_0.pdbqt'))

ドッキングの実行

結合ポケットの情報、タンパク質、リガンド構造の準備が完了したのでドッキングを実行します。ドッキングを行う空間はDoGSiteScorerのWeb版で実行しPyMolでも確認した以下のポケットです。

• 中心座標・・・ [0.63, 18.01, 37.64]
• サイズ・・・ 22.90 → 23にします。

from oddt.docking import AutodockVina

# パラメーターの設定
DockingParams = AutodockVina.autodock_vina(protein=protein, 
                                           size = (23, 23, 23), 
                                           center = [0.63, 18.01, 37.64], 
                                           exhaustiveness = 8, num_modes = 9, energy_range = 3,
                                           n_cpu = 3, 
                                           executable='../../../../autodock_vina_1_1_2_mac/bin/vina') 

引数executableはインストールしたAutodockVinaのlocalのpathを指定します。ここでは相対パスで指定しています。パラメーターの設定が完了したのでリガンドを渡してドッキングを実行します。

DockingResult = DockingParams.dock(ligand)

len(DockingResult)
#  9

9個の結果が生成されました。それぞれODDTのMolオブジェクトで、スコアはdataメソッドに格納されています。to_dict()で辞書に変換することができるので、辞書を経由してDataFrameに変換してみます。

data_list = [mol.data.to_dict() for mol in DockingResult]
data_df = pd.DataFrame(data_list)
data_df = data_df[['vina_affinity', 'vina_rmsd_lb', 'vina_rmsd_ub','vina_rmsd_input', 'vina_rmsd_input_min']]
data_df

	vina_affinity	vina_rmsd_lb	vina_rmsd_ub	vina_rmsd_input	vina_rmsd_input_min
0	-6.8	0.000	0.000	43.380367	43.380367
1	-6.8	2.927	5.787	43.486763	43.486763
2	-6.6	11.173	12.829	37.1132	37.1132
3	-6.6	2.950	4.507	43.876957	43.876957
4	-6.5	3.989	6.204	43.624863	43.624863
5	-6.5	3.524	5.691	44.41204	44.41204
6	-6.4	4.251	6.471	45.132126	45.132126
7	-6.4	7.300	8.536	43.149113	43.149113
8	-6.3	8.683	9.959	42.101715	42.101715

ここではvinaの出力に関係する列のみを取り出しています。vina_affinityとなっている列が親和性で最も良いもので-6.8となっています。DataFrameをcsvファイルに書き出しておきます。またドッキング結果のリガンドの情報もSDFに書き出しておきます。

# csv作成
data_df.to_csv("data/AutoDockVina_DockingResult.csv")

# sdf作成
ODDTsdfWriter = oddt.toolkit.Outputfile('sdf', 'data/ODDT_DockingResult_ligand_0.sdf')
for mol in DockingResult:
    ODDTsdfWriter.write(mol)
ODDTsdfWriter.close()

引数のsdfの部分を変えればmol2等別の形式で出力することも可能です。

ODDTの出力におけるエラー

ODDTは非常に扱いやすいのですが、問題点としてAutoDockVinaでドッキング後のリガンドを出力すると原子、結合の帰属がおかしなことになっているということがあります。まず出力前の構造をJupyterNotebookで表示させたのちに、PDBQTで出力したリガンドをPyMolに表示させてみます。

DockingResult[0]

ドッキング結果の出力ファイルにおけるエラー

このように、骨格全体の座標としては保たれているものの、原子のタイプ、結合関係が壊れてしまっています。残念ながら出力のファイルタイプを変えてもこのエラーは修正されませんでした。PyMolで眺めて遊びたかったのに・・・

コマンドラインでAutoDockVinaを使う

ODDTで惜しいところまで行きましたが、私の能力では望みの構造の取得にまで到達できませんでした。残念・・・と思っていたらAutoDockVinaの使用方法のページ にコマンドラインで実行する方法が書いてありました。

コマンドはすごい苦手、というかさっぱりわからないのですが、ここまできたらやるしかない。

用意するものは３つ

1. タンパク質のPDBQTファイル（ODDTで作成済み。flexible = Falseで作成していないと動かない）
2. リガンドのPDBQTファイル（ODDTで作成済み）
3. パラメーターを書いたconfigurationファイル(config.txt)

です。

config.txtの中身は以下のような感じ・・・

receptor = 3w33_A.pdbqt

center_x = 0.63
center_y = 18.01
center_z = 37.64

size_x = 23
size_y = 23
size_z = 23

num_modes = 9

receptorにタンパク質のpdbqtファイル名、center_x,y,zはグリッドの中心座標、size_x,y,zはグリッドのサイズです。num_modesは最大いくつの結合ポーズを出すかを指定します。

ファイルの準備ができたら

1. ターミナルでディレクトリに移動
2. AutoDockVinaのインストール箇所へのPATHを通す。（export PATH=$PATH:~~~~~）
3. コマンド実行！

vina --config conf.txt --ligand ligand_0.pdbqt --out out.pdbqt --log log.txt

ターミナルに随時結果が表示される

と、出て結合ポーズを含むout.pdbqtとスコアなどの情報を含むlog.txtが出来ました！！早い！

PyMolで結果を眺めてみます。折角なので元の共結晶のリガンドを含むファイルと重ねてみます。

エモい！よくわからないけど！

スコアが良いものから悪いものへと順番に流れています。スコアの良いものは画面の右、悪いものは左側という印象でした。詳しい解析はT11パートCで扱われるはず・・・。

ちなみにAutoDockVinaは複数のリガンドをドッキングするためのシェルクリプトも用意してくださっています。こちらのページ からダウンロード出来ます。

chmodで実行権限を設定して実行すればまとめて処理できます。すごい便利ですね。食わず嫌いせずにシェルスクリプト勉強しないと・・・。

まとめ

以上、今回はドッキングを行ってみました。オンラインサービスは自分で環境を準備しなくてもよく、パワーのある計算機を使えるので素敵ですが、いつ止まるかわからないという不安があります。お手元のPCでも試せると便利かもしれません。心置きなく適当なことできますしね！

ODDTを導入したものの結果をうまく出力できなかったり、と相変わらずグダグダになってしまいましたが、T11パートBのオンラインではできなかったことをボンヤリとオフラインで実行できたので良し、ということにしていただければ幸いです。

色々と間違っていることがあると思うのでご指摘いただければ幸いです。

*1: DogSite文献 J. Chem. Inf. Model. (2010), 50, 11, 2041-2052
Bioinformatics (2012), 28, 15, 2074–2075

*2:ProteinsPlus

*3:J. Comput. Chem. 2010(31)455-61

*4:ドッキングシミュレーションのやり方【AutoDock vina】

*5: タンパク質ファイルの読み込みが遅く、セルが [＊] のまま進まないときはJupyter Notebookの拡張機能Variable Inspectorが原因の可能性があります。こちらの参考記事　【Jupyter Notebook】Jupyter notebookの動作が重い およびIssue　と同じエラーかはわかりませんが、とりあえずnb_extensionsから外すと次のセルが実行できました。

*6:RDKitによる3次元構造の生成

TeachOpenCADD トピック１１（パートB）〜オンラインAPI/サービスを使った構造に基づくCADD〜

こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル１１をもととしておりCC BY 4.0のライセンスに従っています。

トークトリアル１１は論文発表時のTeachOpenCADDには含まれておらず新しく追加されたものです。
４つのノートブックに分かれており、ここではその3つ目を訳出します。
パートBで取り上げられているWebサーバー、サービスのうちドッキングに関するSwissDock、OPALはともに2020年5月現在サービスが提供されていないようで、実行できなくなっています。使えるものではなくなってしまっていますが、基本的な考え方や手法を知る上では役にたつかもしれません。

~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。

github.com

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

以下、訳。

トークトリアル 11 （パート B）

オンラインAPI/サーバーを使った構造に基づくCADD

Developed at AG Volkamer, Charité

Dr. Jaime Rodríguez-Guerra, Dominique Sydow

このトークとリアルの目的

これは「オンラインWebサービス」についてのトークトリアルのパートBです。:

• 11a. キナーゼ阻害剤の候補をKLIFとPubChemで検索
• 11b. 11aで取得した候補化合物をターゲットタンパク質に対してドッキング
• 11c. 結果を評価し既知のデータと比較

入力構造を取得したあとで、ドッキングソフトウェアを使って良いタンパク質-リガンドポーズを見つけます。

学習の目標

理論

• ドッキングの基本
• 利用可能なソフトウェア
• 知られている限界

実践

• 構造の準備
• 結合サイトの推定
• ドッキング計算の実行
• 結果の保存

レファレンス

• Structure‐Based Virtual Screeningの章 Ch. 6.8 in "Applied Chemoinformatics: Achievements and Future Opportunities" (2018)
• ドッキングソフトウェアをベンチマークテストする方法 (J. Med. Chem. (2006), 49, 23, 6789–6801)
• ドッキングソフトウェアについてのレビュー (Biophysical Reviews (2017), 9, 2, 91–102)
• DoGSiteScorer：結合サイトを同定するプログラム。
  • J. Chem. Inf. Model. (2010), 50, 11, 2041-2052
  • Bioinformatics (2012), 28, 15, 2074–2075
• SwissDock：ドッキングプログラム
  • Nucleic Acids Res. (2011), 39, W270-7.
  • J Comput Chem. (2011), 32, 10, 2149-59
• OPAL Webサービス：必要に応じてリモートでソフトウェアを実行可能
  • 文献： Nucleic Acids Res. (2010), 38, W724-31
  • ドキュメンテーション Opal: Simple Web Services Wrappers for Scientific Applications
• AutoDock Vina：ドッキングプログラム (J Comput Chem. (2010), 31, 455–461)

---

理論

ドッキング

Pagadala, Syed & Tuszynski, Biophysical Reviews 2017, 9, 2, 91–102　から引用して改変:

ドッキングは標的タンパク質の結合サイト内における低分子化合物の振る舞いを探索する手法です。適用されたソフトウェアは探索アルゴリズムを実行し、最小エネルギーに収束するまでリガンドの配座を再帰的に評価します。最後に、候補となるポーズをランク付けするために、親和性スコア関数ΔG （合計の単位はkcal/mol)を静電エネルギーとファンデルワールスエネルギーの合計として用います。生物学的な系におけるこれらの特定の相互作用のドライビングフォースは、結合サイトの表面とリガンドあるいは基質の形状と静電的な状態の相補性への志向です。

これらのスコア関数は計算が速くなるように合わされているので、他の分子モデリング手法と比較してしばしば精確性に劣ります。通常、これらは次に基づいています。

• 分子力学法原理
• 知識に基づくポテンシャル
• 形状と幾何学的な相補性

探索空間の次元を減らすために、いくつかの近似がよく適用されます。

• タンパク質の構造はほとんどが剛直なものと考えられますが、探索空間に近接するいくつかの側鎖は限られた配座セットの探索が許可されていることもあります（ロータマー）。
• リガンドはバーチャルスクリーニングの研究では剛直なものとして考えられますが、より詳細な計算においては結合の二面角を自由に探索することがよく許可されます。二つのオプションの折衷案は、とりうる配座のセットを予め定義しておくことです。

既存ソフトウェアの例

• 商用
  • GOLD
  • Schrödinger
  • FlexX
• 無料（あるいはアカデミックフリー）
  • AutoDock
  • AutoDock Vina
  • DOCK
  • OpenEye

知られている限界

次の近似は計算においてアーティファクトを取り込む原因となりうります。

• タンパク質のほとんどを剛直したものとして扱うので、タンパク質-リガンド結合の動的、適応的な性質はほとんど探索されません。このため偽陽性につながります。リガンドが結合ポケットで適切なポーズを見つけたとしても、タンパク質が近接するエネルギー安定な配座を探索できないと、このポーズは保証されません。言い換えれば、リガンドが結合ポケットにとどまるかどうか確認するために短時間の分子動力学（MD）トラジェクトリを計算することが常に推奨されています。
• スコア関数は解く計算コストが低くなければなりません。良いポーズと悪いポーズを区別するのに十分な精度がある一方で、最も良い複数のポーズを並べ替えるには問題があることがあります。例えば、最も人気のあるドッキングプログラムは、実験で明らかとされているポーズを計算で見いだすことができる一方で、このポーズが提案されたポーズのセットの中で最も良いものであることは稀です。
• 計算コストを減らすために、ドッキング手順はタンパク質の一部分（通常、結合ポケットとして知られている場所の周辺）でのみ実行されます。CADDパイプラインにおいて、正しい結合ポケットを選択することはもう一つの課題となります。
• 計算の精確性を高めるために、適宜構造を削減する必要があります。アミノ酸とリガンドのプロトン化状態は正しい状態を得るにはこつがいることがあり、特に互変異性体（の可能性がある）の場合難しくなります。これは誤った結果を得る、もう一つの原因となります。

これらの限界にも関わらず、依然としてドッキング計算は全ての医薬品研究室で人気のある手法で、他の種類の分子シミュレーションとともに使われています。このトークトリアルのパートでは次の方法を学びます。

1. パートAで取得したタンパク質とリガンドの準備（ローカルで実行）
2. もっとも可能性が高い結合ポケットの推定（オンラインで実行）
3. ドッキング計算の実行（オンラインで実行）

---

構造の準備

OPALWebサービスにあるVinaを使います。しかし、あらかじめ構造を準備することが必要となります。オススメのアプローチはAutoDockToolsにある準備のためのスクリプトを使うことです。これらのツールはツール自体でディストリビューションされており、Python2としか互換性がありませんが、私たちはタンパク質とリガンドの準備に必要な部分を含むPython 3で利用可能なフォークを準備しました。このトークトリアルで必要なことを行うには十分のはずです。

訳注(2020/05)
AutoDockToolsのライセンスは、非商用の科学研究で使用する場合は無料ですが、商用で使用するには正式な登録手続きが必要とのことです。正確な情報はソフトウェアのライセンスを確認してください。
訳注ここまで

結合ポケット予測

ドッキング計算は、合理的に狭い探索空間で実施すると最もパフォーマンスがよくなり、この空間は通常、一つの結合ポケットを含みます。最も良い結合ポケットを計算に基づいて予測するために、Protein.plusで無料でオンラインで利用できるDoGSiteScorerを使うことができます。

Proteins.plus DoGSiteScorer

• 役割： いくつかのCADDツールのためのインタラクティブなWenインターフェース
• Webサイト： http://proteins.plus
• API： RESTに基づくものがあります。requestsをそのまま適用するだけで十分です。
• ドキュメンテーション： https://proteins.plus/help/index
• 文献：
  • J. Chem. Inf. Model. (2010), 50, 11, 2041-2052
  • Bioinformatics (2012), 28, 15, 2074–2075

自動化されたタンパク質活性部位の予測は、大規模なタンパク質機能予測と分類、そしてドラッガビリティ（druggability、医薬品の標的可能性）を推測する上で必須となります。ここで我々は、画像処理に端を発するガウシアン差分（Difference of Gaussian, DOG）の方法に基づきタンパク質の結合サイトを予測する、新規な構造に基づく手法、DoGSiteを提案します。既存の手法とは対照的に、DogSiteは予測したポケットをサブポケットへと分割し、活性部位のトポロジーの洗練された表現を明らかにします。DoGSiteはPDBBindとscPDBデータセットの92%以上で精確に結合ポケットを予測し、現在利用可能な手法のうち最もパフォーマンスの良いものと並ぶ性能があります。PDBBindデータセットの63%で、見つかったポケットはより小さなサブポケットへと細分化することができました。予測の87%で、共結晶化されたリガンドは正確に一つのサブポケットに含まれていました。さらに、リガンドとポケットの組み合わせの網羅する範囲を考慮に入れることで、より精確な予測性能の評価を導入しました。90%の事例で、DoGSiteは少なくともリガンドの半分を含むポケットを予測しました。70%の事例で、さらに各ポケット自体の1/4以上が共結晶化されたリガンドで覆われていました。サブポケットを考慮することで、適用範囲が広がり、後者の性能評価において成功確率83%となりました。

ドッキング

オンラインで無料で利用できるWebサービスが２、３あります。SwissDockとOPAL Webサービス（AutoDock Vinaを含みます）です。

SwissDock

• 役割： ドッキング計算を実行します
• Websサイト： http://www.swissdock.ch
• API： SOAPに基づくものがあります。正式なクライアントはありません。suds-communityを使ってください。
• ドキュメンテーション： http://www.swissdock.ch/pages/soap_access
• 文献：
  • Nucleic Acids Res. (2011), 39, W270-7.
  • J Comput Chem. (2011), 32, 10, 2149-59

SwissDockはターゲットタンパク質と低分子化合物の間で生じうる分子間相互作用を予測するWebサービスです。 SwissDockはドッキングソフトウェアEADock DSSに基づいていますが、そのアルゴリズムは次のステップからなります。
1. ボックスの中（ローカルドッキング）あるいは全ての標的のくぼみの近傍（ブラインドドッキング）で、結合モードを数多く生成します。 2. 同時に、グリッド上でCHARMMエネルギーが予測されます。 3. 最も好ましいエネルギーの結合モードがFACTSで評価され、クラスター化されます。 4. 最も好ましいクラスターをオンラインで可視化し、自分のコンピュータにダウンロードすることができます。

OPAL Webサービス

• 役割： サービスとしてのCADD （CADD as a service)
• Webサイト： http://nbcr-222.ucsd.edu/opal2/dashboard
• API： SOAPに基づくものがあります。正式なクライアントはありません。suds-communityを使ってください。
• ドキュメンテーション： http://nbcr-222.ucsd.edu/opal2/dashboard?command=docs (現在オフラインとなっています）
• 文献：
  • 論文； Nucleic Acids Res. (2010), 38, W724-31
  • ドキュメンテーション： Opal: Simple Web Services Wrappers for Scientific Applications

医学生物的なアプリケーションはますます複雑となってきており、しばしば多数のプロセッサとメモリを有する大規模な高性能計算資源を必要とします。アプリケーションのデプロイの複雑さと、クラスターコンピューティング、グリッドそしてクラウドコンピューティングにおける進歩は、医学生物学研究をサポートする新しい方式を必要としています。サービスとしての科学ソフトウェア（Scientific Software as a Service、SaaS）は、簡単な標準化されたWebインターフェースを通して、医学生物学的なアプリケーションに、大規模に実行可能で透明性のあるアクセスを可能とします。この目的のため、私たちは2007年8月にMEMEと名付けられたバイオインフォマティクスのアプリケーションをサポートするための運用Webサーバーを構築しました（ http://ws.nbcr.net ）。以来、このサーバーは成長を続け、AutoDockとAutoDock Vinaによるドッキング分析、PDB2PQRとAPBSを使った静電計算、そしてSMAPを使ったオフターゲット分析を含むようになりました。サーバー上のアプリケーションは全てOPALで動いています。OPALは、単純なXML環境設定ファイルを書くことで、科学的なコードに一切変更を加えることなく、科学アプリケーションをWebサービスとして簡単に使えるようにするツールキットです。OPALにより、我々の全てのアプリケーションにWebフォームに基づくアクセスとプログラムによるアクセスのどちらでもアクセスすることができるようになります。OPALツールキットは現在、国立生物医学計算資源（National Biomedical Computation Resource、NBCR)と、系列下の共同研究、サービスプロジェクトによる多数の人気のあるアプリケーションへのSOAPに基づくWebサービスアクセスをサポートしています。さらに、OPALはプログラムによるアクセスが可能なので、VisionやKepler、Nimrod/KそしてVisTrailsといった様々なワークフローツールを通して我々のアプリケーションにアクセスすることができます。2007年8月中旬から2009年の終わりにかけて、239,814個のジョブの実行に成功しました。2008年から2009年の間に、１日あたりの実行に成功したジョブは205から411へと２倍以上に増加しました。OPALで実現されているサービスモデルは様々なアプリケーションにとって有用です。Webサービスインターフェースによる他のアプリケーションとの相互運用を提供し、そしてアプリケーション開発者が科学的ツールとワークフローの開発に集中することができるようにします。Webサーバーは次で利用可能です（ http://ws.nbcr.net ）。

実践

パートAからファイルを取得

まず最初に、パート Aで取得したファイルを指定する変数をいくつか定義します。

• PROTEIN: 標的タンパク質の構造を含む.mol2ファイルへのパス。リガンドやイオンは含まない
• COMPLEX: 標的タンパク質の構造と元々のリガンドを結合サイトに含む.pdb ファイルへのパス
• SMILES_FILE: PubChemで見つかった全ての類似化合物を表すSMILESを含む.txtファイルへのパス
• smiles: SMILES_FILEに含まれるSMILES文字列のリスト

PROTEIN = "data/protein.mol2"
COMPLEX = "data/complex.pdb"
SMILES_FILE = "data/similar_smiles.txt"
with open(SMILES_FILE) as f:
    smiles = [line.strip() for line in f]

SwissDockを利用する

SwissDockはSOAPインターフェースを使っているので、sudsをインストールする必要があります。

注： SwissDockのサーバーは最近機能していません。代わりに下のOPALヘ進んでください！ !pip install suds-community

from suds.client import Client
import zlib
import string
import requests

def swissdock_client():
    # 現在サーバーが落ちているようです。。。
    # http://swissdock.vital-it.ch/soap/ は503利用不能を返します。
    # 誤ったドメインを指定しているからでしょうか。。。修正しましょうか？
    SWISSDOCK_WSDL_URL = "http://www.swissdock.ch/soap/wsdl"
    r = requests.get("http://www.swissdock.ch/soap/wsdl")
    r.raise_for_status()
    WSDL = r.text.replace("http://swissdock.vital-it.ch/soap/", "http://www.swissdock.ch/soap/")
    with open("data/swissdock.wsdl", "w") as f:
        f.write(WSDL)
    HERE = _dh[0]
    return Client(f"file://{HERE}/data/swissdock.wsdl")


def prepare_protein(client, protein):
    """
    与えられたPDBファイル（中身は文字列）に対してPSFとCRDを返す。
    """
    encoded_protein = zlib.compress(protein.encode('utf-8'))
    job_id = client.service.prepareTarget(target=encoded_protein)
    while True:
        result = client.service.isTargetPrepared(jobID=job_id)
        if result is None:
            raise ValueError("No such a job present")
        if result in (False, "false", 0):
            time.sleep(5)
        else:  # 準備できました！
            break
    protein_files = client.service.getPreparedTarget(job_id)
    if protein_files is None or len(protein_files) != 2:
        raise ValueError("Could not prepare protein!")
    return protein_files
            

def prepare_ligand(client, ligand):
    """
    与えられたMOL2ファイル（中身は文字列）に対してPDB、RTFとPARを返す。
    
    リガンドはあらかじめプロトン化されていなければなりません!
    """
    encoded_ligand = zlib.compress(ligand.encode('utf-8'))
    job_id = client.service.prepareLigand(ligand=encoded_ligand)
    while True:
        result = client.service.isLigandPrepared(jobID=job_id)
        if result is None:
            raise ValueError("No such a job present")
        if result in (False, "false", 0):
            time.sleep(5)
        else:  # 準備できました！
            break
    ligand_files = client.service.getPreparedLigand(job_id)
    if ligand_files is None or len(ligand_files) != 3:
        raise ValueError("Could not prepare ligand!")
    return ligand_files

def dock(client, protein, ligand, name=None):
    protein_psf, protein_crd = prepare_protein(client, protein)
    ligand_pdb, ligand_rtf, ligand_par = prepare_ligand(client, ligand)
    
    if name is None:
        name = "teachopencadd" + ''.join([random.choice(string.ascii_letters) for _ in range(5)])
    job_id = client.service.startDocking(
        protein_psf, protein_crd,
        ligand_pdb,
        [ligand_rtf],
        [ligand_par],
        name)
    if job_id in (None, "None"):
        raise ValueError("Docking job could not be submitted")
    while not client.service.isDockingTerminated(job_id):
        time.sleep(5)
    all_files = client.service.getPredictedDockingAllFiles(job_id)
    with open('docking_results.zip', 'w') as f:
        f.write(all_files)
    target, docked = client.service.getPredictedDocking(job_id)
    client.service.forget(job_id)
    return target, docked

def smiles_to_pdb(s, out='output.pdb'):
    m = Chem.AddHs(Chem.MolFromSmiles(s))
    AllChem.EmbedMolecule(m)
    if out is None:
        return Chem.MolToPDBBlock(m)
    Chem.MolToPDBFile(m, out)

try:
    import Mol2Writer
except ImportError:
    # RDKitのMol2 writer/readersを手に入れるための醜いハック
    import os
    working_dir = os.getcwd()
    os.chdir(_dh[0])
    !wget https://raw.githubusercontent.com/rdkit/rdkit/60081d31f45fa8d5e8cef527589264c57dce7c65/rdkit/Chem/Mol2Writer.py > /dev/null
    os.chdir(working_dir)
    import Mol2Writer

def step_03_swissdock(protein_pdb, ligand_smiles):
    ligand = Chem.AddHs(Chem.MolFromSmiles(ligand_smiles))
    AllChem.EmbedMolecule(ligand)
    ligand_mol2 = Mol2Writer.MolToCommonMol2Block(ligand)
    client = swissdock_client()
    return dock(client, protein_pdb, ligand_mol2)

現在SwissDockサーバーが利用できないので、下のセルは実行しないようにしています。実行したいのであれば、まずタンパク質のmol2ファイルをPDBに変換する必要があります。OpenBabelを使って行うことができます（パートCでインストールします）。ご自由に都合の良いmol_to_pdb関数を定義してください。

# TODO: mol2からpdbに変換すること！
# protein_pdbblock = mol2_to_pdb(PROTEIN, out=None)
step_03_swissdock(protein_pdbblock, smiles[0])
### OPAL Webサービスを使ったドッキングの実行

SwissDockは現在稼働していません。そこで手を変え別のWebサービスを利用しましょう。インターフェースは少々より原始的ですが、稼働するはずです。ですが、タンパク質とリガンドはローカル環境でAutoDockToolsを使って準備しなければなりません。Python3で利用可能なフォークを準備しましたが、テストは十分にできていません。ですが、ここでの目的には十分機能するように見えます。

次を実行してインストールすることができます：

!pip install https://github.com/jaimergp/autodocktools-prepare-py3k/archive/master.zip

ドッキング計算の手順は次のようになっています。

1. AutoDockToolsでタンパク質とリガンドを準備します。（ローカル環境で実施）
2. Proteins.plus' DoGSiteScorerを使って結合ポケットの可能性があるもののうち最も良いものを見つけます。この情報をVinaによる計算の設定を行うために使います。（幾何学的な中心と探索空間のサイズ）
3. OPAL上でVinaの計算を実行します。

import time
import os
from io import StringIO

構造の準備

構造の準備はAutoDockTooksライブラリの適切な部分を実行するだけでできます。

• MolKitで構造のファイルを読み込む
• 準備のための正しい機能を適用する。タンパク質にはAD4ReceptorPreparation を、リガンドにはAD4LigandPreparationを適用します。

準備自体は次のような項目に注意して行います。

• タンパク質とリガンドに水素原子を付加する。
• タンパク質に属さないおかしな残基を取り除く。
• 原子タイプと部分電荷を割り当てる。
• リガンドの回転可能（torsionable）な分岐をみつける。

結果はPDBQTファイルとしてディスクに書き込まれます。

######################
#
# 構造の準備
#
######################

import MolKit
from AutoDockTools.MoleculePreparation import AD4ReceptorPreparation, AD4LigandPreparation

def opal_prepare_protein(protein):
    """
    AutoDockを使うにはPDBQTファイルが必要です。
    """
    mol = MolKit.Read(protein)[0]
    mol.buildBondsByDistance()
    RPO = AD4ReceptorPreparation(mol, outputfilename=protein+'.pdbqt')
    return protein + '.pdbqt'

def opal_prepare_ligand(ligand):
    """
    AutoDockを使うにはPDBQTファイルが必要です。
    """
    mol = MolKit.Read(ligand)[0]
    mol.buildBondsByDistance()
    RPO = AD4LigandPreparation(mol, outputfilename=ligand+'.pdbqt')
    return ligand + '.pdbqt'

結合サイトの予測

Proteins.plusのDoGSiteScorerは、PDBデータベースにあるタンパク質の処理が必要なだけならば、とても簡単に利用できるREST APIを提供しています（dogsite_scorer_submit_with_pdbidを参照してください）。ですが、我々はKLIFSから取得した構造を使っているので、公式のPDBに登録されている構造の位置と向きが、KLIFSのものと一致しているか保証することができません。両者を重ね合わせて得られた変換行列を、取得したポケットに適用することもできますが、標準Webインターフェースで提供されている方法と同様に、単純にPDBファイルをProteins.plusにアップロードするだけの方がより簡単でしょう。

しかし、REST APIはそのようなオプションを提供していないので、我々はこの処理をリバースエンジニアリングする必要があります。適切なHTTPリクエストの場所をみつけるためには、Chromデベロッパーツールのネットワークタブを開き、通常のWebサイトを利用した際にどのHTTPリクエストが実行されているか書き留めておく必要があります。ここで有効なHTTPリクエストを手に入れるためのポイントは、信頼性トークン(Authenticity token)とHTTPヘッダーです。

この調査の結果は統合して整理しdogsite_scorer_submit_with_custom_pdb 関数に反映させています。技術的な詳細について興味があれば、関数のコメントを読んでください。

このアプローチを行うには、結合サイトを予測するためにHTMLコードを構文解析するため、BeautifulSoupも必要となります。

!pip install BeautifulSoup4

この方法を使えば、最も可能性の高い結合ポケットの幾何学的な中心とそのサイズを取得することができます。Vinaの計算を設定するにはどちらの値も必要となります。

######################
#
# 結合ポケット予測
#
######################

from bs4 import BeautifulSoup
import requests

def dogsite_scorer_submit_with_pdbid(pdbid, chain):
    """
    これは公式APIですが、PDBコードしか使えません。
    
    パラメーター
    ----------------
    pdbid : str
        ４文字のPDB　ID
    chain : str
        解析するタンパク質の鎖（chain）
        
    戻り値
    ----------------
    str
        ジョブが実行中か完了したかについて更新情報を受け取るために用いるURL
    """
    # Proteins.plusにジョブを投げる
    r = requests.post("https://proteins.plus/api/dogsite_rest",
        json={
            "dogsite": {
                "pdbCode": pdbid,
                "analysisDetail": "1",
                "bindingSitePredictionGranularity": "1",
                "ligand": "",
                "chain": chain
            }
        },
        headers= {'Content-type': 'application/json', 'Accept': 'application/json'}
    )

    r.raise_for_status()
    # 計算の更新情報を取得するためには"location"を問い合わせる必要があります。
    return r.json()['location']

def dogsite_scorer_submit_with_custom_pdb(pdbfile):
    """
    カスタムしたPDBをアップロードするために、実際のHTMLフロントエンドを模倣する必要があります。
    
    1. HTMLメタヘッダーからCSRF（クロスサイトリクエストフォージェリ）トークンを取得します。
    2. アップロードするためにファイルをPOSTします。
    3. 返ってくるHTMLページはURL IDを含んでおり、次に内部の省略されたジョブIDを取得するために利用できます。
    　　　　（Webインターフェースを使っている時と同様に）Webサーバーがフロントエンドで実行する非同期の呼び出しを真似することで、
        返ってきたHTMLページのURL IDを使ってパブリックなジョブのAPI IDを取得することができます。
    4. パブリックジョブIDを手に入れたらREST APIの利用に切り替えることができます。
    """
    # プロセスの間、途中のcookieを保持するために`session`を使う必要があります。
    session = requests.Session()
    r = session.get("https://proteins.plus/")
    r.raise_for_status()
    # ホームページには我々のリクエストを有効にするのに必要なCSRFトークンが含まれています。
    #　そうでないと安全ではありません！リクエストを行なっている間、これを使わなければいけないので、
    # session HTTPヘッダーの一部として設定しておくのが最善の方法です。
    html = BeautifulSoup(r.text)
    token = html.find('input', {'name': 'authenticity_token'}).attrs['value']
    session.headers['X-CSRF-Token'] = token

    # 1. ファイルをアップロード
    with open(pdbfile, 'rb') as f:
        r = session.post("https://proteins.plus", files={'pdb_file[pathvar]': f})
    r.raise_for_status()
    
    # REST APIがファイルのアップロードをサポートしているのであれば、すでにパブリックIDを手に入れていることになりますが、
    # それまでは応急処置としてこの方法で済ませておく必要があります。
        
    # 2. 内部のlocation IDの取得
    html = BeautifulSoup(r.text)
    pdb_id = html.find('input', {'name': "dogsite[pdbCode]"}).attrs['value']

    # 3. 内部のジョブIDの取得
    session.headers['Referer'] = "https://proteins.plus" + pdb_id
    r = session.post(f"https://proteins.plus/{pdb_id}/dogsites",
            json={"dogsite": {"pdbCode": pdb_id}},
            headers= {'Content-type': 'application/json', 
                      'Accept': 'application/json'}
    )
    r.raise_for_status()
    job_id = r.json()['job_id']
    time.sleep(3)  # サーバーがリクエストを処理できるように、処理を続ける前に少し待ちます。
    
    # 4. パブリックなジョブIDの取得
    while True:
        r = session.get(f"https://proteins.plus/{pdb_id}/dogsites/{job_id}?_={round(time.time())}",
                        headers= {
                            'Accept': 'application/json, text/javascript, */*',
                            'Sec-Fetch-Mode': 'cors',
                            'Sec-Fetch-Site': 'same-origin',
                            # どうもこの下の１行が大事なようです。
                            # これ無しではエラー406が投げられます。
                            'X-Requested-With': 'XMLHttpRequest'}
                       )
        r.raise_for_status()
        if 'Calculation in progress...' in r.text:  # まだ完了していません。
            time.sleep(5)
            continue
        if 'Error during DogSiteScorer calculation' in r.text:  # 形式のおかしなファイル？
            raise ValueError('Could not run DoGSiteScorer!')
        break
    
    results_id = None
    for lines in r.text.splitlines():
        for line in lines.split('\\n'):
            if 'results/dogsite' in line:
                results_id = line.split('/')[3]
                break
    if results_id is None:
        raise ValueError(r.text)
        
    return f"https://proteins.plus/api/dogsite_rest/{results_id}"
    

def dogsite_scorer_guess_binding_site(protein):
    """
    タンパク質の最も可能性が高い結合サイトを取得するためにProteins.plusのDoGSiteScorerを使います。
    """
    if len(protein) == 4:  # PDBコード
        job_location = dogsite_scorer_submit_with_pdbid(protein)
    elif protein.endswith('.pdb'):
        job_location = dogsite_scorer_submit_with_custom_pdb(protein)
    else:
        raise ValueError("`protein` must be a PDB ID or a path to a .pdb file!")
    
    # 計算が完了したか確認します。
    while True:
        result = requests.get(job_location)
        result.raise_for_status()  # 失敗した場合、ここで止まります。
        if result.status_code == 202:  # まだ実行中です。
            time.sleep(5)
            continue
        break
    
    # residuesファイルは幾何学的中心と半径をPDBファイルのコメントとして保持しています。
    # 最初のファイル（residues[0]）が最もスコアの良いポケットです。
    pdb_residues = requests.get(result.json()['residues'][0]).text
    for line in pdb_residues.splitlines():
        line = line.strip()
        if line.startswith('HEADER') and 'Geometric pocket center at' in line:
            fields = line.split()
            center = [float(x) for x in fields[5:8]]
            radius = float(fields[-1])
            break
    return center, radius  # これがVinaの計算に必要となるものです。

OPAL上でVinaを実行

（１）タンパク質とリガンドを準備し、（２）探索空間の推定が完了したら、OPAL Webサーバーに実際の計算をサブミットできます。以下の処理を行います。

1. sudsを使ってSOAPクライアントを初期化する
2. SOAP XMLリクエストに加えられるように、ファイルをbase64文字列 としてエンコードする
3. ジョブリクエストをサブミットし、ジョブIDを受け取る
4. 計算が完了したかどうか確認するため、サーバーにジョブIDを使って問い合わせる
5. requestsで関連するファイルをダウンロードする

上記の手順はOPAL Webサイトにはあまり記述されておらず（実際、ドキュメンテーションは利用できなくなってきます)）、これらのサーバーに依存しているUCSF Chimeraモジュールのソースコードの中身を調査することでわかりました。

通常、計算には５−１５分かかるので、４番目のステップのためにJupyterの便利な機能を使ってリアルタイムに出力ファイルの中身を更新します（iprint()関数を参照してください）。こうすることで、計算がうまく進んでいるかどうか盲目的に信じる代わり進行度合いを確認できます！

######################
#
# 計算の実行
#
######################

from suds.client import Client
from IPython.display import display, clear_output, HTML
from rdkit import Chem
from rdkit.Chem import AllChem

VINA_CONFIG = """
center_x = {center[0]}
center_y = {center[1]}
center_z = {center[2]}
size_x = {size[0]}
size_y = {size[1]}
size_z = {size[2]}
"""

def opal_run_docking(protein, ligand, center, size, stream_output=True):
    """
    OPAL Webサービスに接続しジョブをサブミットする
    """
    client = Client("http://nbcr-222.ucsd.edu/opal2/services/vina_1.1.2?wsdl")
    files = 'receptor.pdbqt', 'ligand.pdbqt', 'vina.conf'
    with open(protein) as f:
        protein_contents = f.read()
    with open(ligand) as f:
        ligand_contents = f.read()
    file_map = [
        {'name': 'receptor.pdbqt',
         'contents': base64ify(protein_contents)},
        {'name': 'ligand.pdbqt',
         'contents': base64ify(ligand_contents)},
        {'name': 'vina.conf',
         'contents': base64ify(VINA_CONFIG.format(center=center, size=size))},
        {'name': 'results.pdbqt',
         'contents': ''},
    ]
    cli_args = "--receptor receptor.pdbqt --ligand ligand.pdbqt --config vina.conf --out results.pdbqt"
    
    response = client.service.launchJob(cli_args, inputFile=file_map)
    job_id = response.jobID
    url = f"http://nbcr-222.ucsd.edu/opal-jobs/{job_id}"
    message = "Waiting for job " + url
    while True:
        r = requests.get(url + "/vina.out")
        try:
            r.raise_for_status()
        except:  # 最初のチェック段階では出力ファイルはまだ出来ていないかもしれません。
            iprint(message)
        else:
            iprint(f"{message}\n{r.text}")
        if client.service.queryStatus(job_id).code == 2:
            time.sleep(10)
            continue
        print('\nFinished!')
        break
        
    output_response = client.service.getOutputs(job_id)
    output_files = {
        'stdout.txt': requests.get(output_response.stdOut).text,
        'stderr.txt': requests.get(output_response.stdErr).text,
    }
    for f in output_response.outputFile:
        if f.name in files:
            continue
        r = requests.get(f.url)
        r.encoding = 'utf-8'
        r.raise_for_status()
        contents = r.text
        output_files[f.name] = contents 
        time.sleep(0.1)
    
    return output_files

######################
#
# ユーティリティー
#
######################

import base64

def base64ify(bytes_or_str):
    """
    Py2k base64encodeの動きを真似する
    """
    if isinstance(bytes_or_str, str):
        input_bytes = bytes_or_str.encode('utf8')
    else:
        input_bytes = bytes_or_str

    output_bytes = base64.urlsafe_b64encode(input_bytes)
    return output_bytes.decode('ascii')

def iprint(s):
    """
    この関数を使って出力をプリントし、その前のものを上書きすることで、
    継続的に更新されているようにみえるようにすることができます:)
    """
    clear_output(wait=True)
    s = s.replace("\n", "<br />")
    display(HTML(f'<pre>{s}</pre>'))

全部まとめる

必要な関数をすべて定義し終わったので、パイプラインを作成することができます。

def step_03_opal(protein, smiles, pdbcomplex):
    """
    タンパク質の構造とSMILES文字列のリストを与え、以下の手順を実施します・
    ステップ：
        1. AutoDock Vinaのためのタンパク質を準備（ローカルで実行） 
        2. DoGSiteScorerを使って最も可能性の高い結合サイトを見つける
        3. 各リガンドについてRDKitを使って3DのPDBファイルを書き出し、
            OPALでAutoDockを実行する
    
    すべて行うのに大体5-15分かかります。
    
    結果は出力ファイルを含む辞書です。我々が主に興味があるのは結果の result['results.pdbqt']です。
    """
    prepared_protein = opal_prepare_protein(protein)
    center, radius = dogsite_scorer_guess_binding_site(pdbcomplex)
    size = [radius] * 3  # Vinaは立方体のボックス以外に対応していますが、単純のため立方体を使います
    for i, smile in enumerate(smiles):
        smiles_to_pdb(smile, f'data/ligand{i}.pdb')
        prepared_ligand = opal_prepare_ligand(f'data/ligand{i}.pdb')
        result = opal_run_docking(prepared_protein, prepared_ligand, center, size)
    return result

実行します！マジックコマンド%timeでどれくらい時間がかかったか測ります。

# `smiles`の最初のリガンドだけを処理します。
%time result = step_03_opal(PROTEIN, smiles[:1], COMPLEX)

訳注(2020/05)
上記のセルはエラーを返しました。OPALがサービスを停止しているのか現在実行できなくなっています。 NBCRと関係するWebサイト、Webサービスが2020年4月30日をもって終了したようです（NBCRリンク ）。
訳注ここまで

出力結果を理解する

resultは辞書で、出力ファイルとその中身のテキストに相当するいくつかのkeyをもちます。我々が主に興味があるのは次です。

• ドッキングされたリガンドを含むresults.pdbqt。複数のモデルをふくむ修正されたPDBファイルです。タンパク質のグリッドを保持しているので、もともとのタンパク質の構造と一緒に各リガンドのモデルを開くだけで大丈夫です。
• vina.outは上記で見たテキストの出力です。テーブルのような情報を見ることができます。

あとで取り出せるようにディスクに保存しておきます。

with open('data/results.pdbqt', 'w') as f:
    f.write(result['results.pdbqt'])
with open('data/vina.out', 'w') as f:
    f.write(result['vina.out'])

ドッキング結果の可視化

計算を実行しファイルをダウンロードしたら、可視化しましょう！やり方はパートCを見てください。

ディスカッション

OPALはドッキング計算を行うVinaを無料で提供していますが、サブミットする前に入力ファイルをローカル環境で準備する必要があります。準備の中には、探索空間の定義も含まれており、通常既知の結合サイトの周囲を指定します。視覚的に推定する代わりに、Proteins.plusのDoGSiteScorerサーバーを使って、最も可能性が高い結合サイトの中心と半径を計算しました。これら２つのサービスは異なるコミュニケーションインターフェースを使っていました。

Protein.PlusはRESTを使っていますが、KLIFSとは異なりswagger.jsonの定義を提供していないので、自分のリクエストを手動で構築しなければなりませんでした。幸い、Webサービスが単純だったので２、３個のリクエストが必要なだけ、、、と思っていました！！現在のAPIはカスタマイズしたPDBのアップロードを許可していないので、通常のブラウザーを使っているふりをするために多少のGET、POSTリクエストを使う必要がありました。正しいリクエストを推測するうえで、Chromデベロッパーツールのネットワークタブはとても便利でした。Webサービスのリバースエンジニアリングが必要になったら、便利なツールの一つです！

OPALはXMLベースの標準SOAPで構築されており、JSONベースのRESTよりもちょと扱いにくですが、sudsを使うことでとても簡単になりました。sudsは利用できるメソッドを教えてくれますが、これらはあまりきちんと説明されていません。幸い、UCSF Chimeraのモジュールにコードの例がいくつか散らばっていました。ジョブIDを手に入れたらOPALにサブミットされたジョブはパブリックなものとなり、ファイルはリアルタイムに更新されるので、サーバーにN秒毎に問い合わせを行い、動的にdisplay()関数で表示されるHTMLを更新することで、計算の出力をライブでプレビューすることができました。この技は自分のマシンで行なっている計算を問い合わせるのにも使うことができます。ファイルの中身を読み、現在の出力を消して新しい出力をHTML()オブジェクトとして表示するだけです（ipirnt()関数を参照してください）。

クイズ

• リモートサーバーでドッキングが成功したかどうかどうすればわかりますか？
• なぜAutoDock Vinaの入力ファイルをローカル環境で準備する必要がありましたか？
• 応用：MCule'sドッキングサーバーをJupyter Notebookから使うための正しいHTTPリクエストを見つけることはできますか？

訳以上

追記

今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。

残念ながらノートブックの実行はできませんでした。オンラインサービスを使うとなるとサービスが停止したり安定しないといった問題もあるんですね。
パートBの内容はアプリケーションをプログラム上で扱うためのコードばかりで私は全くついていけませんでした。文章から理解する能力がないので絵で説明してほしいです・・・。