magattacaのブログ

日付以外誤報

花咲爺の隣の家で 〜植物ホルモンとサリドマイドと分子糊〜

皆さん高校の理科は何を選択されましたか? バネの運動で早々と物理を諦めた私は生物一択でした。

で、(おそらく)全ての生物選択の高校生が直面する問題

「茎どっちに伸びるのか問題」

f:id:magattaca:20200606225129p:plain
Fig.1 チンアナゴではない

そんな乱暴なことしていいんですか!って感じの実験ですよね。

この屈性の背後にあるのが植物ホルモン オーキシンですが、その作用メカニズムにはサリドマイドとの類似性が指摘されているそうです。*1

1. 植物ホルモンの分子実体?

さて、植物ホルモン。「高校生物の教科書に載っているくらいなんだから研究なんてし尽くされているんじゃないの?」と思っていましたが、 意外にも(?)その分子的な実体、メカニズムについては最近になってやっと明らかとなったことが沢山あるようです。

例えば花成ホルモン フロリゲン。こちらは植物の花芽分化を決定づける因子、つまり花を咲かせるホルモンです。

花咲かジジイホルモン(by 某高校教師)」なんて科学的にも産業的にもとても魅力的な研究対象ですが、 1936年にその存在が予言されて以来、ずっと正体は謎のままで「幻の植物ホルモン」とも呼ばれていたそうです。

明らかになったのは2000年代に入ってから。。。

まず、

  • 2007年 フロリゲンの分子実体はタンパク質Hd3a/FTであることが*2*3

ついで、

  • 2011年 フロリゲンの受容体が14-3-3タンパク質であることが*4

明らかとされました。

f:id:magattaca:20200606230714p:plain
Fig. 2 フロリゲン受容体

フロリゲンをめぐる一連の研究には日本の研究の貢献も大きいらしく、 以下の日本語レビューでは一筋縄には行かなかった発見の経緯が書かれていて非常に面白かったです。*5

leading.lifesciencedb.jp

2. オーキシン

さて、フロリゲンの同定が困難であった一つの理由は、他の植物ホルモンが低分子有機化合物であるのに対して、フロリゲンがタンパク質であったことだそうです。

例えば、冒頭取り上げた屈性に関する植物ホルモン、オーキシンの分子実体はインドール-3-酢酸(IAA)です。

f:id:magattaca:20200606230827p:plain
Fig.3 オーキシンの分子構造

分子量の差がすごいですね。こちらは1930年代には単離同定されていたようです。

人尿から単離」というパワーワードは脇において、、、インドール-3-酢酸をみればすぐにトリプトファンとの類似性に気づくと思います。 酵素反応ならちょちょいと簡単に作ってしまいそうですが、意外にもこの生合成経路、明らかになったのは2011年だそうです!*6 *7 *8

f:id:magattaca:20200606231100p:plain
Fig. 4 オーキシンの生合成主要経路

分子サイズが小さければ簡単というわけにはいかないんですね。

3. オーキシンの作用メカニズム

発見から100年たっても活発な研究がつづくオーキシン。ではその作用メカニズムはどのようになっているのでしょうか?

外部刺激の応答に関わる植物ホルモンと聞くと、シグナル伝達物質、例えば

  • 細胞膜に受容体があって下流にシグナルが流れる??
  • ステロイドみたいに核内受容体と結合して転写調節???

といったあたりが思い浮かぶと思います。

・・・ちゃうねんて!!! タンパク質分解誘導するんやって!!! びっくりやわ!!!

f:id:magattaca:20200606231729p:plain
Fig. 5 オーキシンの作用メカニズムはタンパク質分解誘導

2000年代になり、まず

ついで

  • 2007年、オーキシンと受容体、相互作用するペプチドとの複合体の共結晶構造 *11

により、その認識機構が明らかとなりました。

カニズムの詳細は流れは以下の通りです。

まず、オーキシンの作用の最終的なアウトプットはオーキシン応答遺伝子群の転写促進です。(Fig. 5)

この遺伝子群はARF(Auxin Responce Factor)タンパク質群により転写調節されています。一方、ARFタンパク質群はAux/IAAタンパク質が結合することでその機能が抑制されています。 オーキシンはこのリプレッサー Aux/IAAタンパク質の分解を誘導することで、ARFタンパク質群の「抑制を抑制」します。その結果、オーキシン応答遺伝子群の転写が活性化されます。

・・・抑制の抑制 マイナス x マイナス = プラス ってやつですね。

次に、オーキシンによるタンパク質分解誘導のメカニズムですが、ユビキチンープロテアソーム系を利用しています。

f:id:magattaca:20200606232613p:plain
Fig. 6 オーキシンは分子糊としてユビキチン-プロテアソーム系を利用する

共結晶構造でオーキシンの結合しているTIR1はSCFユビキチンーリガーゼ複合体のうち、基質の認識に関わるタンパク質です。 TIR1がAux/IAAタンパク質を認識し、ユビキチン化が進行、ユビキチンラベルを認識したプロテアソームにより分解が行われます。 結晶構造で明らかとなったように、オーキシンはTIR1のくぼみに結合し、認識表面を調節することで、 Aux/IAAタンパク質をTIR1が認識できるようにします。 つまり、オーキシンはタンパク質間の結合を促進する分子糊(Molecular Glue)として働いていることになります。

f:id:magattaca:20200607000252p:plain
Fig. 7 Molecular Glue オーキシンの結合様式

すごい仕組みですね。ずっとメカニズムが不明というのも納得です。結晶構造解析の威力すごいです。*12

さて、植物に関して高校教科書にのっているようなトピックでも、「詳細・分子レベルのメカニズムは最近になってようやくわかったことも多く、まだまだ面白い話がありそうだ」ということが分かりました。

医薬品に関しても同様に、教科書にも載っているし臨床で使われているけどメカニズムはよくわかっていない、実用化されてからメカニズムの研究が進んだという例は色々とあります。

というわけで次のトピック、サリドマイドに移ります。

4. サリドマイドのメカニズム

サリドマイドの臨床上の問題に関しては非専門家の私が書ける内容ではないので、ここでは「催奇形性の問題があるが現在でも多発性骨髄腫の治療に用いられている」ということに留めておきます。

サリドマイド睡眠薬としての利用は1960年前後にまで遡るのですが、そもそもどのように薬効、副作用が生じているのか?そのメカニズムは2010年頃から漸く分かってきたようです。

まず、2010年 サリドマイドの結合するタンパク質がcerebron(CRBN)であることが東京工業大学 半田宏教授らのグループにより報告されました。*13

当時、色々なメディアで取り上げられたので耳にされた方も多いかもしれません。

独自に開発したアフィ二ティクロマトグラフィー用担体FGビーズ(通称、半田ビーズ)を利用し、サリドマイドを固定化したビーズに結合するタンパク質を単離する、という手法で同定しています。  

独自のツール開発とそれを用いた長年の謎を解明する糸口をみつける研究って凄い格好良いですよね!

さて、こうして見つかったタンパク質CRBNですが、分子の機能としては他のタンパク質(CCB1、Cul4A、Roc1)と共にユビキチンリガーゼ複合体を形成します。

2014年には複合体の結晶構造が解かれています。*14

f:id:magattaca:20200606233536p:plain
Fig. 8 サリドマイド結合タンパク質の複合体結晶構造

CRBN表面のくぼみにサリドマイドが結合している様です。オーキシンの時と同じ感じになってきましたね。先の図を使いまわすと・・・

f:id:magattaca:20200606233732p:plain
Fig. 9 サリドマイドは分子糊としてユビキチン-プロテアソーム系を利用する

結晶構造の報告された2014年 Natureの文献中では、CRBNが認識、分解誘導する元々のタンパク質(endogeneous substrate)はMEIS2で、 サリドマイド(thalidomideとその類縁体)存在下では、標的タンパク質がIKZF1 / IKZF3に変化するという提案がなされています。

つまり、サリドマイドは、植物におけるオーキシンのように、

「ユビキチンーリガーゼ複合体に結合して、特定の基質タンパク質との結合を促進する接着剤(分子糊Molecular Glue)として働くことで、その分解を誘導する」

ということになるようです。

5. サリドマイドを出発点としたタンパク質分解誘導薬

さて、ユビキチン-リガーゼ複合体に結合することがわかったサリドマイド。結晶構造では、タンパク質とぶつからずに分子の構造をさらに拡張できそうな方向も示唆されました。

