TeachOpenCADD トピック１１（パートB）〜オンラインAPI/サービスを使った構造に基づくCADD〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル１１をもととしておりCC BY 4.0のライセンスに従っています。
トークトリアル１１は論文発表時のTeachOpenCADDには含まれておらず新しく追加されたものです。
４つのノートブックに分かれており、ここではその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. 結果を評価し既知のデータと比較
入力構造を取得したあとで、ドッキングソフトウェアを使って良いタンパク質-リガンドポーズを見つけます。
学習の目標
理論
- ドッキングの基本
- 利用可能なソフトウェア
- 知られている限界
実践
- 構造の準備
- 結合サイトの推定
- ドッキング計算の実行
- 結果の保存
レファレンス
- Structure‐Based Virtual Screeningの章 Ch. 6.8 in "Applied Chemoinformatics: Achievements and Future Opportunities" (2018)
- ドッキングソフトウェアをベンチマークテストする方法 (J. Med. Chem. (2006), 49, 23, 6789–6801)
- ドッキングソフトウェアについてのレビュー (Biophysical Reviews (2017), 9, 2, 91–102)
- DoGSiteScorer：結合サイトを同定するプログラム。
- SwissDock：ドッキングプログラム
- OPAL Webサービス：必要に応じてリモートでソフトウェアを実行可能
- AutoDock Vina：ドッキングプログラム (J Comput Chem. (2010), 31, 455–461)
理論
ドッキング
Pagadala, Syed & Tuszynski, Biophysical Reviews 2017, 9, 2, 91–102 から引用して改変:
ドッキングは標的タンパク質の結合サイト内における低分子化合物の振る舞いを探索する手法です。適用されたソフトウェアは探索アルゴリズムを実行し、最小エネルギーに収束するまでリガンドの配座を再帰的に評価します。最後に、候補となるポーズをランク付けするために、親和性スコア関数ΔG （合計の単位はkcal/mol)を静電エネルギーとファンデルワールスエネルギーの合計として用います。生物学的な系におけるこれらの特定の相互作用のドライビングフォースは、結合サイトの表面とリガンドあるいは基質の形状と静電的な状態の相補性への志向です。
これらのスコア関数は計算が速くなるように合わされているので、他の分子モデリング手法と比較してしばしば精確性に劣ります。通常、これらは次に基づいています。
- 分子力学法原理
- 知識に基づくポテンシャル
- 形状と幾何学的な相補性
探索空間の次元を減らすために、いくつかの近似がよく適用されます。
- タンパク質の構造はほとんどが剛直なものと考えられますが、探索空間に近接するいくつかの側鎖は限られた配座セットの探索が許可されていることもあります（ロータマー）。
- リガンドはバーチャルスクリーニングの研究では剛直なものとして考えられますが、より詳細な計算においては結合の二面角を自由に探索することがよく許可されます。二つのオプションの折衷案は、とりうる配座のセットを予め定義しておくことです。
既存ソフトウェアの例
- 商用
- GOLD
- Schrödinger
- FlexX
- 無料（あるいはアカデミックフリー）
- AutoDock
- AutoDock Vina
- DOCK
- OpenEye
知られている限界
次の近似は計算においてアーティファクトを取り込む原因となりうります。
- タンパク質のほとんどを剛直したものとして扱うので、タンパク質-リガンド結合の動的、適応的な性質はほとんど探索されません。このため偽陽性につながります。リガンドが結合ポケットで適切なポーズを見つけたとしても、タンパク質が近接するエネルギー安定な配座を探索できないと、このポーズは保証されません。言い換えれば、リガンドが結合ポケットにとどまるかどうか確認するために短時間の分子動力学（MD）トラジェクトリを計算することが常に推奨されています。
- スコア関数は解く計算コストが低くなければなりません。良いポーズと悪いポーズを区別するのに十分な精度がある一方で、最も良い複数のポーズを並べ替えるには問題があることがあります。例えば、最も人気のあるドッキングプログラムは、実験で明らかとされているポーズを計算で見いだすことができる一方で、このポーズが提案されたポーズのセットの中で最も良いものであることは稀です。
- 計算コストを減らすために、ドッキング手順はタンパク質の一部分（通常、結合ポケットとして知られている場所の周辺）でのみ実行されます。CADDパイプラインにおいて、正しい結合ポケットを選択することはもう一つの課題となります。
- 計算の精確性を高めるために、適宜構造を削減する必要があります。アミノ酸とリガンドのプロトン化状態は正しい状態を得るにはこつがいることがあり、特に互変異性体（の可能性がある）の場合難しくなります。これは誤った結果を得る、もう一つの原因となります。
これらの限界にも関わらず、依然としてドッキング計算は全ての医薬品研究室で人気のある手法で、他の種類の分子シミュレーションとともに使われています。このトークトリアルのパートでは次の方法を学びます。
- パートAで取得したタンパク質とリガンドの準備（ローカルで実行）
- もっとも可能性が高い結合ポケットの推定（オンラインで実行）
- ドッキング計算の実行（オンラインで実行）
構造の準備
OPALWebサービスにあるVinaを使います。しかし、あらかじめ構造を準備することが必要となります。オススメのアプローチはAutoDockToolsにある準備のためのスクリプトを使うことです。これらのツールはツール自体でディストリビューションされており、Python2としか互換性がありませんが、私たちはタンパク質とリガンドの準備に必要な部分を含むPython 3で利用可能なフォークを準備しました。このトークトリアルで必要なことを行うには十分のはずです。
訳注(2020/05)
AutoDockToolsのライセンスは、非商用の科学研究で使用する場合は無料ですが、商用で使用するには正式な登録手続きが必要とのことです。正確な情報はソフトウェアのライセンスを確認してください。
訳注ここまで
結合ポケット予測
ドッキング計算は、合理的に狭い探索空間で実施すると最もパフォーマンスがよくなり、この空間は通常、一つの結合ポケットを含みます。最も良い結合ポケットを計算に基づいて予測するために、Protein.plusで無料でオンラインで利用できるDoGSiteScorerを使うことができます。
Proteins.plus DoGSiteScorer
- 役割： いくつかのCADDツールのためのインタラクティブなWenインターフェース
- Webサイト： http://proteins.plus
- API： RESTに基づくものがあります。
requestsをそのまま適用するだけで十分です。 - ドキュメンテーション： https://proteins.plus/help/index
- 文献：
自動化されたタンパク質活性部位の予測は、大規模なタンパク質機能予測と分類、そしてドラッガビリティ（druggability、医薬品の標的可能性）を推測する上で必須となります。ここで我々は、画像処理に端を発するガウシアン差分（Difference of Gaussian, DOG）の方法に基づきタンパク質の結合サイトを予測する、新規な構造に基づく手法、DoGSiteを提案します。既存の手法とは対照的に、DogSiteは予測したポケットをサブポケットへと分割し、活性部位のトポロジーの洗練された表現を明らかにします。DoGSiteはPDBBindとscPDBデータセットの92%以上で精確に結合ポケットを予測し、現在利用可能な手法のうち最もパフォーマンスの良いものと並ぶ性能があります。PDBBindデータセットの63%で、見つかったポケットはより小さなサブポケットへと細分化することができました。予測の87%で、共結晶化されたリガンドは正確に一つのサブポケットに含まれていました。さらに、リガンドとポケットの組み合わせの網羅する範囲を考慮に入れることで、より精確な予測性能の評価を導入しました。90%の事例で、DoGSiteは少なくともリガンドの半分を含むポケットを予測しました。70%の事例で、さらに各ポケット自体の1/4以上が共結晶化されたリガンドで覆われていました。サブポケットを考慮することで、適用範囲が広がり、後者の性能評価において成功確率83%となりました。
ドッキング
オンラインで無料で利用できるWebサービスが２、３あります。SwissDockとOPAL Webサービス（AutoDock Vinaを含みます）です。
SwissDock
- 役割： ドッキング計算を実行します
- Websサイト： http://www.swissdock.ch
- API： SOAPに基づくものがあります。正式なクライアントはありません。
suds-communityを使ってください。 - ドキュメンテーション： http://www.swissdock.ch/pages/soap_access
- 文献：
SwissDockはターゲットタンパク質と低分子化合物の間で生じうる分子間相互作用を予測するWebサービスです。 SwissDockはドッキングソフトウェアEADock DSSに基づいていますが、そのアルゴリズムは次のステップからなります。
1. ボックスの中（ローカルドッキング）あるいは全ての標的のくぼみの近傍（ブラインドドッキング）で、結合モードを数多く生成します。 2. 同時に、グリッド上でCHARMMエネルギーが予測されます。 3. 最も好ましいエネルギーの結合モードがFACTSで評価され、クラスター化されます。 4. 最も好ましいクラスターをオンラインで可視化し、自分のコンピュータにダウンロードすることができます。
OPAL Webサービス
- 役割： サービスとしてのCADD （CADD as a service)
- Webサイト： http://nbcr-222.ucsd.edu/opal2/dashboard
- API： SOAPに基づくものがあります。正式なクライアントはありません。
suds-communityを使ってください。 - ドキュメンテーション： http://nbcr-222.ucsd.edu/opal2/dashboard?command=docs (現在オフラインとなっています）
- 文献：
医学生物的なアプリケーションはますます複雑となってきており、しばしば多数のプロセッサとメモリを有する大規模な高性能計算資源を必要とします。アプリケーションのデプロイの複雑さと、クラスターコンピューティング、グリッドそしてクラウドコンピューティングにおける進歩は、医学生物学研究をサポートする新しい方式を必要としています。サービスとしての科学ソフトウェア（Scientific Software as a Service、SaaS）は、簡単な標準化されたWebインターフェースを通して、医学生物学的なアプリケーションに、大規模に実行可能で透明性のあるアクセスを可能とします。この目的のため、私たちは2007年8月にMEMEと名付けられたバイオインフォマティクスのアプリケーションをサポートするための運用Webサーバーを構築しました（ http://ws.nbcr.net ）。以来、このサーバーは成長を続け、AutoDockとAutoDock Vinaによるドッキング分析、PDB2PQRとAPBSを使った静電計算、そしてSMAPを使ったオフターゲット分析を含むようになりました。サーバー上のアプリケーションは全てOPALで動いています。OPALは、単純なXML環境設定ファイルを書くことで、科学的なコードに一切変更を加えることなく、科学アプリケーションをWebサービスとして簡単に使えるようにするツールキットです。OPALにより、我々の全てのアプリケーションにWebフォームに基づくアクセスとプログラムによるアクセスのどちらでもアクセスすることができるようになります。OPALツールキットは現在、国立生物医学計算資源（National Biomedical Computation Resource、NBCR)と、系列下の共同研究、サービスプロジェクトによる多数の人気のあるアプリケーションへのSOAPに基づくWebサービスアクセスをサポートしています。さらに、OPALはプログラムによるアクセスが可能なので、VisionやKepler、Nimrod/KそしてVisTrailsといった様々なワークフローツールを通して我々のアプリケーションにアクセスすることができます。2007年8月中旬から2009年の終わりにかけて、239,814個のジョブの実行に成功しました。2008年から2009年の間に、１日あたりの実行に成功したジョブは205から411へと２倍以上に増加しました。OPALで実現されているサービスモデルは様々なアプリケーションにとって有用です。Webサービスインターフェースによる他のアプリケーションとの相互運用を提供し、そしてアプリケーション開発者が科学的ツールとワークフローの開発に集中することができるようにします。Webサーバーは次で利用可能です（ http://ws.nbcr.net ）。
実践
パートAからファイルを取得
まず最初に、パート Aで取得したファイルを指定する変数をいくつか定義します。
PROTEIN: 標的タンパク質の構造を含む.mol2ファイルへのパス。リガンドやイオンは含まないCOMPLEX: 標的タンパク質の構造と元々のリガンドを結合サイトに含む.pdbファイルへのパスSMILES_FILE: PubChemで見つかった全ての類似化合物を表すSMILESを含む.txtファイルへのパスsmiles:SMILES_FILEに含まれるSMILES文字列のリスト
PROTEIN = "data/protein.mol2" COMPLEX = "data/complex.pdb" SMILES_FILE = "data/similar_smiles.txt" with open(SMILES_FILE) as f: smiles = [line.strip() for line in f]
SwissDockを利用する
SwissDockはSOAPインターフェースを使っているので、sudsをインストールする必要があります。
注： SwissDockのサーバーは最近機能していません。代わりに下のOPALヘ進んでください！ !pip install suds-community
from suds.client import Client import zlib import string import requests def swissdock_client(): # 現在サーバーが落ちているようです。。。 # http://swissdock.vital-it.ch/soap/ は503利用不能を返します。 # 誤ったドメインを指定しているからでしょうか。。。修正しましょうか？ SWISSDOCK_WSDL_URL = "http://www.swissdock.ch/soap/wsdl" r = requests.get("http://www.swissdock.ch/soap/wsdl") r.raise_for_status() WSDL = r.text.replace("http://swissdock.vital-it.ch/soap/", "http://www.swissdock.ch/soap/") with open("data/swissdock.wsdl", "w") as f: f.write(WSDL) HERE = _dh[0] return Client(f"file://{HERE}/data/swissdock.wsdl") def prepare_protein(client, protein): """ 与えられたPDBファイル（中身は文字列）に対してPSFとCRDを返す。 """ encoded_protein = zlib.compress(protein.encode('utf-8')) job_id = client.service.prepareTarget(target=encoded_protein) while True: result = client.service.isTargetPrepared(jobID=job_id) if result is None: raise ValueError("No such a job present") if result in (False, "false", 0): time.sleep(5) else: # 準備できました！ break protein_files = client.service.getPreparedTarget(job_id) if protein_files is None or len(protein_files) != 2: raise ValueError("Could not prepare protein!") return protein_files def prepare_ligand(client, ligand): """ 与えられたMOL2ファイル（中身は文字列）に対してPDB、RTFとPARを返す。 リガンドはあらかじめプロトン化されていなければなりません! """ encoded_ligand = zlib.compress(ligand.encode('utf-8')) job_id = client.service.prepareLigand(ligand=encoded_ligand) while True: result = client.service.isLigandPrepared(jobID=job_id) if result is None: raise ValueError("No such a job present") if result in (False, "false", 0): time.sleep(5) else: # 準備できました！ break ligand_files = client.service.getPreparedLigand(job_id) if ligand_files is None or len(ligand_files) != 3: raise ValueError("Could not prepare ligand!") return ligand_files def dock(client, protein, ligand, name=None): protein_psf, protein_crd = prepare_protein(client, protein) ligand_pdb, ligand_rtf, ligand_par = prepare_ligand(client, ligand) if name is None: name = "teachopencadd" + ''.join([random.choice(string.ascii_letters) for _ in range(5)]) job_id = client.service.startDocking( protein_psf, protein_crd, ligand_pdb, [ligand_rtf], [ligand_par], name) if job_id in (None, "None"): raise ValueError("Docking job could not be submitted") while not client.service.isDockingTerminated(job_id): time.sleep(5) all_files = client.service.getPredictedDockingAllFiles(job_id) with open('docking_results.zip', 'w') as f: f.write(all_files) target, docked = client.service.getPredictedDocking(job_id) client.service.forget(job_id) return target, docked def smiles_to_pdb(s, out='output.pdb'): m = Chem.AddHs(Chem.MolFromSmiles(s)) AllChem.EmbedMolecule(m) if out is None: return Chem.MolToPDBBlock(m) Chem.MolToPDBFile(m, out)
try: import Mol2Writer except ImportError: # RDKitのMol2 writer/readersを手に入れるための醜いハック import os working_dir = os.getcwd() os.chdir(_dh[0]) !wget https://raw.githubusercontent.com/rdkit/rdkit/60081d31f45fa8d5e8cef527589264c57dce7c65/rdkit/Chem/Mol2Writer.py > /dev/null os.chdir(working_dir) import Mol2Writer
def step_03_swissdock(protein_pdb, ligand_smiles): ligand = Chem.AddHs(Chem.MolFromSmiles(ligand_smiles)) AllChem.EmbedMolecule(ligand) ligand_mol2 = Mol2Writer.MolToCommonMol2Block(ligand) client = swissdock_client() return dock(client, protein_pdb, ligand_mol2)
現在SwissDockサーバーが利用できないので、下のセルは実行しないようにしています。実行したいのであれば、まずタンパク質のmol2ファイルをPDBに変換する必要があります。OpenBabelを使って行うことができます（パートCでインストールします）。ご自由に都合の良い
mol_to_pdb関数を定義してください。
# TODO: mol2からpdbに変換すること！ # protein_pdbblock = mol2_to_pdb(PROTEIN, out=None) step_03_swissdock(protein_pdbblock, smiles[0]) ### OPAL Webサービスを使ったドッキングの実行
SwissDockは現在稼働していません。そこで手を変え別のWebサービスを利用しましょう。インターフェースは少々より原始的ですが、稼働するはずです。ですが、タンパク質とリガンドはローカル環境でAutoDockToolsを使って準備しなければなりません。Python3で利用可能なフォークを準備しましたが、テストは十分にできていません。ですが、ここでの目的には十分機能するように見えます。
次を実行してインストールすることができます：
!pip install https://github.com/jaimergp/autodocktools-prepare-py3k/archive/master.zip
ドッキング計算の手順は次のようになっています。
- AutoDockToolsでタンパク質とリガンドを準備します。（ローカル環境で実施）
- Proteins.plus' DoGSiteScorerを使って結合ポケットの可能性があるもののうち最も良いものを見つけます。この情報をVinaによる計算の設定を行うために使います。（幾何学的な中心と探索空間のサイズ）
- OPAL上でVinaの計算を実行します。
import time import os from io import StringIO
構造の準備
構造の準備はAutoDockTooksライブラリの適切な部分を実行するだけでできます。
MolKitで構造のファイルを読み込む- 準備のための正しい機能を適用する。タンパク質には
AD4ReceptorPreparationを、リガンドにはAD4LigandPreparationを適用します。
準備自体は次のような項目に注意して行います。
- タンパク質とリガンドに水素原子を付加する。
- タンパク質に属さないおかしな残基を取り除く。
- 原子タイプと部分電荷を割り当てる。
- リガンドの回転可能（torsionable）な分岐をみつける。
結果はPDBQTファイルとしてディスクに書き込まれます。
###################### # # 構造の準備 # ###################### import MolKit from AutoDockTools.MoleculePreparation import AD4ReceptorPreparation, AD4LigandPreparation def opal_prepare_protein(protein): """ AutoDockを使うにはPDBQTファイルが必要です。 """ mol = MolKit.Read(protein)[0] mol.buildBondsByDistance() RPO = AD4ReceptorPreparation(mol, outputfilename=protein+'.pdbqt') return protein + '.pdbqt' def opal_prepare_ligand(ligand): """ AutoDockを使うにはPDBQTファイルが必要です。 """ mol = MolKit.Read(ligand)[0] mol.buildBondsByDistance() RPO = AD4LigandPreparation(mol, outputfilename=ligand+'.pdbqt') return ligand + '.pdbqt'
結合サイトの予測
Proteins.plusのDoGSiteScorerは、PDBデータベースにあるタンパク質の処理が必要なだけならば、とても簡単に利用できるREST APIを提供しています（dogsite_scorer_submit_with_pdbidを参照してください）。ですが、我々はKLIFSから取得した構造を使っているので、公式のPDBに登録されている構造の位置と向きが、KLIFSのものと一致しているか保証することができません。両者を重ね合わせて得られた変換行列を、取得したポケットに適用することもできますが、標準Webインターフェースで提供されている方法と同様に、単純にPDBファイルをProteins.plusにアップロードするだけの方がより簡単でしょう。
しかし、REST APIはそのようなオプションを提供していないので、我々はこの処理をリバースエンジニアリングする必要があります。適切なHTTPリクエストの場所をみつけるためには、Chromデベロッパーツールのネットワークタブを開き、通常のWebサイトを利用した際にどのHTTPリクエストが実行されているか書き留めておく必要があります。ここで有効なHTTPリクエストを手に入れるためのポイントは、信頼性トークン(Authenticity token)とHTTPヘッダーです。
この調査の結果は統合して整理しdogsite_scorer_submit_with_custom_pdb 関数に反映させています。技術的な詳細について興味があれば、関数のコメントを読んでください。
このアプローチを行うには、結合サイトを予測するためにHTMLコードを構文解析するため、
BeautifulSoupも必要となります。
!pip install BeautifulSoup4
この方法を使えば、最も可能性の高い結合ポケットの幾何学的な中心とそのサイズを取得することができます。Vinaの計算を設定するにはどちらの値も必要となります。
###################### # # 結合ポケット予測 # ###################### from bs4 import BeautifulSoup import requests def dogsite_scorer_submit_with_pdbid(pdbid, chain): """ これは公式APIですが、PDBコードしか使えません。 パラメーター ---------------- pdbid : str ４文字のPDB ID chain : str 解析するタンパク質の鎖（chain） 戻り値 ---------------- str ジョブが実行中か完了したかについて更新情報を受け取るために用いるURL """ # Proteins.plusにジョブを投げる r = requests.post("https://proteins.plus/api/dogsite_rest", json={ "dogsite": { "pdbCode": pdbid, "analysisDetail": "1", "bindingSitePredictionGranularity": "1", "ligand": "", "chain": chain } }, headers= {'Content-type': 'application/json', 'Accept': 'application/json'} ) r.raise_for_status() # 計算の更新情報を取得するためには"location"を問い合わせる必要があります。 return r.json()['location'] def dogsite_scorer_submit_with_custom_pdb(pdbfile): """ カスタムしたPDBをアップロードするために、実際のHTMLフロントエンドを模倣する必要があります。 1. HTMLメタヘッダーからCSRF（クロスサイトリクエストフォージェリ）トークンを取得します。 2. アップロードするためにファイルをPOSTします。 3. 返ってくるHTMLページはURL IDを含んでおり、次に内部の省略されたジョブIDを取得するために利用できます。 （Webインターフェースを使っている時と同様に）Webサーバーがフロントエンドで実行する非同期の呼び出しを真似することで、 返ってきたHTMLページのURL IDを使ってパブリックなジョブのAPI IDを取得することができます。 4. パブリックジョブIDを手に入れたらREST APIの利用に切り替えることができます。 """ # プロセスの間、途中のcookieを保持するために`session`を使う必要があります。 session = requests.Session() r = session.get("https://proteins.plus/") r.raise_for_status() # ホームページには我々のリクエストを有効にするのに必要なCSRFトークンが含まれています。 # そうでないと安全ではありません！リクエストを行なっている間、これを使わなければいけないので、 # session HTTPヘッダーの一部として設定しておくのが最善の方法です。 html = BeautifulSoup(r.text) token = html.find('input', {'name': 'authenticity_token'}).attrs['value'] session.headers['X-CSRF-Token'] = token # 1. ファイルをアップロード with open(pdbfile, 'rb') as f: r = session.post("https://proteins.plus", files={'pdb_file[pathvar]': f}) r.raise_for_status() # REST APIがファイルのアップロードをサポートしているのであれば、すでにパブリックIDを手に入れていることになりますが、 # それまでは応急処置としてこの方法で済ませておく必要があります。 # 2. 内部のlocation IDの取得 html = BeautifulSoup(r.text) pdb_id = html.find('input', {'name': "dogsite[pdbCode]"}).attrs['value'] # 3. 内部のジョブIDの取得 session.headers['Referer'] = "https://proteins.plus" + pdb_id r = session.post(f"https://proteins.plus/{pdb_id}/dogsites", json={"dogsite": {"pdbCode": pdb_id}}, headers= {'Content-type': 'application/json', 'Accept': 'application/json'} ) r.raise_for_status() job_id = r.json()['job_id'] time.sleep(3) # サーバーがリクエストを処理できるように、処理を続ける前に少し待ちます。 # 4. パブリックなジョブIDの取得 while True: r = session.get(f"https://proteins.plus/{pdb_id}/dogsites/{job_id}?_={round(time.time())}", headers= { 'Accept': 'application/json, text/javascript, */*', 'Sec-Fetch-Mode': 'cors', 'Sec-Fetch-Site': 'same-origin', # どうもこの下の１行が大事なようです。 # これ無しではエラー406が投げられます。 'X-Requested-With': 'XMLHttpRequest'} ) r.raise_for_status() if 'Calculation in progress...' in r.text: # まだ完了していません。 time.sleep(5) continue if 'Error during DogSiteScorer calculation' in r.text: # 形式のおかしなファイル？ raise ValueError('Could not run DoGSiteScorer!') break results_id = None for lines in r.text.splitlines(): for line in lines.split('\\n'): if 'results/dogsite' in line: results_id = line.split('/')[3] break if results_id is None: raise ValueError(r.text) return f"https://proteins.plus/api/dogsite_rest/{results_id}" def dogsite_scorer_guess_binding_site(protein): """ タンパク質の最も可能性が高い結合サイトを取得するためにProteins.plusのDoGSiteScorerを使います。 """ if len(protein) == 4: # PDBコード job_location = dogsite_scorer_submit_with_pdbid(protein) elif protein.endswith('.pdb'): job_location = dogsite_scorer_submit_with_custom_pdb(protein) else: raise ValueError("`protein` must be a PDB ID or a path to a .pdb file!") # 計算が完了したか確認します。 while True: result = requests.get(job_location) result.raise_for_status() # 失敗した場合、ここで止まります。 if result.status_code == 202: # まだ実行中です。 time.sleep(5) continue break # residuesファイルは幾何学的中心と半径をPDBファイルのコメントとして保持しています。 # 最初のファイル（residues[0]）が最もスコアの良いポケットです。 pdb_residues = requests.get(result.json()['residues'][0]).text for line in pdb_residues.splitlines(): line = line.strip() if line.startswith('HEADER') and 'Geometric pocket center at' in line: fields = line.split() center = [float(x) for x in fields[5:8]] radius = float(fields[-1]) break return center, radius # これがVinaの計算に必要となるものです。
OPAL上でVinaを実行
（１）タンパク質とリガンドを準備し、（２）探索空間の推定が完了したら、OPAL Webサーバーに実際の計算をサブミットできます。以下の処理を行います。
sudsを使ってSOAPクライアントを初期化する- SOAP XMLリクエストに加えられるように、ファイルを
base64文字列 としてエンコードする - ジョブリクエストをサブミットし、ジョブIDを受け取る
- 計算が完了したかどうか確認するため、サーバーにジョブIDを使って問い合わせる
requestsで関連するファイルをダウンロードする
上記の手順はOPAL Webサイトにはあまり記述されておらず（実際、ドキュメンテーションは利用できなくなってきます)）、これらのサーバーに依存しているUCSF Chimeraモジュールのソースコードの中身を調査することでわかりました。
通常、計算には５−１５分かかるので、４番目のステップのためにJupyterの便利な機能を使ってリアルタイムに出力ファイルの中身を更新します（iprint()関数を参照してください）。こうすることで、計算がうまく進んでいるかどうか盲目的に信じる代わり進行度合いを確認できます！
###################### # # 計算の実行 # ###################### from suds.client import Client from IPython.display import display, clear_output, HTML from rdkit import Chem from rdkit.Chem import AllChem VINA_CONFIG = """ center_x = {center[0]} center_y = {center[1]} center_z = {center[2]} size_x = {size[0]} size_y = {size[1]} size_z = {size[2]} """ def opal_run_docking(protein, ligand, center, size, stream_output=True): """ OPAL Webサービスに接続しジョブをサブミットする """ client = Client("http://nbcr-222.ucsd.edu/opal2/services/vina_1.1.2?wsdl") files = 'receptor.pdbqt', 'ligand.pdbqt', 'vina.conf' with open(protein) as f: protein_contents = f.read() with open(ligand) as f: ligand_contents = f.read() file_map = [ {'name': 'receptor.pdbqt', 'contents': base64ify(protein_contents)}, {'name': 'ligand.pdbqt', 'contents': base64ify(ligand_contents)}, {'name': 'vina.conf', 'contents': base64ify(VINA_CONFIG.format(center=center, size=size))}, {'name': 'results.pdbqt', 'contents': ''}, ] cli_args = "--receptor receptor.pdbqt --ligand ligand.pdbqt --config vina.conf --out results.pdbqt" response = client.service.launchJob(cli_args, inputFile=file_map) job_id = response.jobID url = f"http://nbcr-222.ucsd.edu/opal-jobs/{job_id}" message = "Waiting for job " + url while True: r = requests.get(url + "/vina.out") try: r.raise_for_status() except: # 最初のチェック段階では出力ファイルはまだ出来ていないかもしれません。 iprint(message) else: iprint(f"{message}\n{r.text}") if client.service.queryStatus(job_id).code == 2: time.sleep(10) continue print('\nFinished!') break output_response = client.service.getOutputs(job_id) output_files = { 'stdout.txt': requests.get(output_response.stdOut).text, 'stderr.txt': requests.get(output_response.stdErr).text, } for f in output_response.outputFile: if f.name in files: continue r = requests.get(f.url) r.encoding = 'utf-8' r.raise_for_status() contents = r.text output_files[f.name] = contents time.sleep(0.1) return output_files ###################### # # ユーティリティー # ###################### import base64 def base64ify(bytes_or_str): """ Py2k base64encodeの動きを真似する """ if isinstance(bytes_or_str, str): input_bytes = bytes_or_str.encode('utf8') else: input_bytes = bytes_or_str output_bytes = base64.urlsafe_b64encode(input_bytes) return output_bytes.decode('ascii') def iprint(s): """ この関数を使って出力をプリントし、その前のものを上書きすることで、 継続的に更新されているようにみえるようにすることができます:) """ clear_output(wait=True) s = s.replace("\n", "<br />") display(HTML(f'<pre>{s}</pre>'))
全部まとめる
必要な関数をすべて定義し終わったので、パイプラインを作成することができます。
def step_03_opal(protein, smiles, pdbcomplex): """ タンパク質の構造とSMILES文字列のリストを与え、以下の手順を実施します・ ステップ： 1. AutoDock Vinaのためのタンパク質を準備（ローカルで実行） 2. DoGSiteScorerを使って最も可能性の高い結合サイトを見つける 3. 各リガンドについてRDKitを使って3DのPDBファイルを書き出し、 OPALでAutoDockを実行する すべて行うのに大体5-15分かかります。 結果は出力ファイルを含む辞書です。我々が主に興味があるのは結果の result['results.pdbqt']です。 """ prepared_protein = opal_prepare_protein(protein) center, radius = dogsite_scorer_guess_binding_site(pdbcomplex) size = [radius] * 3 # Vinaは立方体のボックス以外に対応していますが、単純のため立方体を使います for i, smile in enumerate(smiles): smiles_to_pdb(smile, f'data/ligand{i}.pdb') prepared_ligand = opal_prepare_ligand(f'data/ligand{i}.pdb') result = opal_run_docking(prepared_protein, prepared_ligand, center, size) return result
実行します！マジックコマンド%timeでどれくらい時間がかかったか測ります。
# `smiles`の最初のリガンドだけを処理します。 %time result = step_03_opal(PROTEIN, smiles[:1], COMPLEX)
訳注(2020/05)
上記のセルはエラーを返しました。OPALがサービスを停止しているのか現在実行できなくなっています。
NBCRと関係するWebサイト、Webサービスが2020年4月30日をもって終了したようです（NBCRリンク ）。
訳注ここまで
出力結果を理解する
resultは辞書で、出力ファイルとその中身のテキストに相当するいくつかのkeyをもちます。我々が主に興味があるのは次です。
- ドッキングされたリガンドを含む
results.pdbqt。複数のモデルをふくむ修正されたPDBファイルです。タンパク質のグリッドを保持しているので、もともとのタンパク質の構造と一緒に各リガンドのモデルを開くだけで大丈夫です。 vina.outは上記で見たテキストの出力です。テーブルのような情報を見ることができます。
あとで取り出せるようにディスクに保存しておきます。
with open('data/results.pdbqt', 'w') as f:
f.write(result['results.pdbqt'])
with open('data/vina.out', 'w') as f:
f.write(result['vina.out'])
ドッキング結果の可視化
計算を実行しファイルをダウンロードしたら、可視化しましょう！やり方はパートCを見てください。
ディスカッション
OPALはドッキング計算を行うVinaを無料で提供していますが、サブミットする前に入力ファイルをローカル環境で準備する必要があります。準備の中には、探索空間の定義も含まれており、通常既知の結合サイトの周囲を指定します。視覚的に推定する代わりに、Proteins.plusのDoGSiteScorerサーバーを使って、最も可能性が高い結合サイトの中心と半径を計算しました。これら２つのサービスは異なるコミュニケーションインターフェースを使っていました。
Protein.PlusはRESTを使っていますが、KLIFSとは異なりswagger.jsonの定義を提供していないので、自分のリクエストを手動で構築しなければなりませんでした。幸い、Webサービスが単純だったので２、３個のリクエストが必要なだけ、、、と思っていました！！現在のAPIはカスタマイズしたPDBのアップロードを許可していないので、通常のブラウザーを使っているふりをするために多少のGET、POSTリクエストを使う必要がありました。正しいリクエストを推測するうえで、Chromデベロッパーツールのネットワークタブはとても便利でした。Webサービスのリバースエンジニアリングが必要になったら、便利なツールの一つです！
OPALはXMLベースの標準SOAPで構築されており、JSONベースのRESTよりもちょと扱いにくですが、sudsを使うことでとても簡単になりました。sudsは利用できるメソッドを教えてくれますが、これらはあまりきちんと説明されていません。幸い、UCSF Chimeraのモジュールにコードの例がいくつか散らばっていました。ジョブIDを手に入れたらOPALにサブミットされたジョブはパブリックなものとなり、ファイルはリアルタイムに更新されるので、サーバーにN秒毎に問い合わせを行い、動的にdisplay()関数で表示されるHTMLを更新することで、計算の出力をライブでプレビューすることができました。この技は自分のマシンで行なっている計算を問い合わせるのにも使うことができます。ファイルの中身を読み、現在の出力を消して新しい出力をHTML()オブジェクトとして表示するだけです（ipirnt()関数を参照してください）。
クイズ
- リモートサーバーでドッキングが成功したかどうかどうすればわかりますか？
- なぜAutoDock Vinaの入力ファイルをローカル環境で準備する必要がありましたか？
- 応用：MCule'sドッキングサーバーをJupyter Notebookから使うための正しいHTTPリクエストを見つけることはできますか？
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
残念ながらノートブックの実行はできませんでした。オンラインサービスを使うとなるとサービスが停止したり安定しないといった問題もあるんですね。
パートBの内容はアプリケーションをプログラム上で扱うためのコードばかりで私は全くついていけませんでした。文章から理解する能力がないので絵で説明してほしいです・・・。
