KLIFSで遊んだ話（T11パートAの備忘録）
トークトリアル１１イントロに引き続きトークトリアル１１パートAも耳慣れない言葉が多く内容のフォローが難しかったです。
と、いうわけで今回も素人による適当解説をしますよ！！
KLIFS
今回はKLIFS（Kinase-Ligand Interaction Fingerprints and Structures）について調べてみます。トークトリアルの説明は文字のみで、データの取得に用いただけだったのですが、実際にWebページ にアクセスしてみるととても素敵な感じでした！
KLIFS Webページで遊ぶ
折角なのでT11パートAでコードで実施していた内容をWebページ上でも行なってみます。T11 パートAの内容は以下でした。
- Kinaseファミリー 「EGFR」 を検索
- PDB ID 3w33 を選択
- 描画
まずは検索から。上タブからSearchページに移動します。タンパク質やリガンドだけでなく色々な検索方法が用意されています。興味深いところでは結合ポケットのなかの水分子や、リガンド-キナーゼ相互作用パターンによる検索といったものまであります。目移りしてしまいますが、まずはT11通りにKinaseファミリーで検索します。
Search by classificationにて「Family：EGFR」「Organism：Human」として検索しました。
ヒットした構造（355エントリー）がテーブルとして表示されました。下図に化合物２次元構造がありますが、Orthosteric ligandの列の化合物の名称にマウスオーバーするだけで描画されます。レスポンスも早くてすごい！
またテーブル左上Plot results on Kinomeをクリックするとキノーム（ゲノム上のプロテインキナーゼ一式）の中で、どの位置に相当するか図示されます。
検索結果のテーブルはExcel、csv形式でダウンロードすることも可能で、またDownload Structuresで構造をダウンロードすることもできるようです（スクロールしていくとボタンがあります。）
「PDB ID: 3W33」の詳細を見るためにDetails列のShowをクリックして見ます。こんな感じ
エントリーの詳細な情報に加えて３Dの描画やリガンドの活性値情報もあります。ChEMBLから取得されているようで、EGFR以外のキナーゼに対する活性も合わせてみることができるので、選択性やオフターゲットの雰囲気もつかめそうです。
スクロールするともっと面白い項目もあります。
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の１文字表記。側鎖の動きが大きくDFG-in構造、DFG-out構造とよばれるコンフォメーションの変化がある
などなど。また、阻害剤にも結合するサイト・阻害様式によってType I ~ IV、共有結合型などの分類があるそうです。
こういうの格好いいですよね！共通の情報が整理されて見通し良くなってく感じが分野の発展の歴史を見ているようで楽しいです。
まとめると「医薬品としてキナーゼ阻害剤は大事だが選択性が課題となる。背景には結合サイトにおける保存されたキナーゼの構造的特徴があり、特徴に従い分類、キナーゼー阻害剤相互作用の細かな差異を分析することで課題解決の有用な情報となる。」・・・ってことですかね？
つまり、そういう情報まとめたKLIFSは素敵なリソース！！
KLIFS文献参照
やっとKLIFSの文献に戻ってFigure 2. を引用します。こちらはKLIFSがどのように複合体のアノテーションを行なっているかを説明した図です。
先に見た単語がちらほら見えますね。この図で面白いのはキナーゼーリガンド相互作用をフィンガープリントして表現しているFig.2 A)下部だと思います。各アミノ酸残基について、疎水性相互作用や水素結合といった7種類の相互作用がリガンドとの間にみられるか？というのを、0、１のビットで表記しています（kinase-ligand interaction finger print, IFP）。
KLIFSは触媒ドメインを含む85残基を解析の対象としています。したがって、各複合体のキナーゼーリガンド相互作用の情報は「85残基x相互作用7種=595ビット」のフィンガープリントにエンコーディングされることになります（？あってます？）。これだけ情報が圧縮されてるから「類似の相互作用の複合体は他にあるのか？」といった検索が容易になるんですね。なるほど。
さて、KLIFSでは相互作用のフィンガープリント化と結合ポケットの部位の分類（3 major pockets + 12 subpockets）を行なったことで、各複合体のリガンドの結合ポーズ（binding mode）の分類が可能になったそうです。
また、興味深いことにリガンドに加えて水分子についても解析が行われています。下に引用した図のように、保存された水分子のクラスター位置が見つかっています。こちらもWebページ上で確認、検索することが可能になっています。
以上、非常にざっくりですがKLIFSで行われていることの雰囲気です。分類と解析とが体系化されることで自動化が可能になり、データベースの継続的なアップデートに繋がっているようです。
KLIFSのAPI
KLIFSのWebページについてはだいたいわかったので、T11パートAではAPIを使ってどのように情報が取り出されていたのかを見ていきたいと思います。
取り上げるのは途中で定義された以下の３つの関数です。
_all_kinase_families(): KLIFSデータベースに含まれているキナーゼファミリーの名前を取得_kinases_from_family(family, species="HUMAN"): 指定したキナーゼファミリーに含まれるキナーゼの名前を取得_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 をみてみます。
たとえば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種類と同じものがきちんと入っているようです。
関数_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を引数として与え、以下の３つを返す関数です。、
- キナーゼIDに基づく複合体構造（PDB形式）
- タンパク質のみの構造（mol2形式）
- キナーゼIDに基づく共結晶構造中のリガンドのリスト（SMILESのみを抜き出して出力）
1、２は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
２つ結果が得られました。今回はhumanの情報をとりたいので一つ目のIDを取り出します。
EGFR_ID = EGFR_id_result[0].kinase_ID print(EGFR_ID) # 406
では_protein_and_ligand_structure()関数の１行目、複合体の構造のリストの取得です。ここでは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)
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経由で２次元に戻してから描画 m782 = Chem.MolFromMol2Block(ligand_782) smi782 = Chem.MolToSmiles(m782) m782_2d = Chem.MolFromSmiles(smi782) Draw.MolToImage(m782_2d)
先に確認したように、「structure_id: 782」のPDB_IDは3W33で、KLIFSのWebページで見ていた複合体構造と同じでした。確かにここでは同じリガンドが取り出されているようです。
まとめ
以上、今回はT11 パートAで扱われていたデータベースKLIFSで遊んで見ました。素晴らしいWebページ版が提供されているので、まずはこちらで遊んでどのような情報が得られるのか感触をつかんでからの方がプログラムの内容を理解しやすのではないかな？と思いました。一目で結果を理解しやすいGUIベースのページは、プログラムが意図した通りの動きをしているかを検証する目的でも良さそうです。
プログラムの内容についてREST、SWAGGERは相変わらず定義はよく分かりませんが、ぼんやりと何がしたいかはわかったような気がしないでもないような？仕様を定義ってのは、文法みたいなイメージで良いのでしょうか？？？
素人が適当にしらべたことをだらだら書いているので色々と間違いが多いと思います。ご指摘いただければ幸いです。
*1: キナーゼを標的とした構造生物学および創薬の現状 日本結晶学会誌 59，174-181(2017)
TeachOpenCADD トピック１１（パートA）〜オンラインAPI/サービスを使った構造に基づくCADD〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル１１をもととしておりCC BY 4.0のライセンスに従っています。
トークトリアル１１は論文発表時のTeachOpenCADDには含まれておらず新しく追加されたものです。
４つのノートブックに分かれており、ここではその2つ目を訳出します。
トークトリアル１１で取り上げられている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で類似の化合物を検索します。
学習の目標
理論
- パイプラインと含まれているWebサービスの解説
実践
- パイプラインを構成する関数を書く
- ケーススタディー：EGFR阻害剤
レファレンス
- KLIFS：キナーゼー阻害剤相互作用データベース
- データベースそのものの解説 (J. Med. Chem. (2014), 57, 2, 249-277)
- オンラインサービスの解説 (Nucleic Acids Res. (2016), 44, 6, D365–D371)
- PubChem：低分子化合物のデータベース(Nucleic Acids Res. (2019), 47, D1102-1109)
- NGLView：Notebookのためのインタラクティブな可視化ソフト(Bioinformatics (2018), 34, 1241–124)
理論
ここまでで、さまざまなPythonライブラリをつかってオンラインWebサービスにクエリを投げる方法をみてきました。オンラインWebサービスからパイプラインを構築してみましょう！
パイプラインと含まれているWebサービスの解説
- キナーゼーリガンド相互作用フィンガープリントと構造のデータベース（KLIFS: Kinase-Ligand Interaction Fingerprints and Structures database）はアムステルダムのVU大学医薬品化学分野で開発されたデータベースで、触媒キナーゼドメインのタンパク質構造に関する情報（PDBから収集）と、リガンドとの相互作用に関する情報とを提供しています。このデータベースから精選されたタンパク質構造を取得し、そのリガンドの情報をつかってPubChemやChEMBLといった他のデータベースから類似したリガンドを取得することができます。
- KLIFによるリガンド情報を使って、PubChemで類似化合物を探します。
- タンパク質構造と候補となるリガンドをいくつか手に入れたら、OPALWebサービスで提供されている組み込みのVinaを使ってオンラインでドッキングを実施します。また、proteins.plus' DoGSiteScorerで化合物をドッキングする結合サイトの候補を探します。（パート B）
- 結果を
nglviewで可視化し、相互作用を PLIPを使ってレポートします。（パート C）
このトークトリアル パート(A)では、KLIFSデータベースから入力として用いるキナーゼー阻害剤構造を取得し、候補阻害剤として評価することが可能な類似化合物をPubChemで検索します。関連する出力はディスクのdata/に書き込み、あとで使います。
KLIFS
- 役割： キナーゼーリガンド相互作用プロファイルデータベース
- Webサイト： http://klifs.vu-compmedchem.nl/
- API： RESTベース、Swaggerで利用可能なものがあります。正式なクライアントはありません。
bravadoを使ってください。 - ドキュメンテーション: http://klifs.vu-compmedchem.nl/swagger/
- 文献：
- データベースそのものの解説 (J. Med. Chem. (2014), 57, 2, 249-277)
- オンラインサービスの解説 (Nucleic Acids Res. (2016), 44, 6, D365–D371)
キナーゼーリガンド相互作用フィンガープリントと構造のデータベース（KLIFS: Kinase-Ligand Interaction Fingerprints and Structures database）はアムステルダムのVU大学医薬品化学分野で開発されたデータベースで、触媒キナーゼドメインのタンパク質構造と、キナーゼ阻害剤の阻害様式に焦点をあてています。一貫したシステマティックなプロトコルに基づき、（現在、ヒトとマウスの）全てのキナーゼ構造とキナーゼリガンドの結合様式を直接互いに比較することができます。さらに、結合サイト全体を取り囲む85残基の分類がなされているので、例えば、キナーゼー阻害剤選択性を決める重要な相互作用を見つけるといった目的のため、キナーゼー阻害剤の相互作用パターンを互いに比較することができます。
PubChem
- 役割： 低分子化合物データベース
- Webサイト： https://pubchem.ncbi.nlm.nih.gov/
- API: RESTベースのものが利用可能です。正式なクライアントはありません。
bravadoを使ってください。 - ドキュメンテーション： https://pubchemdocs.ncbi.nlm.nih.gov/pug-rest
- 文献：
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()
preview_smiles(ligands[0]) # Clc1ccc(NC(=O)c2c(NCc3ccncc3)cccc2)cc1
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上の最も類似した１０個の化合物を探します。 return similar_compounds_pubchem(smiles, n_records=10)
取得した化合物をRDKitで描画します。
similar_smiles = query_pubchem_for_similar_compounds(ligands) multi_preview_smiles(*similar_smiles)
ケーススタディー: EGFR
対応するWikipediaの記事から改変して引用します。
上皮成長因子受容体（Epidermal Growth Factor Receptor; EGFR）はErbBファミリー受容体のメンバーで、密接に関係する４つの受容体チロシンキナーゼ、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()
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
(訳注 ADPの誤記だと思います。)
類似したリガンドを試しに見てみましょう。
similar_smiles_egfr = query_pubchem_for_similar_compounds(egfr_ligands) multi_preview_smiles(*similar_smiles_egfr)
トークトリアルの次のパートのためにディスクに結果を書き出しておきましょう！
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 イントロより）
トークトリアル１１はオンラインAPI/サービスを利用しています。T1~10と比べてCADD以前にプログラミングの部分の敷居が高いという印象です。・・・正直コードが何をしているかプログラミングのできない私には分からなかったです（レベル低くてすみません）。
と、いうわけで素人による適当解説をしますよ！！WEBスクレイパーに俺はなる！！
題材
T11の初っ端、イントロから敷居の高い単語が連発でしたが、一番気になった「WkipediaProteinogenic amino acidからのアミノ酸のテーブルを取得する部分」を見ていきます。
行われていることは以下の通りです。
| 手順 | ライブラリ | やってること |
|---|---|---|
| データの取得 | requests | Wikipediaの該当ページにアクセス HTMLデータを取得 |
| データの抽出 | BeautifulSoup | HTMLの文字列から 目的のテーブルの場所を指定して 情報を抜き出す |
| テーブルの再構成 | Pandas | 抜き出した情報をDataFrameで テーブルに再構成して表示 |
では順番に見ていきます。。*1
requestsによるHTMLデータの取得
「入門 Python3 9章」によると、WWWの単純な形としてウェブクライアントが
- ウェブサーバーにHTTP(Hypertext Transfer Protocl)で接続し
- サーバーにURL（Uniform Resource Locator）を要求し、
- HTML(Hypertext Markup Language)を受け取る
という流れになっているそうです。で、Pythonはこれが得意なんだそうな。
クライアントがサーバーに要求(リクエスト)を送って応答（レスポンス）を受け取るウェブのデータの交換の標準プロトコルがHTTPで、そのHTTPメソッドにGETやPOSTというものがあるそうです。
では、Pythonのrequestsを使って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の情報をちゃんと取得できているようです。
BeautifulSoupによる構文解析
取得したHTMLの情報から該当の場所を探し出します。といってもHTMLの書式が分からない私にはさっぱりです。でも大丈夫! ......そう、「Google Chrome」ならね！ *2
「右クリック → 検証」だ！
HTMLをただ表示するだけではありません。HTML上でマウスオーバーするだけで、通常の表示ページでの該当箇所がハイライトされて対応を確認することができます。
ここを頼りに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 > | テーブルのデータセルを作成 |
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要素は見出し以降の５番目（index 4）の位置に相当します。
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
２２の行がdataに格納されました。アミノ酸って20個じゃなかったっけ？と思いましたが、テーブルを見直すとSelenocysteine、Pyrrolysineが含まれていました。こんなアミノ酸知らなかった。。。
あとはデータフレームにするだけです。
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まであるのですが先が思いやられますね！
色々間違っていそうなのでご指摘いただければ幸いです。
ではでは
*1: 参考にさせていただいた記事
図解！PythonでWEB スクレイピングを極めろ！(サンプルコード付きチュートリアル)
図解！Python BeautifulSoupの使い方を徹底解説！(select、find、find_all、インストール、スクレイピングなど)
図解！PythonのRequestsを徹底解説！(インストール・使い方)
Pythonでスクレイピングをする最初の一歩、Webページを丸ごと取得する方法
PythonでWebページを取得できたかどうかのエラーチェックと安全な中止の仕方
*2: とてもわかりやすい記事 Python Webスクレイピング テクニック集「取得できない値は無い」JavaScript対応
*4: テーブルの中では書式の問題で<>とタグにスペースを挟んでいますが、実際は不要です。念のため
*5:初めからtableタグで指定すればよかったのでは？とも思いましたが、同じページにある複数のテーブルでtableタグとその属性が同じになっていたので指定しづらいから使わなかったのだと思います。多分
TeachOpenCADD トピック１１（イントロ）〜オンラインAPI/サービスを使った構造に基づくCADD〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル１１をもととしておりCC BY 4.0のライセンスに従っています。
トークトリアル１１は論文発表時のTeachOpenCADDには含まれておらず新しく追加されたものです。
４つのノートブックに分かれており、ここではその一つ目を訳出します。
トークトリアル１１で取り上げられているWebサーバー、サービスには既に停止しているものもあり、そのまま使えるものではなくなってしまっていますが、基本的な考え方や手法を知る上では役にたつかもしれません。
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下、訳。
トークトリアル 11 （イントロ）
オンラインAPI/サービスを使った構造に基づくCADD
Developed at AG Volkamer, Charité
Dr. Jaime Rodríguez-Guerra, Dominique Sydow
このトークトリアルの目的
Webサービスは、ユーザーのインストールの煩わしさを避けられるので、ソフトウェア利用の上で便利な方法です。Web UIは、通常ワークフローの自動化の可能性を犠牲にして、簡単に利用できるようにします。幸運にも、アクセスを自動化するためのAPI（アプリケーションプログラミングインターフェース / Application Programming Interface）を提供するWebサービスの数は増えています。コンピューター支援医薬品デザインの分野における例としては次のようなものがあります：
- ChEMBL
- RCSB PDB
- KLIFS
- Proteins.plus
- SwissDock
以下のノートブックでは、医薬品デザインの文脈において、オンラインWebサービスをPythonからプラグラム的に使う方法を学びます。最終的な目標は、（ほとんど）すべてのローカル環境での実行なしに、もっぱらWebサービスのみに依るパイプライン全体を構築することです！
注: 簡単のため、このトークトリアルの実践編はイントロダクションと以下の追加の３つのノートブックに分けていてます。
- 11a. キナーゼ阻害剤の候補をKLIFとPubChemで検索
- 11b. 11aで取得した候補化合物をターゲットタンパク質に対してドッキング
- 11c. 結果を評価し既知のデータと比較
学習の目標
理論
Pythonを使ってWebサービスにプログラムでアクセスする方法の種類：
理論
プログラムによるアクセスの種類
プログラムのAPI
現代のWebサービスは、特にデータベースの場合、データへのアクセスのために標準化された方法を提供することができます。通常、これはサイトの特定のURLにアクセスし、機械が読むことができる形式の結果を要求することができることを意味します。
例えば、UniProt はタンパク質に関するあらゆる種類の情報のためのデータベースです。もし（癌に関与している）Srcのような特定のタンパク質を探しているのであれば、セクションごとにコンテンツがよく組織化されたこの美しいWebページにたどり着くでしょう。ですが、もしURLに.fastを追加したならば、タンパク質の配列をFASTA形式で取得することができます。
https://www.uniprot.org/uniprot/P12931 -> https://www.uniprot.org/uniprot/P12931.fasta --- >sp|P12931|SRC_HUMAN Proto-oncogene tyrosine-protein kinase Src OS=Homo sapiens OX=9606 GN=SRC PE=1 SV=3 MGSNKSKPKDASQRRRSLEPAENVHGAGGGAFPASQTPSKPASADGHRGPSAAFAPAAAE PKLFGGFNSSDTVTSPQRAGPLAGGVTTFVALYDYESRTETDLSFKKGERLQIVNNTEGD WWLAHSLSTGQTGYIPSNYVAPSDSIQAEEWYFGKITRRESERLLLNAENPRGTFLVRES ETTKGAYCLSVSDFDNAKGLNVKHYKIRKLDSGGFYITSRTQFNSLQQLVAYYSKHADGL CHRLTTVCPTSKPQTQGLAKDAWEIPRESLRLEVKLGQGCFGEVWMGTWNGTTRVAIKTL KPGTMSPEAFLQEAQVMKKLRHEKLVQLYAVVSEEPIYIVTEYMSKGSLLDFLKGETGKY LRLPQLVDMAAQIASGMAYVERMNYVHRDLRAANILVGENLVCKVADFGLARLIEDNEYT ARQGAKFPIKWTAPEAALYGRFTIKSDVWSFGILLTELTTKGRVPYPGMVNREVLDQVER GYRMPCPPECPESLHDLMCQCWRKEPEERPTFEYLQAFLEDYFTSTEPQYQPGENL
これは特定のURLスキームを介して、プログラムによりWebサービスにアクセスすることを提供する方法です。ですが、各Webサービスは各々独自の仕組みを提案しており、開発者はケースバイケースでスクリプトを実装せざるを得ません。
幸運にも、この種のプログラムによるアクセスを提供する標準化された方法がいくつかあります：
このトークトリアルでは、RESTとSOAP APIだけを扱います。
HTTPベースのRESTful API
この種のプログラムアクセス/規定は、プログラムによるアクセスを要求するクライアント（スクリプト、ライブラリ、プログラム）に対して特定のエントリーポイントを定義します。api.webservice.comといったようなものです。バージョン管理することができるので、サービス提供者はすでに存在している実装を中断せずにスキームを更新することができます（api.webservice.com/v2 をデプロイした時でさえでも、まだapi.webservice.com/v1が機能し続けます）
この種のAPIには通常、プラットフォームで利用可能なアクションを全て解説したうまくまとめられたドキュメンテーションがあります。例えば、GitHub REST API for listing repositoriesをみてください。各引数とオプションが、利用例とともにどのようにドキュメント化されているか、見てみることができるでしょう。
唯一の難しい点は、必要なURLを正しい方法で構築することです：
https://api.github.com/users/volkamerlab/repos
訳注(2020/05)
ノートブックではここにURLの戻り値が書かれていますが、冗長で内容とは関係しないので省略します。GitHubのノートブックを確認してください。
訳注ここまで
これはたまたまJSONフォーマットの辞書でした！jsonライブラリを使って構文解析しPythonの辞書に変換することが簡単にできます。さらに良いことには、それをする必要すらありません。requestsを使うことで、以下の操作が３行でできます。
import requests response = requests.get("https://api.github.com/users/volkamerlab/repos") response.raise_for_status() # この行はエラーが生じる可能性をチェックします。 result = response.json() result
訳注(2020/05)
ノートブックではrequestsで得た結果が出力されていますが、冗長なので省略します。先にURLを使った場合と同じ結果が得られることが確認されています。
訳注ここまで
URLをf文字列(f-string)でパラメーター化すれば、どのユーザーのレポジトリでもリストするように関数を書くことができます。
def repos(username): response = requests.get(f"https://api.github.com/users/{username}/repos") response.raise_for_status() return response.json() repos("volkamerlab")
訳注(2020/05)
ノートブックではここに関数の出力が書かれていますが省略します。
訳注ここまで
APIのクライアントを生成
便利だと思いましたか？まだ終わりではありません！
REAT APIのスキーマはプログラムによりSwagger/OpenAPI definitionsと呼ばれるドキュメントで表現することができます。これによりSwagger/OpenAPIスキーマを実装するREST APIのPythonクライアントを動的に生成することができます。 これはGitHubのための実例です。
もちろん、Pythonで行うためのライブラリもあります。
bravadoagithub
Bravadoを使う場合、
# 以下でImportErrorが出た場合はこれを実行してください !pip install bravado
from bravado.client import SwaggerClient GITHUB_SWAGGER = "https://api.apis.guru/v2/specs/github.com/v3/swagger.json" client = SwaggerClient.from_url(GITHUB_SWAGGER) client # SwaggerClient(https://api.github.com/)
これで、メソッドとして全てのAPIアクションがあるかclientオブジェクトを調べて遊ぶことができます。ケーススタディとしてKLIFSを実例として取り上げます。
ヒント: clientについてこのNotebook上で調べるにはclient?を使ってください。
訳注(2020/05)
KLIFSの実例はトークトリアル11パートAを参照してください。
訳注ここまで
SOAP（Simple Object Access Protovol）のような他の標準的な方法も非常に似た原理で動作します。必要なことはクライアント生成ライブラリによって処理される定義ファイルを提供することだけです。SOAPの場合、定義ファイルは.wsdl(Web Services Description Language)のフォーマットとなっています。最も人気のあるライブラリの一つはsudsです（現在、suds-communityとしてインストール可能です）。
# sudsのインストールが必要であればこのセルを"code"に変換してください !pip install suds-community
import suds.client wsdl = 'http://www.swissdock.ch/soap/wsdl' client = suds.client.Client(wsdl) # Sudsは("?でアクセスできる) __docstring__ を持ちませんが、__str__() メソッドに多くのことが記載されています。 # 言い換えれば、ドキュメンテーションを入手するにはオブジェクトを"print"しなければなりません。 print(client)
訳注（2020/05)
残念ながら現在swisdockのサーバーは機能しておらずHTTP Error 502: Proxy Errorというエラーが表示されます。
訳注ここまで
ドキュメント構文解析
時折運悪く、Webサービスが、機械可読なドキュメントを生成する標準化されたAPIを提供していない場合があるでしょう。代わりに、必要な情報を得るために、通常のWebページを使って、HTMLコードを構文解析しなければなりません。これはウェブスクレイピング（(web)scraping）と呼ばれており、通常、 価値のあるデータを含む正しいHTMLタグとIDを見つけることを含みます（サイドバーやトップメニュー、フッター、広告 etc.といったものを無視します）。
スクレイピングでは、基本的に２つのことを行います。
requestsを使ってWebページにアクセスしHTMLコンテンツを取得します。BeautifulSoupあるいはrequests-htmlで、HTML文字列を構文解析します。
このWikipediaの記事のタンパク質を構成する(proteinogenic)アミノ酸のテーブルを構文解析して見ましょう！
import requests from bs4 import BeautifulSoup import pandas as pd r = requests.get("https://en.wikipedia.org/wiki/Proteinogenic_amino_acid") r.raise_for_status() # ここで正しいステップを推測するためには、手動でHTMLコードを調べる必要があります。 # ヒント： HTMLの定義にたどり着くにはWebページ上で右クリック＋検証（inspect content）を使ってください ;) html = BeautifulSoup(r.text) header = html.find('span', id="General_chemical_properties") table = header.find_all_next()[4] table_body = table.find('tbody') 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) 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 |
ブラウザのリモートコントロールBrowser remote control
数年前のトレンドではJavaScriptを使って動的にHTMLドキュメントを生成するサーバーを構築することがあちこちで行われていました（Wikipediaのように）。言い換えれば、HTMLはサーバに組み込まれクラアント（あなたのブラウザ）に送られます。
しかし、最近のトレンドではアプリケーション全体をJavaScriptフレームワークで全部構築することに向かっています。つまり、HTMLのコンテンツがクライアント上で動的に生成されます。伝統的な構文解析は機能せず、JavaScriptフレームワークをホストする代替のHTMLコードをダウンロードするだけです。これを回避するには、HTMLはクライアント側のJavaScriptエンジンで翻訳されなければなりません。
現在のノートブックではこれは取扱いませんが、興味があれば次のプロジェクトをチェックしてください。
ディスカッション
このイントロダクション理論編では、PythonインタープリターからオンラインのWebサービスにプログラムでアクセスする様々な方法を見てきました。これらのテクニックを活用することで、Jupyter Notebook上で自動化されたパイプラインを構築することができるようになります。このトークトリアルの続きのパートを読んでCADDへの適用例を見てください。
クイズ
- Webサイトの特定の箇所をスクレイピングするために正確なHTMLタグとIDをみつけるためにはどうすれば良いでしょうか？自動化できるでしょうか？
- プログラムによるAPI、あるいは手動で作成したスクレイピングのどちらが良いですか？
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
トークトリアル１１はちょっと私のようなプログラミングができないものには辛いです。何をやっているかよくわかりません。HTMLのタグなどWebスクレイピングについても勉強しないといけなさそうです・・・
とりあえずこれまでの訳を置いておくレポジトリを作成しました。
TeachOpenCADD トピック１０ 〜結合サイトの類似性とオフターゲット予測〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル１０をもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを！
以下、訳。
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
トークトリアル 10
結合サイトの類似性とオフターゲット予測
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Angelika Szengel, Marvis Sydow and Dominique Sydow
Note: このノートブックはセルを順番に実行してください。一度に全てのセルを実行することもできますが、PyMolのイメージが思った通りにならない可能性があります。
このトークとリアルの目的
このトークトリアルでは、タンパク質全体および結合サイトの構造的類似性を使ってオフターゲットを予測します。オフターゲットとは意図せず医薬品のターゲットとなっているタンパク質で、望まぬ副作用につながったり、あるいは、当初の意図とは異なる望ましい医薬品の適応に繋がることがあります（ドラッグリポジショニング）。結合サイトの比較の主要なステップについて議論し、基本的な手法、つまり構造間の幾何学的な差異（２つの構造の平均二乗偏差（root mean square deviation, rmsd））を実装します。
学習の目標
理論
- オフターゲットタンパク質
- コンピューターによるオフターゲット予測：結合サイト比較
- 類似性の簡単な指標としてのペアワイズ平均二乗偏差（Pairwise RMSD）
- チロシンキナーゼ阻害剤、イマチニブ、
実践
- 着目するリガンドを読み込み、可視化する（イマチニブ/STI)
- PDBから全てのタンパク質-STI複合体を取得
- PDB構造を可視化
- PDBの構造をアラインメント（タンパク質全体）
- ペアワイズRMSDの取得（タンパク質全体）
- PDBの構造をアラインメント（結合サイト）
- ペアワイズRMSDの取得（結合サイト）
レファレンス
結合サイト比較：
- 結合サイト比較についてのレビュー： (Curr. Comput. Aided Drug Des. (2008), 4, 209-20) および (J. Med. Chem. (2016), 9, 4121-51)
- PyMolのコマンド
alignについてのドキュメンテーション (PyMolWiki:align) - Wikipedia記事：平均二乗偏差（RMSD）(Wikipedia: RMSD)および構造の重ね合わせ（superposition）(Wikipedia: structural superposition)
- 構造の重ね合わせ (Book chapter: Algorithms, Applications, and Challenges of Protein Structure Alignment in Advances in Protein Chemistry and Structural Biology (2014), 94, 121-75)
イマチニブ（Imatinib）：
- イマチニブについてのレビュー (Nat. Rev. Clin. Oncol. (2016), 13, 431-46)
- イマチニブの選択性の低さ（Promiscuity of imatinib） (J. Biol. (2009), 8, 10.1186/jbiol134)
- Imatinibに関するChEMBLの情報 (ChEMBL: Imatinib)
- ImatinibのPDB上の情報 (PDB: STI)
- イマチニブの副作用 (Drugs.com: Imatinib)
- イマチニブの副作用 (BMC Struct. Biol. (2009), 9, 10.1186/1472-6807-9-7)
理論
オフターゲットタンパク質
オフターゲットは、医薬品あるいは医薬品の代謝物（の一つ）と相互作用するタンパク質で、ターゲットタンパク質として意図されていないものです。オフターゲットによって起こされる分子の応答は、望まぬ副作用につながる可能性があり、あまり害の無いものから非常に有害な副作用まであります。オフターゲットが生じる主な理由は、本来のターゲット（on-target）とオフターゲット（off-target）が互いの結合サイトに類似の構造的モチーフを共に有しており、それゆえ類似のリガンドに結合することができる、ということです。
コンピューターによるオフターゲット予測：結合サイト比較
コンピューター支援による潜在的なオフターゲットの予測の目的は、潜在的な危険性をもつ医療用物質を開発するリスクを最小にすることです。結合サイトの類似性評価にはいくつかのアルゴリズム的アプローチがありますが、常に次の３つの主要なステップで構成されています。
- 結合サイトの符号化（エンコーディング）：結合サイトを様々な記述子テクニックでエンコードし、ターゲットデータベースに格納する
- 結合サイトの比較：様々な類似性指標を使って、クエリの結合サイトをターゲットデータベースと比較する
- ターゲットランキング：適切なスコア化に基づきターゲットを順位づけする
様々な類似性指標と既存のツールについてのより詳細な情報は、結合サイト比較に関する次の素晴らしいレビュー２本を参照して下さい(Curr. Comput. Aided Drug Des. (2008), 4, 209-20 そしてJ. Med. Chem. (2016), 9, 4121-51)。
Figure 1: 結合サイト比較手法の主な工程（Dominique Sydowによる図）。
類似性の簡単な指標としてのペアワイズ平均二乗偏差
類似性のスコア化のための簡単かつ直接的な手法は平均二乗偏差（RMSD）の計算値を使うことです。RMSDはアラインメントをとった２つの構造の原子間の距離の２乗の平均の平方根となります (Wikipedia: RMSD)。
２つの構造間で比較するそれぞれの原子を見つけるために、まず、配列に基づく方法あるいは配列とは独立なアラインメントアルゴリズムによって、２つの構造のアラインメントを取る必要があります (Book chapter: Algorithms, Applications, and Challenges of Protein Structure Alignment in Advances in Protein Chemistry and Structural Biology (2014), 94, 121-75)。
チロシンキナーゼ阻害剤、イマチニブ
キナーゼはATPからタンパク質へとリン酸を移す酵素で、それによりシグナル伝達や代謝、タンパク質の制御といった様々な細胞のプロセスをコントロールします。 これらのキナーゼが（ゲノムの変異のせいで）恒常的に活性化していると、制御プロセスを歪ませ癌の原因となることがあります。抗がん剤の例としてイマチニブがあげられます(Nat. Rev. Clin. Oncol. (2016), 13, 431-46)。癌の治療に使われる低分子チロシンキナーゼ阻害剤で、より具体的には慢性骨髄性白血病（Chronic Myeloid Leukaemia、CML)）と消化管間質腫瘍（GastroIntestinal Stromal Tumour、GIST）の治療に使われます。
イマチニブの選択性は完全なものでは無く、その主要なターゲット以外のチロシンキナーゼも標的とすることが示されています。このことはドラッグリポジショニングに使われており、つまり、イマチニブは異なる癌種に対しても承認されています (J. Biol. (2009), 8, 10.1186/jbiol134)。しかしながら、アレルギー反応や感染、出血、頭痛といった望まぬ副作用も示します (Drugs.com: Imatinib)。
実践
以下では、イマチニブに結合するPDB構造を取得しフィルタリングします。イマチニブの結合したタンパク質の（構造の解かれたタンパク質との）構造類似性を調べます。類似性の指標は（簡単な類似性指標として）ペアワイズRMSD計算を使います。潜在的なオフターゲットをみつけるための最初の評価として、このシンプルな手法が利用できることを示します。
import os import pprint import pickle import glob import time from rdkit import Chem from rdkit.Chem.Draw import IPythonConsole from rdkit.Chem import Draw, AllChem from rdkit.Chem import PyMol IPythonConsole.ipython_useSVG=True import pandas as pd import numpy as np import random import matplotlib.pyplot as plt import seaborn as sns from pypdb import * from biopandas.pdb import PandasPdb
対象のリガンドを読み込み可視化する（イマチニブ（Imatinib/STI））
イマチニブ（よく使われる略語： STI）のSMILES形式はChEMBLデータベース(ChEMBL: Imatinib) あるいはPDBデータベース(PDB: STI)から、略語「STI」で取得することができます。ここでは単純に手作業でChemical Component Summaryテーブルの「Isomeric SMILES」エントリーから文字列をコピーし、リガンドを読み込みます。
訳注(2020/05)
イマチニブに関するPDBのChemical Component Summarのリンクはこちら(RCSB PDB: STI )です。
訳注ここまで
sti = Chem.MolFromSmiles('CN1CCN(Cc2ccc(cc2)C(=O)Nc2ccc(C)c(Nc3nccc(n3)-c3cccnc3)c2)CC1')
Draw.MolToImage(sti)
STIの３次元構造を調べるために、オープンソースツールのPyMolを使います（トークトリアルT8のイントロダクションを参照してください）。PyMolでSTIを眺める前に、３次元座標を計算する必要があります。
まず、化合物に水素原子を追加しますが、水素原子はSMILES形式でいつも明示的に記述されているとは限りません。次に、距離幾何学（ディスタンスジオメトリー法、Distance Geometry）を使って化合物の初期座標を取得し、UFF（Universal Force Field）力場を使って化合物構造の最適化を行います。
sti_mol = Chem.AddHs(sti) sti_mol
AllChem.EmbedMolecule(sti_mol) AllChem.UFFOptimizeMolecule(sti_mol) sti_mol
これで、PyMolに進む準備ができました！
ターミナルからPyMolを開きます（ここではJupyter notebookからコントロールしています）。
os.popen('pymol -R')
objPMV = PyMol.MolViewer()経由でPyMolとJupyter notebookを接続するために、PyMolが完全に立ち上がるまで待つ必要があります。
# エラーへの対処方法：PyMolが読み込まれるまでまつ nrTry = 0 # 現在の試行回数 ttw = 10 # 待ち時間（秒単位） while nrTry < ttw: # PyMolが読み込まれオブジェクトが保存されるまで試行する nrTry += 1 try: objPMV = PyMol.MolViewer() # PyMolオブジェクトの保存 break # PyMolが読み込まれたらループを止める except ConnectionRefusedError: # PyMolがまだ読み込まれていない場合の例外処理 time.sleep(1) # 待つ... if nrTry == ttw: # ttw試行の後：エラーメッセージを表示 print("Error: PyMol did not start correctly.\n" + "Try again and/or check if PyMol is installed completely.")
STIをPyMolに読み込み、このJupyter notebookにPyMolフレームを表示させます。
# PyMolを初期化（すでにPyMolにオブジェクトを読み込んでいるといけないので初期化する） objPMV.server.do("reinitialize") # PyMolに化合物を読み込み表示 objPMV.ShowMol(sti_mol) # sticksで表示 objPMV.server.do('cmd.show("sticks","molecule")') # PyMolの背景を白色に設定 objPMV.server.do("bg_color white") # 現在のPyMolフレームのレイトレースイメージを作成 objPMV.server.do("ray 900, 900") objPMV.GetPNG(h=300)
PDBから全てのタンパク質-STI複合体を取得
PDBのような公開データベースでイマチニブ/STIを調べ、ターゲットタンパク質として報告されているタンパク質を検索することができます。PDBでは、イマチニブは通常STIとして省略されています。以下では、両方の名称を検索し結果を統合します。
PDBを検索
まず、対象のリガンド（STI）に結合するタンパク質をPDBから全て取得します。
search_dict = make_query('STI') # PDBでリガンドSTIに結合するタンパク質を検索 found_pbd_ids = do_search(search_dict) print(found_pbd_ids) print("\nNumber of structures connected with STI in the PDB: " + str(len(found_pbd_ids))) # ['1AVU', '1AVW', '1AVX', '1BA7', '1FPU', '1IEP', '1M52', '1OPJ', '1R8N', '1R8O', '1T45', '1T46', '1XBB', '2BEA', '2BEB', '2DRE', '2ESU', '2ET2', '2HYY', '2OIQ', '2PL0', '3EA7', '3EA8', '3FW1', '3GVU', '3HEC', '3I2A', '3I2X', '3K5V', '3M3V', '3MS9', '3MSS', '3OEZ', '3PYY', '3QQM', '3S8J', '3S8K', '4BKJ', '4CSV', '4GCN', '4GCO', '4H9W', '4HA2', '4R7I', '4TLP', '5FNW', '5FNX', '5FZU', '5FZY', '5FZZ', '5G00', '5MQT', '6HD4', '6HD6', '6I0I', '6JOL', '6KTN', '6NPE', '6NPU', '6NPV', '6OTU', '6WE5'] # Number of structures connected with STI in the PDB: 62
検索対象のリガンドに使う名称（ここではイマチニブ（imatinib））によって検索結果が変わることがあることに注意してください。
search_dict2 = make_query('Imatinib') # PDBでリガンド Imatinib結合するタンパク質を検索 found_pbd_ids2 = do_search(search_dict2) print(found_pbd_ids2) print("\nNumber of structures connected with Imatinib in the PDB: " + str(len(found_pbd_ids2))) # ['1IEP', '1M52', '1OPJ', '1T46', '1XBB', '2F4J', '2GQG', '2HYY', '2OIQ', '2PL0', '2XYN', '3EL7', '3EL8', '3FW1', '3G0E', '3G0F', '3G6G', '3G6H', '3GVU', '3HEC', '3HEG', '3K5V', '3MS9', '3MSS', '3OEZ', '3OF0', '3PYY', '3QLF', '3QLG', '4BKJ', '4CSV', '4R7I', '5MQT', '6HD4', '6HD6', '6JOL', '6KTM', '6KTN', '6NPE', '6NPU', '6NPV'] # Number of structures connected with Imatinib in the PDB: 41
両方の検索結果を統合し、重複しない（ユニークな）エントリーだけを残します。
pdb_ids = list(set(found_pbd_ids + found_pbd_ids)) print("Number of structures connected with STI/Imatinib in the PDB: " + str(len(pdb_ids))) # Number of structures connected with STI/Imatinib in the PDB: 62
訳注(2020/05)
set()で２つのリストの和集合をとっています（重複を削除）。set型となるため再度list()でリスト型にしています。
訳注ここまで
PDBデータセットをフィルタリング
以下のクライテリアに従ってデータセットをフィルタリングするために、pypdbのget_entity_info関数を使ってPDB構造のメタ情報を取得します。
- 実験手法でフィルタリング (
xray) - 分解能でフィルタリング (3 Å以下)
- （簡単にするため）1本鎖のPDB構造のみを残す（single chain）
- Imatinibに結合した構造のみを残す（例、Imatinibに関連づけられているが、結合はしていないPDB構造が除かれる）
- 2019年以前に登録されたPDB IDのみを残す（トークトリアル公開時点のリソースにデータセットをあわせる）
PDBを検索する方法についての情報はトークトリアル 8を参照してください。
# PDBからメタ情報を取得 entity_info = [] for i in pdb_ids: entity_info.append(get_entity_info(i)) entity_info[0]
{'Method': {'@name': 'xray'},
'Entity': {'@id': '1',
'@type': 'protein',
'Chain': [{'@id': 'A'}, {'@id': 'B'}]},
'structureId': '6I0I',
'bioAssemblies': '1',
'release_date': 'Wed Sep 04 00:00:00 PDT 2019',
'resolution': '1.45'}
# リストをデータフレームに変換 entity_info_pd = pd.DataFrame(entity_info) a = [int(i.split(" ")[5]) for i in entity_info_pd["release_date"].tolist()] entity_info_pd.head()
| Method | Entity | structureId | bioAssemblies | release_date | resolution | |
|---|---|---|---|---|---|---|
| 0 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': [{'@... | 6I0I | 1 | Wed Sep 04 00:00:00 PDT 2019 | 1.45 |
| 1 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 5FZU | 1 | Wed Jul 06 00:00:00 PDT 2016 | 2.43 |
| 2 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 2PL0 | 1 | Tue Oct 09 00:00:00 PDT 2007 | 2.80 |
| 3 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 1R8N | 1 | Tue May 25 00:00:00 PDT 2004 | 1.75 |
| 4 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': [{'@... | 3M3V | 1 | Wed Mar 23 00:00:00 PDT 2011 | 2.70 |
# 1. 実験手法でフィルタリング entity_info_pd = entity_info_pd[entity_info_pd["Method"] == {'@name': 'xray'}] # 2. 分解能でフィルタリング entity_info_pd = entity_info_pd[entity_info_pd["resolution"].astype(float) <= 3.0] # 3. 1本鎖の構造のみを保持（簡単化のため） entity_info_pd = entity_info_pd[[type(i) == dict for i in entity_info_pd["Entity"]]] entity_info_pd = entity_info_pd[[type(i["Chain"]) == dict for i in entity_info_pd["Entity"]]] pdb_ids = entity_info_pd["structureId"].tolist() print("Number of structures after filtering: " + str(len(pdb_ids))) # Number of structures after filtering: 24
entity_info_pd.head()
| Method | Entity | structureId | bioAssemblies | release_date | resolution | |
|---|---|---|---|---|---|---|
| 1 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 5FZU | 1 | Wed Jul 06 00:00:00 PDT 2016 | 2.43 |
| 2 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 2PL0 | 1 | Tue Oct 09 00:00:00 PDT 2007 | 2.80 |
| 3 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 1R8N | 1 | Tue May 25 00:00:00 PDT 2004 | 1.75 |
| 5 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 1T45 | 1 | Tue Jun 15 00:00:00 PDT 2004 | 1.90 |
| 6 | {'@name': 'xray'} | {'@id': '1', '@type': 'protein', 'Chain': {'@i... | 1AVU | 1 | Wed Oct 28 00:00:00 PST 1998 | 2.30 |
訳注
Chainに関する情報はカラム「Entity」の辞書型の要素のうち、「Chain」をkeyとするvalueとして格納されています。
Chainの数に応じて、
* Chainが1つの時 → 「'Chain': {'@id': 'A'}」 → valueのtypeはdict
* Chainが2つの時 → 「'Chain': [{'@id': 'A'}, {'@id': 'B'}]」 → valueのtypeはlist
のようになっており、複数のChainがある場合はvalueのtypeがlist（辞書を要素とするリスト）になります。
従って、フィルタリング３行程目では「type(i["Chain"]) == dict」 とすることで「Chainが1つのもののみ」を残していると思われます。
訳注終わり
以下では、BioPandasというパッケージを使用します。BioPandasにはpandasデータフレームに生物学的高分子（PDBからmol2ファイルまで）の分子構造を読み込むために役立つ関数があります。PDBファイルを使った作業を容易にするためにPandasPdbオブジェクトを使います。
# 4. Imatinibの結合した構造だけを保持 def check_if_ligand_present(pdb_id, ligand_name): ppdb = PandasPdb().fetch_pdb(pdb_id) # PDB（原子情報、座標）を取得 return sum(ppdb.df["HETATM"]["residue_name"] == ligand_name) > 0 # エントリー STIがあるかどうかｂ確認 entity_info_pd = entity_info_pd[[check_if_ligand_present(pdb_id, "STI") for pdb_id in pdb_ids]] # 関数を適用 pdb_ids = entity_info_pd["structureId"].tolist() print("Number of structures after filtering: " + str(len(pdb_ids))) # Number of structures after filtering: 9
訳注(2020/05)
BioPandasのドキュメンテーションはこちらです。
PandasPdbオブジェクトのアトリビュートPnadasPdb.dfDataFrameの辞書となっており、PDBファイルにpandasデータフレームのようにアクセスできます（BioPandas Tutorial)。ここではPDBファイルのHETATMレコードにアクセスし残基名(residue_name)でリガンドの名称を取得しています。
また、check_if_ligand_presentは引数ligand_nameが一つでもあればTrue、なければFalseを返す関数です。上記のセルではリスト内包表記で「STI有り無しのboolリスト」を取得し、これを使ってTrueとなる行のみをDataFrameから抽出しています。
訳注ここまで
# 5. 2019年以前に登録されたPDB DIのみを保持 entity_info_pd = entity_info_pd[[int(i.split()[5]) < 2019 for i in entity_info_pd["release_date"].tolist()]] pdb_ids = entity_info_pd["structureId"].tolist() print("Number of structures after filtering: " + str(len(pdb_ids))) # Number of structures after filtering: 8
訳注(2020/05)
「release date」は「例：Wed Jul 06 00:00:00 PDT 2016」という形式なので、半角スペースでsplitした6番目（インデックス5)が「年」になります。
訳注ここまで
# 6. 手作業で目視で調べたあとで3GVUを除く（STIを２つ含む：以下の自動化ワークフローに適していない） pdb_ids.remove("3GVU")
pdb_ids # random.shuffle(pdb_ids) # IDの順序を変えたい時に使用してください # ['2PL0', '3FW1', '3HEC', '4R7I', '1XBB', '4CSV', '1T46']
フィルタリングしたPDB IDの保存
さらに解析を進めるため、フィルタリングしたデータセットのPDB IDを保存します（後ほど、このファイルに従ってPDB IDを処理するPyMolのPythonスクリプトを使います）。
pickle.dump(pdb_ids, open("../data/T10/pdb_ids.p", "wb"))
PDB構造の可視化
まず、タンパク質データセットの3次元構造を目視で調べるためPyMolに全構造を読み込みます。
このJupyter notebook上でのPyMolフレームの固定化したイメージ形式での可視化に加えて、PyMolアプリケーション上で直接３次元構造を眺め、インタラクティブに検証することをお勧めします。PyMolアプリケーションはこのトークトリアルと並行して立ち上り、操作されているはずです。
# PyMolの再初期化 objPMV.server.do("reinitialize") # タンパク質ファイルの読み込み for pdb_id in pdb_ids: cmd = "fetch " + pdb_id objPMV.server.do(cmd) # タンパク質をcartoonとして表示（使用しているpymolのversionに応じて必要な可能性があります） objPMV.server.do('cmd.show("cartoon","all")') # オブジェクトを隠す objPMV.server.do("hide polar_contacts") # 背景を白色に設定 objPMV.server.do("bg_color white") # 水分子とイオンを削除 objPMV.server.do("remove solvent") # 中心に置き、ズーム objPMV.server.do("center all") objPMV.server.do("zoom all") objPMV.server.do("ray 400,400") # Jupyter noteboomにPyMolフレームを表示 objPMV.GetPNG(h=500)
このイメージは美しくカラフルでくるくると巻いていますが、まだ有益な情報を与えてくれません。次のステップではこの構造を互いに重ね合わせます（アラインメント）。
PDBの構造をアラインメント（タンパク質全体）
PyMolはさまざまなレベルの配列類似性と構造類似性に適した、異なるアラインメント手法を提供しています：
superコマンドはかなりの構造類似性（decent structural similarity）のあるタンパク質に適しています（配列には依存しない）。alignコマンドはかなりの配列類似性（decent sequence similarity）のあるタンパク質に適しています（配列に依存）。cealignコマンドはほとんどあるいは全く配列類似性がない（little to no sequence similarity）タンパク質に対して頑健な手法です（配列依存性は中間）
このトークトリアルではalignコマンドを選択し、（配列に基づいて）RMSDを最小化することで構造を重ね合わせます。
注：このアプローチは類似の配列を持つ構造に有利なように解析されます（alignコマンドはかなりの配列類似性をもつタンパク質のペアでパフォーマンスがより良くなります）。配列類似性が低いタンパク質の比較にはsuperあるいはcealignコマンドの方が選択肢として良いかもしれません。（タンパク質ペアの配列や構造類似性がわからない）自動化されたワークフローにの場合、３つの全ての手法に基づいてRMSDを計算し、最も結果の良いものをさらなる解析のために残すことが解決策となるかもしれません。
まず、pdb_idsリストの全ての構造を一番最初の構造に対してアラインメントをとります。
immobile_pdb_id = pdb_ids[0] for mobile_pdb_id in pdb_ids[1:]: objPMV.server.do("align " + mobile_pdb_id + ", " + immobile_pdb_id) objPMV.server.do("center all") objPMV.server.do("zoom all") objPMV.server.do("ray 400,400") objPMV.GetPNG(h=500)
タンパク質のうちの一つ、3FW1は他のタンパク質と比較して重なりが悪くなっています。重なりの良いタンパク質を表示させるために、このタンパク質を隠します。
objPMV.HideObject("3FW1") objPMV.server.do("center all") objPMV.server.do("zoom all") objPMV.server.do("ray 400,400") objPMV.GetPNG(h=500)
多くの部分（例、ヘリックス）では構造の重なりの程度は高いですが、他の部分では重なりは低いか悪くなっています。
ペアワイズRMSDの取得（タンパク質全体）
（PyMolとPyMol Pythonスクリプトでは可能ですが）このJupyter notebook上ではPyMolからRMSD値を取得する関数を見つけられなかったので、PyMolのPythonスクリプトを呼びます。PyMol Pythonスクリプトでは全てのタンパク質のペアについて重ね合わせと最終的なRMSDのリファインを行います。結果としてえらえれるRMSD値をさらなる解析のために保存しておきます。
まず、このタスクのために使う、PyMolのPythonスクリプトを見ます。
f = open("./pymol_scripts/pymol_align_proteins.py") file_content = f.readlines() for i in file_content: print(i, end="") """ ''' This script aligns all pairs of input proteins. ''' # 訳者実行環境の都合上pymolのpathを追加するために以下2行追記 import sys sys.path.append('/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages') # 追記終わり。以降オリジナルスクリプト import pymol import pickle import pandas as pd import sys import numpy as np # Launch PyMol pymol.finish_launching() # Load PDB IDs pdb_ids = pickle.load(open("../data/T10/pdb_ids.p", "rb")) # Create DataFrame to store align results align_df = pd.DataFrame(index=pdb_ids, columns=pdb_ids) def align_structures(m1, m2, align_df): # Fetch PDB structures pymol.cmd.fetch(m1, "i") pymol.cmd.fetch(m2, "m") # Align proteins aln = pymol.cmd.align(mobile="m", target="i", quiet=1, object="{}_{}".format(m1, m2)) # Save align results to DataFrame align_df.loc[m1, m2] = aln # Save sequence alignment to file pymol.cmd.save(filename="../data/T10/alignment/alignment_proteins_" + "{}_{}".format(m1, m2) + ".aln", selection="{}_{}".format(m1, m2)) # Reinitialize PyMol for next structure pair to be loaded pymol.cmd.reinitialize() return align_df # Iterate over all structure pairs for n, immobile in enumerate(pdb_ids): for mobile in pdb_ids[n:]: align_df = align_structures(mobile, immobile, align_df) align_df = align_structures(immobile, mobile, align_df) # Quit PyMol pymol.cmd.quit() # Save results to file pickle.dump(align_df, open("../data/T10/align_df_proteins.p", "wb")) """
os.popen("python ./pymol_scripts/pymol_align_proteins.py")
訳注
アラインメント結果が保存されていない場合、Pythonスクリプトのエラーの可能性があります。
私の実行環境ではターミナル上でModuleNotFoundError: No module named 'pymol' が確認できたため、PyMolのPathを追加したところ実行することができました。
訳注終わり
このPyMol Pythonスクリプトは全てのペアワイズRMSDの結果をファイルに保存します。それでは結果をJupyter notebookに読み込みます。
align_df_proteins = pickle.load(open("../data/T10/align_df_proteins.p", "rb"))
align_df_proteinsは各ペアワイズの比較についてPyMolのalignコマンドの戻り値、つまり7項目のタプル、を含むデータフレームです：
- リファイン後のRMSD
- リファイン後の重ね合わされた原子の数
- リファインメントサイクルの試行数
- リファイン前のRMSD
- リファイン前の重ね合わされた原子の数
- アラインメントスコアの生の値
- 重ね合わされた残基の数
１つのタンパク質ペアを例に、このデータ構造に慣れ親しんでみます。
print("Structures in the data set:") cols = align_df_proteins.columns rows = align_df_proteins.index print(cols) print(rows) example_col = cols[0] example_row = rows[2] print("\nExample align return values for the {}-{} pair: ".format(example_col, example_row)) print(align_df_proteins.loc[example_col, example_row]) print("\nRMSD value (in Angström) after refinement for the {}-{} pair: ".format(example_col, example_row)) print(align_df_proteins.loc[example_col, example_row][0]) """ Structures in the data set: Index(['1XBB', '2PL0', '4R7I', '3HEC', '4CSV', '1T46', '3FW1'], dtype='object') Index(['1XBB', '2PL0', '4R7I', '3HEC', '4CSV', '1T46', '3FW1'], dtype='object') Example align return values for the 1XBB-4R7I pair: (1.0288203954696655, 1167, 5, 6.308946132659912, 1539, 404.5, 245) RMSD value (in Angström) after refinement for the 1XBB-4R7I pair: 1.0288203954696655 """
さて、全てのペアについてリファイン後のRMSD（あるいはその他のalignの戻り値のどれでも）を取り出し、さらなる解析のためPandsデータフレームとします。
従って、pymol_align_resultsの全てのペアに関するalign関数の戻り値、および、知りたいpositionの戻り値の位置を入力とする関数を定義します。
def extract_align_info(align_df, position): if position in [0, 3, 5]: return align_df.applymap(lambda x: np.float(x[position])) if position in [1, 2, 4, 6]: return align_df.applymap(lambda x: np.int(x[position])) else: print("Position not available.")
例えば、タンパク質の各ペアについて、重ね合わされた残基の数をチェックすることができます（alignの戻り値の最後の位置）。
extract_align_info(align_df_proteins, 6)
| 1XBB | 2PL0 | 4R7I | 3HEC | 4CSV | 1T46 | 3FW1 | |
|---|---|---|---|---|---|---|---|
| 1XBB | 292 | 272 | 245 | 159 | 249 | 258 | 81 |
| 2PL0 | 272 | 290 | 267 | 207 | 269 | 261 | 146 |
| 4R7I | 245 | 267 | 344 | 186 | 263 | 308 | 26 |
| 3HEC | 159 | 207 | 186 | 350 | 243 | 196 | 90 |
| 4CSV | 249 | 269 | 263 | 243 | 276 | 269 | 47 |
| 1T46 | 258 | 261 | 308 | 196 | 269 | 316 | 101 |
| 3FW1 | 81 | 146 | 26 | 90 | 47 | 101 | 238 |
ここでは、ほとんどのタンパク質のペアについて、残基の大部分が重ね合わせできていることがわかります（ここではEGFR配列の長さは270から350の範囲をとっています）。しかしながら、3FW1の配列のアラインメントは低くなっています。
参考までに：全てのペアの配列アラインメントは../data/T10/alignment/でチェックできます。例を下に示します。
example_alignment_path = glob.glob("../data/T10/alignment/*")[3] f = open(example_alignment_path) file_content = f.readlines() for i in file_content: print(i, end="") """ CLUSTAL i MALEEIRPKEVYLDRKLLTLED---------------------------------KELGSGNF m MKKGHHHHHHGQKPKYQVRWKIIESYEGNSYTFIDPTQLPYNEKWEFPRNNLQFGKTLGAGAF *.**.*.* i GTVKKGYYQM----KKVVKTVAVKILKNEANDPALKDELLAEANVMQQL-DNPYIVRMIGICE m GKVVEATAFGLGKEDAVLK-VAVKMLKSTAHADE-KEALMSELKIMSHLGQHENIVNLLGACT *.*.. . ****.** *..*..*...*..* ....**...*.* i AESWMLVM-EMAELGPLNKYLQQNR-----------------------HVKDKNIIELVHQVS m HGGPVLVITEYCTYGDLLNFLRRKAEAMLGPSLSPGQDPEGLDKEDGRPLELRDLLHFSSQVA *....*.*...*... ............**. i MGMKYLEESNFVHRDLAARNVLLVTQHYAKISDFGLSKALRADENYYKAQTHGKWPVKWYAPE m QGMAFLASKNCIHRDVAARNVLLTNGHVAKIGDFGLARDIMNDSNYIVKGNA-RLPVKWMAPE .**..*...*..***.*******. .***.* ..****.*** i CINYYKFSSKSDVWSFGVLMWEAFSYGQKPYRGM-KGSEVTAMLEKGERMGCPAGCPREMYDL m SIFDSVYTVQSDVWSYGILLWEIFSLGLNPYPGILVNSKFYKLVKDGYQMAQPAFAPKNIYSI .*........*****.*.*.**.**.*..**.* . .......*..*..**..*...*.. i MNLCWTYDVENRPGFAAVELRLRNYYYDVVNEGHHHHHH m MQACWALEPTHRPTFQQITSFLQEQAQEDRR-------- *..**......**.*... """
次のステップでは、リファイン後のRMSD値のデータフレームを作成します（alignの戻り値の最初の位置）。
rmsd = extract_align_info(align_df_proteins, 0)
rmsd
| 1XBB | 2PL0 | 4R7I | 3HEC | 4CSV | 1T46 | 3FW1 | |
|---|---|---|---|---|---|---|---|
| 1XBB | 0.000000 | 1.549436 | 1.028820 | 2.560788 | 0.897684 | 1.132188 | 15.648476 |
| 2PL0 | 1.549436 | 0.000000 | 1.161164 | 1.894447 | 0.871177 | 1.021322 | 19.819407 |
| 4R7I | 1.028820 | 1.161164 | 0.000000 | 2.022787 | 0.736076 | 0.687878 | 7.503986 |
| 3HEC | 2.560788 | 1.894447 | 2.022787 | 0.000000 | 1.683216 | 1.912085 | 14.695195 |
| 4CSV | 0.897684 | 0.871177 | 0.736076 | 1.683216 | 0.000000 | 0.671700 | 16.296894 |
| 1T46 | 1.132188 | 1.021322 | 0.687878 | 1.912085 | 0.671700 | 0.000000 | 15.937148 |
| 3FW1 | 15.648476 | 19.819407 | 7.503986 | 14.695195 | 16.296894 | 15.937148 | 0.000000 |
この改善されたRMSDの結果をヒートマップとして可視化します。
sns.heatmap(rmsd, linewidths=1, annot=True, cbar_kws={"label": "RMSD ($\AA$)"}, cmap="Blues") plt.show()
リファイン後のRMSDに基づいてタンパク質の類似性を見るためにヒートマップをクラスタリングします。
def plot_clustermap(rmsd, title): g = sns.clustermap(rmsd, linewidths=1, annot=True, cbar_kws={"label": "RMSD ($\AA$)"}, cmap="Blues") plt.setp(g.ax_heatmap.get_xticklabels(), rotation=0) plt.setp(g.ax_heatmap.get_yticklabels(), rotation=0) sns.set(font_scale=1.5) # プロットの保存ーテキストボックスを含めるためにbbox_inchesを使います： # https://stackoverflow.com/questions/44642082/text-or-legend-cut-from-matplotlib-figure-on-savefig?rq=1 plt.savefig("../data/T10/bsc_{}.png".format(title), dpi=300, bbox_inches="tight", transparent=True) plt.show()
plot_clustermap(rmsd, "protein")
RMSDの比較により、一つのタンパク質、つまり3FW1が他の全てのタンパク質とは異なっていることが示されました（既にアラインメントを３次元で目視で検証した結果と、重ね合わされた残基の数に基づいておこなった議論のとおりです）。
タンパク質は触媒する化学反応によって、いわゆるEC（Enzyme Commission）番号で分類されます。ここでは、このEC番号をタンパク質が属する酵素グループを確認するために使います。
# PDBからPDB IDのEC番号を取得 pdb_all_info = [get_all_info(pdb_id) for pdb_id in pdb_ids] ec_numbers = [i["polymer"]["enzClass"]["@ec"] for i in pdb_all_info] target_set = {"pdb_id": pdb_ids, "ec_number": ec_numbers} target_set = pd.DataFrame(target_set) target_set
| pdb_id | ec_number | |
|---|---|---|
| 0 | 2PL0 | 2.7.10.2 |
| 1 | 3FW1 | 1.10.5.1 |
| 2 | 3HEC | 2.7.11.24 |
| 3 | 4R7I | 2.7.10.1 |
| 4 | 1XBB | 2.7.10.2 |
| 5 | 4CSV | 2.7.10.2 |
| 6 | 1T46 | 2.7.10.1 |
3FW1、つまりヒトキノン還元酵素２（human quinone reductase 2、NQO2）はEC class 1、つまり酸化還元酵素（oxidoreductases）に属している一方で、他のタンパク質は全てEC class 2.7、つまりリン酸転移酵素（phosphorus transferases）に属しています。EC class 2.7はチロシンキナーゼ（EC 2.7.10.2）を含んでおり、これはイマチニブの目標とするターゲットです。3FW1は「医薬品デザインと慢性骨髄性白血病（chronic myelogenous leukemia）患者の治療に潜在的に関与しうる」オフターゲットとして報告がなされています(BMC Struct. Biol. (2009), 9, 10.1186/1472-6807-9-7)。
PDBの構造をアラインメント（結合サイト）
ここまでアラインメントとRMSDのリファインメントにタンパク質全体の構造を使ってきました。ですが、リガンドはタンパク質の結合サイトだけに結合し、従ってタンパク質全体よりもむしろ結合サイトの類似性がオフターゲットを予測するうえでより確からしい基準とされています。
リガンド原子のいずれかの10 Å 以内となる残基全てを選択することで、タンパク質の結合サイトを定義します。これらの結合サイトの残基をアラインメントに用い、 Cɑ 原子（タンパク質主鎖）だけをRMSDのリファインメントに用います。ここでは、pdb_idsリストの全ての構造を一番最初の構造に対してアラインメントをとります。
# PyMolの再初期化 objPMV.server.do("reinitialize") # タンパク質ファイルの読み込み for pdb_id in pdb_ids: cmd = "fetch " + pdb_id objPMV.server.do(cmd) # cartoonとしてタンパク質を表示（使用しているpymolのversionに応じて必要な可能性があります） #objPMV.server.do('cmd.show("cartoon","all")') # オブジェクトを隠す objPMV.server.do("hide polar_contacts") # 背景を白色に設定 objPMV.server.do("bg_color white") # 水分子とイオンを除去 objPMV.server.do("remove solvent") # 結合サイトのアラインメント immobile_pdb_id = pdb_ids[0] for mobile_pdb_id in pdb_ids[1:]: # STIのいずれかの原子と一定の半径以内にある原子を選択し、選択を残基全体に拡張 objPMV.server.do("select mobile_bs, byres " + mobile_pdb_id + " within 10 of (" + mobile_pdb_id + " and resn STI)") objPMV.server.do("select immobile_bs, byres " + immobile_pdb_id + " within 10 of (" + immobile_pdb_id + " and resn STI)") # アラインメントを実行 objPMV.server.do("align mobile_bs, immobile_bs") # 中心に置き、ズーム objPMV.server.do("center all") objPMV.server.do("zoom all") objPMV.server.do("ray 400,400") # Jupyter notebookにPyMolフレームを表示 objPMV.GetPNG(h=500)
ペアワイズRMSDの取得（結合サイト）
タンパク質構造全体のアラインメントとRMSDのリファインメントで先に示したのと同様に、ここでは全てのタンパク質の結合サイトのペアについてアラインメントとRMSDを計算するPyMol Pythonスクリプトを実行します。
使用するPyMolコマンドはPyMol Pythonスクリプトの中で説明されています。
f = open("./pymol_scripts/pymol_align_bindingsites.py") file_content = f.readlines() for i in file_content: print(i, end="") """ ''' This script aligns all pairs of input proteins (binding sites). ''' # pymolのpathを追加するために以下2行追記 import sys sys.path.append('/usr/local/Cellar/pymol/2.3.0/libexec/lib/python3.7/site-packages') # 追記終わり import pymol import pickle import pandas as pd import sys import numpy as np # Launch PyMol pymol.finish_launching() # Get the parameter (radius from ligand) radius = sys.argv[1] # Load PDB IDs pdb_ids = pickle.load(open("../data/T10/pdb_ids.p", "rb")) # Create DataFrame to store align results align_df = pd.DataFrame(index=pdb_ids, columns=pdb_ids) def align_structures(m1, m2, align_df): # Fetch PDB structures pymol.cmd.fetch(m1, "i") pymol.cmd.fetch(m2, "m") # Select binding site pymol.cmd.select("i_bs", "byres i within " + radius + " of (i and resn STI)") pymol.cmd.select("m_bs", "byres m within " + radius + " of (m and resn STI)") # Align binding sites aln = pymol.cmd.align(mobile="m_bs", target="i_bs", quiet=1) # Save align results to DataFrame align_df.loc[m1, m2] = aln # Reinitialize PyMol for next structure pair to be loaded pymol.cmd.reinitialize() return align_df # Iterate over all structure pairs for n, immobile in enumerate(pdb_ids): for mobile in pdb_ids[n:]: align_df = align_structures(mobile, immobile, align_df) align_df = align_structures(immobile, mobile, align_df) # Quit PyMol pymol.cmd.quit() # Save results to file pickle.dump(align_df, open("../data/T10/align_df_bindingsites.p", "wb")) """
PyMol Pythonスクリプトをターミナル経由で実行します（結合サイトの原子を定義するためリガンドータンパク質半径を入力とます）。
# PyMolのalign関数を追加って結合サイトの比較を実施 os.popen("python ./pymol_scripts/pymol_align_bindingsites.py 10")
結合サイト比較のためにalignデータフレームを読み込みます。
align_df_bindingsites = pickle.load(open("../data/T10/align_df_bindingsites.p", "rb"))
データフレームからRMSD値を抜き出します。
rmsd_bs = extract_align_info(align_df_bindingsites, 0)
extract_align_info(align_df_bindingsites, 6)
| 2PL0 | 3FW1 | 3HEC | 4R7I | 1XBB | 4CSV | 1T46 | |
|---|---|---|---|---|---|---|---|
| 2PL0 | 87 | 6 | 67 | 81 | 51 | 79 | 85 |
| 3FW1 | 6 | 38 | 4 | 6 | 7 | 8 | 6 |
| 3HEC | 67 | 4 | 83 | 68 | 48 | 67 | 69 |
| 4R7I | 81 | 6 | 68 | 90 | 59 | 83 | 86 |
| 1XBB | 51 | 7 | 48 | 59 | 65 | 57 | 55 |
| 4CSV | 79 | 8 | 67 | 83 | 57 | 83 | 79 |
| 1T46 | 85 | 6 | 69 | 86 | 55 | 79 | 98 |
RMSDの結果のクラスタリングしたヒートマップを表示します。
# ペアワイズRMSD値をクラスタリングしたヒートマップとして表示 plot_clustermap(rmsd_bs, "bs")
アラインメントをとった結合サイトのRMSD値はデータセット（EC number 2.7）における3FW1（EC number 1.10.5.1）の非類似性を示していて、PyMol上の目視による検証でもSTIがタンパク質表面上に結合していることがわかります。1XBB-4CSVと1XBB-3HECのペアも類似性が低いことを示している一方で、データセットの残りはRMSD値が低くなっています。
ここで計算したRMSD値は残基の選択（結合サイトの定義）と、あらかじめ与えられた配列のアラインメントの質に依存します。
# ディレクトリを綺麗にする（ダウンロードしたPDBとPyMolアラインメントファイルを除去） os.popen("rm ./*.cif") os.popen("rm ../data/T10/alignment/*.aln")
ディスカッション
このトークトリアルでは、イマチニブに結合するタンパク質のセットにおける類似性と非類似性を評価するために、タンパク質全体と結合サイトについて、配列のアラインメントとそれに続くRMSDのリファインメントを行いました。しかしながら、イマチニブのオフターゲットの予測のためには、解かれた構造の巨大なデータベース（PDB）と、イマチニブの目的のターゲット（チロシンキナーゼ）の結合サイトとを比較する必要があります。この場合、類似性の低い配列の比較もおこなうことになるので、より洗練された検索を可能とするために、配列に依存しないアラインメントアルゴリズムを使う手法や、結合サイトの物理化学的特性を含む手法といった、より洗練された手法が求められます。
クイズ
- 医薬品のオンターゲット（on-target）とオフターゲット（off-target）という用語を説明してください。
- なぜ結合サイト類似性が、クエリタンパク質のオフターゲットを見つけるのに使えるか説明してください？
- オフターゲットを予測する上で(i)タンパク質全体と(ii)タンパク質結合サイトのRMSD値がどう役にたつか議論してください。
- 結合サイトの情報にについて代わりとなるアプローチについて考えてください（結合サイトの比較のためにどうやって結合サイトをエンコードするか？）。
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回はオフターゲットの予測を目的としたタンパク質の構造比較でした。内容としては計算リソースの課題もあり、同一リガンドに結合することが確からしい（共結晶構造が取得されるレベルの情報がある）タンパク質の比較（RMSD）でした。検証したエントリーの中にEC番号がキナーゼですらないタンパク質も含まれており、こんなタンパク質も共結晶が取られているんだ！ということにびっくりしました。これを実際にPDBの構造全体に行おうとしたらとんでもなく大変なのではないでしょうか？？
TeachOpenCADDを開発されたVolkamer研究室は研究テーマとして「オフターゲット予測」をあげていらっしゃいますので、おそらく効率的な検索手法の開発をされているのだと思います。不勉強で研究内容の中身まで踏み込めておりません。すみません。
昨今では配列からのタンパク質構造の予測についてDeep Learningを用いた手法（AlphaFold ）の精度向上も話題になっています。配列の1Dアラインメントのみから2D、3Dの構造予測を踏まえてオフターゲットを予測する、なんてこともできるようになるのかもしれませんね（？）。
TeachOpenCADD トピック９ 〜リガンドベースファーマコフォア〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのREADMEをもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを！
以下、訳。
- トークトリアル 9
- リガンドベースファーマコフォア（Ligand-based pharmacophores）
- このトークトリアルの目的
- 学習の目標
- レファレンス
- 理論
- 実践
- nglviewで可視化
- ディスカッション
- クイズ
- 追記
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
トークトリアル 9
リガンドベースファーマコフォア（Ligand-based pharmacophores）
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Pratik Dhakal, Florian Gusewski, Jaime Rodríguez-Guerra and Dominique Sydow
注: このノートブックのセルは一つずつ順番に実行してください。全てのセルを一度に実行することも可能ですが、NGLviewのいくつかは意図した通りに動かないかもしれません。
このトークトリアルの目的
このトークトリアルでは、以前のトークトリアルで選択し、重ね合わせを行った既知のEGFRリガンドを用いて、各リガンドのドナー、アクセプター、そして疎水性ファーマコフォアフィーチャーを見つけます。次に、これらのフィーチャーをクラスター化し、組み合わせファーマコフォア（ensemble pharmacophore)を定義します。これは既知のEGFRリガンドのセットの特性を表し、バーチャルクリーニングによる新奇なEGFRリガンドの探索に用いることができます。
学習の目標
理論
- ファーマコフォアモデリング
- 構造ベース（structure-based）とリガンドベース（ligand-based）ファーマコフォアモデリング
- ファーマコフォアを用いたバーチャルスクリーニング
- クラスタリング：k平均法（k-means）
実践
- 先のトークトリアルであらかじめ重ね合わせたリガンドを取得
- NGLViewでリガンドを表示
- ファーマコフォアフィーチャー（pharmacophore feature）を抽出
- 全リガンドのファーマコフォアフィーチャーを表示
- 水素結合ドナー
- 水素結合アクセプター
- 疎水性コンタクト
- フィーチャータイプ毎にフィーチャーの座標を収集
- 組み合わせファーマコフォアの作成
- クラスターの表示
- 水素結合ドナー
- 水素結合アクセプター
- 疎水性コンタクト
- 組み合わせファーマコフォアの表示
レファレンス
- IUPACにおけるファーマコフォアの定義 (Pure & Appl. Chem (1998), 70, 1129-43)
- LigandScoutにおける３Dファーマコフォア (J. Chem. Inf. Model. (2005), 45, 160-9)
- 書籍の章： Pharmacophore Perception and Applications (Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 259-82)
- 書籍の章： Structure-Based Virtual Screening (Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 313-31).
- Monty Kierとファーマコフォアの概念の起源 (Internet Electron. J. Mol. Des. (2007), 6, 271-9)
- Nik StieflによるRDKitを使ったファーマコフォアモデリングのデモンストレーション (RDKit UGM 2016 on GitHub)
理論
ファーマコフォア
コンピューター支援医薬品デザイン（computer-aided drug design）において、ファーマコフォアを用いた医薬品とターゲット分子との相互作用の表現はよく確立された手法です。ファーマコフォアという用語は1998年IUPACの作業部会によって定義されました。即ち、「ファーマコフォアは、特定の生物学的標的の構造との最適な分子間相互作用を確実なものとし、その生物学的応答を引き起こす（あるいは阻害する）ために必要な、立体的、電子的特徴の組み合わせ」です。 (Pure & Appl. Chem. (1998), 70, 1129-43)
言い換えれば、ファーマコフォアはいくつかのファーマコフォアフィーチャーで構成されていて、ファーマコフォアフィーチャーは研究対象のターゲット分子に結合することが観察されているリガンドの重要な立体的、物理化学的特性を表現します。そのような物理化学的特性（フィーチャータイプとも呼ばれます）としては水素結合ドナー/アクセプター、疎水性/芳香族性相互作用、あるいは、正/負電荷を帯びた官能基があり、立体的特性はこれらのフィーチャーの３次元配置によって定義されます。
構造ベースファーマコフォアモデリングとリガンドベースファーマコフォアモデリング
ファーマコフォアモデリングには、生物学的な問いと、手に入るデータソースに応じて、２つの主要なアプローチが使われます。即ち、構造ベースファーマコフォアモデリングとリガンドベースファーマコフォアモデリングです。
構造ベースファーマコフォアモデル(Structure-based pharmacophore models) はタンパク質-リガンド複合体から作成されます。フィーチャーはタンパク質とリガンドの間で観察される相互作用によって定義され、リガンドの結合に関与することが示されているこれらのリガンド部位だけがバーチャルスクリーニングに使われることを保証します。 しかしながら、タンパク質-リガンド複合体の構造は全ての標的タンパク質で手に入る訳ではありません。この場合、例えばドッキングによってリガンドを標的の結合サイトにモデリングすることで複合体の構造を生成するか、あるいは、タンパク質とリガンドが相互作用する可能性がある場所を見つけるために、標的タンパク質の結合サイトだけを使ってファーコフォアモデリングを行います。
リガンドベースファーマコフォアモデル(Ligand-based pharmacophore models) は研究対象の標的分子に結合することが知られている一連のリガンドに基づきます。これらのリガンドに共通してみられる化学的特徴でファーマコフォアモデルを作成します。この手法は複数の既知リガンドがある標的タンパク質で、タンパク質ーリガンド複合体構造が無い場合に使われます。このトークトリアルでは、既知のEGFRリガンドセットを使ってリガンドベースのファーマコフォアモデリングを行います。
ファーマコフォアモデリングに関してさらに知りたい場合は、 (Pharmacophore Perception and Applications: Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 259-82) and (J. Chem. Inf. Model. (2005), 45, 160-9)をお勧めします。
Figure 1: タンパク質-リガンド相互作用を表す構造ベースのファーマコフォア（Dominique Sydowによる図）
ファーマコフォアを用いたバーチャルスクリーニング
先にトークトリアル 4で説明したように、バーチャルクリーニング（Virtual Screening, VS）は、（クエリによって代表される）研究対象のターゲット分子に結合する可能性が最も高い（ライブラリ中の）低分子を見つけるために、巨大な化合物ライブラリに対してクエリ（例えばこのトークトリアル 9ではファーマコフォアモデル、トークトリアル 4ではクエリ化合物）のスクリーニングを実施することです。ファーマコフォアベースのバーチャルスクリーニングでは、化合物ライブラリがファーマコフォアモデルに対して化合物毎に突き合わせられ、最も整合性の良かった結果にもとづきランク付けされます(Structure-Based Virtual Screening: Applied Chemoinformatics, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, (2018), 1, 313-31)。
クラスタリング：k平均法
このトークトリアルでは、いくつかのリガンドベースファーマコフォアのフィーチャーポイントをクラスタリングすることで組み合わせファーマコフォアを作成します。クラスタリングのアルゴリズムとして、データセットをk個のクラスターに分類するために使われる、k平均法クラスタリングを使います：
- k個の異なる重心（centroid）を選択し、データセットの各ポイントを最も近い重心に割り当てる
- 現在のクラスターに基づいて新しい重心を計算し、データセットの各ポイントが新しい重心のうち最も近いものに割り当てられる
- 重心が安定化するまでこの過程を繰り返す
実践
import os, glob # RDKit from rdkit import RDConfig, Chem, Geometry, DistanceGeometry from rdkit.Chem import ChemicalFeatures, rdDistGeom, Draw, rdMolTransforms, AllChem from rdkit.Chem.Draw import IPythonConsole, DrawingOptions from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib from rdkit.Numerics import rdAlignment IPythonConsole.ipython_useSVG=True import collections import pandas as pd import math from sklearn import datasets, cluster import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from collections import Counter # For handling the labels import operator
先のトークトリアルのあらかじめ重ね合わせたリガンドを取得
先のトークトリアルで重ね合わせたリガンドを取得します。
まず、全リガンドPDBファイルへのファイルパスを取得します。
mol_files = list(glob.glob("../data/T8/*_lig.pdb")) mol_files # ['../data/T8/5UG8_lig.pdb', # '../data/T8/3POZ_lig.pdb', # '../data/T8/5UG9_lig.pdb', # '../data/T8/5HG8_lig.pdb']
pdb_ids = [os.path.basename(i).split("_")[0] for i in mol_files] pdb_ids # ['5UG8', '3POZ', '5UG9', '5HG8']
訳注(2020/05)
まず、Python標準ライブラリのglobモジュールを使うことでワイルドカード*を使って一連のファイルのパスを取得しています。
次にos.pathモジュールのbasename()で、パス文字列からファイル名を取得し、splitで分割することでPDB IDだけにしています。
訳注ここまで
次に、RDKitを使ってこれらのPDBファイルから全てのリガンドを読み込みます。
mols = [] for mol_file in mol_files: mol = Chem.MolFromPDBFile(mol_file, removeHs=False) if mol is None: print(mol_file, 'could not be read') else: Chem.SanitizeMol(mol) print(Chem.MolToSmiles(mol)) mols.append(mol) rangeMols = range(1, len(mols)+1) print('Number of molecules: ', len(mols)) # CCC(O)N[C@@H]1CN(C2NC(NC3CNN(C)C3)C3NCN(C(C)C)C3N2)C[C@H]1F # CC(C)(O)CC(O)NCCN1CCC2NCNC(NC3CCC(OC4CCCC(C(F)(F)F)C4)C(Cl)C3)[C@H]21 # CCC(O)N[C@@H]1CN(C2NC(NC3CN(C)NC3OC)C3NCN(C(C)C)C3N2)C[C@H]1F # CCC(O)NC1CCCC(OC2NC(NC3CNN(C)C3)NC3NCCC32)C1 # Number of molecules: 4
Draw.MolsToGridImage(mols, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])
ここで問題に行き当たりました：PDBファイルからリガンドを読み込むと、RDKitはリガンドに対して例えば芳香環の情報を割り当てません。RDKitの関数AssignBondOrdersFromTemplateを使うと、リファレンス分子に基づいて分子の結合に情報を割り当てます。リファレンス分子として例えば、我々の場合化合物のSMILESパターンを使います。
もっと情報が知りたい場合は(RDKitディスカッションの「PDBの非タンパク質分子の芳香族が認識されない」)と (RDKitドキュメンデーションのAssignBondOrdersFromTemplate)をチェックしてください。
# PDBリガンド構造のSMILESを読み込む ligs = pd.read_csv("../data/T8/PDB_top_ligands.csv", sep="\t") # pdb_idsと同じ順番でSMILESを取得 ligs_smiles = [ligs[ligs["@structureId"]==pdb_id]["smiles"].values[0] for pdb_id in pdb_ids] # SMILESからRDKit Molオブジェクトを生成 refmols = [Chem.MolFromSmiles(smiles) for smiles in ligs_smiles] # SMILESパターン（refmols）に基づき化合物（mols)に結合次数を割り当てる mols = [AllChem.AssignBondOrdersFromTemplate(refmol, mol) for refmol, mol in zip(refmols, mols)]
Draw.MolsToGridImage(mols, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])
２次元で化合物を見ることもできます（この例では元々の座標を保持しておくために化合物をコピーしておきます）。
mols_2D = [] for mol in mols: tmp=Chem.Mol(mol) AllChem.Compute2DCoords(tmp) mols_2D.append(tmp) Draw.MolsToGridImage(mols_2D, molsPerRow=4, legends=["From PDB ID: "+i for i in pdb_ids])
nglviewで可視化
ImportErrorあるいはModuleNotFoundErrorの例外がでた場合は、下のセルのタイプをCodeに切り替えて実行してください。
!pip install nglview !nglview install
import nglview as nv
import time def show_ligands(rdkit_mols): v = nv.NGLWidget() for mol in rdkit_mols: c = v.add_component(Chem.MolToPDBBlock(mol), ext="pdb") time.sleep(0.1) c.clear() c.add_ball_and_stick(multipleBond=True) return v
v = show_ligands(mols) v
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
ファーマコフォアフィーチャーの抽出
上で述べたように、このトークトリアルの目的はリガンドセットからリガンドベースのファーマコフォアの組み合わせを作成することです。まず最初に、リガンド毎にファーマコフォアフィーチャーを取り出す必要があります。そこで、フィーチャーファクトリーを（デフォルトのフィーチャーの定義で）読み込みます。
RDKitドキュメンテーションの「chemical features and pharmacophores」も参照してください。.
ffact = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'))
RDKitに実装されているファーマコフォアフィーチャーを見てみましょう。
list(ffact.GetFeatureDefs().keys()) """ ['Donor.SingleAtomDonor', 'Acceptor.SingleAtomAcceptor', 'NegIonizable.AcidicGroup', 'PosIonizable.BasicGroup', 'PosIonizable.PosN', 'PosIonizable.Imidazole', 'PosIonizable.Guanidine', 'ZnBinder.ZnBinder1', 'ZnBinder.ZnBinder2', 'ZnBinder.ZnBinder3', 'ZnBinder.ZnBinder4', 'ZnBinder.ZnBinder5', 'ZnBinder.ZnBinder6', 'Aromatic.Arom4', 'Aromatic.Arom5', 'Aromatic.Arom6', 'Aromatic.Arom7', 'Aromatic.Arom8', 'Hydrophobe.ThreeWayAttach', 'Hydrophobe.ChainTwoWayAttach', 'LumpedHydrophobe.Nitro2', 'LumpedHydrophobe.RH6_6', 'LumpedHydrophobe.RH5_5', 'LumpedHydrophobe.RH4_4', 'LumpedHydrophobe.RH3_3', 'LumpedHydrophobe.tButyl', 'LumpedHydrophobe.iPropyl'] """
一例として、参考例の化合物の全てのフィーチャーを取得してみます。
m1 = mols[0] feats = ffact.GetFeaturesForMol(m1) print('Number of features found:',len(feats)) # Number of features found: 14
フィーチャーのタイプ（RDKitではfamilyと呼ばれています）はGetFamily()で取得することができます。
feats[0].GetFamily() # 'Donor'
参考例の化合物のフィーチャーのタイプについてその頻度を取得します。
feats_freq = collections.Counter([x.GetFamily() for x in feats]) feats_freq """ Counter({'Donor': 2, 'Acceptor': 6, 'PosIonizable': 1, 'Aromatic': 3, 'Hydrophobe': 1, 'LumpedHydrophobe': 1}) """
訳注(2020/05)
Python標準ライブラリcollectionsのCounterクラスにより、与えられたリストの要素の出現頻度について辞書型オブジェクトが生成されています。keyが要素、値が出現頻度で、ここではフィーチャータイプと頻度の辞書になっています。
訳注ここまで
上記の関数をリガンドセットの全化合物に適用します。化合物ごとのフィーチャータイプの頻度をデータフレームとして表示します。
# 化合物ごとのフィーチャータイプの頻度を取得 mols_feats_freq = [] for i in mols: feats = [x.GetFamily() for x in ffact.GetFeaturesForMol(i)] feats_freq = collections.Counter(feats) mols_feats_freq.append(feats_freq) # データをデータフレームとして表示 p = pd.DataFrame(mols_feats_freq, index=["m"+str(i) for i in range(1, len(mols)+1)]).fillna(0).astype(int) p.transpose()
| m1 | m2 | m3 | m4 | |
|---|---|---|---|---|
| Donor | 2 | 3 | 2 | 4 |
| Acceptor | 6 | 5 | 7 | 5 |
| PosIonizable | 1 | 0 | 1 | 0 |
| Aromatic | 3 | 4 | 3 | 4 |
| Hydrophobe | 1 | 3 | 1 | 2 |
| LumpedHydrophobe | 1 | 2 | 1 | 1 |
この先、このトークとリアルでは次のフィーチャータイプのみに焦点をあてます。すなわち水素結合アクセプター（acceptors）、水素結合ドナー（donors）、そして疎水性コンタクト（hydrophobics）です。
フィーチャータイプごと、そして化合物ごとに、RDKitのフィーチャーオブジェクトを取得します。
acceptors = [] donors = [] hydrophobics = [] for i in mols: acceptors.append(ffact.GetFeaturesForMol(i, includeOnly='Acceptor')) donors.append(ffact.GetFeaturesForMol(i, includeOnly='Donor')) hydrophobics.append(ffact.GetFeaturesForMol(i, includeOnly='Hydrophobe')) features = {"donors": donors, "acceptors": acceptors, "hydrophobics": hydrophobics}
全てのリガンドのファーマコフォアフィーチャーを表示
ファーマコフォアフィーチャーのタイプは大抵の場合、定義された色で表示されます。例えば、水素結合ドナー、水素結合アクセプターそして、疎水性コンタクトはそれぞれ緑、赤、黄色でよく表されます。
feature_colors = {"donors": (0,0.9,0), # Green
"acceptors": (0.9,0,0), # Red
"hydrophobics": (1,0.9,0)} # Yellow
def visualize_features(molecules, feature_type, features, color="yellow", sphere_radius=0.5): print("Number of", feature_type, "in all ligands:", sum([len(i) for i in features])) v = show_ligands(molecules) for i, feature in enumerate(features, 1): for feat in feature: loc = list(feat.GetPos()) label = f"{feature_type}_{i}" v.shape.add_sphere(loc, color, sphere_radius, label) return v
考察中のフィーチャータイプのフィーチャーを可視化するためにこの関数を使います。
水素結合ドナー
feature_type = "donors" v = visualize_features(mols, feature_type, features[feature_type], feature_colors[feature_type]) v # Number of donors in all ligands: 11
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
水素結合アクセプター
feature_type = "acceptors" v = visualize_features(mols, feature_type, features[feature_type], feature_colors[feature_type]) v # Number of acceptors in all ligands: 23
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
疎水性コンタクト
feature_type = "hydrophobics" v = visualize_features(mols, feature_type, features[feature_type], feature_colors[feature_type]) v # Number of hydrophobics in all ligands: 7
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
フィーチャータイプ毎にフィーチャーの座標を集める
（フィーチャータイプ毎に）フィーチャーをクラスタリングしたいので、（フィーチャータイプ毎に）フィーチャーのすべての座標を集めます。
features_coord = {"donors": [list(item.GetPos()) for sublist in features["donors"] for item in sublist],
"acceptors": [list(item.GetPos()) for sublist in features["acceptors"] for item in sublist],
"hydrophobics": [list(item.GetPos()) for sublist in features["hydrophobics"] for item in sublist]}
これで全てのフィーチャーの座標が手に入りました。例えばacceptorフィーチャーの全ての座標は次のようになります。
features_coord["acceptors"] """ [[-13.788, 14.818, -27.097], [-12.118, 16.261, -28.026], [-11.376, 12.959, -29.238], [-16.74, 10.887, -25.841], [-16.005, 17.279, -22.5], [-15.388, 19.42, -27.168], [-11.052, 12.941, -29.208], [-12.995, 19.311, -27.816], [-10.269, 15.094, -29.72], [-8.077, 20.488, -31.958], [-12.686, 20.459, -24.958], [-13.863, 14.66, -27.134], [-12.175, 16.118, -27.994], [-11.411, 12.856, -29.266], [-17.026, 10.809, -26.396], [-15.601, 9.901, -28.086], [-16.044, 17.165, -22.506], [-15.445, 19.229, -27.215], [-13.162, 14.493, -28.185], [-11.603, 12.928, -29.136], [-12.896, 16.874, -28.352], [-16.001, 17.24, -22.411], [-16.523, 11.38, -26.098]] """
組み合わせファーマコフォアの作成
組み合わせファーマコフォアを作成するために、k平均法クラスタリングをつかってフィーチャータイプ毎にフィーチャーをクラスタリングします。
k平均法クラスタリングの静的パラメーターを設定
kq: このパラメーターで、フィーチャーポイントの数に応じてフィーチャータイプ毎のクラスター数 k を決定します。つまり、フィーチャータイプ毎に：
k = number_of_features / kq
となります。
# k quotient (kq) はk平均法でkを決めるために使われます： k = number of feature points / kq # kq は全てのクラスターについてk（フィーチャークラスター）が少なくとも1クラスターで、4-5クラスターよりも大きくならないように選択するべきです。 kq = 7
クラスター選別のための静的パラメーターを設定
min_cluster_size: リガンドの集合のなかのほとんどの化合物がもっているフィーチャーを含む可能性があるクラスターだけを残します。したがって、この変数は、我々のリガンドの集合の化合物数の75%に設定します。
top_cluster_number: このパラメーターで、最も大きいクラスターだけを選択します。
# クラスタリングの閾値：数値=閾値のパーセンテージ min_cluster_size = int(len(mols) * 0.75) # 上位のフィーチャーだけを表示 top_cluster_number = 4
k平均法クラスタリングとクラスター選択の関数を定義
k平均法クラスタリングから導いたクラスターの中心を計算する関数を定義します。
def clustering(feature_coord, kd): ''' この関数は入力のフィーチャーの座標のk平均法クラスタリングを計算します ''' # パラメーターkを、"k quotient"でフィーチャーの数を割った値として定義 k = math.ceil(len(feature_coord) / kq) k = 2 if k == 1 else k # このトークトリアルの例のにあうように疎水性コンタクトのkを調整 print('Clustering: \nVariable k in k-means: %d of %d points\n'%(k, len(feature_coord))) # k平均法を初期化 k_means = cluster.KMeans(n_clusters=k) # k平均法クラスタリングを計算 k_means.fit(feature_coord) # クラスターを返す return k_means
クラスターをサイズで並べ替え、最も大きいクラスターのインデックスのリストを出力する関数を定義します。
def get_clusters(k_means, min_cluster_size, top_cluster_number): ''' この関数は入力されたk平均法クラスタリングの情報を取得します： * 各フィーチャーのクラスターのラベルを取得 * クラスターのサイズを数え、クラスターサイズによってクラスターのインデックスを並べ替える * サイズに基づきクラスターを選択 * 選択したクラスターのインデックスを返す ''' # サイズによってクラスターを並べ替え、最大のもののみを表示 feature_labels = k_means.labels_ print('Cluster labels for all features: \n%s\n'% feature_labels) feature_labels_count = Counter(feature_labels) print('Cluster label counter: \n%s\n'% feature_labels_count) feature_labels_count = sorted(feature_labels_count.items(), key=operator.itemgetter(1), reverse=True) print('Sorted cluster label counters: \n%s\n'% feature_labels_count) # 閾値よりも大きなクラスターで、最大のものの番号を取得（選択されたクラスター） cluster_indices_sel = [] for cluster_index, cluster_size in feature_labels_count: # feature_labels_count = list of (cluster_index, cluster_size) if cluster_size >= min_cluster_size and top_cluster_number > 0: cluster_indices_sel.append(cluster_index) top_cluster_number -= 1 print('Cluster indices of selected clusters: \n%s\n'% cluster_indices_sel) return cluster_indices_sel
訳注(2020/05)
feature_labels_countは、collectionsのCounterでクラスターのラベルを集計（要素の出現頻度の辞書を作成）した後に並べ替えています。counterクラスは辞書型で順番が保証されないので、items()とsorted()を使うことでソートされたタプル(key, value)のリストを得ています。辞書のvalueでソートしたいので、sorted()の引数のkeyに、operator.itemgetter()を使用しています。（itemgetter(0)ではなく）itemgetter(1)とすることで辞書のvalueが使われることになります。
訳注ここまで
フィーチャーをクラスタリング
各フィーチャータイプについて、定義したclustering関数を使ってk平均法クラスタリングを実施します。
k_means = {"donors": clustering(features_coord["donors"], kq),
"acceptors": clustering(features_coord["acceptors"], kq),
"hydrophobics": clustering(features_coord["hydrophobics"], kq)}
"""
Clustering:
Variable k in k-means: 2 of 11 points
Clustering:
Variable k in k-means: 4 of 23 points
Clustering:
Variable k in k-means: 2 of 7 points
"""
適切なクラスターの選択
各フィーチャータイプについて、定義したget_clusters関数をつかって適切なクラスターを選択します。
print("Hydrogen bond donors\n") cluster_indices_sel_don = get_clusters(k_means["donors"], min_cluster_size, top_cluster_number) """ Hydrogen bond donors Cluster labels for all features: [1 0 0 0 1 1 0 1 1 0 1] Cluster label counter: Counter({1: 6, 0: 5}) Sorted cluster label counters: [(1, 6), (0, 5)] Cluster indices of selected clusters: [1, 0] """
print("Hydrogen bond acceptors\n") cluster_indices_sel_acc = get_clusters(k_means["acceptors"], min_cluster_size, top_cluster_number) """ Hydrogen bond acceptors Cluster labels for all features: [0 0 0 2 1 1 0 1 0 3 1 0 0 0 2 2 1 1 0 0 0 1 2] Cluster label counter: Counter({0: 11, 1: 7, 2: 4, 3: 1}) Sorted cluster label counters: [(0, 11), (1, 7), (2, 4), (3, 1)] Cluster indices of selected clusters: [0, 1, 2] """
print("Hydrophobic contacts\n") cluster_indices_sel_h = get_clusters(k_means["hydrophobics"], min_cluster_size, top_cluster_number) """ Hydrophobic contacts Cluster labels for all features: [1 0 0 1 1 0 1] Cluster label counter: Counter({1: 4, 0: 3}) Sorted cluster label counters: [(1, 4), (0, 3)] Cluster indices of selected clusters: [1, 0] """
cluster_indices_sel = {"donors": cluster_indices_sel_don,
"acceptors": cluster_indices_sel_acc,
"hydrophobics": cluster_indices_sel_h}
選択したクラスターの座標を取得
def get_selected_cluster_center_coords(k_means, cluster_indices_sel, feature_type): ''' この関数は選択したクラスターのクラスター中心の座標を（インデックスで）取得します。 ''' # ある特定のフィーチャータイプについてクラスター中心を取得 cluster_centers = k_means[feature_type].cluster_centers_ # （インデックスにより要素を選択するために）リストに型変換したのち、PandasのSeriesに型変換する cluster_centers = pd.Series(cluster_centers.tolist()) # 選択したクラスターのインデックスでクラスター中心を選択 cluster_centers_sel = cluster_centers[cluster_indices_sel[feature_type]] # リストに型変換しリストを返す return list(cluster_centers_sel)
cluster_centers_sel = {"donors": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "donors"),
"acceptors": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "acceptors"),
"hydrophobics": get_selected_cluster_center_coords(k_means, cluster_indices_sel, "hydrophobics")}
cluster_centers_sel["acceptors"] """ [[-12.155727272727272, 14.545636363636364, -28.48690909090909], [-14.937714285714286, 18.586142857142857, -24.93914285714286], [-16.4725, 10.744250000000001, -26.605249999999998]] """
クラスターを表示
フィーチャータイプ毎に、全ての化合物と全てのフィチャーポイントと一緒にクラスター中心を可視化します。
def visualize_clusters(molecules, feature_type, features, clusters, color="yellow", feature_radius=0.5, cluster_radius=1): v = visualize_features(molecules, feature_type, features, color=color, sphere_radius=feature_radius) for i, center in enumerate(clusters, 1): v.shape.add_sphere(list(center), color, cluster_radius, f"cluster_{feature_type}_{i}") return v
水素結合ドナー
feature_type = "donors" v = visualize_clusters(mols, feature_type, features[feature_type], cluster_centers_sel[feature_type], feature_colors[feature_type]) v # Number of donors in all ligands: 11
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
水素結合アクセプター
feature_type = "acceptors" v = visualize_clusters(mols, feature_type, features[feature_type], cluster_centers_sel[feature_type], feature_colors[feature_type]) v # Number of acceptors in all ligands: 23
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
疎水性コンタクト
feature_type = "hydrophobics" v = visualize_clusters(mols, feature_type, features[feature_type], cluster_centers_sel[feature_type], feature_colors[feature_type]) v # Number of hydrophobics in all ligands: 7
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
組み合わせファーマコフォアの表示
この最後のステップでは、クラスタリングしたファーマコフォアのフィーチャー（即ち、水素結合ドナー、アクセプターと疎水性コンタクト）を組み合わせて、一つの組み合わせファーマコフォアを作成します。この組み合わせファーマコフォアは４つの選択したリガンドのファーマコフォアの特徴を表します。
v = show_ligands(mols) # クラスターの読み込み for feature_type in cluster_indices_sel.keys(): centers = cluster_centers_sel[feature_type] for i, loc in enumerate(centers): sphere_radius = 1 feature_color = feature_colors[feature_type] label = f"{feature_type}_c{i}" v.shape.add_sphere(loc, feature_color, sphere_radius, label) v
訳注(2020/05)
Jupyter Notebook上ではここに3Dの描画があらわれ、インタラクティブな観察が可能となります。
訳注ここまで
v.render_image(),
v._display_image()
ディスカッション
このトークトリアルでは、EGFRに結合することが知られているリガンドのセットを、予め重ね合わせたものを用いて、組み合わせファーマコフォアモデルを作成しました。これで、既知のリガンドに観測された立体的、物理化学的特徴を示し、EGFR結合サイトに結合するかもしれないような新奇な低分子を見つけることを目的として、低分子化合物の巨大なライブラリのバーチャルスクリーニングを行うために、この組み合わせファーマコフォアモデルを使うことができます。
スクリーニングの前に、通常、ファーアマコフォアモデルはさらに最適化されます。例えば、スクリーニングに用いるフィーチャーの特徴の数を減らすために、生物学的な知識（ある相互作用は重要である一方、他の相互作用はそうではないという報告されているかもしれません）、あるいは化学的専門知識に基づいて、いくつかのフィーチャーは除外されるかもしれません。
このトークトリアルではバーチャルスクリーニングは扱いませんが、RDKitを使ってファーマコフォアモデリングとバーチャルスクリーニングをデモンストレーションしている、Nik Stieflによる素晴らしいチュートリアルに言及しておきます。(RDKit UGM 2016 on GitHub).
ファーマコフォアフィーチャーのクラスタリングを行うためにk平均法クラスタリングを用いました。このクラスタリング手法は、使用者が予めクラスターの数を定義する必要があるという欠点があります。通常、クラスターの数は、クラスタリングの実施前（あるいはクラスターの精度を高めていく過程）に、ポイントの分布を目視で調査して決めるので、ファーマコフォアを自動的に作成するには妨げとなります。この解決策として、密度に基づくクラスタリング手法（とk平均法クラスタリングの組み合わせも）が挙げられます。
クイズ
- ファーマコフォアフィーチャーとファーマコフォアを説明してください。
- 構造ベースのファーマコフォアモデリングとリガンドベースのファーマコフォアモデリングの違いを説明してください。
- 組み合わせファーマコフォアを導く方法を説明してください。
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回のトピックに関連して日本語で読める記事として株式会社理論創薬研究所 吉森篤史氏による以下の記事をお勧めします。 RDKitを利用したファーマコフォアの定義とその使用方法について解説がされていてとてもわかりやすいです。
- SAR News No.19 (Oct. 2010) 「Chemoinformatics Toolkits を用いた創薬システム開発における ラピッドプロトタイピング」
欲を言えば作ったファーマコフォアモデルによるバーチャルスクリーニングの実施まで行なってほしかったですが、計算リソースとデータセットの準備だけでも大変そうですね。
勝手にいろいろと引用させていただいておりますので問題があれば削除します。
TeachOpenCADD トピック８ 〜タンパク質データの取得〜
こちらはTeachOpenCADDの試訳です。TeachOpenCADD GitHubレポジトリのトークトリアル８をもととしておりCC BY 4.0のライセンスに従っています。
GitHubとはてなのMarkdownの取り扱いの違いのせいかフォーマットの不具合がありますがお目こぼしを！
~~ 追記 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下のページに翻訳したノートブックを保存しています。よろしければご利用ください。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
以下、訳。
トークトリアル 8
タンパク質データの取得: Protein Data Bank (PDB)
Developed in the CADD seminars 2017 and 2018, AG Volkamer, Charité/FU Berlin
Anja Georgi, Majid Vafadar and Dominique Sydow
このトークトリアルの目的
このトークトリアルでは、次のトークトリアルのための土台を準備します。次のトークトリアルではEGFRに対するリガンドベースの組み合わせファーマコフォア（ligand-based ensemble pharmacophore）を作成します。したがって準備として(i) PDBデータベースからEGFRの全てのPDB IDを取得し、(ii) X線結晶構造解析による構造で、最も質の良い５つのタンパク質ーリガンド複合体構造を取得、そして、(iii) ５つの構造を互いに3Dで重ね合わせ（アラインメント）、次のトークトリアルで使用するためにリガンドを抜き出し保存します。
学習の目標
理論
実践
- 検索するクエリタンパク質の選択
- クエリタンパク質のPDBエントリーに関する統計値を取得
- クエリタンパク質の全PDB IDを取得
- PDBエントリーのメタ情報を取得
- PDBエントリーのメタ情報をフィルタリングし並べ替え
- 上位の構造からリガンドのメタ情報を取得
- 上位リガンド化合物を描画
- タンパク質ーリガンド IDのペアを作成
- PDB構造ファイルの取得
- PDB構造の重ね合わせ
レファレンス
- Protein Data Bank (PDB website)
- PyPDB pythonパッケージ (Bioinformatics (2016), 32, 159-60)
- PyPDB pythonパッケージドキュメンテーション (PyPDB website)
- PyMol selection algebra (PyMolWiki: selection algebra)
理論
Protein Data Bank (PDB)
Protein Data Bank (PDB) は構造生物学に関する情報のデータベースの中で最も包括的なものの一つで、構造ゲノム科学と医薬品デザインといった構造生物学の分野で重要な情報源です。(PDB Webサイト)
構造のデータは、X線結晶構造解析（最も良く使われている手法）、核磁気共鳴（NMR）そしてクライオ電子顕微鏡（cryo-EM）といった構造決定手法によって取得されています。各エントリーについて、データベースに含まれている情報には(i) タンパク質、リガンド、補因子、水分子、そしてイオンといった原子と、原子同士をつなぐ結合の３次元座標、そして(ii) PDB ID、著者、登録日、用いられた構造決定手法や構造の分解能といった構造情報に関するメタ情報があります。
構造の分解能（resolution）は、集められたデータの質の指標で、単位はÅ（オングストローム）です。この値が小さいほど、構造の質はより良いものとなります。
PDB Webサイトはタンパク質の構造（と、手に入る場合はリガンドとの相互作用）の３次元描画と、構造の質の指標（メトリクス）とを提供しています。例として、上皮増殖因子（EGFR）のPDBエントリー PDB ID 3UG5 で見ることができます。
Figure 1: 上皮増殖因子（EGFR）の例としてPDB ID 3UG5のタンパク質構造（灰色）と相互作用するリガンド（緑色）を描画（Dominique Sydowによる図）。
PyPDB
PyPDBはPDBのためのPythonプログラミングインターフェースで、Python3でのみ機能します。 このパッケージを利用することで、バイオインフォマティクスのワークフローにPDB自動検索を組み込み、既存の検索結果に基づいて多数の検索をおこなうプロセスを簡略化することができます。 また、PDBエントリーの情報について高度な検索を行うことも可能になります。PDBは現在RESTful APIを使用しており、標準的なHTMLボキャブラリーを使用して情報を取得することができます。PyPDBはこれらのオブジェクトをXML文字列に変換します。 (Bioinformatics (2016), 32, 159-60)
提供されている機能（関数）のリストはPyPDBドキュメンテーションのWebサイトで見ることができます(PyPDB website)。
実践
# 必要なライブラリのインポート from pypdb import * from pymol import * from rdkit import Chem from rdkit.Chem.Draw import IPythonConsole from rdkit.Chem import Draw from rdkit.Chem import PandasTools IPythonConsole.ipython_useSVG=True import pprint import glob import pandas as pd from array import array import numpy as np import collections import matplotlib.pyplot as plt %matplotlib inline
検索するクエリタンパク質の選択
このトークトリアルではEGFRをクエリタンパク質として使います。EGFRのUniProt IDはP00533で、以降のPDBデータベース検索でこのIDを使います。
クエリタンパク質のPDBエントリーに関する統計値を取得
まず最初の質問です。「EGFRについて、毎年いくつのエントリーがPDBに登録されていて、全体ではいくつになるのでしょうか？」
PDB Webサイト で検索ターム「P00533」で検索することができます。2018年10月現在、PDBでの検索結果は179でした。
訳注(2020/05)
2020年5月現在、検索結果は207でした。
訳注ここまで
pypdbを使うと、PDBデータベースから全てのEGFR構造の登録日（deposition date）を取得することができます。find_dates関数のパラメーターmax_resultsを設定するために、登録されている構造の数が必要となります。
# 注: max_results パラメーターのデフォルトは100でEGFRのエントリーに対して少なすぎます。 # max_results > EGFR構造の最大の数 となるとエラーが出ます。 # 従って実行前に結果がいくつあるかチェックしました。(#179) # 訳注: 2020年5月時点では207で、以下のコードも合わせて変更しています。207ではStopIterationとなったので205としています。 # このデータベースクエリには時間がかかるかもしれません（数分程度） # all_dates = find_dates("P00533", max_results=179) all_dates = find_dates("P00533", max_results=205)
print("Number of EGFR structures found: " + str(len(all_dates))) # Number of EGFR structures found: 205
# 登録日の例 all_dates[:3] # ['2002-03-28', '2002-06-17', '2002-06-17']
登録日から登録年の情報を抽出し、各年毎の登録数のヒストグラムを計算します。
# 登録年の抽出 all_dates = np.asarray(all_dates) all_years = np.asarray([int(depdate[:4]) for depdate in all_dates]) # ヒストグラムの計算 bins = max(all_years)-min(all_years) # ビンの数 = 登録年の範囲 subs_v_time = np.histogram(all_years, bins) # 全エントリー（2018年を除く）のプロット # 訳注: 最新年を除く設定なので用いたデータにより除かれる年は異なります dates, num_entries = subs_v_time[1][:-1], subs_v_time[0] # ヒストグラムの表示 fig = plt.figure() ax = plt.subplot(111) ax.fill_between(dates, 0, num_entries) ax.set_ylabel("New entries per year") ax.set_xlabel("Year") ax.set_title("PDB entries for EGFR") plt.show()
クエリタンパク質の全PDB IDを取得
次に、pypdb関数のmake_queryとdo_searchを使って、クエリタンパク質EGFRについての全てのPDB構造を取得します。
search_dict = make_query("P00533") # max_resultsが180以上の場合タイムアウトになる可能性があります。 found_pdb_ids = do_search(search_dict) print("PDB IDs found for query: ") print(found_pdb_ids) print("\nNumber of structures: " + str(len(found_pdb_ids))) # PDB IDs found for query: # ['1IVO', '1M14', '1M17', '1MOX', '1XKK', '1YY9', '1Z9I', '2EB2', '2EB3', '2GS2', '2GS7', '2ITN', '2ITO', '2ITP', '2ITQ', '2ITT', '2ITU', '2ITV', '2ITW', '2ITX', '2ITY', '2ITZ', '2J5E', '2J5F', '2J6M', '2JIT', '2JIU', '2JIV', '2KS1', '2M0B', '2M20', '2N5S', '2RF9', '2RFD', '2RFE', '2RGP', '3B2U', '3B2V', '3BEL', '3BUO', '3C09', '3GOP', '3GT8', '3IKA', '3LZB', '3NJP', '3OB2', '3OP0', '3P0Y', '3PFV', '3POZ', '3QWQ', '3UG1', '3UG2', '3VJN', '3VJO', '3VRP', '3VRR', '3W2O', '3W2P', '3W2Q', '3W2R', '3W2S', '3W32', '3W33', '4G5J', '4G5P', '4HJO', '4I1Z', '4I20', '4I21', '4I22', '4I23', '4I24', '4JQ7', '4JQ8', '4JR3', '4JRV', '4KRL', '4KRM', '4KRO', '4KRP', '4LI5', '4LL0', '4LQM', '4LRM', '4R3P', '4R3R', '4R5S', '4RIW', '4RIX', '4RIY', '4RJ4', '4RJ5', '4RJ6', '4RJ7', '4RJ8', '4TKS', '4UIP', '4UV7', '4WD5', '4WKQ', '4WRG', '4ZAU', '4ZJV', '4ZSE', '5C8K', '5C8M', '5C8N', '5CAL', '5CAN', '5CAO', '5CAP', '5CAQ', '5CAS', '5CAU', '5CAV', '5CNN', '5CNO', '5CZH', '5CZI', '5D41', '5EDP', '5EDQ', '5EDR', '5EM5', '5EM6', '5EM7', '5EM8', '5FED', '5FEE', '5FEQ', '5GMP', '5GNK', '5GTY', '5GTZ', '5HCX', '5HCY', '5HCZ', '5HG5', '5HG7', '5HG8', '5HG9', '5HIB', '5HIC', '5J9Y', '5J9Z', '5JEB', '5LV6', '5SX4', '5SX5', '5U8L', '5UG8', '5UG9', '5UGA', '5UGB', '5UGC', '5UWD', '5WB7', '5WB8', '5X26', '5X27', '5X28', '5X2A', '5X2C', '5X2F', '5X2K', '5XDK', '5XDL', '5XGM', '5XGN', '5XWD', '5Y25', '5Y9T', '5YU9', '5ZTO', '5ZWJ', '6ARU', '6B3S', '6D8E', '6DUK', '6JRJ', '6JRK', '6JRX', '6JWL', '6JX0', '6JX4', '6JXT', '6JZ0', '6P1D', '6P1L', '6P8Q', '6S89', '6S8A', '6S9B', '6S9C', '6S9D', '6V5N', '6V5P', '6V66', '6V6K', '6V6O', '6VH4', '6VHN', '6VHP'] # Number of structures: 205
訳注(05/2020)
RCSB PDBのWebページで「P00533」を検索するとヒットしたエントリーは207でしたが、ここでは205個のエントリーがヒットしています。「2.5.2節：クエリタンパク質のPDBエントリーに関する統計値を取得」においてmax_results=207とすると途中でとまってしまいましたが、205では実行可能でした。タイムアウトが原因か、もしくはデータベースの問題かはわかりませんでした。
訳注ここまで
PDBエントリーのメタ情報を取得
describe_pdbを使って構造のメタ情報を取得します。メタ情報は構造毎に辞書として格納されます。
注：ここではPDB構造のメタ情報を取得するだけで、構造そのもの（３次元座標）はまだ取得しません。
# このデータベースクエリには少し時間がかかるかもしれません。 pdbs = [] for i in found_pdb_ids: pdbs.append(describe_pdb(i)) pdbs[0]
{'relatedPDB': {'@pdbId': '1JL9',
'@details': '1JL9 contains dymeric human EGF molecules.'},
'structureId': '1IVO',
'title': 'Crystal Structure of the Complex of Human Epidermal Growth Factor and Receptor Extracellular Domains.',
'pubmedId': '12297050',
'expMethod': 'X-RAY DIFFRACTION',
'resolution': '3.30',
'keywords': 'TRANSFERASE/SIGNALING PROTEIN',
'nr_entities': '2',
'nr_residues': '1350',
'nr_atoms': '8813',
'deposition_date': '2002-03-28',
'release_date': '2002-10-16',
'last_modification_date': '2011-07-13',
'structure_authors': 'Ogiso, H., Ishitani, R., Nureki, O., Fukai, S., Yamanaka, M., Kim, J.H., Saito, K., Shirouzu, M., Yokoyama, S., RIKEN Structural Genomics/Proteomics Initiative (RSGI)',
'citation_authors': 'Ogiso, H., Ishitani, R., Nureki, O., Fukai, S., Yamanaka, M., Kim, J.H., Saito, K., Inoue, M., Shirouzu, M., Yokoyama, S.',
'status': 'CURRENT'}
PDBエントリーのメタ情報をフィルタリングし並べ替え
取得した情報を関連するPDB構造をフィルタリングして絞り込むために使いたいので、データセットを辞書からより扱いやすいデータフレームに変換します。
pdbs = pd.DataFrame(pdbs) pdbs.head()
| relatedPDB | structureId | title | pubmedId | expMethod | resolution | keywords | nr_entities | nr_residues | nr_atoms | deposition_date | release_date | last_modification_date | structure_authors | citation_authors | status | pubmedCentralId | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | {'@pdbId': '1JL9', '@details': '1JL9 contains ... | 1IVO | Crystal Structure of the Complex of Human Epid... | 12297050 | X-RAY DIFFRACTION | 3.30 | TRANSFERASE/SIGNALING PROTEIN | 2 | 1350 | 8813 | 2002-03-28 | 2002-10-16 | 2011-07-13 | Ogiso, H., Ishitani, R., Nureki, O., Fukai, S.... | Ogiso, H., Ishitani, R., Nureki, O., Fukai, S.... | CURRENT | NaN |
| 1 | {'@pdbId': '1M17', '@details': 'Epidermal Grow... | 1M14 | Tyrosine Kinase Domain from Epidermal Growth F... | 12196540 | X-RAY DIFFRACTION | 2.60 | TRANSFERASE | 1 | 333 | 2452 | 2002-06-17 | 2002-09-04 | 2011-07-13 | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | CURRENT | NaN |
| 2 | {'@pdbId': '1M14', '@details': 'Apo-form Epide... | 1M17 | Epidermal Growth Factor Receptor tyrosine kina... | 12196540 | X-RAY DIFFRACTION | 2.60 | TRANSFERASE | 1 | 333 | 2540 | 2002-06-17 | 2002-09-04 | 2011-07-13 | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | Stamos, J., Sliwkowski, M.X., Eigenbrot, C. | CURRENT | NaN |
| 3 | [{'@pdbId': '1IGR', '@details': '1IGR contains... | 1MOX | Crystal Structure of Human Epidermal Growth Fa... | 12297049 | X-RAY DIFFRACTION | 2.50 | transferase/growth factor | 2 | 1102 | 8607 | 2002-09-10 | 2003-09-10 | 2011-07-13 | Garrett, T.P.J., McKern, N.M., Lou, M., Ellema... | Garrett, T.P.J., McKern, N.M., Lou, M., Ellema... | CURRENT | NaN |
| 4 | NaN | 1XKK | EGFR kinase domain complexed with a quinazolin... | 15374980 | X-RAY DIFFRACTION | 2.40 | TRANSFERASE | 1 | 352 | 2299 | 2004-09-29 | 2004-12-07 | 2011-07-13 | Wood, E.R., Truesdale, A.T., McDonald, O.B., Y... | Wood, E.R., Truesdale, A.T., McDonald, O.B., Y... | CURRENT | NaN |
print("Number of PDB structures for EGFR: " + str(len(pdbs))) # Number of PDB structures for EGFR: 205
データセットを以下のクライテリアに沿ってフィルタリングしていきます。
1. 実験手法：X線回折（X-ray diffraction）
X線回折（X-RAY DIFFRACTION）によって解かれた構造のみを残します。X線回折は最もよく使われている構造決定手法です。
pdbs = pdbs[pdbs.expMethod =="X-RAY DIFFRACTION"] print("Number of PDB structures for EGFR from X-ray: " + str(len(pdbs))) # Number of PDB structures for EGFR from X-ray: 199
2. 構造分解能
分解能が3Å（Angström）以下の構造のみを残します。分解能の値が小さいほど、構造の質が良くなります（＝質が良いとは、原子に割り当てられた３次元座標の確からしさが高くなるということです）。3Åよりも下の場合、原子レベルの配置を決めることができ、したがって構造に基づく医薬品デザイン（structure-based drug design）に関連する構造を選択するための閾値として使われることがよくあります。
pdbs_resolution = [float(i) for i in pdbs.resolution.tolist()] pdbs = pdbs[[i <= 3.0 for i in pdbs_resolution]] print("Number of PDB structures for EGFR from X-ray with resolution <= 3.0 Angström: " + str(len(pdbs))) # Number of PDB structures for EGFR from X-ray with resolution <= 3.0 Angström: 164
データセットを構造の分解能で並べ替えます。
pdbs = pdbs.sort_values(["resolution"], ascending=True, na_position='last')
（分解能でソートした）上位のPDB構造を確認します。
pdbs.head()[["structureId", "resolution"]]
| structureId | resolution | |
|---|---|---|
| 153 | 5UG9 | 1.33 |
| 141 | 5HG8 | 1.42 |
| 152 | 5UG8 | 1.46 |
| 50 | 3POZ | 1.50 |
| 56 | 3VRP | 1.52 |
3. リガンドの結合した構造
次のトークトリアルでリガンドベースの組み合わせファーマコフォア（ensemble ligand-based pharmacophore）を作成するので、データフレームから結合するリガンドを含まないPDB構造を取り除きます。pypdb関数のget_ligandsを使うことでPDB構造に含まれるリガンドの確認と取得ができます。PDBでリガンドと分類されているものにはリガンド、補因子だけでなく溶媒とイオンも含まれています。リガンドの結合した構造だけに絞り込むため、(i) リガンドに分類されている物を含まないすべての構造を取り除き、そして(ii) 分子量（MW）100Da（ダルトン）以上のリガンドを含まない構造をすべて取り除きます。分子量100Daとするのは多くの溶媒とイオンの分子量がそれよりも小さいからです。注：この方法は簡単ですが、全ての溶媒とイオンを除去する方法ではありません。
# データフレームからPDB IDを取得 # pdb_ids = pdbs["structureId"].get_values().tolist() # 訳注: オリジナル資料の上記コードはpandas version 0.25.0以降では使用できなくなっています。 pdb_ids = pdbs["structureId"].to_numpy().tolist()
# 構造の除去 # (i) リガンドを持たないもの # (ii) 分子量(MW)が100 Da(Dalton)以上のリガンドを持たないもの mw_cutoff = 100.0 # 分子量のカットオフ値（単位Da） # このデータベースクエリはしばらく時間がかかるかもしれません removed_pdb_ids = [] for i in pdb_ids: ligand_dict = get_ligands(i) # (i) リガンドがない場合に構造を取り除く if ligand_dict["ligandInfo"] is None: pdb_ids.remove(i) # リストからリガンドのないPDB IDを取り除く removed_pdb_ids.append(i) # リガンドのないPDB IDを格納 # (ii) mw_cutoff以上の分子量のリガンドを一つも持たない構造の除去 else: # リガンド情報の取得 ligs = ligand_dict["ligandInfo"]["ligand"] # 専門的な事項：リガンドが一つの時、辞書型からリスト型へキャストします（続けてリスト内包表記処理するため） if type(ligs) == dict: ligs = [ligs] # リガンドと分類されているもの毎に分子量を取得 mw_list = [float(i["@molecularWeight"]) for i in ligs] # 分子量（MW）100Da（ダルトン）以上のリガンドを含まない構造を除去 if sum([mw > mw_cutoff for mw in mw_list]) == 0: pdb_ids.remove(i) # リストからリガンドのないPDB IDを取り除く removed_pdb_ids.append(i) # リガンドのないPDB IDを格納 print("PDB structures without a ligand (removed from our data set):") print(removed_pdb_ids) # PDB structures without a ligand (removed from our data set): # ['2EB2', '1M14', '2GS2', '2RFE', '4I1Z']
print("Number of structures with ligand: " + str(len(pdb_ids))) # Number of structures with ligand: 159
上位の構造からリガンドのメタ情報を取得
次のトークトリアルでは、最も分解能の高い上位のtop_num構造からリガンドベースの組み合わせファーマコフォアを作成します。
top_num = 4 # 上位の構造の数 pdb_ids = pdb_ids[0:top_num] pdb_ids # ['5UG9', '5HG8', '5UG8', '3POZ']
get_ligandsを使って上位top_num個のリガンドについてPDB情報を取得し、csv ファイルとして保存します（リガンド毎の辞書）。
構造が複数のリガンドを含んでいる場合、最も大きなリガンドを選びます。注：これは単純な方法で、タンパク質の結合サイトに結合するリガンドを選択するうえで完璧な方法ではありません。このアプローチではタンパク質に結合する補因子を選んでしまう可能性もあります。ですので、以降のステップで使用する前に、PyMolで自動的に選択した上位のリガンドをチェックしてください。
ligands_list = [] for i in pdb_ids: ligands = get_ligands(i)["ligandInfo"]["ligand"] # 専門的な事項：リガンドが一つの時、辞書型からリスト型へキャストします（続けてリスト内包表記処理するため） if type(ligands) == dict: ligands = [ligands] weight = 0 this_lig = {} # 幾つかリガンドが含まれている場合、最大のものを取得 for lig in ligands: if float(lig["@molecularWeight"]) > weight: this_lig = lig weight = float(lig["@molecularWeight"]) ligands_list.append(this_lig) # フォーマットをデータフレームに変更 ligs = pd.DataFrame(ligands_list) ligs
| @structureId | @chemicalID | @type | @molecularWeight | chemicalName | formula | InChI | InChIKey | smiles | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 5UG9 | 8AM | non-polymer | 445.494 | N-[(3R,4R)-4-fluoro-1-{6-[(3-methoxy-1-methyl-... | C20 H28 F N9 O2 | InChI=1S/C20H28FN9O2/c1-6-15(31)23-13-9-29(7-1... | MJLFLAORJNTNDV-CHWSQXEVSA-N | CCC(=O)N[C@@H]1CN(C[C@H]1F)c2nc(c3c(n2)n(cn3)C... |
| 1 | 5HG8 | 634 | non-polymer | 377.4 | N-[3-({2-[(1-methyl-1H-pyrazol-4-yl)amino]-7H-... | C19 H19 N7 O2 | InChI=1S/C19H19N7O2/c1-3-16(27)22-12-5-4-6-14(... | YWNHZBNRKJYHTR-UHFFFAOYSA-N | CCC(=O)Nc1cccc(c1)Oc2c3cc[nH]c3nc(n2)Nc4cnn(c4)C |
| 2 | 5UG8 | 8BP | non-polymer | 415.468 | N-[(3R,4R)-4-fluoro-1-{6-[(1-methyl-1H-pyrazol... | C19 H26 F N9 O | InChI=1S/C19H26FN9O/c1-5-15(30)24-14-9-28(8-13... | CGULPICMFDDQRH-ZIAGYGMSSA-N | CCC(=O)N[C@@H]1CN(C[C@H]1F)c2nc(c3c(n2)n(cn3)C... |
| 3 | 3POZ | 03P | non-polymer | 547.957 | N-{2-[4-({3-chloro-4-[3-(trifluoromethyl)pheno... | C26 H25 Cl F3 N5 O3 | InChI=1S/C26H25ClF3N5O3/c1-25(2,37)14-22(36)31... | ZYQXEVJIFYIBHZ-UHFFFAOYSA-N | CC(C)(CC(=O)NCCn1ccc2c1c(ncn2)Nc3ccc(c(c3)Cl)O... |
ligs.to_csv('../data/T8/PDB_top_ligands.csv', header=True, index=False, sep='\t')
上位リガンド分子の描画
PandasTools.AddMoleculeColumnToFrame(ligs, 'smiles') Draw.MolsToGridImage(mols=list(ligs.ROMol), legends=list(ligs['@chemicalID']+', '+ligs['@structureId']), molsPerRow=top_num)
タンパク質ーリガンド IDペアの作成
pairs = collections.OrderedDict() for idx, row in ligs.iterrows(): pairs[str(row['@structureId'])] = str(row['@chemicalID']) print(pairs) # OrderedDict([('5UG9', '8AM'), ('5HG8', '634'), ('5UG8', '8BP'), ('3POZ', '03P')])
訳注(2020/05)
OrderedDictは順序つきの辞書です。ここではタンパク質のPDB IDとリガンドのIDをそれぞれDataFrameの@structureId列と@chemicalID列から取得し、keyとvalueとしています。
訳注ここまで
PDB構造ファイルの取得
では、pypdbのget_pdb_file関数を使ってPDBからPDB構造ファイルを取得します。すなわちタンパク質、リガンド（そして可能な場合には補因子や水分子、そしてイオンといった他の原子あるいは分子）の３次元座標を取得します。
利用可能なファイルフォーマットはpdbとcifで、タンパク質（とリガンド、補因子、水分子、イオン）の原子の３次元座標と、原子間の結合に関する情報とが格納されています。ここではpdbファイルを使用します。
# PDBファイルを取ってきてローカルに保存する for prot, lig in pairs.items(): pdb_file = get_pdb_file(prot, filetype='pdb', compression=False) with open('../data/T8/'+ prot + '.pdb', 'w') as f: f.write(pdb_file)
PDB構造の重ね合わせ
次のトークトリアルでリガンドベースの組み合わせファーマコフォアを作成したいので、全ての構造を互いに３次元で重ね合わせする（アラインメントを取る）必要があります。このタスクのため、分子描画プログラムPyMolを使用します。PyMolはJupyter notebookの中でも使うことができます。PyMolは２つの構造間の原子の距離が最小になるように、一度に２つの構造の重ね合わせをします。
まず最初にコマンドラインからPyMolを立ち上げます（quietモード、すなわちGUIが開かないモード）。
# PyMolをquietモードで立ち上げる pymol.pymol_argv = ['pymol','-qc'] pymol.finish_launching()
重ね合わせのために、リファレンスとなる構造のファイル（動かさないPDB、'target'）を選択し、それに対してもう一方の（'query'）構造ファイルを、PyMolのコマンドcmd.align(query, target)を使って重ね合わせます。全ての cmd. コマンドはPyMolのコマンドです。
重ね合わせた構造を新しい座標とともにpdb ファイルとして保存します。また、構造ファイルからリガンドを取り出して別のpdbファイルとして保存し、次のトークトリアルで使用します。
# アラインメントのログをファイルに保存 f = open("../data/T8/alignments.log", "w") # アラインメントの間、動かさない構造と動かす構造とを区別するための変数 immobile_pdb = True refAlignTarget='non' refAlignQuery='non' # 最初のタンパク質に対してタンパク質のアラインメントを取る for prot, lig in pairs.items(): # 動かさない構造（アラインメントのリファレンス構造） if immobile_pdb: target = prot f.write('Immobile target: ' + prot + '\n') # pdbファイルの読み込み（タンパク質とリガンドの複合体） targetFile = cmd.load('../data/T8/' + target + '.pdb') # 精密化したアラインメントのための名前を格納 refAlignTarget='('+target+' within 5 of resn '+lig+')' # 複合体をpdbファイルとして保存 cmd.save('../data/T8/' + target + '_algn.pdb', selection=target) # 選択した名前のリガンドだけを選択 ligObj = cmd.select('ligand', target + ' and resn ' + lig) # 選択したものをpdbファイルとして保存 cmd.save('../data/T8/' + target + '_lig.pdb', selection='ligand', format='pdb') # リガンドの選択を削除 cmd.delete(ligObj) # targetが選択されました。 immobile_pdb = False # 動かす構造（リファレンス構造に対して重ね合わせる） else: query = prot f.write('-- align %s to %s \n' %(query, target)) # pdbファイルの読み込み（タンパク質のリガンドの複合体） queryFile = cmd.load('../data/T8/' + query + '.pdb') # 結合サイトに焦点をあてて構造（タンパク質）を重ね合わせる refAlignQuery= '('+query+' within 5 of resn '+lig+')' values = cmd.align(refAlignQuery, refAlignTarget) # 構造の重ね合わせができない場合（即ち、RMSD > 5A の場合）、スキップする if values[0] > 5: f.write('--- bad alignment: skip structure\n') else: # 複合体をpdbファイルとして保存 cmd.save('../data/T8/' + query + '_algn.pdb', selection=query) # リガンドのみを選択 ligObj = cmd.select('ligand', query + ' and resn ' + lig) # 選択したものをpdbファイルとして保存 cmd.save('../data/T8/' + query + '_lig.pdb', selection='ligand', format='pdb') # リガンドの選択を削除 cmd.delete(ligObj) # "query"の選択を削除 cmd.delete(query) # "target"の選択を削除 cmd.delete(target) f.close()
訳注(2020/05)
Pymol Wikiを参照して上のセルの補足をします。
アラインメントのコマンドcmd.alignに渡す引数のrefAlignTargetとrefAlignQueryはどちらも( PDB_ID within 5 resn of Ligand_ID)という文字列になる変数です。これは「タンパク質（PDB_ID）のうち、残基名(resn, residune name)で指定されたリガンド（Ligand_ID）の5Å以内に入る部分」を指し示すようです。リガンドの周囲を指定してアラインメントをとっているので「結合サイトに焦点をあてて重ね合わせる」ことになるようです。
また、cmd.alignの戻り値は7つあり、一番最初（index 0）が「精密化した後のRMSD」となっているため、if values[0] > 5という条件文になっているようです（PyMol Wiki: Align）。
訳注ここまで
PyMolアプリケーションを終了します。
# PyMolを終了
pymol.cmd.quit()
全てのリガンドpdbファイルが存在していることを確認します。もしファイルが欠けていたなら、手作業でPyMOLでタンパク質ーリガンド複合体構造をチェックしてください。
mol_files = [] for file in glob.glob("../data/T8/*_lig.pdb"): mol_files.append(file) mol_files # ['../data/T8/5UG8_lig.pdb', # '../data/T8/3POZ_lig.pdb', # '../data/T8/5UG9_lig.pdb', # '../data/T8/5HG8_lig.pdb']
ディスカッション
このトークトリアルでは、PDBからタンパク質とリガンドのメタ情報と構造情報を取得する方法を学びました。X線構造だけを残し、分解能とリガンドがあるかどうかによって絞りこみました。最終的な目標は、重ね合わせたリガンドのセットを使って、次のトークトリアルでリガンドベースの組み合わせファーマコフォアを作成することです。
ファーマコフォアモデリングのためにリガンドについての情報量を多くするため、PDB構造を分解能でフィルタリングするだけでなく、リガンドの多様性をチェックすること（類似性による分子クラスタリングについてはトークトリアル 5を参照）と、リガンドの生理活性をチェックする（即ち、活性があるリガンドだけを含める）ことをお勧めします。
クイズ
- Protein Data Bankの含むデータの種類を要約してください。
- 構造の分解能が何を表し、このトークトリアルではどうやって、また、なぜ分解能についてフィルタリングしたのかを説明してください。
- 構造の重ね合わせ（アラインメントを取ること）とはどういう意味か説明し、このトークとリアルで実施したアラインメントについて議論してください。
訳以上
追記
今回も訳出、用語の選択に誤りがあるかもしれません。ご指摘いただければ幸いです。
今回、 取得したPDBの構造を絞り込むうえで分解能が用いられていましたが、登録されている構造の正確さについては色々と注意が必要だそうです。
@torusengoku先生の記事： X線結晶構造中の低分子や水の妥当性評価 で驚きの事例を紹介してくださっています。
PDBからの複数の構造について、KNIMEを活用する方法を 或る化みす途のブログ さんが紹介してくださっています。
* Dear Pymol Beginner x KNIME -KNIMEを使ってPDBファイルをまとめてDLしてみたり-
PyMol 機能ありすぎてよくわからないのですが、とりあえず pymol-book が素晴らしいので読んでみようと思います（知らなかった・・・）。