で、その後展開としてはどうなるか?というと

  • ① 置換基をぶら下げて選択性をあげる or 分解されるタンパク質基質を変える(Cerebron Modulator)
  • ② リンカーを介して別のタンパク質結合部位を追加することで基質を変える(PROTACDegronimid

といったタンパク質分解誘導薬の研究に派生していったようです。

以下のブログ記事が非常に分かりやすくてオススメです。というかこちらを読めば私のこの記事を読む必要はありません。*15

luckprepopp.com

さておき、一応有機化学をかじった身ととしてはどんな構造に展開されて行ったのかがやっぱり気になります。②のPROTACは大きくてしんどいので、①の方を見てみたいと思います。

6. Cerebron Modulator

上でご紹介した記事に書かれていましたが、Molecular Glueとしてのサリドマイド誘導体の研究はCelgene社*16が活発にすすめているそうです。

承認済みのサリドマイド誘導体と一緒にCelgene社による化合物を描いてみました。

f:id:magattaca:20200606234854p:plain
Fig. 10 サリドマイド誘導体の構造

サリドマイドとCRBNの複合体結晶構造では、サリドマイドの2つの環構造のうち芳香側が外側に向かって出ていました。 新しい化合物はどれも芳香環側の修飾となっているようですね。

とても興味深いことにCC-885については分解するタンパク質基質の変化が報告されています。GPST1というタンパク質の分解が誘導され、急性骨髄性白血病の細胞に対して強い増殖抑制作用を示したとのことです。*17

複合体の結晶構造をみてみましょう。

f:id:magattaca:20200606235207p:plain
Fig. 11 分子糊CC-855の複合体構造

CC-855がCRBNとGSPT1の両方に接触し、まさに分子糊、Molecular Glueとして働いているようです。すごい。

長年の謎であった標的分子の解明とそれに基づいた新しい科学の展開。ワクワクします。

サリドマイドのメカニズムやCelebron Modulatorについては東京医科大学 伊藤先生、半田先生が日本語で解説してくださっていてとても面白いのでオススメです。*18 *19

www.jstage.jst.go.jp

first.lifesciencedb.jp

7. オーキシンの応用、オーキシンデグロン法

さて、ここまでタンパク質分解を誘導する分子糊という「植物ホルモン オーキシンと医薬品 サリドマイドのメカニズムの類似性」をみてきました。

また、研究の発展として「サリドマイドのタンパク質分解誘導薬としての展開」をみました。

では、話を戻してサリドマイドよりも先に分子メカニズムが明らかとされたオーキシンではどのような応用がなされているのでしょうか?

国立遺伝学研究所 鐘巻研究室では、 元々は植物における作用であるオーキシンの分解誘導を酵母や哺乳動物細胞でもつかえるようにした技術、 オーキシンデグロン(Auxin-Inducible Degron: AID)法を確立されたそうです。

研究室HPの説明動画(こちらのページの下の方)が 非常に分かりやすいので興味を持たれた方はご覧になると良いと思います。・・・私では到底説明できない。

ざっくり雰囲気でいうと、「分解したいタンパク質にラベルをつけておくと、オーキシンを添加することによって分解が始まるようにできる」そうです。

つまり、「コンディショナルノックアウトのコンディションがオーキシン」ってこと、、、でしょうか???

上の説明動画中では蛍光タンパク質での例が紹介されていましたが、みるみる細胞の光が消えていく様子はなかなかに衝撃的です。

  • レスポンススピードが速い(半減期30分!)

  • オフターゲットの可能性が低い

ということがツールとしての特徴とのことです。

なんと言うかそもそも植物のメカニズムが動物細胞で使えることが衝撃でした。 それだけユビキチン-プロテアソーム系の保存度が高い、重要な機能ということなのでしょうか?・・・よく分かりません。

8. 花咲爺の隣の家で

最後に、冒頭、フロリゲンにまで話をもどしましょう。

幻の植物ホルモンフロリゲンはタンパク質で花の形成に関わるものでした。
私が高校生だったころにはまだ同定されておらず、生物教師は「これを見つければ花咲か爺さんになれる。大金持ちだ!」といっていました。

ですが、2020年の性格の悪い私はこう考えてしまいます。

「フロリゲン タンパク質と結合する化合物を見つけよう!そいつとオーキシンをいい感じにつなげば、フロリゲン分解PROTACのできあがり!すべての桜を散らせてやるぜ!」

・・・おしまい!!

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

*1:こちらのレビューを読んだのが切っ掛けです。フリーで読めてとても面白い。 Stanton BZ, Chory EJ, Crabtree GR. Chemically induced proximity in biology and medicine. Science. 2018;359(6380):eaao5902. doi:10.1126/science.aao5902

*2:Science 2007(316)1030

*3:Science 2007(316)1033

*4: Nature 2011(476)332

*5:ライフサイエンス領域融合レビュー「花成ホルモン”フロリゲン”の構造と機能」 辻 寛之・田岡健一郎・島本 功 領域融合レビュー, 2, e004 (2013) DOI: 10.7875/leading.author.2.e004

*6:PNAS2011(108)18512

*7:理化学研究所プレスリリース

*8: 独立行政法人 理化学研究所 環境報告書2012 環境研究紹介① 「オーキシン」生合成主要経路を解明

*9:Nature2005(435)441

*10: Nature2005(435)446

*11: Nature2007(446)640

*12:タンパク質分解を誘導するというメカニズム自体が機能解明が難しかった一因でもあったそうな。 細胞内から関係するタンパク質をとってこようにも、相手が分解されていなくなるので取るに取れないみたいな(・・・リファレンスを散逸したので噂話程度で聞き流してください)。

*13:Science2010(327)1345

*14:Nature2014(512)49

*15:本当に。ここまで読んでくださった方、すみません。オススメ記事 Luck Is What Happens When Preparation Meets Opportunity サリドマイドをはじめとした免疫調節薬 (IMiDs)の作用機序:タンパク質分解誘導薬

*16: 今はBristol-Myers Squibb社?

*17:Nature2016(535)252

*18: 抗ガン作用を持つ新たなセレブロモジュレーターの開発 ファルマシア 2017(53)328

*19:急性骨髄性白血病に対し効果のある新たなセレブロモジュレーターの開発 ライフサイエンス新着論文レビュー

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

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

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

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

以下、訳。

トークトリアル 11 (パートC)

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

Developed at AG Volkamer, Charité

Dr. Jaime Rodríguez-Guerra, Dominique Sydow

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

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

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

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

学習の目標

理論

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

実践

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

レファレンス


理論

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

リガンドの結合は、リガンドとタンパク質ポケットの表面との間、あるいはタンパク質ータンパク質境界面における非共有結合性相互作用によって主に支配されています。このプロセスは、静電的相補性と形状相補性、誘導適合(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Å以内の距離にある2つの原子として定義されています。最も頻度の高いトップ100のリガンドータンパク質原子ペアは7つの相互作用のタイプに分類することができました(下の図を参照してください)。最も頻度高くみられた相互作用の中には、疎水性相互作用(hydophobic contact)、水素結合そしてπ-スタッキングといった、リガンドデザインにおいてよく知られており、広く使われている相互作用がありました。それらに続いて、弱い水素結合、塩橋、アミドスタッキング、そしてカチオンーπ相互作用といったものがありました。

f:id:magattaca:20200523143124p:plain

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

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

f:id:magattaca:20200523143205p:plain

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の出力には3つの列があります。

  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レイアウトを使って水平に配置された2つの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応用の説明

最初に、JavaScriptPythonのコードを混ぜていることに気づいたかもしれません!これは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)
以下は画像ですが、ノートブック上ではインタラクティブなビューワーが表示されます。
訳注ここまで

f:id:magattaca:20200523143924p:plain

viewer.render_image(),
(Image(value=b'', width='99%'),)
# 静的な出力
viewer._display_image()

f:id:magattaca:20200523143719p:plain

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')]}}

この辞書は2つのレベルからなります。最初のレベルは見つかった結合サイトです。一つのリガンドだけをドッキングしているので、一つだけ見つかることが予想されます。各結合サイトについて、8個のリストを含むもう一つの辞書があり、リストは各特定の相互作用のためのものです。各リストは最初の列に列名を持ち、つづいて(該当するものがあれば)データが格納されています。これを踏まえて、次のように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.DataFrameLIGCOOPROTCOO)の間に円柱(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)
以下は画像ですが、ノートブック上ではインタラクティブなビューワーが表示されます。
訳注ここまで

f:id:magattaca:20200523144541p:plain

app.center.render_image(),
(Image(value=b'', width='99%'),)
# 静的なプレビュー
app.center._display_image()

f:id:magattaca:20200523144613p:plain

ディスカッション

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

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

クイズ

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

訳以上

追記

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

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

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

RMSDについてsymmetry correctionを行ったか否かの2種類の記述があります。以下は、私の理解です。
リガンドの各原子を区別してラベルしている場合、分子に対称性があると、化学的には同じ構造でも原子のラベルとしては異なるという可能性があります(ex. ベンゼンの6つの炭素原子にラベルをつける場合)。ドッキングの結果の比較として、同じリガンドの異なる結合ポーズを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 OSPython 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ファイルには原子の座標のみで、結合関係の情報が入っていないようでしたので、その辺りが影響しているのかもしれません。

終わりに

これでトークトリアル11含め全てのトークトリアル終了となります。長々とお付き合いいただきありがとうございました!GitHubにノート入れておきましたのでよろしければご利用ください。

github.com

私は、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

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

f:id:magattaca:20200517172822p:plain
ProteinPlusトップページ

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

f:id:magattaca:20200517172921p:plain

素敵・・・。

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

f:id:magattaca:20200517173006p:plain
DoGSiteScorerの設定

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

f:id:magattaca:20200517173114p:plain
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で確認してみます。

f:id:magattaca:20200517173246p:plain
ポケットの可視化

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

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

  • 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)。ProteinPlusREST 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なのでうまくいっているようです。ドキュメンテーションには各ステータスコードの場合に得られる結果も記載されています。

f:id:magattaca:20200517173419p:plain
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

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

f:id:magattaca:20200517173513p:plain
GETの結果

このうちresiduesの一番最初(インデックス0)が最もスコアの良いポケットだそうです。Webで結果をダウンロードした時と同じですね。
residuesはPDBファイル形式となっています。今欲しい情報は結合ポケットの中心の座標と半径で、こちらはHEADERGeometric 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

インストールはpipcondaでできます。

# 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で3次元構造におこします。AllChem.EmbedMolecule()で3次元にしますが、性能が良いという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

# 3次元構造を生成してリストに格納
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)

f:id:magattaca:20200517174259p:plain

立体感でてますね!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]

f:id:magattaca:20200517174608p:plain
ドッキング結果の出力ファイルにおけるエラー

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

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

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

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

用意するものは3つ

  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

f:id:magattaca:20200517174802p:plain
ターミナルに随時結果が表示される

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

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

f:id:magattaca:20200517175038g:plain

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

スコアが良いものから悪いものへと順番に流れています。スコアの良いものは画面の右、悪いものは左側という印象でした。詳しい解析は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 トピック11(パートB)〜オンラインAPI/サービスを使った構造に基づくCADD〜

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

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

以下、訳。

トークトリアル 11 (パート B)

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

Developed at AG Volkamer, Charité

Dr. Jaime Rodríguez-Guerra, Dominique Sydow

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

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

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

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

学習の目標

理論

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

実践

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

レファレンス


理論

ドッキング

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

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

ドッキング

オンラインで無料で利用できるWebサービスが2、3あります。SwissDockとOPAL Webサービス(AutoDock Vinaを含みます)です。

SwissDock

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

OPAL Webサービス

医学生物的なアプリケーションはますます複雑となってきており、しばしば多数のプロセッサとメモリを有する大規模な高性能計算資源を必要とします。アプリケーションのデプロイの複雑さと、クラスターコンピューティング、グリッドそしてクラウドコンピューティングにおける進歩は、医学生物学研究をサポートする新しい方式を必要としています。サービスとしての科学ソフトウェア(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年の間に、1日あたりの実行に成功したジョブは205から411へと2倍以上に増加しました。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
        4文字の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',
                            # どうもこの下の1行が大事なようです。
                            # これ無しではエラー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を実行

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

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

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

通常、計算には5−15分かかるので、4番目のステップのために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サーバーを使って、最も可能性が高い結合サイトの中心と半径を計算しました。これら2つのサービスは異なるコミュニケーションインターフェースを使っていました。

Protein.PlusはRESTを使っていますが、KLIFSとは異なりswagger.jsonの定義を提供していないので、自分のリクエストを手動で構築しなければなりませんでした。幸い、Webサービスが単純だったので2、3個のリクエストが必要なだけ、、、と思っていました!!現在の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の内容はアプリケーションをプログラム上で扱うためのコードばかりで私は全くついていけませんでした。文章から理解する能力がないので絵で説明してほしいです・・・。

KLIFSで遊んだ話(T11パートAの備忘録)

トークトリアル11イントロに引き続きトークトリアル11パートAも耳慣れない言葉が多く内容のフォローが難しかったです。

と、いうわけで今回も素人による適当解説をしますよ!!

KLIFS

今回はKLIFS(Kinase-Ligand Interaction Fingerprints and Structures)について調べてみます。トークトリアルの説明は文字のみで、データの取得に用いただけだったのですが、実際にWebページ にアクセスしてみるととても素敵な感じでした!

f:id:magattaca:20200510211647p:plain
KLIFSトップページ

KLIFS Webページで遊ぶ

折角なのでT11パートAでコードで実施していた内容をWebページ上でも行なってみます。T11 パートAの内容は以下でした。

  1. Kinaseファミリー 「EGFR」 を検索
  2. PDB ID 3w33 を選択
  3. 描画

まずは検索から。上タブからSearchページに移動します。タンパク質やリガンドだけでなく色々な検索方法が用意されています。興味深いところでは結合ポケットのなかの水分子や、リガンド-キナーゼ相互作用パターンによる検索といったものまであります。目移りしてしまいますが、まずはT11通りにKinaseファミリーで検索します。

Search by classificationにて「FamilyEGFR」「OrganismHuman」として検索しました。

f:id:magattaca:20200510192718p:plain

ヒットした構造(355エントリー)がテーブルとして表示されました。下図に化合物2次元構造がありますが、Orthosteric ligandの列の化合物の名称にマウスオーバーするだけで描画されます。レスポンスも早くてすごい!

f:id:magattaca:20200510192743p:plain

またテーブル左上Plot results on Kinomeをクリックするとキノーム(ゲノム上のプロテインキナーゼ一式)の中で、どの位置に相当するか図示されます。

f:id:magattaca:20200510192821p:plain

検索結果のテーブルはExcelcsv形式でダウンロードすることも可能で、またDownload Structuresで構造をダウンロードすることもできるようです(スクロールしていくとボタンがあります。)

PDB ID: 3W33」の詳細を見るためにDetails列のShowをクリックして見ます。こんな感じ

f:id:magattaca:20200510192907p:plain

エントリーの詳細な情報に加えて3Dの描画やリガンドの活性値情報もあります。ChEMBLから取得されているようで、EGFR以外のキナーゼに対する活性も合わせてみることができるので、選択性やオフターゲットの雰囲気もつかめそうです。

スクロールするともっと面白い項目もあります。

f:id:magattaca:20200510193025p:plain

Kinase-ligand interaction patternは、リガンドがタンパク質のどの残基とどのような相互作用をしているか一目でわかるようになっています。私は結晶構造の描画をみても「格好いい!! で、どこが大事なんですかね???」となってしまうので、こういう見所(?)テーブルがあるのは非常に良いなと思います。また、一番下にひっそりとInteraction pattern searchなるものがあり、類似の相互作用パターンをもつ複合体をKLIFSから探してくれます。エモい・・・

KLIFSの背景

KLIFSのWebページの充実ぶりすごいですね。ちょっと背景を知りたいので論文を参照してみたいと思います。オンラインサービスに関する文献 Nucleic Acids Res. (2016), 44, 6, D365–D371 はオープンアクセスでした。

T11 パートAにも説明がありましたが、KLIFSはPDBからキナーゼの結晶構造を抜き出し、リガンドとの相互作用に関して分析処理を加えたデータベースです。特徴としてキナーゼの触媒ドメインに着目していることが挙げられます。

内部に踏み込む前にまずはキナーゼの構造について一般的な内容から

キナーゼの構造の特徴

こちらの日本語レビュー「 キナーゼを標的とした構造生物学および創薬の現状」が非常に分かりやすかったです。 *1 ざっくりと理解したことを紹介します。

キナーゼは基質をリン酸化するタンパク質ですが、リン酸基のソースはATPです。つまり共通して触媒キナーゼドメインの結合サイトにATPが結合します。
同じものが結合するならポケットも似ているだろ、ということでキナーゼ触媒ドメインはお互いにある程度似た構造をしています。
多くのキナーゼ阻害剤はこのATP結合領域に結合します。なので、異なるキナーゼ間で構造が似ていると、標的以外のキナーゼにも結合して阻害する可能性あります。
OpenTeachCADDで度々出てくる選択性の問題ですね。

よく保存された構造をもつキナーゼの結合ポケットですが、その特徴的な部分には名前がつけられています。例えば・・・

  • Hinge : ATPのアデニンと相互作用する残基部位
  • Gate Keeper : ATP結合ポケットの一番奥の残基部位。ATP結合サイトの立体構造への影響大
  • DFG motif : Asp-Phe-Glyの1文字表記。側鎖の動きが大きくDFG-in構造、DFG-out構造とよばれるコンフォメーションの変化がある

などなど。また、阻害剤にも結合するサイト・阻害様式によってType I ~ IV、共有結合型などの分類があるそうです。

こういうの格好いいですよね!共通の情報が整理されて見通し良くなってく感じが分野の発展の歴史を見ているようで楽しいです。

まとめると「医薬品としてキナーゼ阻害剤は大事だが選択性が課題となる。背景には結合サイトにおける保存されたキナーゼの構造的特徴があり、特徴に従い分類、キナーゼー阻害剤相互作用の細かな差異を分析することで課題解決の有用な情報となる。」・・・ってことですかね?

つまり、そういう情報まとめたKLIFSは素敵なリソース!!

KLIFS文献参照

やっとKLIFSの文献に戻ってFigure 2. を引用します。こちらはKLIFSがどのように複合体のアノテーションを行なっているかを説明した図です。

f:id:magattaca:20200510193619p:plain
Nucleic Acids Res. (2016), 44, 6, D365–D371 Fig.2 Annotation structural kinase-ligand interaction data in KLIFS.より引用

先に見た単語がちらほら見えますね。この図で面白いのはキナーゼーリガンド相互作用をフィンガープリントして表現しているFig.2 A)下部だと思います。各アミノ酸残基について、疎水性相互作用や水素結合といった7種類の相互作用がリガンドとの間にみられるか?というのを、0、1のビットで表記しています(kinase-ligand interaction finger print, IFP)。

KLIFSは触媒ドメインを含む85残基を解析の対象としています。したがって、各複合体のキナーゼーリガンド相互作用の情報は「85残基x相互作用7種=595ビット」のフィンガープリントにエンコーディングされることになります(?あってます?)。これだけ情報が圧縮されてるから「類似の相互作用の複合体は他にあるのか?」といった検索が容易になるんですね。なるほど。

さて、KLIFSでは相互作用のフィンガープリント化と結合ポケットの部位の分類(3 major pockets + 12 subpockets)を行なったことで、各複合体のリガンドの結合ポーズ(binding mode)の分類が可能になったそうです。
また、興味深いことにリガンドに加えて水分子についても解析が行われています。下に引用した図のように、保存された水分子のクラスター位置が見つかっています。こちらもWebページ上で確認、検索することが可能になっています。

f:id:magattaca:20200510193833p:plain
Nucleic Acids Res. (2016), 44, 6, D365–D371 Supplementary Figure S2 A)より改変して引用

以上、非常にざっくりですがKLIFSで行われていることの雰囲気です。分類と解析とが体系化されることで自動化が可能になり、データベースの継続的なアップデートに繋がっているようです。

KLIFSのAPI

KLIFSのWebページについてはだいたいわかったので、T11パートAではAPIを使ってどのように情報が取り出されていたのかを見ていきたいと思います。
取り上げるのは途中で定義された以下の3つの関数です。

  1. _all_kinase_families(): KLIFSデータベースに含まれているキナーゼファミリーの名前を取得
  2. _kinases_from_family(family, species="HUMAN"): 指定したキナーゼファミリーに含まれるキナーゼの名前を取得
  3. _protein_and_ligand_structure(*kinase_ids): 指定したKLIFSのキナーゼIDに紐づいたリガンド-キナーゼ複合体構造 (PDB)、タンパク質のみの構造(MOL2)、リガンド(SMILESリスト)を取得

REST?SWAGGER?

コードの中身を見る前にまずはよくわからない単語たち・・・

  • REST(Representational State Transfer): Web APIの仕様を決める上での基本的な考え方。リソースのURIと、アクセスするメソッド(HTTPのGET、POSTなどの)の組み合わせの設計方針。*2
  • SWAGGER: REST APIを定義、仕様を管理するためのフレームワーク。OpenAPI仕様をもとにした実装

・・・わからない。KLIFSのAPI をみてみます。

f:id:magattaca:20200510194315p:plain

たとえばKinase groupsについて情報が欲しければ、リソースのURI(Uniform Resource Identifier)「~/kinase_groups」に対して、GETメソッドを使えば良いということのようですね。
また、今回のトークトリアルでは、このSwaggerサービスを使うためのクライアントの生成をライブラリBravadoを使っておこなっています。Bravadoのドキュメンテーションこちらです。

bravadoによるSwagger用のクライアントの作成

それではbravadoを使って、KLIFS APIのURLからクライアントを作成します。URLからわかるようにjson形式となっているようです。

from bravado.client import SwaggerClient

# URLを定義
KLIFS_API_DEFINITIONS = "http://klifs.vu-compmedchem.nl/swagger/swagger.json"

# SwaggerClient : A client for accessing a Swagger-documented RESTful service
KLIFS_CLIENT = SwaggerClient.from_url(KLIFS_API_DEFINITIONS, config={'validate_responses': False})

作成したSwaggerClientの属性をdir()で確認して見ます。

print(dir(KLIFS_CLIENT))

#  ['Information', 'Interactions', 'Ligands', 'Structures']

これらが利用可能なリソースで、KLIFSのAPIに書かれている4種類と同じものがきちんと入っているようです。

f:id:magattaca:20200510194400p:plain

関数_all_kinase_families()の中身

こちらはKLIFSからKinaseファミリーを取得する関数です。こちらが何をしているか順番に見ていきます。

def _all_kinase_families():
    return KLIFS_CLIENT.Information.get_kinase_families().response().result

まず、リソースとしてinformationを選択し、そのkinase_familyに対してGETで情報を取得します(get_kinase_families())。Bravadoのドキュメンテーションによると、メソッドを実行することで戻り値としてHttPFutureを得るそうです(Futures and responses)。

確認してみます。

print(KLIFS_CLIENT.Information.get_kinase_families())

#  <bravado.http_future.HttpFuture object at 0x117469ef0>

確かに戻り値はHttpFutureオブジェクトとなっています。ここからレスポンスを得るにはHttpFuture.responce()とすればよく、この戻り値としてBravadoResponseインスタンスが得られます。Swaggerの結果を得るためにはこのインスタンスに対してBravadoResponse.resultとすれば良いそうです。

順番にやってみます。

print("レスポンス:", KLIFS_CLIENT.Information.get_kinase_families().response())

#  レスポンス: <bravado.response.BravadoResponse object at 0x117da7048>
print("type:", type(KLIFS_CLIENT.Information.get_kinase_families().response().result))
print("result", KLIFS_CLIENT.Information.get_kinase_families().response().result)

#  type: <class 'list'>
#  result ['A6', 'ABC1', 'AKT', 'ALK', 'AUR', 'Abl', 'Ack', 'Akt', 'Alk', 'Alpha', 'Aur', 'Axl', 'BCR', 'BRD', 'BUB', 'Bud32', 'CAMK-Unique', 'CAMK1', 'CAMK2', 'CAMKK', 'CAMKL', 'CASK', 'CCK4', 'CDC7', 'CDK', 'CDKL', 'CK1', 'CK2', 'CLK', 'Csk', 'DAPK', 'DCAMKL', 'DDR', 'DMPK', 'DYRK', 'EGFR', 'Eph', 'FAK', 'FAST', 'FGFR', 'Fer', 'G11', 'GRK', 'GSK', 'H11', 'Haspin', 'IKK', 'IRAK', 'IRE', 'InsR', 'JakA', 'JakB', 'KIS', 'LISK', 'LRRK', 'Lmr', 'MAPK', 'MAPKAPK', 'MAST', 'MLCK', 'MLK', 'MOS', 'Met', 'Musk', 'NAK', 'NDR', 'NEK', 'NKF1', 'NKF2', 'NKF3', 'NKF4', 'NKF5', 'NRBP', 'Other-Unique', 'PAN3', 'PDGFR', 'PDHK', 'PDK1', 'PEK', 'PHK', 'PIK', 'PIKK', 'PIM', 'PIP', 'PIPK', 'PKA', 'PKC', 'PKD', 'PKG', 'PKN', 'PLK', 'PSK', 'RAD53', 'RAF', 'RCK', 'RGC', 'RIO', 'RIPK', 'RSK', 'RSKL', 'RSKR', 'RSKb', 'Ret', 'Ror', 'Ryk', 'SCY1', 'SGK', 'SRPK', 'STE-Unique', 'STE11', 'STE20', 'STE7', 'STKR', 'Sev', 'SgK071', 'SgK493', 'SgK495', 'SgK496', 'Slob', 'Src', 'Syk', 'TAF1', 'TBCK', 'TIF1', 'TK-Unique', 'TKL-Unique', 'TLK', 'TOPK', 'TSSK', 'TTBK', 'TTK', 'Tec', 'Tie', 'Trbl', 'Trio', 'Trk', 'ULK', 'VEGFR', 'VPS15', 'VRK', 'WEE', 'WNK', 'Wnk', 'YANK']

無事BravadoResponseオブジェクトの生成と、結果として全キナーゼファミリーを要素としてもつリストを得ることができました。

Bravadoの動作の流れをまとめると・・・

  • SwaggerClient作成 > リソース選択 > get_URIで情報取得 > HttpFuture > response()でレスポンス取得 > BravadoResponse > resultで結果取得

・・・長い。

以上が定義した_all_kinase_families()関数になります。

関数_kinases_from_family(family, species="HUMAN")の中身

こちらは指定したGETメソッドの引数にキナーゼファミリーの名前と種(デフォルトの引数はHUMANに設定)を与えてキナーゼの名前を取得する関数です。
操作自体は先の例と同様です。

def _kinases_from_family(family, species="HUMAN"):
    return KLIFS_CLIENT.Information.get_kinase_names(kinase_family=family, species=species).response().result

関数_protein_and_ligand_structure()の中身

こちらはキナーゼのIDを引数として与え、以下の3つを返す関数です。、

  1. キナーゼIDに基づく複合体構造(PDB形式)
  2. タンパク質のみの構造(mol2形式)
  3. キナーゼIDに基づく共結晶構造中のリガンドのリスト(SMILESのみを抜き出して出力)

1、2はStructuresリソースから、3はLigandsリソースから得ています。

def _protein_and_ligand_structure(*kinase_ids):
    structures = KLIFS_CLIENT.Structures.get_structures_list(kinase_ID=kinase_ids).response().result
    molcomplex = KLIFS_CLIENT.Structures.get_structure_get_pdb_complex(structure_ID=structures[0].structure_ID).response().result
    protein = KLIFS_CLIENT.Structures.get_structure_get_protein(structure_ID=structures[0].structure_ID).response().result
    ligands = KLIFS_CLIENT.Ligands.get_ligands_list(kinase_ID=kinase_ids).response().result
    print(f"Chosen KLIFS entry with PDB ID {structures[0].pdb} with chain {structures[0].chain} and alternate model {structures[0].alt}")
    return molcomplex, protein, [ligand.SMILES for ligand in ligands]

EGFRを例として順番に見ていきます。準備としてまず例とするEGFRのkinase_IDを取得します。

EGFR_id_result = KLIFS_CLIENT.Information.get_kinase_ID(kinase_name="EGFR").response().result
print(EGFR_id_result)
print("number of result:", len(EGFR_id_result))

#  [KinaseInformation(HGNC='EGFR', family='EGFR', full_name='epidermal growth factor receptor', group='TK', iuphar=1797, kinase_ID=406, kinase_class='', name='EGFR', pocket='KVLGSGAFGTVYKVAIKELEILDEAYVMASVDPHVCRLLGIQLITQLMPFGCLLDYVREYLEDRRLVHRDLAARNVLVITDFGLA', species='Human', uniprot='P00533'), KinaseInformation(HGNC='Egfr', family='EGFR', full_name='epidermal growth factor receptor', group='TK', iuphar=0, kinase_ID=663, kinase_class='', name='EGFR', pocket='KVLGSGAFGTVYKVAIKELEILDEAYVMASVDPHVCRLLGIQLITQLMPYGCLLDYVREYLEDRRLVHRDLAARNVLVITDFGLA', species='Mouse', uniprot='Q01279')]
# number of result: 2

2つ結果が得られました。今回はhumanの情報をとりたいので一つ目のIDを取り出します。

EGFR_ID = EGFR_id_result[0].kinase_ID
print(EGFR_ID)

#  406

では_protein_and_ligand_structure()関数の1行目、複合体の構造のリストの取得です。ここではkinase_IDをリストで与えてやる必要があるので[406]として渡します。

EGFR_structures = KLIFS_CLIENT.Structures.get_structures_list(kinase_ID=[406]).response().result
print(len(EGFR_structures))    

#  324

全部で324個の構造情報が得られました。次の工程で複合体のPDBファイルを得るには、構造情報のうちstructure_idが必要となります。一つ目の構造情報で確認します。

EGFR_structure_id = EGFR_structures[0].structure_ID
print(EGFR_structure_id)

#  782

このstructure_idで複合体を取得します。今回は引数はintのままで大丈夫です。

EGFR_molcomplex = KLIFS_CLIENT.Structures.get_structure_get_pdb_complex(structure_ID=782).response().result
print(type(EGFR_molcomplex))
#  str

冗長なので省略しますがprintするとPDB形式になっていることが確認できます。ちなみに中身を見るとPDB_ID: 3w33でKLIFSのWebページで眺めていたものと同じでした。

同様にして同じstructure_idでタンパク質単体も取り出します。

EGFR_protein = KLIFS_CLIENT.Structures.get_structure_get_protein(structure_ID=782).response().result
print(type(EGFR_protein))
# str

こちらも省略しますがprintするとmol2形式になっていることが確認できます。

最後にリガンドのリストを取得します。EGFRのkinase_IDを使いますが、ここでもリスト[406]にして渡します。

EGFR_ligands = KLIFS_CLIENT.Ligands.get_ligands_list(kinase_ID=[406]).response().result
print(len(EGFR_ligands))

#  94

94個のリガンド情報を取得できました。ここからSMILESを取り出しますが、参考までにどのような情報があるか見てみます。

print(type(EGFR_ligands[0]))
print(dir(EGFR_ligands[0]))

# <class 'abc.ligandDetails'>
#  ['InChIKey', 'Name', 'PDB-code', 'SMILES', 'ligand_ID']

SMILES以外にもInChIKeyやNameなどの情報があるようです。

EGFR_SMILES = []
for lig in EGFR_ligands:
    smi = lig.SMILES
    EGFR_SMILES.append(smi)

print("number of ligands: ", len(EGFR_SMILES))
print("Example SMILES: ", EGFR_SMILES[0])

# number of ligands:  94
# Example SMILES:  P(=O)(OP(=O)(O)O)(OC[C@H]1O[C@@H](N2c3ncnc(N)c3N=C2)[C@H](O)[C@@H]1O)O

94個全ての構造のSMILESを取り出すことができました!これで_protein_and_ligand_structure()関数の全行程の確認ができました。

T11パートAのKLIFSに関する部分は、ここまで見てきたコードを踏まえて見返すとなんとなくやりたいことがわかりました。

最後に取り出したリガンドの構造を描画して見ましょう。

from rdkit import Chem
from rdkit.Chem import Draw

m = Chem.MolFromSmiles(EGFR_SMILES[0])
Draw.MolToImage(m)

f:id:magattaca:20200510194746p:plain

ADPでした。

T11パートAの結果を振り返って

以上、KLIFSのWebページ版とAPIとを比較してきました。T11パートAの結果と見比べることで、取り出されたリガンド構造に関して気づいた点があります。

T11パートAの最後、ケーススタディー: EGFRでEGFRのリガンドを取り出した時、同時に取り出したPDB共結晶構造中のリガンドが選ばれたと思っていました。実行時に得られたPDB_ID3w33とリガンドADPでした。ですが、上で見たようにKLIFSのWebページで同じPDB_ID3w33を確認するとorthosteric ligndはADPではない化合物でした。

この差が不思議だったのですが、APIを順番に見ていくことでT11 パートAの最後はkinase_IDで指定したリガンドを全て取り出し、その最初の一つを表示させているだけ、ということがわかりました。

もし同時にPDBのリガンドを指定して取得したいのなら、同じstructure_IDを使って以下のようにすれば良さそうです。

# mol2形式でリガンドを取得

ligand_782 = KLIFS_CLIENT.Structures.get_structure_get_ligand(structure_ID=782).response().result

# Mol2の中身は3次元で見づらいのでSMILES経由で2次元に戻してから描画

m782 = Chem.MolFromMol2Block(ligand_782)
smi782 = Chem.MolToSmiles(m782)
m782_2d = Chem.MolFromSmiles(smi782)
Draw.MolToImage(m782_2d)

f:id:magattaca:20200510194808p:plain

先に確認したように、「structure_id: 782」のPDB_IDは3W33で、KLIFSのWebページで見ていた複合体構造と同じでした。確かにここでは同じリガンドが取り出されているようです。

まとめ

以上、今回はT11 パートAで扱われていたデータベースKLIFSで遊んで見ました。素晴らしいWebページ版が提供されているので、まずはこちらで遊んでどのような情報が得られるのか感触をつかんでからの方がプログラムの内容を理解しやすのではないかな?と思いました。一目で結果を理解しやすいGUIベースのページは、プログラムが意図した通りの動きをしているかを検証する目的でも良さそうです。

プログラムの内容についてREST、SWAGGERは相変わらず定義はよく分かりませんが、ぼんやりと何がしたいかはわかったような気がしないでもないような?仕様を定義ってのは、文法みたいなイメージで良いのでしょうか???

素人が適当にしらべたことをだらだら書いているので色々と間違いが多いと思います。ご指摘いただければ幸いです。

TeachOpenCADD トピック11(パートA)〜オンラインAPI/サービスを使った構造に基づくCADD〜

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

トークトリアル11は論文発表時のTeachOpenCADDには含まれておらず新しく追加されたものです。
4つのノートブックに分かれており、ここではその2つ目を訳出します。
トークトリアル11で取り上げられているWebサーバー、サービスには既に停止しているものもあり、そのまま使えるものではなくなってしまっていますが、基本的な考え方や手法を知る上では役にたつかもしれません。

以下、訳。

トークトリアル 11 (パート A)

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

Developed at AG Volkamer, Charité

Dr. Jaime Rodríguez-Guerra, Dominique Sydow

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

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

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

KLIFデータベースにクエリを投げ、キナーゼの構造と阻害剤を一つ取得します。それから、PubChemで類似の化合物を検索します。

学習の目標

理論

実践

レファレンス


理論

ここまでで、さまざまなPythonライブラリをつかってオンラインWebサービスにクエリを投げる方法をみてきました。オンラインWebサービスからパイプラインを構築してみましょう!

パイプラインと含まれているWebサービスの解説

  1. キナーゼーリガンド相互作用フィンガープリントと構造のデータベース(KLIFS: Kinase-Ligand Interaction Fingerprints and Structures database)はアムステルダムのVU大学医薬品化学分野で開発されたデータベースで、触媒キナーゼドメインのタンパク質構造に関する情報(PDBから収集)と、リガンドとの相互作用に関する情報とを提供しています。このデータベースから精選されたタンパク質構造を取得し、そのリガンドの情報をつかってPubChemやChEMBLといった他のデータベースから類似したリガンドを取得することができます。
  2. KLIFによるリガンド情報を使って、PubChemで類似化合物を探します。
  3. タンパク質構造と候補となるリガンドをいくつか手に入れたら、OPALWebサービスで提供されている組み込みのVinaを使ってオンラインでドッキングを実施します。また、proteins.plus' DoGSiteScorerで化合物をドッキングする結合サイトの候補を探します。(パート B)
  4. 結果をnglviewで可視化し、相互作用を PLIPを使ってレポートします。(パート C)

このトークトリアル パート(A)では、KLIFSデータベースから入力として用いるキナーゼー阻害剤構造を取得し、候補阻害剤として評価することが可能な類似化合物をPubChemで検索します。関連する出力はディスクのdata/に書き込み、あとで使います。

KLIFS

キナーゼーリガンド相互作用フィンガープリントと構造のデータベース(KLIFS: Kinase-Ligand Interaction Fingerprints and Structures database)はアムステルダムのVU大学医薬品化学分野で開発されたデータベースで、触媒キナーゼドメインのタンパク質構造と、キナーゼ阻害剤の阻害様式に焦点をあてています。一貫したシステマティックなプロトコルに基づき、(現在、ヒトとマウスの)全てのキナーゼ構造とキナーゼリガンドの結合様式を直接互いに比較することができます。さらに、結合サイト全体を取り囲む85残基の分類がなされているので、例えば、キナーゼー阻害剤選択性を決める重要な相互作用を見つけるといった目的のため、キナーゼー阻害剤の相互作用パターンを互いに比較することができます。

PubChem

PubChemはアメリ国立衛生研究所(NIH、National Institutes of Health)のオープン化学データベースです。「オープン」とはつまり、自分の科学的なデータをPubChemに入れることができ、他の利用者がそのデータを使用するかしれないということを意味します。2004年の立ち上げ以来、PubChemは科学者、学生そして一般の人々のための化学に関する重要な情報源となってきました。毎月、我々のWebサイトとプログラムによるサービスは世界中の数百万の利用者にデータを提供しています。


実践

パイプラインの構築

KLIFSから情報を取得する

KLIFSデータベースから、ランダムなキナーゼファミリーのランダムなキナーゼ(mol2形式)と、対応するリガンド(SMILES)を選びます。また、パートBで結合ポケットを計算するために、タンパク質ーリガンド複合体のPDB構造を取得します。

from bravado.client import SwaggerClient
KLIFS_API_DEFINITIONS = "http://klifs.vu-compmedchem.nl/swagger/swagger.json"
KLIFS_CLIENT = SwaggerClient.from_url(KLIFS_API_DEFINITIONS, config={'validate_responses': False})
def _all_kinase_families():
    return KLIFS_CLIENT.Information.get_kinase_families().response().result

def _kinases_from_family(family, species="HUMAN"):
    return KLIFS_CLIENT.Information.get_kinase_names(kinase_family=family, species=species).response().result

def _protein_and_ligand_structure(*kinase_ids):
    structures = KLIFS_CLIENT.Structures.get_structures_list(kinase_ID=kinase_ids).response().result
    molcomplex = KLIFS_CLIENT.Structures.get_structure_get_pdb_complex(structure_ID=structures[0].structure_ID).response().result
    protein = KLIFS_CLIENT.Structures.get_structure_get_protein(structure_ID=structures[0].structure_ID).response().result
    ligands = KLIFS_CLIENT.Ligands.get_ligands_list(kinase_ID=kinase_ids).response().result
    print(f"Chosen KLIFS entry with PDB ID {structures[0].pdb} with chain {structures[0].chain} and alternate model {structures[0].alt}")
    return molcomplex, protein, [ligand.SMILES for ligand in ligands]
import random
import time

def random_kinase_structure():
    """
    ランダムなキナーゼファミリーのランダムなキナーゼを取得する
    """
    attempts = 20
    families = _all_kinase_families()
    while attempts:  # キナーゼIDの中には利用可能な構造がないものもあります
        family = random.choice(families)
        kinase = random.choice(_kinases_from_family(family))
        try:
            molcomplex, protein, ligands = _protein_and_ligand_structure(kinase.kinase_ID)
        except:
            attempts -= 1
            time.sleep(1)
        else:                   
            print("Chosen", kinase.name, "kinase with ID", kinase.kinase_ID, "from family", family)
            return molcomplex, protein, ligands
    print("Could not find a valid kinase. Try again!")
    return None, None, None


def kinase_structure_from_family(family):
    """
    キナーゼファミリーの名前(`_all_kinase_families()`を確認してください)が与えられると、ランダムな構造を取得します。
    """
    attempts = 20
    while attempts:  # キナーゼIDの中には構造のリストで見つけ流ことができないものもあります・・・。
        kinase = random.choice(_kinases_from_family(family))
        try:
            molcomplex, protein, ligands = _protein_and_ligand_structure(kinase.kinase_ID)
        except:
            attempts -= 1
            time.sleep(1)
        else:                   
            print("Chosen", kinase.name, "kinase with ID", kinase.kinase_ID, "from family", family)
            return molcomplex, protein, ligands
    print("Could not find a valid kinase. Try again!")
    return None, None, None

うまくいくか確認してみましょう!例えば、ランダムに選んだキナーゼファミリーの中からランダムにキナーゼが欲しいとします。その場合、次のようにrandom_kinase_structure()を使うことができます。

molcomplex, protein, ligands = random_kinase_structure()

#  Chosen KLIFS entry with PDB ID 3hng with chain A and alternate model A
#  Chosen FLT1 kinase with ID 483 from family VEGFR
ligands

#  ['Clc1ccc(NC(=O)c2c(NCc3ccncc3)cccc2)cc1']

訳注(2020/05)
上記は実行の度に結果が変わります。オリジナルのノートブックでは以下を得ています。

  • PDB ID 6cot、PKGファミリー(ID 42)、PRKG1キナーゼ
  • リガンド: ['P(=O)(OP(=O)(O)NP(=O)(O)O)(OC[C@H]1OC@@Hc3N=C2)C@H[C@@H]1O)O', 'CCCOc1ccc(c(c1C(=O)c2ccc(cc2)C(=O)N[C@@H]3CNC[C@H]3NC(=O)c4ccc5c(c4)cn[nH]5)F)OC']

以降の本文はこの出力を踏まえたものであることをお知りおきください。
訳注ここまで

試しにタンパク質をnglviewで、リガンドをrdkitで見てみましょう。まだnglviewをインストールしていなければ以下のセルを実行してください。

!pip install nglview  
!nglview install  
!nglview enable  
import nglview as nv
from tempfile import NamedTemporaryFile
from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole, MolsToGridImage # 分子を描画するのに必要です。

def preview_molecule_contents(contents, ext="mol2"):
    # これは一時的なファイルで、自動的に消去されます。
    v = nv.NGLWidget()
    v.add_component(contents, ext=ext)
    return v

def preview_smiles(smiles):
    print(smiles)
    return Chem.MolFromSmiles(smiles)

def multi_preview_smiles(*smiles):
    legends = [f"{s[:30]}..." for s in smiles]  # テキストが重なるのを避けるためにSMILES文字列を短くします。
    molecules = [Chem.MolFromSmiles(s) for s in smiles]
    return MolsToGridImage(molecules, molsPerRow=3, subImgSize=(300, 300), maxMols=len(molecules),
                           legends=legends, useSVG=True)
v = preview_molecule_contents(protein)
v
# NGLWidget()
v.render_image(),
#  (Image(value=b'', width='99%'),)
v._display_image()

f:id:magattaca:20200510084440p:plain

preview_smiles(ligands[0])

# Clc1ccc(NC(=O)c2c(NCc3ccncc3)cccc2)cc1

f:id:magattaca:20200510084506p:plain

PubChemを検索

上で取得したSMILES文字列を使って、PubChemで類似の化合物を検索しましょう。

ヒント:ChEMBLデータベースを使って似たような操作を行う方法はTalktorial T1 を参照してください。

import requests
def similar_compounds_pubchem(smiles, threshold=75, n_records=10):
    # 類似化合物をPubChemで検索する
    url = f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/similarity/smiles/{smiles}/JSON?Threshold={threshold}&MaxRecords={n_records}"
    r = requests.get(url)
    r.raise_for_status()
    key = r.json()['Waiting']['ListKey']
    # 応答は非同期なので、処理が終わったかどうか25秒間毎秒問い合わせます。
    attempts = 25
    while attempts:
        r = requests.get(f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/listkey/{key}/cids/JSON")
        r.raise_for_status()
        response = r.json()
        if 'IdentifierList' in response:
            cids = response['IdentifierList']['CID']
            break
        attempts -= 1
        time.sleep(1)
    else:
        raise IOError("Could not find matches for " + smiles)
    # 化合物IDを戻り値として得られますが、SMILESが必要です。
    url = f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/{','.join(map(str, cids))}/property/CanonicalSMILES/JSON"
    r = requests.get(url)
    r.raise_for_status()
    return [item['CanonicalSMILES'] for item in r.json()['PropertyTable']['Properties']]
def query_pubchem_for_similar_compounds(ligands):
    # 現在のキナーゼの最初のリガンドを取得します。
    smiles = ligands[0]
    # PubChem上の最も類似した10個の化合物を探します。
    return similar_compounds_pubchem(smiles, n_records=10)

取得した化合物をRDKitで描画します。

similar_smiles = query_pubchem_for_similar_compounds(ligands)
multi_preview_smiles(*similar_smiles)

f:id:magattaca:20200510084531p:plain

ケーススタディー: EGFR

対応するWikipediaの記事から改変して引用します。

上皮成長因子受容体(Epidermal Growth Factor Receptor; EGFR)はErbBファミリー受容体のメンバーで、密接に関係する4つの受容体チロシンキナーゼ、EGFR (ErbB-1)、HER2/neu (ErbB-2)、Her3 (ErbB-3)そしてHer4 (ErbB-4)のサブファミリーです。多くのがん種で、EGFRの発現あるいは活性に影響を与える変異が、発がんに繋がります。 ヒトにおいてEGFRと他の受容体チロシンキナーゼのシグナルの欠損はアルツハイマー病といった疾患と関連づけられており、一方で過剰発現はさまざまな種類の腫瘍の発達と関連づけられています。受容体細胞外ドメインのEGFR結合サイトの阻害、あるいは細胞内のチロシンキナーゼ活性の阻害によるEGFRシグナルの遮断により、EGFRを発現する腫瘍の成長を阻害したり患者の症状を改善することができます。

つまり、私たちはEGFRファミリーのメンバーをターゲットとすることができる阻害剤の候補を探すことに興味があります。さあ、上記と同じステップを繰り返しましょう!但し今回はこの特定のキナーゼファミリーをターッゲトとします。

egfr_molcomplex, egfr_protein, egfr_ligands = kinase_structure_from_family('EGFR')

# Chosen KLIFS entry with PDB ID 3w33 with chain A and alternate model A
# Chosen EGFR kinase with ID 406 from family EGFR

訳注(2020/05)
今回はたまたまオリジナルのノートブックと同じ結果が得られましたが、EGFRについては複数の構造が登録されているため、実行の度に取得されるPDB IDは異なります。
訳注ここまで

タンパク質を試しに見てみましょう。

v = preview_molecule_contents(egfr_protein)
v
# NGLWidget()
v.render_image(),
#  (Image(value=b'', width='99%'),)
v._display_image()

f:id:magattaca:20200510084623p:plain

EGFRのリガンド(ATP)を試しに見てみましょう。

preview_smiles(egfr_ligands[0])

#  P(=O)(OP(=O)(O)O)(OC[C@H]1O[C@@H](N2c3ncnc(N)c3N=C2)[C@H](O)[C@@H]1O)O

f:id:magattaca:20200510084651p:plain

(訳注 ADPの誤記だと思います。)

類似したリガンドを試しに見てみましょう。

similar_smiles_egfr = query_pubchem_for_similar_compounds(egfr_ligands)
multi_preview_smiles(*similar_smiles_egfr)

f:id:magattaca:20200510084709p:plain

トークトリアルの次のパートのためにディスクに結果を書き出しておきましょう!

import os
os.makedirs('data', exist_ok=True)
with open('data/similar_smiles.txt', 'w') as f:
    f.write('\n'.join(similar_smiles_egfr))
with open('data/protein.mol2', 'w') as f:
    f.write(egfr_protein)
with open('data/complex.pdb', 'w') as f:
    f.write(egfr_molcomplex)

ディスカッション

このノートブックでは利用可能な技術に応じて異なるWebサービスにアクセスし使用する方法を学びました。きちんとした説明のあるAPIから、実際のWebブラウザーを真似しようとする手作りのスクレイパーまで取り上げました。

ここまで、KLIFSデータベースでEGFRファミリーのメンバーを検索し、キナーゼの構造と(当然のことながら)その内在性リガンドであるATPを取得しました。次に、PubChemから75%以上の類似度をもつ化合物を取得しました。次のパートではこれらの中から一つを選んで構造モデリングを実施します。


クイズ

  • "EGFR"ファミリーについてKLIFは幾つのキナーゼを提供しているでしょうか?
  • 類似性検索をより厳密にすることはできますか?
  • 提案された候補のなかにすでに承認された阻害剤はあるでしょうか?(ヒント:PKIDBスクレイプして、SMILESのリストを確認することができます)

訳以上

追記

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

今回でてきたPubChemにはPythonからアクセスするためのライブラリPubChemPyもあります。ドキュメントをざっとみる限り類似性検索の機能はなさそうでした。

PubChemPyの使い方と注意点は「化学の新しいカタチ」さんが以下の記事でわかりやすく解説してくださっていますのでご参照いただければと思います。

Webスクレイピング備忘録(T11 イントロより)

トークトリアル11はオンラインAPI/サービスを利用しています。T1~10と比べてCADD以前にプログラミングの部分の敷居が高いという印象です。・・・正直コードが何をしているかプログラミングのできない私には分からなかったです(レベル低くてすみません)。

と、いうわけで素人による適当解説をしますよ!!WEBスクレイパーに俺はなる!!

題材

T11の初っ端、イントロから敷居の高い単語が連発でしたが、一番気になった「WkipediaProteinogenic amino acidからのアミノ酸のテーブルを取得する部分」を見ていきます。

f:id:magattaca:20200509235206p:plain

行われていることは以下の通りです。

手順 ライブラリ やってること
データの取得 requests Wikipediaの該当ページにアクセス
HTMLデータを取得
データの抽出 BeautifulSoup HTMLの文字列から
目的のテーブルの場所を指定して
情報を抜き出す
テーブルの再構成 Pandas 抜き出した情報をDataFrameで
テーブルに再構成して表示

では順番に見ていきます。。*1

requestsによるHTMLデータの取得

「入門 Python3 9章」によると、WWWの単純な形としてウェブクライアントが

  1. ウェブサーバーにHTTP(Hypertext Transfer Protocl)で接続し
  2. サーバーにURL(Uniform Resource Locator)を要求し、
  3. HTML(Hypertext Markup Language)を受け取る

という流れになっているそうです。で、Pythonはこれが得意なんだそうな。

クライアントがサーバーに要求(リクエスト)を送って応答(レスポンス)を受け取るウェブのデータの交換の標準プロトコルHTTPで、そのHTTPメソッドにGETPOSTというものがあるそうです。

では、Pythonrequestsを使ってGETリクエストをおくり、指定したURLのHTML情報を受け取ります。(requests.get(URL)

import requests

r = requests.get("https://en.wikipedia.org/wiki/Proteinogenic_amino_acid")

requests.get関数の戻り値として得たResponseオブジェクトにはリクエストがうまくいったかどうかの情報(ステータス)も含まれています。HTTPステータスコードという3桁の数字で表されており、成功した場合は「2XX」、失敗すると「4XX」や「5XX」といったコードになります(Xも数字)。

今回のリクエストはうまくいったのか?Responseオブジェクトのstatus_code属性を調べます。

print(r.status_code)
#  200

200なので上手くいっているようです。URLをわざと間違えて失敗してみます。

r_error = requests.get("https://en.wikipedia.org/wiki/shippai")

print(r_error.status_code)
# 404

404 Not Found : 未検出 ! 存在しないページを要求したのでちゃんと失敗しました。4XXはクライアント側の誤りで、サーバー側がおかしい時は5XXというステータスコードになるそうです。

リクエストが失敗した時に備えて、エラーを認識してスクリプトを停止する(例外処理)メソッドも用意されています(raise_for_status())。こちらもT11で使われていました。

r_error.raise_for_status()

"""
    ---------------------------------------------------------------------------

    HTTPError                                 Traceback (most recent call last)

    <ipython-input-13-af32409c08e5> in <module>
    ----> 1 r_error.raise_for_status()
    ~~~ 省略~~~

        939 
        940         if http_error_msg:
    --> 941             raise HTTPError(http_error_msg, response=self)
        942 
        943     def close(self):

    HTTPError: 404 Client Error: Not Found for url: https://en.wikipedia.org/wiki/Shippai
"""

「URLが見つからない」というエラーが出ました。親切ですね!

では上手くいった場合どのような情報が取得できているのでしょうか? Resposeオブジェクトのtext属性で確認できるそうです。長いので最初の300文字取り出します。

print(r.text[:300])

"""
    <!DOCTYPE html>
    <html class="client-nojs" lang="en" dir="ltr">
    <head>
    <meta charset="UTF-8"/>
    <title>Proteinogenic amino acid - Wikipedia</title>
    <script>document.documentElement.className="client-js";RLCONF={"wgBreakFrames":!1,"wgSeparatorTransformTable":["",""],"wgDigitTransformTable":["",""],"wg
"""

Wikipediaの該当ページProteinogenic amino acidソースを表示させてみると上の出力と一致しました。HTMLの情報をちゃんと取得できているようです。

f:id:magattaca:20200509235557p:plain

BeautifulSoupによる構文解析

取得したHTMLの情報から該当の場所を探し出します。といってもHTMLの書式が分からない私にはさっぱりです。でも大丈夫! ......そう、「Google Chrome」ならね! *2

「右クリック → 検証」だ!

f:id:magattaca:20200509235709p:plain

HTMLをただ表示するだけではありません。HTML上でマウスオーバーするだけで、通常の表示ページでの該当箇所がハイライトされて対応を確認することができます。

f:id:magattaca:20200509235744p:plain

ここを頼りにHTML構文解析を行なって該当のテーブルの情報を抜き出していきます。関係するHTMLのタグを先に確認しておきましょう。HTMLクイックリファレンスがとても便利なページでした。*3

初歩の初歩で恐縮ですがタグはHTMLの目印で、「<●●>(開始タグ)~<\/●●>(終了タグ)」という形で目印をつけたい部分を囲んで指定します。(タグも知らずに偉そうに解説かいてごめんなさい。)*4

タグ 内容
< h1 > ~ < h6 > 見出し。hはheadingの略で、
数字1~6は見出しの上位、下位を示す(1が最上位)
< span > ひとかたまりの範囲として定義。
囲んだ範囲にスタイルシートを適用(cf. DIV)
style属性、id属性、class属性など
< table > テーブルを作成
< thead > テーブルのヘッダ行を定義
< tbody > テーブルのボディ部分を定義
< tr > テーブルの一行を定義
< th > テーブルの見出しセルを作成
< td > テーブルのデータセルを作成

f:id:magattaca:20200510000139p:plain

HTMLの書式もだいたいわかったのでBeautifulSoupを使ってデータを抽出していきます。先に見たようにrequestsで取得したHTMLテキストはResponseオブジェクトの.text属性でした。これをBeautifulSoupに渡してやります。

from bs4 import BeautifulSoup

html = BeautifulSoup(r.text)

BeautifulSoupのHTML解析では、HTMLタグをメソッドfind()find_all()を使って検索し、該当箇所を検索、指定します。find(タグ)は引数に一致する最初の一つfind_all(タグ)は一致する全ての要素を取得します。

まずは、目的のテーブル(General chemical properties)の見出しを探します。span要素で、id属性「id="General_chemical_properties"」 となっていました。一つしかないのでfind()で探します。

header = html.find('span', id="General_chemical_properties")

print(type(header))
print(header)

# <class 'bs4.element.Tag'>
# <span class="mw-headline" id="General_chemical_properties">General chemical properties</span>

望みのタグの箇所を指定できているようです。テーブルはこの見出しの下にあるので、順に辿っていきます。find_all_next()を使うことで以降の要素を全て取得することができ、さらにインデックスで要素を指定できます。

下図の通りtable要素は見出し以降の5番目(index 4)の位置に相当します。

f:id:magattaca:20200510000333p:plain

table = header.find_all_next()[4]

冗長なので省略しますが、print(table)として確認すると< table > ~ < /table >で囲まれる要素が取得されていることがわかります。*5

取り出したい情報はテーブルからヘッダを除いたボディ部分なのでtbodyタグを選択し抽出します。

table_body = table.find('tbody')

あとはここから取り出した情報をテーブルとして再構成すればOKです。

データの抽出とデータフレームの作成

DataFrameを作成する準備として、リストのリストを作成します。外側のリストはテーブルの各行を要素とするリストで、内部のリストは各行の各セルを要素とするリストです。

テーブルの各行はHTMLのtrタグを用いることで識別でき、各セルはtdにより識別できます。各行(row)、各セル(cell)についてforループを回すことでデータを順番に取り出していきます。forループでは、データ全体を格納するリストdataの中に、各行毎に空のリスト[]を追加したのち、各セルの中身を該当の行のリスト(リストdataの中で一番新しい最後の要素data[-1])に追加していきます。

data = []
for row in table_body.find_all('tr'):
    cells = row.find_all('td')
    if cells:
        data.append([])
    for cell in cells:
        cell_content = cell.text.strip()
        try:  # 可能であればfloatに変換します
            cell_content = float(cell_content)
        except ValueError:
            pass
        data[-1].append(cell_content)
print(len(data))

#  22

22の行がdataに格納されました。アミノ酸って20個じゃなかったっけ?と思いましたが、テーブルを見直すとSelenocysteinePyrrolysineが含まれていました。こんなアミノ酸知らなかった。。。

あとはデータフレームにするだけです。

import pandas as pd

pd.DataFrame.from_records(data)
0 1 2 3 4 5
0 A Ala 89.09404 6.01 2.35 9.87
1 C Cys 121.15404 5.05 1.92 10.7
2 D Asp 133.10384 2.85 1.99 9.9
3 E Glu 147.13074 3.15 2.1 9.47
4 F Phe 165.19184 5.49 2.2 9.31
5 G Gly 75.06714 6.06 2.35 9.78
6 H His 155.15634 7.6 1.8 9.33
7 I Ile 131.17464 6.05 2.32 9.76
8 K Lys 146.18934 9.6 2.16 9.06
9 L Leu 131.17464 6.01 2.33 9.74
10 M Met 149.20784 5.74 2.13 9.28
11 N Asn 132.11904 5.41 2.14 8.72
12 O Pyl 255.31000 ? ? ?
13 P Pro 115.13194 6.3 1.95 10.64
14 Q Gln 146.14594 5.65 2.17 9.13
15 R Arg 174.20274 10.76 1.82 8.99
16 S Ser 105.09344 5.68 2.19 9.21
17 T Thr 119.12034 5.6 2.09 9.1
18 U Sec 168.05300 5.47 1.91 10
19 V Val 117.14784 6 2.39 9.74
20 W Trp 204.22844 5.89 2.46 9.41
21 Y Tyr 181.19124 5.64 2.2 9.21

できました!

最初からPandasを使う

ちなみに最初からPandasを使って直接HTMLのtableを取得することも可能だそうです。*6

url = "https://en.wikipedia.org/wiki/Proteinogenic_amino_acid"

dfs = pd.read_html(url)

print("ページのテーブルの数: ", len(dfs))

print("最初のテーブルを出力")

dfs[0]

# ページのテーブルの数:  9
# 最初のテーブルを出力
Amino acid Short Abbrev. Avg. mass (Da) pI pK1(α-COOH) pK2(α-+NH3)
0 Alanine A Ala 89.09404 6.01 2.35 9.87
1 Cysteine C Cys 121.15404 5.05 1.92 10.70
2 Aspartic acid D Asp 133.10384 2.85 1.99 9.90
3 Glutamic acid E Glu 147.13074 3.15 2.10 9.47
4 Phenylalanine F Phe 165.19184 5.49 2.20 9.31
5 Glycine G Gly 75.06714 6.06 2.35 9.78
6 Histidine H His 155.15634 7.60 1.80 9.33
7 Isoleucine I Ile 131.17464 6.05 2.32 9.76
8 Lysine K Lys 146.18934 9.60 2.16 9.06
9 Leucine L Leu 131.17464 6.01 2.33 9.74
10 Methionine M Met 149.20784 5.74 2.13 9.28
11 Asparagine N Asn 132.11904 5.41 2.14 8.72
12 Pyrrolysine O Pyl 255.31000 ? ? ?
13 Proline P Pro 115.13194 6.30 1.95 10.64
14 Glutamine Q Gln 146.14594 5.65 2.17 9.13
15 Arginine R Arg 174.20274 10.76 1.82 8.99
16 Serine S Ser 105.09344 5.68 2.19 9.21
17 Threonine T Thr 119.12034 5.60 2.09 9.10
18 Selenocysteine U Sec 168.05300 5.47 1.91 10
19 Valine V Val 117.14784 6.00 2.39 9.74
20 Tryptophan W Trp 204.22844 5.89 2.46 9.41
21 Tyrosine Y Tyr 181.19124 5.64 2.20 9.21

・・・めっちゃ楽やん。

まとめ

以上、Webスクレイピングでした。ウェブのやりとりの仕組みからHTMLの書式までど素人には辛い内容ですね。T11の残り、本体パートA~Cまであるのですが先が思いやられますね!

色々間違っていそうなのでご指摘いただければ幸いです。

ではでは