GW読書感想文

GW10連休、皆様いかがお過ごしでしょうか？充実したお休みでしたでしょうか？ GW期間中もお仕事だった皆様はお疲れ様です。ありがとうございました。

さて、とりあえず創薬レイドバトルがひと段落ついたので関連記事をGithubにまとめなおしました。レポジトリはこちらです。ノートブック形式ですが中身はMarkdownというよくわからないファイルですがご笑覧いただければ幸いです。

で、やはりどうにも自分には基本的な知識が足りていない（何がわからないかもわからない状態）ということが分かったので、とりあえず専門用語を拾うべくいくつか本を読んでみました。以下、とりとめもない備忘録です。

1冊目：コンピュータシステムの理論と実装
2冊目: StanとRでベイズ統計モデリング、3冊目: Pythonで体験するベイズ推論
４冊目：はじめてのパターン認識
５冊目：カーネル多変量解析
６冊目：岩波データサイエンス vol.5 スパースモデリングと多変量データ解析
まとめ

1冊目：コンピュータシステムの理論と実装

コンピュータシステムの理論と実装 ―モダンなコンピュータの作り方

作者: Noam Nisan,Shimon Schocken,斎藤康毅
出版社/メーカー: オライリージャパン
発売日: 2015/03/25
メディア: 単行本（ソフトカバー）
この商品を含むブログ (4件) を見る

あれ？そもそもパソコンってどうやって動いているんだろう？？みたいなレベルなので、とりあえずコンピューターサイエンスっぽい本を読んでみました。

最も基本的なユニット（電子ゲート）からはじめて、最終的にコンピュータシステムをつくるに至る過程を、一つずつ段階を踏みながら実践していこう！という内容です。「イントロダクション」に引用されている言葉がこの本のスタンスを最も明確に示していると思うので引用します。

「人に重要な影響を与える唯一の学びは、自己発見・自己本位による学び、つまり、経験によって理解された真実のみである」 by カール・ロジャーズ（心理学者）

実際に作り上げるという経験を通して「コンピュータの動く仕組み」を理解しよう、というのが本書の目的ですが、では具体的にどうすればよいのか？著者は「抽象と実装」を繰り返すことで可能になると述べています。複雑なプロジェクトを成し遂げるためには、「抽象化」によりモジュールに分割し、それを組み立て直すというアプローチを取るべきだというのがその主張です。

「どのようにおこなうか（how）」という詳細を無視して、「この要素が何をするのか（what）」だけを考えて抽象化すること。この段階が重要かつ実例に触れて学ぶべきものなので、「抽象化」は著者が解説・提供し、その実装の経験を読者が積む、という構成になっています。そのために本書では解くべき課題とともに、導入の容易な実行環境、手順が提供されています。

特に私のような初学者にとって親切だと感じた点は、本書全体としては「ボトムアップでコンピュータをつくりあげる」という構成になっているものの、最初のイントロダクションで「これから何をおこなおうとしているのか？」トップダウン的にプロジェクト全体の俯瞰を提供してくれている点です。読み進めていくうちに何度も「あれ？いま何やってんすかね？？」と迷子になりかけてましたが、その度にイントロに戻ることで全体の中での位置付けを確認することができました。

と、まるで理解したかのように書いていますが私は１章の論理回路で実装をあきらめました（知的体力が無さすぎる・・・）

とりあえず本書の中の図を切りはりすると以下のような話です。たぶん・・・

f:id:magattaca:20190506232548p:plain

NANDからスタートして作って行く話が、パーセプトロンの組み合わせと重み変えていろんな関数表現できるよ、みたいなDeepなお話に繋がっている気がしてちょっと賢くなった気分。（違うか・・・）

2冊目: StanとRでベイズ統計モデリング、3冊目: Pythonで体験するベイズ推論

StanとRでベイズ統計モデリング (Wonderful R)

作者: 松浦健太郎,石田基広
出版社/メーカー: 共立出版
発売日: 2016/10/25
メディア: 単行本
この商品を含むブログ (10件) を見る

Pythonで体験するベイズ推論 PyMCによるMCMC入門

作者: キャメロンデビッドソン=ピロン,玉木徹
出版社/メーカー: 森北出版
発売日: 2017/04/06
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

お正月に「データ解析のための統計モデリング入門」という本をよみました（緑本？）。統計モデリングについて、簡単なモデルから徐々に複雑なモデルへと、非常にわかりやすく解説されており、各所でオススメされている素敵な本なのですが、後半、ベイズ統計モデルの話になったところで私は理解があやふやになりました。

モデルが複雑化し、パラメータが増えるに従ってパラメータの推定が困難になる。（なるほど）
多数のパラメータの推定方法としてMCMCサンプリングがある。（・・・ほう）
MCMCサンプリングで得られる定常分布はパラメータの確率分布（・・・はあ）
だが頻度主義では「真の値」が１つあるという考えだから分布はあり得ない（・・・それはまた）
でもベイズ統計モデルは確率分布としてあつかう枠組み（・・・ふーん）
だから、今考えている統計モデルはベイズ統計モデルだったと見直そう（・・・・・・えっ？）

っていう感じの流れでした。４番目くらいで頭がついていかなくなりました。（３段階が思考の限界なのです。）　　

しかしベイズという名前はあちこちで見かけるので、やはりお勉強しなければならないのではないか、ということで２冊ほどパラパラめくってみました。

今回読んだ２冊はいずれも実践的な面を重視しており、ベイズモデリングのためのソフトウェアとしてそれぞれ「RとStan」、「PythonとPyMC」が使われています。コード例が多数記載されていて、特に「グラフの描画」がこまめに挿入されている点が実際の解析過程をなぞっているようで面白いと思いました。

が、しかし、ここまで多数の例を懇切丁寧に紹介されても私の理解力ではしっくりこない。ぼんやりとネットをうろうろしていたところ、以下の解説記事を見つけました。

作って遊ぶ機械学習　ベイズ推論による機械学習の基本

下の本の著者の方が書かれているブログのようです。（著書は未読です。すみません。ちらっと見たら数式が多そうで手がでませんでした。）

機械学習スタートアップシリーズベイズ推論による機械学習入門 (KS情報科学専門書)

作者: 須山敦志,杉山将
出版社/メーカー: 講談社
発売日: 2017/10/21
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

上記のブログ記事の「ベイズ学習における基本的な流れ」を読み、やっと書籍で行われていたコードの流れが少しわかるようになりました。引用させていただきます。

尤度関数（観測データの表現方法）を定義する
事前分布（パラメータの散らばり具合）を設定する。
ベイズの定理を使って事後分布を更新する。（=学習）
事後分布から予測分布を計算し、未知の値の確率分布を求める。（=予測）

こちらを見るまでは、一つ一つのコードを見ているうちに全体として何をやっているか追えなくなる、ということを繰り返していました。*1

流れがわかるようになったところで、改めてベイズ推論についてです。
考え方の導入としては③「Pythonで体験するベイズ推論」の「第１章　ベイズ推論の考え方」が「気持ち」を理解する上でわかりやすかったです。

曰く、ベイズ主義とは「信念」であり、さらに「人によって信念が異なっていることも定義が許している」点が、我々が日常で行なっている推論と類似しており自然な考え方、とのこと。（何を言っているかよくわからないと思うので気になった方は書籍の１章の最初だけでもご覧になっていただければと思います。）

②「StanとRでベイズ統計モデリング」２章のとても簡潔にまとめられた特徴と合わせると以下のような感じです。

f:id:magattaca:20190506233025p:plain

もう一つのキーワードMCMCですが、こちらは事後分布を求めるところで出てきます。複雑なモデルではデータの得られる確率（p（D）)を求める計算の難易度が高いため、事後分布（p(θ|D))に比例する分布 p(D|θ)xp(θ)から乱数サンプルを多数発生させて事後分布の代わりとするそうです。*2

以上のような導入を踏まえて、読者の理解が進むよう、書籍②、③ともに豊富な事例を提供してくださっています。

「実際にモデルを作成するにははどのような分布を選択すれば良いのか？」
「MCMCの短所として、時間がかかったり収束しないということがあるけど、そんな時どうすればよいのか？」

といった話題に踏み込んで解説してくださっていますので、「頭のいい人はそう言う考え方でアプローチするのか」と非常に興味深かったです。

といってもほとんど理解できていないのでもう一回読んだ方が良さそう・・・。

４冊目：はじめてのパターン認識

はじめてのパターン認識

作者: 平井有三
出版社/メーカー: 森北出版
発売日: 2012/07/31
メディア: 単行本（ソフトカバー）
購入: 1人クリック: 7回
この商品を含むブログ (5件) を見る

「機械学習　教科書」で検索すると必ず出てくる「はじパタ」こと「はじめてのパターン認識」・・・。

重い腰をあげて漸く購入したのですが、私には難しすぎました。本当に皆さんここからはじめているのでしょうか？分厚すぎず教科書として無駄がない良い本なのだろうなあ、ということは印象ですが、はじめっからハードル高すぎでは・・・？？？*3

個人的には　第５章　k最近傍法で

ボロノイ図がでてきた（「CASTpの理論で見たやつだ！」）ところと
次元の呪いの説明（「これがサクサクメロンパンの呪いか！」）

がピークでした。

ただ、分からないなりに読み進めていく中で、雰囲気として「統計モデリングのお話と機械学習ではちょっと重点のおいている場所が違う」というのを感じました。

先に読んだベイズ関連の書籍では、背後にあるモデルを意識していることを前提としていました。一方で「パターン認識」では「今あるデータを分析し、新しいデータの予測を当てること」に重点があり、なんなら「現実の理由がわからないならそれでもいいんじゃない」くらいの雰囲気を感じました。

これも数式のもつ意味合いを理解できない私の力不足のせいでしょう。いつか理解できる日が来るのでしょうか？？？

５冊目：カーネル多変量解析

カーネル多変量解析―非線形データ解析の新しい展開 (シリーズ確率と情報の科学)

作者: 赤穂昭太郎
出版社/メーカー: 岩波書店
発売日: 2008/11/27
メディア: 単行本
購入: 7人クリック: 180回
この商品を含むブログ (32件) を見る

はじパタ「第８章 サポートベクトルマシン」で「カーネルトリック」という言葉を目にし、「なんか知らんけど響きが格好いい！！」となったので購入しました。（安直）

結論としてこちらの本、門外漢にとってはとても楽しかったのです。なぜかというと全７章のうち、最初の４章が「カーネル法すごいぞ！」「こんなことできてしまうぞ！」という「ご利益」の説明に重点がおかれており、証明や理由の込み入った部分については適宜飛ばしつつでも読める、というように書かれているからです。

もちろん残りの後半の章で、前半で触れられなかった数学的な話が展開されています（・・・つまり私は後半を飛ばした）。

第１章ではまず現代の多変量解析においてカーネルがどのような役割を果たすか、その位置づけが語られます。

そもそもカーネル法とは、「カーネル関数の重み付きの和で表したモデルを正則化付きで最適化するデータ解析手法」と言えるそうです。
カーネル法のご利益は以下のような感じです。

線形モデルと異なり入力データとして非線形なもの（文字列やグラフ構造でもOK）を取ることができる
「カーネル関数の重み付きの和」を扱うので決めるべきパラメータ自体には線形性がある
サンプル数（入力データ）と同じだけの自由度があるため複雑な関数を表現できる ⇄　汎化能力が問題となる（正則化で解決を試みる）
カーネルの計算とその後の処理を分離できるので、計算をモジュール化できる。

次いで第２章では、カーネル法で使われるカーネル関数はどのようなもので、どうやって導入することができるのか、複数の観点から定義が試みられています。
注目の「カーネルトリック」は、「正定値性からの導入」で取り上げられています。

曰く、「正定値性をカーネル関数の定義とすれば」、「内積計算が不要」となり、「計算量をぐっと減らすことができる」そうな・・・。　　この一つ前に語られている「特徴抽出からの導入」ではカーネル関数が「特徴ベクトルΦ(x)どうしの内積」として定義されていました。
これに対して「正定値性を満たすことがわかっている関数」を利用すれば、内積計算で導く必要なくカーネル関数を用意できることになるそうです。（・・・？？？）*4

ここまでの例示は回帰の話題が取り扱われているのですが、第３章、第４章と、より複雑な問題へ話題が移ります。
第３章では固有値問題の例として「カーネル主成分分析」や「カーネル正準相関分析」が、
第４章では凸計画問題の例として「サポートベクトルマシン」や「スパース性」について触れられています。

私はそもそも「カーネル」が接頭語としてついていない「主成分分析」「正準相関分析」もあやふやな段階だったので、「そんなお話もあるのね」程度にしかわかりませんでしたが、だいたい第３章の手順としては、繰り返し

リプレゼンター定理　→　正則化項の導入
ラグランジュの未定乗数法　→　最適化問題

が行われていました。

第４章 SVMの解説はとても面白く、「はじパタ」ではいまいち理解できなかった部分についてもう少し追うことができました。
まず、そもそもの出発点として、クラス分類の識別関数の誤差に対するペナルティの考慮からスタートしています。

問題意識として、

損失関数に２乗誤差を用いると外れ値の影響が大きくなる（ロバスト性）
サンプル数が多すぎる場合計算量が大きくなりすぎる（スパース性）

という点があげられています。これに対して

「最小２乗ではない物」を使いたい　＋　「凸関数」をつかいたい、と言うことで損失関数に「区分線形損失」を導入
これを解いた解は、サンプルの一部分しか使わないため（スパース）、計算量が節約できる

という解説がなされています。より具体的には、

カーネル関数と２次の正則化（リプレゼンター定理）により凸二次計画問題（QP）へ
制約つきの２次関数の最小問題 → ラグランジュの未定乗数法・KKT条件
相対問題による計算の単純化

といった流れです。この中で、

KKT条件の相補性条件から解のスパース性（サポートベクトルが全体のサンプルに比べてスパースになる）

が導かれるという説明がなされます。書き写しておきながら全然理解していませんが、未来の自分のためのメモです。ご容赦ください。

「それでは、具体的なカーネル関数の中身はどのようなものなのか？」「どう設計すれば良いのか？」という話題は第５章以降に展開されていますが、またしても私の能力では追いきれませんでした。

６冊目：岩波データサイエンス vol.5 スパースモデリングと多変量データ解析

岩波データサイエンス Vol.5

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2017/02/16
メディア: 単行本（ソフトカバー）
この商品を含むブログ (3件) を見る

GW読書感想文、最後ははじめての「岩波データサイエンス」です。これはもう単純に「スパース性」の話題を知りたかった、というだけです。

結論・・・「えるわん」

なんだかすごくL1でlassoでした。めっちゃ０になるらしいです。

明治大学　データ化学工学研究室　金子先生のページに、「Lassoの回帰係数がなぜ0になりやすいのか？」について解説が記載されていました。「係数が0になりやすい」というのを目にし、吝嗇家としては「もったいない・・・」となりましたが、見方を変えて「特徴選択」という観点からは、「0になることにこそ価値がある」という話になるようです。

その発想はなかった・・・。

フィンガープリントを計算して見た際に、「立っているビットの数が少ない」という印象をうけ、「なんだかスパースを考えないといけないのでは！！！」というのが興味をもったきっかけだったのですが、とりあえずL1正則化項を入れれば解決、ということでしょうか？？？

まとめ

以上、いつにもましてとりとめもない備忘録でした。理解できないとはわかっていましたが、ここまでわからないとは・・・。そもそももっと腰を落ち着けて時間をかけて読むべき本ばかりなのでしょうが、如何せん計画を立てて勉強できない自分の性格をなんとかしないといけなさそうです。

カリキュラムがある大学はいいなあ。

*1:薄々感じていましたがどうやら私は物事の全体像、大きな流れを把握する能力がとても低いようです！まさしく先の「コンピュータシステムの理論と実装」で言及されていた「抽象化」の能力が欠けてるってやつですね！

*2:※ p(θ|D) = p(D|θ) x p(θ) / p(D) を想定

*3:そもそも「添字のある式は読めない」のに機械学習の理論を勉強しようと言うのが間違ってる？？？

*4:例としてはガウスカーネルが挙げられています。SVMみかけるRBFカーネルもこれのことでしょうか？？

2019-04-07

スコア関数を変えてドッキング結果を再解析する話

さて、前回の記事でドッキングによるSBVSもなんとか終わったので後は化合物IDのリストを提出するだけ・・・、というところですがとても気になる文献を知ってしまいました。

Performance of machine-learning scoring functions in structure-based virtual screening. Wójcikowski, M., Ballester, P.J., and Siedlecki, P., Sci.Rep. 2017(7)46710 pubmed

こちらの文献は前回使用したODDT(Open Drug Dicovery Toolkit)と同一の筆頭著者によるもので、どうやらドッキングのスコアを機械学習を用いて改善したという内容のようです。しかもODDT同様に、このスコア関数（RF-ScoreVS）を利用できるとのこと！！

前回の記事で参照した「タンパク質計算科学基礎と創薬への応用」*1 において、ドッキングのスコア関数は依然として精度に課題がある（実験値との乖離がある）といった趣旨の記述がありました。AutoDock VinaはAutoDockよりも精度が改善している、とはいえ10年ほど前のソフトウェアです。より改善したものがあるのであれば使ってみたい！

（・・・というのは建前でMachine Learningって言って格好つけたかった）

1. RF-Score-VS （Sci.Rep. 2017(7)46710）文献紹介
2. RF-Score-VS を使ってみる
3. アッセイ候補提出ファイルの作成
4. まとめ
5. 総括

1. RF-Score-VS （Sci.Rep. 2017(7)46710）文献紹介

まずは文献を読んでみようと思います。Scientific ReportsはOpen Accessなのが良いですね。(CC BY 4.0)
すでに記憶があやふやなのでスコア関数について少し復習します。

1-1. スコア関数の種類

スコア関数にはどのようなものがあるか？
先の記事で参照した日本化学会情報化学部会誌の文献（タンパク質-リガンドドッキングの現状と課題 *2 によると、大きく３つ

①力場に基づいたスコア関数（force-field-based scoring function）
②経験的なスコア関数 (empirical scoring function)
③ 統計に基づいたスコア関数(knowledge-based scoring function)

があり、さらに近年の動きとして、

④複数のスコア関数を利用するコンセンサススコア(concensus scoring)
⑤機械学習に基づいたスコア関数(machine-learning based scoring function)

が開発されている、とのことでした。

「④ コンセンサススコア」について「タンパク質計算科学」の説明を参照します。
「コンセンサス（consensus、合意）」をとる目的は、単一の手法の結果のみに頼らず、複数のスクリーニング結果を統合することでヒット率の向上、当たり外れの低減を目指すことです。結果を統合する方法として、各化合物に対して複数の手法でスコアを得たのち

積集合: 複数の手法で共通に選ばれた化合物をヒットとする
和集合: 複数の手法で少なくとも一つに選ばれた化合物をヒットとする
rank by scoring: 各手法のスコアを「0~1」の範囲にスケールした上で和をとり、新しいスコアとする
rank by rank: 各手法における化合物の順位（ランク）の平均によってランキングする

といった手法が例に挙げられていました。

では、「⑤ 機械学習に基づいたスコア関数」はどのようなものか？　論文の中身に進みます。

1-2. 既存のスコア関数の問題点

関数の構造上の問題

そもそもなぜ機械学習の導入にいたったか？論文中では従来から用いられているスコア関数を"classical scoring functions (SFs)" （おそらく「①力場」や「②経験的」、「③統計」によるスコア関数のこと？）としており、結合親和性の予測精度の向上に限界が見られていることの理由として、以下の問題点を指摘しています。

コンフォメーションのエントロピーをうまく説明しない
溶媒和によるエネルギーの寄与をうまく説明しない

また、従来のスコア関数は線形回帰モデル(linear regression model)を下敷きにしているため、

大量のタンパク質構造や結合に関する情報を吸収してモデルに取り込むことに限界がある(universal SFの限界)
標的のタンパク質（protein family）に焦点を合わせてスコア関数のパラメーターをチューニングしたくても、回帰モデルが変更できる形では提供されていない

といった問題を指摘しています。・・・難しい。

よくわからないので別の文献のFigureを引用します。（WIREs Comput. Mol. Sci.2015(5)405）*3

f:id:magattaca:20190407180233p:plain

Fig 中に示されているようにClassocal SFでは線型結合の式ですが、Machine-learnig SFでは関数となっています。そもそもモデル化の構造が異なっているということで、異なるパフォーマンスが期待できるよね！ということのようです。

モデル作成上の問題

さらに、classical SFのモデル作成時の問題点として、トレーニングセットとテストセットに重複(overlap)があるということが挙げられています。機械学習であればcross validationにより、より厳密なバリデーションの取り扱いが可能になります。

また、classical SFは結晶構造の情報を元に作成されていますが、これは結合するというプラスのベクトルの事例のみの情報を用いて構築していることになります。一方でスコア関数を適用するVirtual Screeningの目的は、スクリーニング結果として得られた「実際に結合する正例(プラスのベクトル/ +ve）」と「結合しない負例（マイナスのベクトル/ -ve）」が入り混じった状態である一連の結合ポーズの中から、正しいtrue binderを見つけ出すということです。つまり、Virtual Screeningでは負例（inactiveのリガンド / -ve）の比率の方が高いのに、その情報が関数作成段階に考慮されていないと指摘されていました。

（なるほど確かに言われてみればそう・・・なのでしょうか？勉強不足のため十分に理解できていません。）

ということで、機械学習を導入することで、これらの問題点を解決して行きましょう！　

1-3. 著者らによる機械学習スコア関数（RF-Score シリーズ）

文献中では機械学習スコア関数の例としてRF-Score、NNScore、SFCscoreといった名前が挙げられています。このうち、RF-Scoreは今回の文献の著者の一人、De. Pedro Ballester（CRCM）が筆頭著者として2010年に報告しているものです*4。RFはRandom Forestの略のようですが、こちらをversion 1 （RF-Score v1）として、RF-Score v2*5が報告されています。特にRF-Score v3の文献はタイトルが「Improving AutoDock Vina Using Random Forest: The Growing Accuracy of Binding Affinity Prediction by the Effective Exploitation of Larger Data Sets.」、とAutoDock Vinaを使った以上とても気になる文献です（アクセスできなかった・・・）。

これらのスコア関数の提案に対する議論、「overfittingしているのでは？」「新しいターゲットにも適用できるの？(applicablity)」に対して、非常に大きなデータセットを使って学習を実施し、virtual screeningのパフォーマンスが上げられることを示した、というのが2017年のSci. Rep.の概要となります。

1-4. RF-Scoreの仕組み

では、ベースとなったRF-Scoreはどのようなものなのか？ちょうど良いスライド資料があったので参照していきたいと思います*6。（元文献を読む能力がなくてすみません）

RF-Scoreではランダムフォレストを使ってスコア関数をモデル化していますが、元々の発想としては"classical SF"のうちの「3. 統計に基づいたスコア関数(knowledge-based scoring function)」を改善する、ということから始まっているようです。

統計的スコア関数の流れは以下のようになっています。

①多数の結晶構造（PDBのタンパク質-リガンド複合体構造(P-L complex)）を比較
②原子（あるいは官能基/フラグメント）にタイプを割り当てる
③中心の原子と周囲の原子のtype-type間の距離のヒストグラムを作成する（存在確率の空間分布）
④存在確率（P）を結合自由エネルギー（ΔF）に変換する（ΔF = -k_BT lnP）
⑤全てのP-Lアトムペアのエネルギーを合計して全体のエネルギーとする。

上記の「④ヒストグラムのエネルギー関数への変換」では"reverse Boltzmann"の手法が使われています。ここで、「温度Tの平衡状態においてタンパク質とリガンドの各原子は独立な粒子である」という仮定が入りますが、これは仮定としては粗すぎる、というのが問題となるようです。
理由として

分子中の原子と原子は結合している（Molecular connectivity）から、原子と原子の距離は独立とはいえない
排除体積の影響（Excluded volume）
平衡状態の仮定に物理的根拠がない
温度(T)による構造変化は小さく、ボルツマン分布により暗に考慮されているほどではない

といったことが挙げられています。以上をふまえて、著者らは「ヒストグラム（距離の分布 / distance distribution）と結合自由エネルギーとを関係付ける」ために、Boltzmannの公式のような公式を仮定するのではなく、「構造情報と結合親和性の関係性を機械学習を使って学習すれば良い」と考えたのがRF-Scoreの出発点となっているそうです。

上記の考えに基づき、以下の操作を行って作成したモデルがRF-Score (v1) となります。

①PDBの構造情報(PDBbind Benchmark)と質の高い結合親和性データ（K_d, K_i）
②タンパク質とリガンドの、閾値以下(d_cutoff)の距離にある原子をみつける
③アトムタイプペアをフィーチャーとして出現率を求める
④各タンパク質-リガンド複合体をフィーチャーのベクトルとして表す
⑤ベクトルを入力（説明変数）、結合親和性を出力（目的変数）としてランダムフォレストによる機械学習を実施

加えてさらに、スコア関数の精度向上にむけて公共データと様々な製薬企業の企業内データを用いてスコア関数の精度向上を目指したコンソーシアム（Scoring Function Consrtium）のデータに基づきフィーチャーの選択を改良した（RF-Score v2 / Data-driven feature selection）や、AutoDock Vinaのフィーチャーを取り入れる(RF-Score v3)、といった改良を加えてモデルの精度を高めてきた、というのが流れとなります。

1-5. RF-Score-VSでの取り組み

先に述べたようにRF-Score-VSでは、学習のデータセットに工夫を加えることでVirtual Screeningにおけるスコアの改善が目指されています。

データセットの工夫

まず、データセットとしてDUD-E(A Database of Useful Decoys: Enganced)データセットが使われています。DUD-E*7はスクリーニングソフトの性能テスト（ベンチマーク）として利用されているDUDの改良版であり、DUDは市販化合物の３次元データベースであるZINCのデータを既知の標的タンパク質ごとに、その既知の活性化合物（active）と、類似した化合物（inacitive）との組み合わせとして整理したデータセットであるとのことです*8。

先述のように、RF-Scoreの学習セットとしてはPDBの情報が用いられていました。これに対してRF-Score-VSの利用している学習セットには

より大きなデータセット
結晶構造ではなく、active リガンド(+ve)とinactive リガンド(-ve)のドッキングポーズ

を利用するという違いがあります。結晶構造はactiveリガンドのポーズしか含まないので、こちらを学習セットとして用いるよりも、activeとinactive両方の結合ポーズを学習データとして用いる方がVirtual Screeningに用いるスコア関数としてはより実際の使用状況に類似しているだろう、というのが背景にあるようです。

学習データ分割方法の工夫

また、興味深いことに学習データセットのcross-validationを行うにあたって、

①標的タンパク質ごと（per-terget）
②訓練データ、テストデータともに全ての種類のタンパクを含む（horizontal split）
③訓練データとテストデータで、標的タンパク質が重複しないようにした上で分割する（vertival split）

という３種類の分割方法が用いられています。

論文中の図を引用します。

f:id:magattaca:20190407181127p:plain

②horizontal splitは「すでに活性化合物（active ligand）が知られている標的についてvirtual screeningを行う」という状況を、③vertical splitは「既知の活性化合物がない状態でvirutal screeningを行う」という状況を模倣する、というのを目的としているそうです。

Virtual Screeningが行われる状況の場合分けを想定して、訓練データとテストデータの分割方法を分けるというのは、なるほど〜という感じでとても面白いと思いました（・・・よく使われる方法なのでしょうか？）。

1-6. RF-Score-VSの性能

評価指標

評価指標（metrics）としてはEnrichment factor(EF)というものが用いられています。

$$ EF_{x\%} = \frac{\frac{ActiveMolecules(x\%)}{Molecules(x\%)}}{\frac{ActiveMolecules(100\%)}{Molecules(100\%)}} $$

ランク付けした化合物の上位x%を取り出して、その中に本当の活性化合物が含まれている割合を求めます。この割合と、全化合物（x = 100）中に含まれる活性化合物の割合の比をとったものがEFとなります。Virtual Screeningの性能指標としてよく用いられるものには、EF以外にROC（Receiver Operating Characteristc）曲線のAUC（Area Under the Curve）があるそうです。
EFを用いた理由として、Virtual Screeningを行って選んだ化合物を実際にアッセイするとしたら、上位数パーセントを選んで行うことになる。その数パーセントに"当たり"が含まれていることが大事なので、その状況をより反映しているEFの方が適している、といった説明がなされていました。

確かに創薬レイドバトルも「スコアでランキングをつけて上位10個をアッセイ候補として提出」ということになっていました。この論文とても勉強になる。

性能

ではAutoDock Vinaをドッキングに用い、得られた結合ポーズを各種スコア関数でスコアリングして比較した結果を引用します。横軸が Enrichment Factor 1% （EF_1%）で、縦軸にスコア関数が並んでいます。

f:id:magattaca:20190407181224p:plain

分割方法として② horizontal splitを用いた、RF-Score-VSに置いて非常に良い結果が得られており、従来のclassical SFと比較して高いEF1%となっています。① per-targetは平均値としては、性能が向上しているものの、標的によって性能のばらつきが大きいという結果になっています。また、③ vertical splitは ② horizontalよりも性能は落ちていますが、classical SFよりも性能が向上していると記載されていました。

AutoDock Vinaのスコア関数を改善したという「RF-Score v3」が上記の結果では「AutoDock Vina Score」と比べて良くなっているように見えないのが気になる点ですが、評価対象のデータセットの質が変わった（結晶構造　→　ドッキングポーズ）ことが理由なのかもしれません。

他のドッキングソフト(Docking engine)を使った比較結果、DUD-Eとは異なるデータセット DEKOIS 2.0 に対して適用した際の結果についても記載されていましたので、ご興味のある方は論文を参照していただければと思います。

要するに、とにかく良い結果が出るっぽい！！！　使いましょう！

2. RF-Score-VS を使ってみる

2-1. ドッキングスコアを再計算

とても有り難いことにRF-Score-VSに関して、データセット（ http://github.com/oddt/rfscorevs ）、学習済みのモデル（http://github.com/oddt/rfscorevs_binary ）ともにGitHubで公開されています。

各OS向けの学習済みのbinaryが用意されているので、ready-to-use RF-Score-VSをローカルに落として来ました。
こちらのbinaryにはRD-SCORE-VS v2 が入っている、とのことです。

使い方は簡単で、ターミナルでダウンロードして来たフォルダ（rd-score-vs_v1）に移動します。

以下でヘルプを参照できます。

$ ./rf-score-vs --help

ドッキングに使ったタンパク質のpdbファイル（今回は5nix_chainAB.pdbという名前）と、AutoDock Vinaのドッキングで得たリガンド構造を含むsdf（ODDT_Dock_All_Result_mod.sdfという名前）を同じフォルダに入れてコマンドを実行します。

./rf-score-vs_v1 --receptor 5nix_chainAB.pdb ODDT_Dock_All_Result_mod.sdf -O ODDT_All_ligands_rescored.sdf

引数 --receptor に続いてpdb、sdfを入力し、 -O　で指定したファイル名で、再計算されたスコア（フィールド名: RFScoreVS_v2）が書き込まれたsdfが出力されました。

今回は全てのドッキングポーズ（50化合物ポーズ合計449個）についてスコアを再計算しています。では結果を確認していきましょう。

2-2. スコアを確認

RDKitでSDFを読み込んでおきます。

from rdkit import Chem
from rdkit.Chem import AllChem, PandasTools 
import pandas as pd

#PandasToolsで読み込み、スコアのTypeをfloatに変換する
Rescored_df = PandasTools.LoadSDF('./ODDT_All_ligands_rescored.sdf')
Rescored_df['RFScoreVS_v2'] = Rescored_df['RFScoreVS_v2'].astype(float)
Rescored_df['vina_affinity'] = Rescored_df['vina_affinity'].astype(float)

f:id:magattaca:20190407181604p:plain

AutoDock Vinaのスコアは結合自由エネルギ-の予測値（kcal/mol）で与えられるのに対し、RF-Score-vsのスコア再計算結果は pK_d/i の単位で与えられます。
学習の段階でDUD-Eの化合物の活性を pK_d/i = 6 よりも弱い（小さい）場合にinactiveとしています。また、トレーニングセットのDecoy（構造類似の不活性化合物?）には一律に pK_d/i = 5.95 を設定しています。従って予測値が6よりも大きければ活性あり、となりそうです。

スコアの分布をみてみましょう。ヒストグラムを描きます。

import matplotlib.pyplot as plt
%matplotlib inline
plt.hist(Rescored_df['RFScoreVS_v2'], bins=20, range=(5.8, 6.8))

f:id:magattaca:20190407181652p:plain

要約統計量を確認します。

Rescored_df['RFScoreVS_v2'].describe()
"""
count    449.000000
mean       6.005943
std        0.075275
min        5.950000
25%        5.960443
50%        5.977587
75%        6.023838
max        6.438919
Name: RFScoreVS_v2, dtype: float64
"""

最小値は 5.95で、inactiveのスコア設定にあった値が出ています。第1四分位数 5.96、第3四分位数 6.02と、3/4近くの結合ポーズが pK_d/i=6以下となっています。活性値の判断基準（6以上）からすると、活性なしと判別されたポーズが 3/4 近いということになりそうです。（・・・この解釈で良いのだろうか？？）

AutoDock Vinaのスコアと一緒にプロットしてみます。

x = Rescored_df['RFScoreVS_v2']
y = Rescored_df['vina_affinity']

plt.scatter(x,y)
plt.title("RF-Score-VS vs AutoDock Vina")
plt.xlabel("RF-Score-VS")
plt.ylabel("AutoDock Vina")
plt.grid(True)

f:id:magattaca:20190407181810p:plain

AutoDock Vinaのスコアは小さいほど良く、RF-Score-VSは大きいほど良い結果、となりますので左上から右下に向かう対角線方向に動いて入れば、予測結果に相関があることになりそうです。

図を見るとRF-Score-VS （x = 5.5, inactive）の線上に様々なy値の点が分布ししており、今回の再スコアリングではかなり異なるスコアリング結果が得られていそうです。特に興味深い結果なのが、y軸で最小値に近い点（前回記事でスコアが良いとして選択していた点）が、再スコアではかなりスコアが悪いという結果になることです。こんなに違う結果になるとは。。。

2-3. 化合物ごとのベストポーズを取得

上のプロットは全499個の結合ポーズ全てなので、ここから各化合物ごとに最も良いスコアとなっているポーズのみを取り出したいと思います。元々の化合物のIDを「original_id」というフィールドに格納しているので、これを使ってgroupbyすれば良さそうです。

# 各化合物ごとのポーズ数を確認してみる
Rescored_df.groupby('original_id').size()
"""
original_id
PB27242793     9
PB300184784    9
PB57131569     9
PB90021090     9
・・・以下略
"""

len(Rescored_df.groupby('original_id').size())
#50

original_idでグループ化したものから、RFScoreVS_v2の値が最大となっているものを取り出して新しいデータフレームとします。*9

id_group = Rescored_df.groupby('original_id')
best_scores = Rescored_df.loc[id_group['RFScoreVS_v2'].idxmax(), :]

# 化合物数を行数で確認
best_scores.shape
# (50, 18)

再度、二つのスコアでプロットしてみます。（コード省略）

f:id:magattaca:20190407182119p:plain

だいぶんスッキリした見た目になりました。要約統計量を確認します。

best_scores['RFScoreVS_v2'].describe()
"""
count    50.000000
mean      6.079934
std       0.111637
min       5.967962
25%       6.004503
50%       6.054512
75%       6.096116
max       6.438919
Name: RFScoreVS_v2, dtype: float64
"""

今回は第１四分位数で6.00、第２四分位数で6.05と、各化合物のベストスコアを抜き出したことで安心できる分布になりました。これなら希望が持てるかもしれません・・・

プロットに合わせて構造を確認してみたいと思います。@iwatobipen先生の記事Plot Chemical space with d3js based libraryを参考にさせていただき、d3jsを使います。

import mpld3  
from mpld3 import plugins
mpld3.enable_notebook()

#描画を作成するためのモジュール
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

#svgを格納するための空のリスト
svgs_list = []

for i in range(len(best_scores)):
    view = rdMolDraw2D.MolDraw2DSVG(300,300)
    
    ROMol_col = best_scores.columns.get_loc('ROMol')
    mol = best_scores.iloc[i, ROMol_col]
    
    prepro_mol = rdMolDraw2D.PrepareMolForDrawing(mol)
    
    view.DrawMolecule(prepro_mol)
    view.FinishDrawing()
    
    svg = view.GetDrawingText()
    svgs_list.append(svg.replace('svg:,',''))

x = best_scores['RFScoreVS_v2']
y = best_scores['vina_affinity']

fig, ax =plt.subplots()
ax.set_xlabel("RF-Score-VS")
ax.set_ylabel("AutoDock Vina")
points = ax.scatter(x, y)

tooltip = plugins.PointHTMLTooltip(points, svgs_list)
plugins.connect(fig, tooltip)

f:id:magattaca:20190407182752g:plain

gifではわかりにくいですが、スポットにカーソルを合わせると構造式が表示されます。

2-4. RF-Score VSにおけるスコアTop 10の取得

Top 10を取り出してSDFとして書き出しておきます。

#スコアでソートし、上から10個を取り出す
best_scores_10 = best_scores.sort_values('RFScoreVS_v2', ascending=False).iloc[:10]
PandasTools.WriteSDF(best_scores_10, 'RF-ScoreVS_Top10.sdf', properties=best_scores.columns)

描画します。前回と同じなのでコードは省略しますが、レジェンドのRFがRF-Score VSのスコアを、ADがAutoDock Vinaのスコアを表します。

f:id:magattaca:20190407183106p:plain

AutoDock Vinaのスコアで選択したTop 10の結果と並べてみます。

f:id:magattaca:20190407183147p:plain

両スコア関数で選択される化合物の違いを見てみてると、ぼんやりとですがRF-Score-VSではリガンドの中心付近にヘテロ芳香族環が位置していたり、AutoDock Vinaによる選択と比較して、分子全体が丸まっているように見えるものなどもあります。

3. アッセイ候補提出ファイルの作成

新しいスコア関数 RD-Score VSによる再スコアリングと結果の取り出しが完了しました。AutoDock Vinaのスコアリングとそれぞれ違いがあって面白いので、各々Top 5を取り出して合わせて10化合物とし、アッセイ候補としたいと思います。（ダイバーシティー！）

提出フォーマットはタブ区切り、ということなのでcsvモジュールを使ってファイルを作ります。まずはリストを作成。

top_10_list = []

# RFScoreVSを基準に5個取り出し
RF_sup = Chem.SDMolSupplier('./RF-ScoreVS_Top10.sdf')

for i in range(5):
    temp_list =[]
    mol = RF_sup[i]

    rank_num = i + 1
    compound_id = mol.GetProp('original_id')
    score = mol.GetProp('RFScoreVS_v2')
    
    temp_list = [rank_num, compound_id,score]
    top_10_list.append(temp_list)

# AutoDock Vinaを基準に5個取り出し
ADVina_sup = Chem.SDMolSupplier('./SBVS_Top10.sdf')

for i in range(5):
    temp_list =[]
    mol = ADVina_sup[i]

    rank_num = i + 6
    compound_id = mol.GetProp('original_id')
    score = mol.GetProp('vina_affinity')
    
    temp_list = [rank_num, compound_id,score]
    top_10_list.append(temp_list)

10化合物についてIDとスコアの組み合わせのリストのリストが作成できたので、タブ区切りのファイルを書き出します。

import csv
with open('./top10_magattaca.txt', 'w') as f:
    writer = csv.writer(f, delimiter='\t')
    writer.writerows(top_10_list)

念のため読み込んで確認

with open('./top10_magattaca.txt') as f:
    print(f.read())
"""
1 Z393761442  6.438919
2 Z243276918  6.404471
3 Z200290200  6.386205
4 Z368242106  6.385228
5 Z278994950  6.207442
6 Z237879862  -10.9
7 Z19455181   -10.9
8 Z872606870  -10.8
9 Z2788907736 -10.4
10    PB90021090  -10.4
"""

無事目的のファイルが作成できていそうです。果たしてこの中に活性のあるものはあるでしょうか？？？

同様に、500化合物のリストも同じファイル形式としておきます。スコアとしてはLBVSを行なったときの類似性評価のTanimoto係数によるスコアを記載しておきます。　（色々こねくり回して494個と500個に少し足りませんがご容赦いただけると信じて）

Top500_sup = Chem.SDMolSupplier('./SCR_compounds_SimScore.sdf')
top_500_list = []
rank_num = 0

for mol in Top500_sup:
    rank_num +=1
    compound_id = mol.GetProp('original_id')
    score = mol.GetProp('Similarity_Score')
    temp_list = [rank_num, compound_id,score]
    top_500_list.append(temp_list)

完了！　これを提出フォームに登録すればOK！

4. まとめ

今回はODDTに関して調べているうちに行き当たった論文、Sci.Rep. 2017(7)46710 に魅かれて機械学習を用いたドッキングスコア関数について調べてみました。実際に適用して見るとAutoDock Vinaによるスコアリングと異なる傾向が見え、とても面白い結果だと感じました。文献中ではまだまだ修正・改善が必要だ、と書かれていましたが、既存のスコア関数のモデルの問題点など、非常に学ぶことが多い論文だと思いました。

前回の記事でドッキングを行なった際には、とりあえず動いてデータが出れば良し！ということで進めていきましたが、AutoDock Vinaが少し古いソフトウェアなのかな？これだけで良いのかな？という懸念がありました。ドッキングそのものは時間がかかるので、何度も条件を変えて試すのは大変ですが、結合ポーズが出ていればスコア関数を修正して評価精度をあげることができる、というのを学べて良かったです。（新しいっぽいことをやって格好つけたかった）

とりあえずデフォルト+α のことが少しできたので満足！

5. 総括

提出も完了したので、最後にこれまで行なってきた処理手順をまとめます。

1. 前処理（脱塩、構造標準化）
1. 指標による絞り込み（主にRule of 5の適用）
1. 部分構造による絞り込み（オルト置換ビフェニル構造）
1. ファーマコフォアスクリーニング（ファーマコフォアを変えて２回実施） ➡︎ 化合物 Top 500 (厳密には494個)
1. LBVS（既知活性化合物に対するフィンガープリントベースの類似性評価による順位づけ） ➡︎ 化合物 Top 50を次工程へ
1. SBVS (AutoDock Vina + RF-Score VS) 　　　　 ➡︎ 化合物 Top 10をアッセイ希望に選択

以上の６ステップです。書いてしまうとこれだけ、という感じもありますが、各工程でびっくりするくらいに行き詰まったので予想以上に時間がかかってしまいました。（当初、３月上旬提出だったのに、ずるずる遅れてしまいました。遅れている上に余計な寄り道ばっかりしてしまいましたし・・・。締め切りを４月中に延ばしてくださった創薬ちゃんに感謝です。）

展望も戦略もないままはじめて右往左往していましたが、行き詰まるたびにTwitterで解決方法や、方針、文献諸々を一からご教示いただき、なんともありがたい気持ちでいっぱいです。この場を借りて皆様にお礼申し上げます。今後ともよろしくお願い申し上げます。（図々しくてすみません）

あとは、選んだ化合物の中に活性を持つものがあることを祈って！果報は寝て待ちましょう。
今回の記事も間違いがたくさんあると思います。ご指摘いただければ幸いです。

*1:タンパク質計算科学基礎と創薬への応用　神谷成敏・肥後順一・福西快文・中村春樹　著　共立出版 2009年初版１刷

*2:Uehara, S., Tanaka, S., CICSJ Bulletin 2016(34)10

*3:Machine-learning scoring functions to improve structure-based binding affinity prediction and virtual screening Ain, Q.U., Aleksandrova, A., Roessler, F.D., and Ballester P.J. WIREs Comput. Mol. Sci.2015(5)405

*4:A machine learning approach to predicting protein-ligand binding affinity with applications to molecular docking. Ballester, P.J. and Mitchell., J.B.O., Bioinformatics. 2010(26)1169

*5:Does a more precise chemical description of protein-ligand complexes lead to more accurate prediction of binding affinity? Ballester, P.J., Schreyer, A., Blundell. T.L., J. Chem. Inf. Model. 2014(54)944)、RF-Score v3 ((Li, H., Leung, KS., Wong, M.H., Ballester P.J., Mol. Inform. 2015(34)115

*6:ML to improve scoring functions binding and virtual screening
RF-Score: a Machine Learning Scoring Function for Protein-Ligand Binding Affinities

*7:Mysinger,M.M.,Carchia,M.,Irwin,J.J.,Shoichet,B.K.,J.Med.Chem.2012(55)6582）

*8:「タンパク質計算科学」の記述参考

*9:pandasのgroupbyを使ってグループ内で最大値（最小値）を持つ行を取得する

2019-04-04

SBVSしようとする話

これまで以下の手順で化合物の絞り込みを行ってきました。

記述子といった指標等を用いた絞り込み
活性化合物に頻出の部分構造を用いた絞り込み
共結晶構造に基づくファーマコフォアによるスクリーニング
フィンガープリントを用いた活性化合物群への類似性に基づくスコア化

創薬レイドバトル2018の提出リストは

候補化合物 ID リスト (1): 500 化合物
候補化合物 ID リスト (2): (1) から更に 10 化合物に絞った優先順位付きリスト

です。すでに上記絞り込みにおいて、手順3.で残り約500個、手順4.でスコアによる順位づけも完了したため、リストはできた！！！・・・と言いたいところですが、スコア上位10化合物にいまいち自信が持てません。

残る課題としては、タンパク質側の情報を明示的には取り扱っていない、ということでしょうか？　手順3、においてファーマコフォア作成時に共結晶構造を参考にはしたものの、情報の一部を抽出したもので間接的な取り扱いです。やはり最後は直接対決！！
・・・ということで、ドッキング（SBVS?）に挑戦しようと思います。
絞り込んだスクリーニング化合物群の中に、無事 PD-L1 のお眼鏡に叶う候補が含まれていることを祈ってシミュレーションしていきましょう。

1. AutoDock
2. ドッキング？
3. ドッキング実施
4. SBVSにむけたスクリーニング化合物の準備
5. AutoDock vinaとODDTでSBVS
6. まとめ

1. AutoDock

早速RDKitをimportして、、、、といきたいところですが、RDKitにはドッキングの機能はないようなので別の方法を探さないといけません。今回は日本語で解説してくださっている記事、資料が検索したら出てきたという安直な理由でAutoDock vinaを使いたいと思います。*1

AutoDockはThe Scripps Research InstituteのThe Olson Laboratoryにより開発されているLigand-Protein Dcokingソフトウェアだそうです。説明によると、結合自由エネルギーの予測と遺伝的アルゴリズムに基づく探索により、タンパク質とフレキシブルなリガンド間でドッキングを実施します。今回使用しようとしてるAutoDock vinaはAutoDockを高速化、精度を高めた改良版です*2。また、ドッキングの実施や結果解析を手助けするGUIツールとしてAutoDockToolsというソフトウェア開発も行われています（こちらは他のツールとまとめてMGL Toolsに含まれています）。

これらを使えばいい感じにタンパク質にリガンドが結合するか否かわかるはず！・・・でも、そもそもドッキングって具体的に何を計算しているのか？？

2. ドッキング？

2-1. ドッキングについて

ドッキングに関して、日本化学会情報化学部会誌（CICSJ、現在はケモインフォマティクス部会？）にちょうど良い文献（タンパク質-リガンドドッキングの現状と課題 *3 ）があったので参照します。

同文献によると、ドッキング計算（protein-ligand docking）は

タンパク質の立体構造を顕に扱うSBDD(structure-based drug design)で特に重要な役割を果たす
計算機により標的タンパク質とリガンドの立体的な結合構造と結合親和性を予測する手法

だそうです。これを利用した仮想スクリーニング（virtual screening, VS）がSBVS（の一つ？）ということになります。

原子レベルで相互作用機序を理解できるため、合理的な医薬品設計の手助けとなる一方、タンパク質のような非常に多数の原子を含む構造を扱うのは計算コストが高すぎるという問題があるため、簡易化したモデルが使われます。
「簡易化」 ≒ 「タンパク質、リガンドのもつ立体構造の柔軟性をどの程度まで考慮に入れるか？」で、柔軟さを許容すればするほど計算コストが高くなります。複雑さの段階に応じて以下の４つのモデルが紹介されていました。

Rigid Docking: 鍵と鍵穴モデル(Lock-and-Key Model)に従い、タンパク、リガンド共に剛体として扱う
Flexible Docking: リガンド側のみ柔軟性を考慮（タンパク質は剛体, Flexible Ligand Model）
Flexible Side-Chains Docking: リガンドに加えて結合ポケットのアミノ酸側鎖の自由度も考慮（Paritally Flexible Protein Model）
Flexible Protein Docking: タンパク、リガンド共に柔軟性を考慮（Induced-Fit Model）

このうち、標準的なものは「3. フレキシブルドッキング」で、AutoDockもここに含まれています。

ドッキングの目的である「結合構造、結合親和性の予測」にあたって、重要となるのが効率的な探査アルゴリズムと高精度なスコア関数の２つです。前者は安定な結合構造を予想、探索するもので、これにより得られた結合構造がスコア関数により順位づけされます。文献中ではどちらも様々な手法の例が挙げられていましたが、AutoDockの場合、探査アルゴリズムとして遺伝的アルゴリズム、スコア関数として力場に基づいたスコア関数が使われています。

2-2. 遺伝的アルゴリズム？

AutoDockの遺伝的アルゴリズムについて情報処理学会(IPSJ) 第75回全国大会論文集（2013年 p.901-902 電子図書館リンク）に記載があったので引用します。*4

AutoDockでは、最適解を探索するアルゴリズムとして遺伝的アルゴリズム (GA: Genetic Algorithm) を用いている。遺伝的アルゴリズムとはデータを遺伝子として表現し、複数の個体を適応度によって選択し、交叉、突然変異などの遺伝的操作を行いながら最適解を探索するアルゴリズムである。

また、愛媛大学　村上研究室のページ（遺伝的アルゴリズム）によると、特徴として

解空間構造が不明
決定的な優れた解法が発見されておらず、
全探索が不可能と考えられるほど広大な解空間を持つ問題に有効

という点があげられるそうです。確かにタンパク質と低分子の結合構造という、無数の組み合わせがありそうな問題を解く上で有効そうです。

先のIPSJ論文集に記載されているGAの流れを少し書き換えると以下のような雰囲気です。

ランダムで初期個体として複数の結合構造を発生（最初の世代）（原子座標：遺伝子、リガンド:個体、染色体 ）
スコア関数に基づき各結合構造のエネルギースコアを計算（適応度の評価）
スコアの良い構造を適応しているとして選択(selection)
選択した構造の間で部分構造を入れ替える（交差(crossover)）
さらに一定の確率で部分構造を変化させる（突然変異(mutation)）（局所最適(local minimum)に収束するのを防いで、大局的最適(global minimum)に近づけるため？）
5.までで残った構造を新しい世代として、再度2.→5.のステップを行う
最大世代数 or 最大評価回数に達したところで反復計算を終了
最終世代の結合構造の中で、最安定エネルギーのものを最適解として出力

生命科学の用語をアルゴリズムに用いるというのは面白いですね。これがミームというやつか・・・（ちがうか）

2-3. スコア関数

探索方法の流れが把握できたので、エネルギーを計算するスコア関数について眺めてみます。先ほど省略しまししたが、CICSJの文献ではスコア関数として３種類が挙げられていました。（AutoDockのスコア関数は1に分類）

力場に基づいたスコア関数（force-field-based scoring function）
経験的なスコア関数 (empirical scoring function)
統計に基づいたスコア関数(knowledge-based scoring function)

1.はファンデルワールス(vdW)相互作用や静電相互作用といった物理化学的な相互作用の計算を積み上げて構成されており、比較的構成が簡単になる一方で、あまりにもモデル化されすぎているため精度があまり高くないという問題があるそうです。そこで、1.の関数に類似したものをもちいつつ、現実の実験データ（既知のタンパク質-リガンド結合自由エネルギー）を取り込んでパラメータをフィッティングさせた経験的なスコア関数 2.があり、こちらの方がΔGの再現性が上がる傾向にあるそうです。1.、2.に対して、3.はさらに実験データに重みを置いており、既知の結晶構造を統計的に処理することで得た確率分布からスコア関数が作成されているそうです。

雰囲気的には「1. → 2. → 3.」の順に「理論 → 実験」の比重が変化している感じです。*5

では、AutoDockのスコア関数ではどのような物理化学的相互作用が考慮されているのか？
以下のような式となっているそうです。

ΔG = ΔE_vdw + ΔE_{electrostatic} + ΔE_H-bond + ΔE_desolvation + ΔS_conf

各項は以下に対応しているとのことです。

ΔE_vdw : ファンデルワールス相互作用（Lennard-Jones (12,6) ポテンシャル）
ΔE_{electrostatic}: 静電相互作用（距離依存性誘電率 *6 ）
ΔE_H-bond : 水素結合による相互作用 (Lennard-Jones (12,10) ポテンシャル/ Goodford directionality)
ΔE_desolvation : 脱溶媒効果
ΔS_conf : リガンドのエントロピー変化の項

上記のうち、最初の３つの項は物理化学的項目としてよく見かけるものですが、最後の２つドッキングならではという気もします。
リガンドがタンパク質に結合するにあたって、溶媒（水）に囲まれた環境から、溶媒が外れ（脱溶媒）タンパク質の結合サイト（より疎水性の高い環境）に移行するという過程をたどると考えます。この時のエネルギー変化を考慮するためのものがΔE_desolvationです。また、溶液中では自由な立体構造（配座、コンフォメーション）をとっているリガンドですが、結合に際してがポケットに合うように固定化されこの自由度が失われます。この分のエントロピー変化を考慮するものがΔS_conf となっているようです。

このあたり字面を追ってもよくわからなかったのですが、ワシントン大学 Jay Ponder Labの講義資料Tutorial & Introduction to AutoDockでなんとなくイメージができました。*7

AutoDock vinaはAutoDockの改良版、とのことなのでスコア関数も変わっているかもしれませんが基本的には同じだと思います。（調べてなくてすみません・・・）

ぼんやりとドッキングの仕組みがわかってきたので実際の作業に移りましょう。

3. ドッキング実施

3-1. AutoDock vina の処理流れ

計算化学.com ドッキングシミュレーションのやり方(AutoDock vina）　にはAutoDock Vinaに必要なファイルとして「タンパクの可動部位」「タンパクの非可動部位」「リガンド」の３つの「PDBQTファイル」と、「inputファイル」の４つが必要と記載されていました。

PDBQTファイルはpdbファイルと異なり、原子の座標に加えて部分電荷(parital charge)とAutoDockのアトムタイプを含むファイル形式となっているそうです*8。より具体的にはPDBファイルに水素原子の負荷、部分電荷（ex. Gastiger部分電荷）の情報を付加してやれば良い、ということのようです。

上記に従ってタンパク質とリガンドのファイルを用意し、以下の流れでドッキングを行います。

step 1. Grid Box（ドッキングシミュレーションを行う探索範囲）の作成タンパク質の結合ポケットを範囲として指定します
（Boxの中心座標 + Boxのサイズ情報）
step 2.ドッキングを実施
step 3.結合親和性（エネルギー）、結合構造を得る

では、具体的にやっていきましょう.

3-2. 使用するタンパク質の選択

実際にどのデータを使ってドッキングを行うか？
今回は共結晶構造があるのでこちらのタンパク質情報を使います。複数報告されていますが、このうち「PDBid: 5NIX」を用いたいと思います。理由は以前の記事で比較した際に、この複合体はリガンドが大きくため、タンパク質側の結合ポケットが広がっている可能性が高い、と考えられるからです。

なぜポケットが大きい方が良いか？「タンパク質計算科学基礎と創薬への応用」の記述を引用します。(p. 217)*9

・・・注意すべきは、タンパク質-化合物ドッキングの精度は通常２Å以上の粗い精度であって、タンパク質と化合物の原子はたいていの場合、vdW衝突しているということである。・・・（中略）・・・したがって、vdWポテンシャルの衝突を減らすことが、ドッキングソフトの重要な技術となっている。

なるほど！それならば最初から空間が大きいもの選んどいた方が良いよね、、、っていう浅い考えです。（この本、積読してないで早く読めばよかった・・・）

3-3. Pymolでタンパク質を前処理

ドッキングに向けて「PDBid: 5NIX」の構造をPymolでちょっと処理しておきます。

まず、この構造は４量体（Chain A、B、C、D）となっているので、簡便のためChain C、Dを削除しておきます。

PyMOL> remove (chain C,D)
 Remove: eliminated 2093 atoms in model "5nix".

ついでに水も除きます。

PyMOL>remove resn hoh
 Remove: eliminated 133 atoms in model "5nix".

remove solventでも良いらしいです。参考

Sequenceを表示させて確認すると、EDO（エチレングリコール）というものもあったので、不要そうなので除去しておきます。

PyMOL>remove resn EDO
 Remove: eliminated 4 atoms in model "5nix".

以上の操作でだいぶスッキリした見た目になりました。

f:id:magattaca:20190403224018p:plain

3-4. Grid Box設定のための座標確認

ドッキングの中心（作成するGridの中心）を決めるために、座標を取得します。今回はファーマコフォアスクリーニングを実施済みの化合物群をドッキングに使用予定のため、ファーマコフォア作成時に利用したリガンド原子部位の中心をGrid中心とします。

PyMOL>pseudoatom center, sele
 Warning: 'center' is a reserved keyword, appending underscore
 ObjMol: created center_/PSDO/P/PSD`1 /PS1

f:id:magattaca:20190402232114p:plain

安易にcenterという名前をつけたら怒られました・・・。このpeudoatomの座標は以下です。

PyMOL>xyz = cmd.get_coords('center_', 1)
PyMOL>print xyz
[[ -5.486217  11.170054 -24.066675]]

全体を5nix_model.pdb、タンパク質のみ（Chain A,B）を5nix_chainAB.pdb、リガンドのみをmol2形式5nix_lig.mol2としてそれぞれ出力しておきました。

3-5. テストドッキング

PDBの共結晶構造から、余分なものを除いたタンパク質のみのPDBファイルと、リガンドのみのmol2ファイルが用意できたのでこれらを用いてドッキングの流れを確認してみます。

参考にしていたページでは準備にAutoDock Tools（MGL Tools）を使うとのことでしたが、私のMac上ではエラーが出てうまく動かなかったので、代わりにUCSF ChimeraのpluginとしてAutoDock vinaを使用します。（Chimeraはアカデミックフリーですが商用利用はライセンスが必要だそうです）*10

3-5-1. タンパク質/リガンドの前処理

AutoDock vinaに合わせて、先に用意したタンパク質PDBファイルを処理します。
ChimeraでPDBファイルを読み込んだのち、Tools の Surface/Binding Analysis から DockPrep　を選択すると新しい画面が開きます。

f:id:magattaca:20190402232521p:plain

この画面の指示に従って、

溶媒の除去 (Delete Solvent)
水素原子の付加 (Add hydrogens)
電荷の付加 (Add charges)

といった必要な処理を選択し、Mol2ファイルで書き出します。

f:id:magattaca:20190402232635p:plain

リガンドのファイルも同様に水素原子と電荷の付加の処理をしました。カルボキシル基には水素原子が付加されておらず、アミノ基には付加されているので、それぞれ負電荷、正電荷を帯びるように処理されているようです。（下図、黄色丸部位）

f:id:magattaca:20190402234105p:plain

3-5-2. ドッキングのジョブを投げる準備

一度「close session」した後、処理したタンパク質とリガンドのmol2ファイルを読み込みます。
Tools の Surface/Binding Analysis から AutoDock Vinaを選択し、開くサブ画面に設定を書き込んでいきます。

f:id:magattaca:20190402234444p:plain — 設定の画面例

上の図で行なっている設定は、

output fileの名前（今回は5nix_redockとした）をつける
タンパク質（receptor）とリガンドを指定
(mol2ファイルを読み込んでいるのに名前は処理前のもので表示された・・・？）
Gridの指定
AutoDock vinaを行う場所の設定（web？ or local？）

Gridの設定は、「Receptor search volume options」の部分で行なっており、Grid中心（Center）を先にpymolで確認した座標を参考に設定、Grid Boxの大きさを size x:20 y:20 z:20としてます。リガンドがBoxからはみ出していたのでGrid中心を少しX軸方向にcenterの座標の位置をずらし、微調整をしてあります（x:-7, y:11, z:-24）。

AutoDock vinaを実行する場所はpublic web serviceか、あるいはlocalを選択できます。今回はlocalで実施するためvinaを入れたpathを指定します。

以上で準備完了、あとは実行のみ

3-5-3. ドッキングの実行と結果

実行・・・Macが発熱します。

結果・・・４つの構造が残ったようです。

mode	affinity (kcal/mol)	RMSD l.b.	RMSD u.b.
1	-11.5	0.000	0.000
2	-11.5	1.155	1.585
3	-11.1	1.280	1.927
4	-10.8	3.173	6.337

pdbqtファイルがいくつか出力されているのでpymolで構造を眺めてみます。 pymolで出力5nix_redock.pdbqtを読み込むと4つの構造をstateとして含む状態で読み込まれます。少しみづらいので「Action → state → split」とすることで一つずつに分けます。

共結晶構造中のリガンド（赤色）とともにそれぞれ描画してみます。

f:id:magattaca:20190402235052p:plain

エネルギー的に近いmode 1、2、3に対してmode　4ではかなり結合ポーズが変化しています。ちょうどファーマコフォアには含めなかったベンジル側鎖がフリップしたような構造となっています。より安定なmode 1、2 はほぼ元々のリガンドに重なっており、今回の結果ではかなり再現性の良い結果となった、と言えそうです。（初期構造もグリッドも共結晶構造をもとに設定しているので当然？かもしれませんが。。。）

スコアに関してですが、香川大学のページの記載（Autodock4の場合）によると、

クラスターの平均エネルギーの差が約2.5 kcal/mol以下の場合，AutoDock force fieldの標準偏差以内であり，エネルギーからどちらが「正しい」かをいうことは出来ない

とのことだそうです。今回この標準偏差よりも安定化しているので、良い結合ポーズ・・・ということで良いのでしょうか？再度「タンパク質計算科学基礎と創薬への応用」の記述を引用します。(p. 225)

・・・通常、ΔGの再現は悪く、2.5~3kcal/mol程度の精度が出ればいいほうであるが、ヒット化合物のΔGは4~8kcal/mol程度、最適化された医薬品で10kcal/mol程度なので、ΔGの計算のみでヒット・リード化合物の活性地を評価するのは困難である。・・・

初版が10年ほど前の書籍なので、現在のソフトウェアの精度はこの時よりも上がっていると思われますが、AutoDock vinaの報告が J.Comp.Chem. 2010(31)455、なので今回の目安としては有効なのではないかと思います。

4. SBVSにむけたスクリーニング化合物の準備

テスト構造での検証、手順の流れがわかってきましたのでようやく本題のスクリーニングに移りたいと思います。スクリーニング対象化合物としては、とりあえずLBVSのスコアTop 50としたいと思います。（私のPCの性能的にこれくらいが限界、、、）

4-1. ファーマコフォアを利用したコンフォマーの発生

まずはリガンドの立体構造が必要なので初期構造を準備します。折角なのでRDKitを利用して、これまでのスクリーニングの結果を応用したいと思います。

手元にあるスクリーニング対象化合物は、活性化合物群に対する類似性をベースにスコアリングを実施したものです。したがって、共結晶構造中のリガンドとある程度の類似性があると思われます。そこで、安定かつ共結晶構造中のリガンドにできるだけ近い配座を用意したいと思います。

具体的には

step 1. ファーマコフォアをembed
step 2. ETKDGでembedした構造を最適化
step 3. 一つの化合物に複数embedできた場合、Feature-mapで最も共結晶のリガンドに近い構造を代表として選出する。

という流れです。

まずはLBVSのスコアを付与したSDFを読み込み、スコアでソートしてTop50化合物を取り出しました。そしてファーマコフォアの設定には、以前ファーマコフォアサーチで用いたものと同じ値を使用し、設定を行いました。この辺りの作業は繰り返しなので省略します。

前回との違いは、このファーマコフォアサーチの関数においてさらに立体構造の最適化のステップを加えることです。以前ファーマコフォアサーチに使った関数はembedできた数のみを返すようにしていましたが、これをconfomerのリストを返すように書き換えます。また、この関数では不斉の情報を取り除いていますが、幸い（？）今回の化合物リストに不斉の情報はなさそうなのでそのまま進めます。

def PConfs(mol):
    #立体の情報を除去
    smi = Chem.MolToSmiles(mol, isomericSmiles=False)
    mol_re = Chem.MolFromSmiles(smi)
    AllChem.Compute2DCoords(mol_re)
    
    #フィーチャーをそもそも持っているか？
    match, mList = EmbedLib.MatchPharmacophoreToMol(mol_re, Ffactory, pcophore)
    
    if match == True:
        #距離の検証
        #距離行列の取得
        bounds = rdDistGeom.GetMoleculeBoundsMatrix(mol_re)

        #ファーマコフォアのマッチング 
        pList = EmbedLib.GetAllPharmacophoreMatches(mList,bounds,pcophore) 
        
        #pListのそれぞれについてFeatureとマッチする原子のidを取得する
        num_match = len(pList)
        phMatches = []
        for i in range(num_match): 
            num_feature = len(pList[i])
            
            phMatch = []
            
            for j in range(num_feature):
                phMatch.append(pList[i][j].GetAtomIds())
                
            phMatches.append(phMatch)
        
        #原子のidとファーマコフォアをもとに3次元構造を埋め込む    
        confs = []

        for phMatch in phMatches:
            bm,embeds,nFail =EmbedLib.EmbedPharmacophore(mol_re, phMatch, pcophore,count=5, silent=1)
            for embed in embeds:
                #水素原子を付加
                embedH = AllChem.AddHs(embed)
                #構造最適化
                AllChem.EmbedMolecule(embedH, AllChem.ETKDG())
                
                #元々の分子のプロパティを追加する作業
                PropNames = list(mol.GetPropNames())
                prop_dict = {}
                #プロパティの値を取り出し
                for n in PropNames:
                    vl = mol.GetProp(n)
                    prop_dict[n] = vl
                #confomerに付与
                for k, v in prop_dict.items():
                    embedH.SetProp(k, v)
                
                #リストに追加
                confs.append(embedH)         
        return confs            
    else:
        return 0

スクリーニングに用いる50個の化合物の一番最初の分子を取り出し、関数の挙動を確認します。

test_m = top50_mols[0]
test_confs = PConfs(test_m)
len(test_confs)
# 4

コンフォマーが４つ生成されました。py3Dmolを使った３Dの描画結果を以下に貼っておきます。

f:id:magattaca:20190403000118p:plain

4-2. Feature Mapによるコンフォマーの選別

次に、得られたconfomerの中から共結晶中のリガンドコンフォマーに最も類似したものをFeature Mapを使って選択します。RDKitのブログ記事Using Feature Mapsを参考にしました。元の文献を読めていないので理解できていないのですが、複数のコンフォマーの中から良いものを選ぶことができそうな感じです。また、比較対象のリガンドとして、ファーマコフォアを満たすのに必要最小限の構造をもつように見える、「PDBid:5N2F」のリガンド（8HW）を選択しました。

RDKitのブログ記事をそのまま関数にしていきます。

from rdkit.Chem.FeatMaps import FeatMaps
from rdkit.Chem import rdMolAlign

def FeatMapScorer(mol, tmp):
    #比較対象(tmp)との間でアラインメントをとる
    o3d = rdMolAlign.GetO3A(mol,tmp)
    o3d.Align()
    
    #feature mapのパラメーターの辞書を作成（デフォルトの値を使用）
    fmParams = {}
    for k in Ffactory.GetFeatureFamilies():
        fparams = FeatMaps.FeatMapParams()
        fmParams[k] = fparams
    
    #featureのうちより一般的なものだけを残すための設定
    keep = ('Donor','Acceptor','NegIonizable','PosIonizable','Aromatic')
    featLists = []
    for m in [mol, tmp]:
        #featureをまず全部取得
        rawFeats = fdef.GetFeaturesForMol(m)
        #keepにある欲しいものだけリストに追加 
        featLists.append([f for f in rawFeats if f.GetFamily() in keep])
    
    #FeatMap オブジェクトを作成
    fms = [FeatMaps.FeatMap(feats=x, weights=[1]*len(x), params=fmParams) for x in featLists]
    
    #feature mapに対するスコアを取得
    score = fms[0].ScoreFeats(featLists[1]) / min(fms[0].GetNumFeatures(), len(featLists[1]))
    
    #scoreを返す
    return score

先のテスト化合物で得た４つのコンフォマーについて、Feat Mapのスコアを求めてみます。

test_scores = []

for conf in test_confs:
    score = FeatMapScorer(conf, temp_8HW)
    test_scores.append(score)

print(test_scores)
#[0.25512226843031627, 0.17970526398716216, 0.05047724449395227, 0.025276456816820638]

スコアが最大となるコンフォマーのみを取り出したいので、numpyのargmaxを利用します。

import numpy as np

test_score_array = np.array(test_scores)
test_best_conf_index = np.argmax(test_score_array)
test_best_conf = test_confs[test_best_conf_index]

上記test_best_confを描画し、無事立体構造が選別できていることが確認できたので、本番、50個の化合物群全てに適用します。

Top50_3D_list = []

for mol in top50_mols:
    confs = PConfs(mol)
    
    FMapScores = []
    
    for conf in confs:
        score = FeatMapScorer(conf, temp_8HW)
        FMapScores.append(score)
    
    ScoreArray = np.array(FMapScores)
    BestConfIndex = np.argmax(ScoreArray)
    BestConf = confs[BestConfIndex]
    
    Top50_3D_list.append(BestConf)

Top50_3D.sdfとして書き出しておきました。

4-3. PyMolでコンフォマーを眺める

本当にテンプレート構造に対して重なりの良いコンフォマーが選ばれているのか、PyMolで確認してみます。テンプレートにした元のリガンド（8HW）と、生成した50化合物のコンフォマーを含むSDFをPyMolで読み込みます。複数の化合物を含むSDFは、1つの化合物を1つのstateとするファイルとして読み込まれているようです。右下再生ボタンを押すと、順番に各stateを表示するムービとして眺めることができました。

f:id:magattaca:20190403230400g:plain

いくつかははみ出しているように見えますが、全体的に良い重なりとなっているように見えます。

試しにタンパク質も同時に表示させて、違う角度から眺めてみます。

f:id:magattaca:20190403230828g:plain

アミノ酸残基の側鎖を表示させていない状態でも、化合物がタンパク質にぶつかっている（重なっている）様子が見て取れます。やはりリガンドのコンフォマーに似せるだけでは、最良なポーズ、とは行かないようです。AutoDock Vinaによるドッキングは、この化合物たちの中から上手く結合ポケットに収まるようなポーズを見つけることができるでしょうか？

5. AutoDock vinaとODDTでSBVS

スクリーニング対象化合物群の３次元構造の生成が完了したので、これを用いてPD-L1に対するSBVSを実施したいと思います。
・・・しかし一つ一つChimeraで実行するのは大変。。。コマンドもイマイチ理解できないし・・・
と思っていたところ、またしても便利なソフトウェアがありました。その名もOpen Drug Dicovery Toolkit!!（略称 ODDT）。さまざまなオープンソースのケモインフォマティクスのツールがあって便利になったものの、それぞれのinput、outputが異なっていて面倒、、、という問題を解決するため、一つのパッケージに統合してもっと便利にしていきましょう！という趣旨のようです。

元の文献はこちら *11。ドキュメントも充実しています（リンク）。@iwatobipen先生のブログ記事でも使用方法が紹介されていました（Open drug discovery toolkit for python）。

早速使っていきましょう。

5-1. ODDTでchimeraでの操作を追試

いきなりスクリーニング本番を行なって失敗するとショックなので、まずは先にUCSF chimeraで行なったテストドッキングと同じことができるか、練習してみましょう。

まずはODDTをimport...

import oddt
from oddt import toolkit

ドッキングに使用する相手のタンパク質のpdbファイル（先にchimeraで読み込んだものと同じもの）を読み込みます。

protein = next(oddt.toolkit.readfile('pdb', './5nix_docking/5nix_chainAB.pdb'))
protein.protein = True

うまく読み込め、タンパク質であると認識してくれているようです。使い方はヘルプでも参照可能です。

help(oddt.toolkit.readfile)
"""
 Help on function readfile in module oddt.toolkits.ob:
readfile(format, filename, opt=None, lazy=False)
"""

次にリガンド（テストドッキングでは共結晶構造のリガンド mol2ファイル）を読み込みます。

ligand = next(oddt.toolkit.readfile('mol2', './5nix_docking/5nix_lig.mol2'))

Autodock vinaを使うためモジュールを読み込みます。

from oddt.docking import AutodockVina

AutodockVina.autodock_vina()で、ドッキングパラメータを指定します。引数はそれぞれ、

protein : ドッキングに使うタンパク質
size : Grid Boxのサイズ、前回同様に (x, y, z) = (20, 20, 20)
center : Grid Boxの中心、前回と同じ座標 (x, y, z) = (-7, 11, -24)
exhaustiveness、num_modes、energy_range: すべてchimeraのAutodockVina pluginと同様
executable: AutodockVinaのlocalのpathを指定　（相対pathの深さから、私のファイル管理がグチャグチャなことがバレますね・・・）

という感じです。

test_docking = AutodockVina.autodock_vina(protein=protein, size=(20, 20,20), center=(-7, 11, -24),
                                         exhaustiveness=8, num_modes=9, energy_range=3, n_cpu=3,
                                         executable='../../../../autodock_vina_1_1_2_mac/bin/vina')

パラメータを設定してインスタンス（？）ができたら、dock()でドッキングを実行です。引数にligandを指定してやります。

dockin_reslut = test_docking.dock(ligand)

4,5分で計算終了しました（chimeraの時と同じくらい）。水素の付加や部分電荷の付加といった作業を行なっていませんが、勝手にそのあたりの前処理もしてくれるみたいです。

len(dockin_reslut)
# 4

生成された結合ポーズは４つ。リガンドの構造のリストとなっており、それぞれの構造に結合エネルギーなどの情報が付与されているようです。ODDTのドキュメントに従って情報を抜き出してみます。

type(dockin_reslut)
# list
dockin_reslut[0].data.keys()
"""
['vina_affinity',
 'vina_rmsd_lb',
 'vina_rmsd_ub',
 'vina_rmsd_input',
 'vina_rmsd_input_min']
"""

親和性やrmsdといった情報があるようです。一つエネルギーを確認してみます。

dockin_reslut[0].data['vina_affinity']
# '-11.4'

Chimeraで計算した時の最小値（最安定）と同じ値となっていそうです。
4つの結果をPandasToolsのデータフレームで読み込むと以下のようになります。

f:id:magattaca:20190403001907p:plain

残りの３つのエネルギーはChimeraの時と異なるようですが、おそらく初期構造をランダムで発生させているため、初期構造に依存して異なる結果が出ているのだと思います。

ODDTは使用方法がpybel(OpenBabel?)に近い書き方となっているそうです（使ったことがないのでわかりません・・・）。ドッキングの結果で得られるリガンドのタイプを確認してみます。

type(dockin_reslut[0])
#oddt.toolkits.ob.Molecule

このobはOpenBanelの略のようです。

SDFに書き出してみます。RDKitと少し異なり、Outputfile()を用意し、構造を一つずつforループで回してwrite()していきます。

TestSDFile = oddt.toolkit.Outputfile('sdf', 'ODDT_docktest.sdf')
for mol in dockin_reslut:
    TestSDFile.write(mol)
TestSDFile.close()

一連の流れが確認でき、妥当な結果（近い結果）が得られていそうなので、実際のスクリーニングに進みます。

5-2. ODDTでスクリーニング本番

以下を実施、、、

protein_5nix = next(oddt.toolkit.readfile('pdb', './5nix_docking/5nix_chainAB.pdb'))

DockingParams = AutodockVina.autodock_vina(protein=protein_5nix, 
                                           size=(20, 20,20), center=(-7, 11, -24),
                                           exhaustiveness=8, num_modes=9, energy_range=3,
                                           n_cpu=4,
                                           executable='../../../../autodock_vina_1_1_2_mac/bin/vina')

Dock_resluts = []

for mol in oddt.toolkit.readfile('sdf', './Top50_3D.sdf'):
    dock_res = DockingParams.dock(mol)
    Dock_resluts.append(dock_res)

3時間ほどMacが騒音を出した後、停止しました・・・

len(Dock_resluts)
# 50

50個全て終わってました。

5-3. 全ドッキングポーズを出力・確認

全50個の化合物について、各化合物で生成されたドッキングポーズを全て格納したリストを作成します。

All_Result = []

for poses in Dock_resluts:
    for pose in poses:
        All_Result.append(pose)

len(All_Result)
# 449

449個のポーズが得られました。プロパティが保持されているかどうか確認してみます。

# keyを確認
Dock_resluts[0][0].data.keys()
"""
['ID',  'MW',  'MolLogP',  'NumHAcceptors',  'NumHDonors',  'NumRotatableBonds',  'PScore',  'PScore2', 'Similarity_Score', 'TPSA', 'original_id',  'vina_affinity', 'vina_rmsd_lb', 'vina_rmsd_ub', 'vina_rmsd_input', 'vina_rmsd_input_min']
"""
#valueを確認
Dock_resluts[0][0].data.values()
"""
['', '457.48600000000027', '3.5305000000000017', '5', '2', '5', '10', '4', '20.832658', '96.970000000000013', 'Z18052460', '-8.0', '0.000', '0.000', '9.41251', '9.408608']
"""

データはキチンと取れていそうなので、間違えて消してしまわないうちに書き出しておきます。一度sdfで書き出した後、RDKitで読み込めるかどうか確認してみます。

#ODDTでSDF書き出し
DockAllResultOutput = oddt.toolkit.Outputfile('sdf', 'ODDT_Dock_All_Result.sdf')
for mol in All_Result:
    DockAllResultOutput.write(mol)
DockAllResultOutput.close()

#RDKitのPandasToolsで読み込み
All_Dock_df = PandasTools.LoadSDF('./ODDT_Dock_All_Result.sdf')
"""
RDKit ERROR: [22:07:38] Explicit valence for atom # 2 C, 5, is greater than permitted
RDKit ERROR: [22:07:38] ERROR: Could not sanitize molecule ending on line 42625
RDKit ERROR: [22:07:38] ERROR: Explicit valence for atom # 2 C, 5, is greater than permitted
・・・以下エラーが続く
"""

エラーがたくさんでました。試しにMarvin viewでひらいてみると、開くことができました。エラーが出た構造を確認してみるとヘテロ環（芳香族）の扱いにおいて結合次数の割り当てに失敗し、炭素に5つ結合が伸びている、といったような不具合が生じているようです。

f:id:magattaca:20190403003131p:plain

テキストエディタでぽちぽちsdfの結合部分を修正し、再度読み込みます。

All_Dock_df = PandasTools.LoadSDF('./ODDT_Dock_All_Result_mod.sdf')
All_Dock_df.head()

f:id:magattaca:20190403003257p:plain

今度はエラーを出さずに読み込めました。

SDFの書き出し、構造の確認が完了したので結果のドッキングポーズをpymolで眺めてみます。先と同様にタンパク質ともともとの共結晶構造中のリガンド、そしてドッキング結果のSDFを読み込み、ムービーを再生します。

f:id:magattaca:20190403233016g:plain

おお！！バーチャルスクリーニング感あふれてる！！

全結合ポーズの中では、結合ポケットの端（タンパク質と溶媒の境界近く）に位置するものも多く見られ、結合ポケット奥深く（共結晶構造リガンドと重なるような位置）まで入り込んでいるものは少ないように見えます。先にみた「タンパク質計算科学」のコメント通り、結合ポケット奥深くに入るポーズは、アミノ酸残基側鎖とのぶつかり合い(vdW相互作用)で弾かれてしまい、なかなか安定なポーズを見つけづらいのでしょうか？　あるいはただ単に探索空間であるGrid Boxの、座標位置、サイズの設定が良くなかったのかもしれません。

5-4. 各化合物のベスト構造のみを出力・確認

つぎに、各化合物について、複数の結合ポーズの中で最安定のもののみを一つずつ取り出してみます。結合親和性(affinity)の値で取り出すと、値が重複した場合に取り出す構造にブレが生じるので、rmsd値が最小のものを取り出していきます。（AutoDockVinaの結果では最安定の結合ポーズと比較したrmsd値が計算されて出力されています。）

Best_Mode_Result = []

for poses in Dock_resluts:
    rmsd_tmp = []
    for pose in poses:
        rmsd = pose.data['vina_rmsd_ub']
        rmsd_tmp.append(rmsd)
    rmsd_arr = np.array(rmsd_tmp)
    best_idx = np.argmin(rmsd_arr)
    best_pose = poses[best_idx]
    Best_Mode_Result.append(best_pose)

こちらについても一度sdfで書き出した後、エラーの構造を修正して、RDKitで読み込みます。

#ODDTでSDF出力
DockBestModeResultOutput = oddt.toolkit.Outputfile('sdf', 'ODDT_Dock_BestMode_Result.sdf')
for mol in Best_Mode_Result:
    DockBestModeResultOutput.write(mol)
DockBestModeResultOutput.close()

# 構造を修正ののち、名前を変えてSDFを保存。RDKitで再度読み込み
BM_df = PandasTools.LoadSDF('./ODDT_Dock_BestMode_Result_mod.sdf')
BM_df['vina_affinity'] = BM_df['vina_affinity'].astype(float)

スコア順にソートして、ドッキングスコアの良い10個を取り出します。

#スコア順にソート
BM_df_st = BM_df.sort_values('vina_affinity', ascending=True)
#Top10を抜き出す
Dock_Top_10 = BM_df_st.iloc[:10]

これらドッキングスコア　ベスト10の化合物群についてもPyMolのムービーで眺めて見ましょう。

f:id:magattaca:20190403234232g:plain

おお！なんだか思っていたよりも良い感じ！！

スコアの良い結合ポーズでは、結合ポケットの奥深くまで化合物が入り込んでおり、反発を抑えつつ上手くタンパク質と相互作用ができていそうな雰囲気がただよっています。

意外な点としては、複数の化合物でビフェニル部分構造が、元の共結晶構造の該当の部分構造とは重ならない形のポーズが、最もスコアが良い結合構造として選ばれていることです。これらでは疎水性のビフェニルが親水性の溶媒接触面側に位置してしまっているようにも見え、ビフェニル部分構造を基準に選んだことが結合親和性に対しては逆効果に働いてしまっているようにも思えます。部分構造の結合位置を指定してドッキングする、といった細かい設定を行っていけば、より元のリガンドに近いような結合ポーズを出せるのかもしれません。

では、これらの化合物がどのような構造式のものなのか？　LBVSの時の類似性スコアとともに10個の化合物を描画してみます。

SBVS_scores = list(BM_df_st['vina_affinity'][:10])
LBVS_scores = list(BM_df_st['Similarity_Score'][:10])

score_legends = []
for s, l in zip(SBVS_scores, LBVS_scores):
    score = 'SB: ' + str(s) + '/ LB: ' + str(round(float(l), 1))
    score_legends.append(score)

Draw.MolsToGridImage(BM_df_st['ROMol'][:10], legends=score_legends)

f:id:magattaca:20190403004218p:plain

リガンドベースで類似性のスコアのみで評価した時とは異なる分子が取り出されているようです。両者を並べてみます。

f:id:magattaca:20190403005401p:plain

丸で囲った部分が今回のドッキング（SBVS）と前回のLBVSで共通して選ばれている構造でTop 10のうち、３つとなりました。類似性のみのLBVSとはまた異なった特徴を持つ化合物を選んでくることができているようです。

6. まとめ

今回は、最後の難関ドッキングによるSBVSを行って見ました。色々とソフトウェアを行ったり来たりしてかなり見通しの悪い感じの記事になってしまいましたが、なんとか最終的にそれっぽい結果を出すことができました。

LBVSで絞り込んだ50個から最終的に10個を選び、LBVSのスコアとは異なるTop 10化合物となりました。SBVSの特徴として、既知の活性化合物とは構造的に異なるものを見つけ出す可能性がある、との記載がありましたので、その意味では異なる結果が出たことでドッキングをした甲斐があった！・・・と思いたい。

（異なる構造が選びたければ先に部分構造や類似性で絞り込むのは順番的に良くないのでは？といった、つらいツッコミはご容赦いただければ幸いです。）

得られた結合ポーズは、当初想定していたのとは異なる結合様式もあった、など、まだまだ調整・工夫の余地がたくさんありそうです。ほとんどデフォルトのコピー&ペーストですし。。。とはいえ、あまりに”ヒトの目”でみて出したいポーズを狙いすぎても折角ドッキングソフトを使っている意味が薄まってしまいますので、今回はこの結果で良し！！　

ということで一通りそれっぽいことができたので満足。

付け焼き刃で調べながら行なったので、間違い等多数あると思います。ご容赦、ご指摘いただければ幸いです。

P.S. 進捗悪すぎてすでに提出期限オーバー感があるのですが・・・今からの提出でも許してもらえるでしょうか？？？

*1:

参考記事一覧:
1. 計算化学.com ドッキングシミュレーションのやり方(AutoDock vina）
2. 計算化学.com AutoDock のスコアファンクションはどれくらい正確か？CASF-2013 ベンチマーク
3. 香川大学　ケミカルバイオロジー研究室 Docking simulation
4. Slide Share AutoDock_vina_japanese_ver.3.0

*2:AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization and multithreading Trott, O., and Olson, A.J., J. Comput. Chem. 2010(31)455

*3:Uehara, S., Tanaka, S., CICSJ Bulletin 2016(34)10

*4:「GPU を用いたエネルギー計算の並列化によるドッキングシミュレーションの高速化」東京工業大学　佐々木崇浩宇田川拓郎関嶋政和

*5:タンパク質計算科学基礎と創薬への応用　神谷成敏・肥後順一・福西快文・中村春樹　著　共立出版 2009年初版１刷の記述も参考にしています。

*6:Mehler, E.L. and Solmajer, T. "Electrostatic effects in proteins: comparison of dielectric and charge models" Protein Engineering 1991(4)903

*7:Jay Ponder MOLECULAR MODELING (Chemistry 478) Spring Term 2014
Lecture 20 Small Molecular & Protein Docking
Lecture 21 Tutorial & Introduction to AutoDock
あるいは
Marc A. Marti-Renom / Structural Genomics Unit / Bioinformatics Department / Prince Felipe Resarch Center (CIPF)
Docking of small molecules. AutoDock.
などを参考にしました。

*8:What is the format of a PDBQT file?

*9:タンパク質計算科学基礎と創薬への応用　神谷成敏・肥後順一・福西快文・中村春樹　著　共立出版 2009年初版１刷

*10:Chimeraの使用にあたって以下のファイルを参考にしました。
* 複合体構造モデリング
 東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム
寺田透
Autodock tutorial VINA with UCSF Chimera
José R. Valverde
CNB/CSIC

*11:Wójcikowski, M., Zielenkiewicz, P., and Siedlecki, P. Open Drug Discovery Toolkit (ODDT): a new open-source player in the drug discovery field Journal of Cheminformatics 2015(7)26

2019-03-24

LBVSしようとする話

さて、ファーマコフォアを用いた化合物の絞り込みも終わりました・・・と言いたい所ですが、またしても間違いに気づきました。ファーマコフォアサーチに用いたsdfは、オルト位置換基で絞り込む前の、ビフェニルで絞り込んだ段階のものでした！！！

順番は入れ替わってしまいますが、ファーマコフォアスクリーニング後の化合物群に対して部分構造（オルト位置換基）で絞り込みをかけても結果は変わらない、、、はず。。。ということで、改めて部分構造のフィルターをかけます。

from rdkit import Chem
from rdkit.Chem import AllChem

ortho_biphenyl= Chem.MolFromSmarts('c1ccc(c(*)c1)c1ccccc1')
P_Matched2_suppl = Chem.ForwardSDMolSupplier('./P_Matched2.sdf')
P_Matched2_mols = [mol for mol in P_Matched2_suppl if mol is not None]

print(len(P_Matched2_mols))
#997

#オルト置換ありのリスト
ortho_cpds = []
#オルト置換なしのリスト
w_o_ortho = []

for mol in P_Matched2_mols:
    if mol.HasSubstructMatch(ortho_biphenyl):
        ortho_cpds.append(mol)
    else:
        w_o_ortho.append(mol)

print('Number of ortho substituted compounds', len(ortho_cpds))
# Number of ortho substituted compounds 494
print('Number of ortho unsubstituted compounds', len(w_o_ortho))
# Number of ortho unsubstituted compounds 503

ファーマコフォアスクリーニングで残った化合物（997個）のうち、オルト位に置換基を有するものは約半数の494個でした。こちらを「P_Matched_ortho.sdf」という名前で保存しておきました。

さて、ようやく今回の本題に、、、

これまで既知の活性化合物の情報をもとに

部分構造（オルト置換ビフェニル）の有無
ファーマコフォア（Aromatic x2、Acceptor x1）を満たすか否か

と、すこしずつ構造を抽象化（？）しながら数を絞り込んできたことになります。次のステップとして、さらに抽象化した基準で、スコアとして評価する、ということを目指したいと思います。

具体的には

フィンガープリントの使用（複数種類）
Tanimoto係数による評価

を行います。複数のフィンガープリントを使用するのは、前回の記事で参照した DeNA 望月氏の資料で単一のものを用いるよりも、統合した方が良い結果が出た、との記載があったからです。

1. フィンガープリントによるスコア化
2. 上位化合物の検証
3. まとめ

1. フィンガープリントによるスコア化

1-1. スコア化手順

具体的な手順としては以下とします。

活性化合物群からファーマコフォア基準を満たすものを選ぶ
フィンガープリントを計算
複数の活性化合物で共通するビットを重要なビットとして抜き出す
3.で抜き出したビットのベクトルと比較した類似度をスコアとする
複数のフィンフガープリントにおけるスコアを統合

1-2. step 1. 活性化合物の選別

まずは活性化合物群から、前回スクリーニングに用いたファーマコフォア基準を満たすものを選別します。すべての化合物を用いないのは、すでにライブラリ側の絞り込みをおこなっているため、条件を等しくした方がファーマコフォアに明示的に取り込めなかった重要な情報を反映させることができるのではないか、と考えたからです。

（ちょっと自分でも何を言っているかよく分からないですが、そんな気がしました。）

同じことの繰り返しなので詳細は省きますが、共闘用シェアデータから抜き出した、環状ペプチドを除く低分子 39個のうち、24個の分子で２段階のファーマコフォアスクリーニングの条件を満たしました。

この化合物群に、共結晶構造から抜き出したリガンド6個を加えて、合計30個をフィンガープリントを計算する対象の活性化合物群とします。(fp_molsというリストに格納しました)

1-3. step 2. フィンガープリントを計算

つづいてフィンガープリントを計算します。以下の６つを使います（Morganを2種類使用）。

MACCS keys RDKit documentation / Python API
RDKit RDKit documentation / Python API / Python API
Morgan(Circular) RDKit documentation / Python API
- ECFP-like (Morgan(radius=2))
- FCFP-like (Morgan(radius=2, ,useFeatures=True))
AtomPairs RDKit documentation / Python API
Avalon Python API

これ以外に、Topological Torsions （ RDKit documentation / Python API )も候補にあったのですが、ビット数が多すぎて異音、発熱とともにjupyter notebookのカーネルが落ちたので省きました。

計算にあたっては以下の記事、資料を参考とさせていただきました。

参考1: 化合物をベクトルにして比較しプロットする *1
参考2: AI創薬のためのケモインフォマティクス入門 6章 *2
参考3: 化学の新しいカタチ RDKitでフィンガープリントを使った分子類似性の判定 *3

maccs_fps = [AllChem.GetMACCSKeysFingerprint(mol) for mol in fp_mols]
rdkit_fps = [AllChem.RDKFingerprint(mol) for mol in fp_mols]
ECFP_fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048) for mol in fp_mols]
FCFP_fps = [AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048, useFeatures=True) for mol in fp_mols]
Apair_fps = [AllChem.GetAtomPairFingerprint(mol) for mol in fp_mols]

from rdkit.Avalon import pyAvalonTools
avalon_fps = [pyAvalonTools.GetAvalonFP(mol) for mol in fp_mols]

1-4. step 3. 共通するビットの抜き出し

ついで各フィンガープリントから複数の活性化合物で共通するビットを抜き出します。

まずはフィンガープリントをnumpyのarrayに変換します。(上記、参考記事1 をご参照ください)
MACCS keysを使って、試し実験を行います。

import numpy as np

#関数作成（コピペ）
def fp2arr(fp):
    from rdkit import DataStructs
    arr = np.zeros((1,))
    DataStructs.ConvertToNumpyArray(fp, arr)
    return arr

#MACCS keysについてアレイを作成、行列化する
maccs_fpMtx = np.array([fp2arr(fp) for fp in maccs_fps])

この行列を使って、以前の記事でPLIFを作成した時と同様の可視化を行おうと思います。

import pandas as pd

#行列からDataFrameへ
maccs_df = pd.DataFrame(maccs_fpMtx)

#各列の合計値をもとめ、DataFrameに追加
sums = maccs_df.sum()
maccs_df = maccs_df.append(sums, ignore_index=True)

#Heatmapで可視化
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(20, 10))
sns.heatmap(maccs_df)

f:id:magattaca:20190324164344p:plain

少々見づらいですが、縦軸上から「0 ~ 29」が化合物(30個)、一番下の行「30」が各列の合計値です。横軸はMACCS keyの各ビット(１６７個)です。数に応じて着色されていて、一番下の行をみるとMACCS keyのビットの番号が大きいもので、特に多数の化合物でビットが立っている傾向があることがわかります。

予想以上にわかりにくいので、最後の行を抜き出して棒グラフにします。

plt.figure(figsize=(20, 10))
sums.plot(kind='bar', color='k', alpha=0.7)

f:id:magattaca:20190324165011p:plain

要約統計量を眺めて見ます。

sums.describe()
"""
count    167.000000
mean      10.377246
std       11.409474
min        0.000000
25%        0.000000
50%        3.000000
75%       21.000000
max       30.000000
dtype: float64
"""

各ビットは30個の化合物のうち、平均して10個でビットがたっているようです。10個以上の活性化合物でビットが立っていることを基準として、重要なビット、と捉えることにしてみます。

maccs_common_bit = []

for x in sums:
    if x >= 10:
        maccs_common_bit.append(1)
    else:
        maccs_common_bit.append(0)

print(sum(maccs_common_bit))
#75

MACCS keysフィンガープリントのうち、活性化合物群で重要と考えられるビットのベクトルが手に入りました。75 個 (167 個中)から構成されているようです。

上記処理を関数としますが、いちいちデータフレームにするのは意味がなさそうなので合計値と、ビットのアレイをそれぞれ作ることとします。まずはMACCS keysを例に手順を確認してみます。

maccs_size = len(fp2arr(maccs_fps[0]))
#合計値を格納するnumpyのarrayを用意
sum_arr = np.zeros(maccs_size)

#numpy行列の足し算
for i in maccs_fps:
    arr = fp2arr(i)
    sum_arr += arr

print(sum_arr.mean())
#10.377245508982035

先ほどと同じ平均値を返しました。DataFrameを介さずにすみそうです。
平均値の小数点以下を切り捨てた整数部分を閾値とするため、mathのfloorを使います

import math
maccs_threshold = math.floor(sum_arr.mean())
print(maccs_threshold)
#10

閾値を求めるところまで確認できたので関数とします。フィンガープリントのサイズと、閾値、各ビットで30個の化合物で立っている数、そして平均以上の数立っているビットのベクトルの４つを戻り値とします。

def common_bit(fps):
    fp_size = len(fp2arr(fps[0]))
    sumarr = np.zeros(fp_size)
    for i in fps:
        arr = fp2arr(i)
        sumarr += arr
    mean = sumarr.mean()
    threshold = math.floor(mean)
    
    common_bit = []
    for x in sumarr:
        if x >= threshold:
            common_bit.append(1)
        else:
            common_bit.append(0)
    return fp_size, threshold, sumarr, common_bit

フィンガープリント6種に対して適用します。

fingerprints = [maccs_fps, rdkit_fps, ECFP_fps, FCFP_fps, Apair_fps, avalon_fps]
fp_names = ['maccs_fps', 'rdkit_fps', 'ECFP_fps', 'FCFP_fps', 'Apair_fps',  'avalon_fps']

#各返り値を格納するためのリストを用意
fp_size_list = []
threshold_list = []
sumarr_list = []
common_bit_list = []

#計算実行
for fps in fingerprints:
    sz, th, ar, cb = common_bit(fps)
    fp_size_list.append(sz)
    threshold_list.append(th)
    sumarr_list.append(ar)
    common_bit_list.append(cb)

各フィンガープリントのサイズと閾値を確認してみます。

for name, size, threshold in zip(fp_names, fp_size_list, threshold_list):
    print(name, ':', size, ':', threshold)

"""
maccs_fps : 167 : 10
rdkit_fps : 2048 : 14
ECFP_fps : 2048 : 0
FCFP_fps : 2048 : 0
Apair_fps : 8388608 : 0
avalon_fps : 512 : 10
"""

ECFP、FCFP、Atom_pairでは閾値が0となってしまっています。ビットベクトルのサイズの大きさの割に立っているビットが少なく（スパース？）、平均値がとても小さな値となってしまっているのでしょうか？？？

これでは閾値として不適合なので、フィンガープリント全体ではなく、合計値がゼロではないビットのみについての平均値を代わりに用いたいと思います。

MACCS keysで実験してみます。

#np.nonzero()は戻り値がインデックスであることに注意
non_zero_sum = sum_arr[np.nonzero(sum_arr)]
print(non_zero_sum.mean())
#16.349056603773583

今度は閾値が16という結果となりました。MACCS keyにおいて非ゼロのビットのみを考慮すると、30個の化合物群のうち、平均して半数で同一のビットが立っているということになります。（この解釈で正しいでしょうか？自信がありません。）

上記の新しい閾値をとりいれた関数に書き換えます。

def common_bit2(fps):
    fp_size = len(fp2arr(fps[0]))
    sumarr = np.zeros(fp_size)
    for i in fps:
        arr = fp2arr(i)
        sumarr += arr
    NonZeroArr = sumarr[np.nonzero(sumarr)]
    NonZeroMean = NonZeroArr.mean()
    threshold = math.floor(NonZeroMean)
    
    common_bit = []
    for x in sumarr:
        if x >= threshold:
            common_bit.append(1)
        else:
            common_bit.append(0)
    return fp_size, threshold, sumarr, common_bit

実行の部分は同じなので省略し、結果のみ記載します。

for name, size, threshold in zip(fp_names, fp_size_list2, threshold_list2):
    print(name, ':', size, ':', threshold)
"""
maccs_fps : 167 : 16
rdkit_fps : 2048 : 14
ECFP_fps : 2048 : 4
FCFP_fps : 2048 : 5
Apair_fps : 8388608 : 14
avalon_fps : 512 : 11
"""

全て閾値が0ではなくなりました！スパース性を解決した！
レイドバトルは非ゼロサムゲーム！（知らんけど）

くどいですが再度MACCS keysの値を確認します。

print('MACCS NonZero: ', len(non_zero_sum))
#MACCS NonZero:  106
print('MACCS base bits: ', sum(common_bit_list2[0]))
#MACCS base bits:  61

MACCS keyに関しては「非ゼロのビットが106個で、このうち61個で閾値よりを満たす」、となりました。先の棒グラフを見てもそこまでおかしな値とはなっていなさそうです。

取り急ぎ、手に入れたビットベクトルを基準としてもちいることとし、辞書に格納しておきます。

base_BitVect_dict={}

for k, v in zip(fp_names, common_bit_list2):
    base_BitVect_dict[k] = v

1-5. step 4. 類似度の指標を作成

基準となるビットベクトルが用意できたので、このベクトルに対してスクリーニング対象の化合物のビットベクトルの類似度を計算したいと思います。

タニモト係数を使いたいのですが、numpyのarrayとしており、RDKitのオブジェクトではなくなってしまったのでDataStructs.TanimotoSimilarity()を使えるかわかりません。

Tanimoto類似度の説目を見る限り、ビットが立っている要素について積集合を和集合で割れば良いということのようですから、比較対象のビットベクトル同士を足し合わせて、「2の要素数」を「1 or 2の要素数」で割ってしまえば良さそうです。

まずはファーマコフォアで絞り込んだスクリーニング対象の化合物（オルト位置換も考慮済みのもの）をSDFから読み込み、scr_molsというリストに格納しました。

一番最初の分子をテスト化合物として取り出し、各フィンガープリントを計算したのち、アレイに変換し、辞書型として格納します。

test_mol = scr_mols[0]

#関数を作成
def array_dict(mol):
    maccsFP_arr = fp2arr(AllChem.GetMACCSKeysFingerprint(mol))
    rdkitFP_arr = fp2arr(AllChem.RDKFingerprint(mol))
    ECFP_arr = fp2arr(AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048))
    FCFP_arr = fp2arr(AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048, useFeatures=True))
    Apair_arr = fp2arr(AllChem.GetAtomPairFingerprint(mol))    
    avalon_arr = fp2arr(pyAvalonTools.GetAvalonFP(mol))
    
    tmp_dict = {}
    fp_names = ['maccs_fps', 'rdkit_fps', 'ECFP_fps', 'FCFP_fps', 'Apair_fps',  'avalon_fps']
    fp_arrs = [maccsFP_arr, rdkitFP_arr, ECFP_arr, FCFP_arr, Apair_arr, avalon_arr]
    for k, v in zip(fp_names, fp_arrs):
        tmp_dict[k] = v
    
    return tmp_dict

#テスト分子で実行
test_array_dict = array_dict(test_mol)

MACCS key に関して、テスト分子と、先に作成した基準ビットベクトルでTanimoto係数を計算してみます。

#それぞれのビットベクトル（array）を読み込む
base_maccsAR = base_BitVect_dict['maccs_fps']
test_maccsAR = test_array_dict['maccs_fps']

#array同士を足し算
maccsAR_sum = base_maccsAR + test_maccsAR

#各値の数を確認
test_count_1 = np.sum(maccsAR_sum == 1)
print('1 :', test_count_1)
# 1 : 32
test_count_2 = np.sum(maccsAR_sum == 2)
print('2 :', test_count_2)
# 2 : 36

#Tanimoto類似度を計算
test_Tanimoto = test_count_2 / (test_count_1 + test_count_2)
print('Tamimoto :', test_Tanimoto)
#Tamimoto : 0.5294117647058824

うまくいったようです。関数化します。

def TaniSimCal(arr1, arr2):
    arr_sum = arr1 + arr2
    count1 = np.sum(arr_sum == 1)
    count2 = np.sum(arr_sum == 2)
    Tanimoto = count2 / (count1 + count2)
    return Tanimoto

スクリーニング対象の化合物について、6種類のフィンガープリントのそれぞれのTanimoto類似度を計算します。

Scores = [MACCS_Scores, RDkit_Scores, ECFP_Scores, FCFP_Scores, Apair_Scores, Avalon_Scores]

MACCS_Scores = []
RDkit_Scores = []
ECFP_Scores = []
FCFP_Scores = []
Apair_Scores = []
Avalon_Scores = []

for mol in scr_mols:
    #スクリーニング対象のFPのアレイを作成（辞書）
    fp_arr_dict = array_dict(mol)
    
    #基準の各ビットベクトルに対するTanimoto類似度を計算
    for fp, sc in zip(fp_names, Scores):
        base_arr = base_BitVect_dict[fp]
        scr_arr = fp_arr_dict[fp]
        
        TaniScore = TaniSimCal(scr_arr, base_arr)
        
        #各フィンガープリントのスコアリストに追加
        sc.append(TaniScore)

スクリーニング対象４９４化合物の、６種類のフィンガープリントについて、活性化合物群から導いた基準ビットベクトルに対する類似度をプロットしてみました。

for i in range(6):
    plt.subplot(2, 3, i+1)
    plt.hist(Scores[i], bins=20, range=(0, 1))
    plt.title(fp_names[i])

f:id:magattaca:20190324173142p:plain

1-6. step 5. 統合したスコアの作成

上記で求めた値を統合して、各化合物の統合スコアとしたいのですが、このままではフィンガープリントごとにスコアの分布のばらつきが大きいです。Tanimoto類似度にこのような処理をするのが適切かわかりませんが、各フィンガープリントにおいてスコアを標準化（平均 0, 分散 1）したのち、和をとった値を統合したスコアとしたいと思います*4。

import statistics

Scores_STD = []
for l in Scores:
    l_mean = statistics.mean(l)
    l_stdev = statistics.stdev(l)
    s = [(i-l_mean) / l_stdev for i in l]
    Scores_STD.append(s)

標準化したスコアを先と同様にプロットします。

for i in range(6):
    plt.subplot(2, 3, i+1)
    plt.hist(Scores_STD[i], bins=20)
    plt.title(fp_names[i])

f:id:magattaca:20190324173742p:plain

標準化が完了したので、各スコアを足しあわせて統合スコアとします。なんとなく0以上の値にしたいので最小値分シフトさせます。

Integrated_SC_arr = np.zeros(494)
for i in Scores_STD:
    Integrated_SC_arr += np.array(i)

#最小値を求める
int_min = min(Integrated_SC_arr)
#最小値の絶対値分シフト
Integrated_SC = [i + abs(int_min) for i in Integrated_SC_arr]

#グラフ化
plt.hist(Integrated_SC, bins=20)

f:id:magattaca:20190324174003p:plain

スコア化が完了したので、SetDoublePropを使って、Molオブジェクトにプロパティ（Similarity_Score）として持たせ、SDFで出力しました（SCR_compounds_SimScore.sdf）。

2. 上位化合物の検証

スコアの付与とデータの書き出しが完了しました。今回設定したスコアで、上位となった化合物はどのようなものか、確認して見たいと思います。PandasToolsでDataFrameに読み込みます。

from rdkit.Chem import PandasTools

SIM_df = PandasTools.LoadSDF('./SCR_compounds_SimScore.sdf')

スコアの順に並べ替えたいので、typeをstrからfloatに変換しておきます。

SIM_df['Similarity_Score'] = SIM_df['Similarity_Score'].astype(float)

新しく追加したスコアをつかって降順にソートします。

SIM_df_st = SIM_df.sort_values('Similarity_Score', ascending=False)

スコアの良いものから１０個取り出してスコアととともに描画します。

from rdkit.Chem import Draw

#Molオブジェクトの取り出し
Top_10 = SIM_df_st.iloc[:10, 8]
#スコアの取り出しとstrへの変換
Top_10_scores = SIM_df_st['Similarity_Score'][:10]
legends = [str(i) for i in Top_10_scores]

Draw.MolsToGridImage(Top_10, legends = legends)

f:id:magattaca:20190324174512p:plain

かなりそれらしい化合物が選択できてきているのではないでしょうか？

少し想定外だったのは、オルト位に置換基の入ったビフェニル構造でスクリーニングした際に、置換基の大きさ（原子数？）をしていなかったためか、オルト位からそのまま分子が伸長したような化合物が上位にヒットしてきているということです。

例えばこれまで何度かみてきた下の構造のように、オルト位置換基をもちつつ、他の位置（メタ位など）からさらに分子が伸長している、というのがぼんやりとした理想だったのですが、、、

f:id:magattaca:20190324174635p:plain

3. まとめ

今回は、リガンドをベースにしたスクリーニングを行うにあたって、これまでよりも抽象化し、フィンガープリントを使って見ました。また、使用にあたって

複数の活性化合物の情報の取り組み
複数の種類のフィンガープリントの組み合わせ

というのを試みて見ました。
上位ランク化合物は、なんとなく頭の中で「こういう化合物がひっかかったらいいな」と思っていた構造とはことなりました。適切な情報にまで落とし込むことと、それを表現することがまだまでできていないようです。

ひとまずLBVSっぽいということにしたい！　

ところで、創薬レイドバトルの候補化合物絞り込みトップ500と、トップ10の提出みたいなのですが、すでに500個ない。どうしよう・・・

それっぽいことをやって見たいというだけで、スコアの標準化や統合（アンサンブル？）などかなり適当なことをしています。間違っている点等あればご指摘いただければ幸いです。

*1:化合物をベクトルにして比較しプロットする https://qiita.com/Mochimasa/items/f1b60246ece7da46f6a9

*2: AI創薬のためのケモインフォマティクス入門 https://github.com/Mishima-syk/py4chemoinformatics/blob/

*3: 化学の新しいカタチ「RDKitでフィンガープリントを使った分子類似性の判定」https://future-chem.com/rdkit-fingerprint/

*4: Pythonで正規化・標準化(リスト、NumPy配列、pandas.DataFrame) https://note.nkmk.me/python-list-ndarray-dataframe-normalize-standardize/

2019-03-23

リガンドを重ね合わせてファーマコフォアを作る話

前回の記事でファーマコフォアスクリーニングを実施し、約2500個の化合物までリストを絞り込みました。思った以上に減らなかったので、さらなる絞り込みのため、ファーマコフォアを作成した流れを振り返りたいと思います。

リガンド-タンパク共結晶構造を取得（PDB, 6構造）
PLIPを使って相互作用を解析
RDKitでPLIF作成
複数の構造で相互作用に関与する残基に着目
4.の残基と近接するリガンドの部分構造を選定（フィーチャー、３つ）
5.のフィーチャーの位置関係からファーマコフォアを作成

このうち、Step 5、Step 6の部分では、共結晶構造１つ（PDB id: 5NIX）をもとに作成をおこなっています。課題として「フィーチャー３つの選択基準がかなり恣意的だった」、ということが挙げられると思います。そこで今回はもう少し、根拠（？）をもってフィーチャーを選択できるか、残りの共結晶構造５つを眺めたいと思います。

1. 今回の流れ
2. Pymolでアラインメント
3. リガンドの座標を出力
4. 重なりの比較
5. ファーマコフォアを再作成
6. ファーマコフォアサーチ
- 6-1. サーチ
- 6-2. 確認
7. まとめ

1. 今回の流れ

といっても、ただ単に眺めるだけでは仕方ないので、今回は「複数の構造を重ね合わせて眺める」という点に主眼を置きたいと思います。

モチベーションは２つ・・・

AI創薬のためのケモインフォマティクス入門 *1 6章、SBDDの項目で異なる化合物の複合体結晶構造の重ね合わせが紹介されていたこと
ファーマコフォアについてのよもやま話 *2 において、結晶構造以外のアプローチとしてリガンド重ね合わせのアプローチが紹介されていたこと

です。折角、異なるリガンドを含む共結晶構造が複数あるのだから、タンパクごと重ね合わせて、そこに含まれているリガンドの重なりを比較してみれば、リガンドのみの重ね合わせよりもより実際の結合のポーズを意識したヒントが得られるかもしれません。

ということで、まずは

Pymolで複合体構造をアラインメント
アラインメント後のリガンドの座標を出力

をやってみます。

2. Pymolでアラインメント

用いる複合体構造を再掲します。①、②、③のグループは以前の記事でPLIFをベースに分けてみたものです。

Group	PDB entry	Chain	リガンド	リガンドID	リガンド分子量	文献
①	5J89	A/B/C/D	BMS-202	6GX	419.52	Oncotarget 2016(7)30323
	5N2D	A/B/C/D	BMS-37	8J8	448.55	J. Med. Chem. 2017(60)5857
②	5J8O	A/B	BMS-8	6GZ	494.42	Oncotarget 2016(7)30323
	5N2F	A/B	BMS-200	8HW	497.49	J. Med. Chem. 2017(60)5857
③	5NIU	A/B/C/D	BMS-1001	8YZ	594.65	Oncotarget 2017(8)72167
	5NIX	A/B/C/D	BMS-1166	8YQ	643.13	Oncotarget 2017(8)72167

pymolの設定などに関しては或る化みす途さんのブログ*3を参考にさせてきただきました。

具体的には、以下3ステップです。

PDBファイル6つをpymolで読み込む
chain C/Dを除いてA/Bのみにする
「remove (chain C,D)」と打ち込む
アラインメントをとる
例)「align 5j89, 5nix」
意) 5j89を5nixに重ね合わせる

そのままの重ね合わせではうまくいかなかったため、ステップ2（chain A/Bのみを残す）を行なっています。全て5nixに重ね合わせました。

次いで、下記の設定で描画を行いました。

set cartoon_transparency, 0.8
リガンドを中心に持ってくる（Command+左クリック）
水を消す（object panelのallでH(hide)からwatersを選択）

以下のような見た目となりました。比較は後回しにして、座標の出力へと進みます。

f:id:magattaca:20190323123633p:plain

3. リガンドの座標を出力

ついで、重ね合わせからリガンドの座標のみを取り出します。以下のように一つずつpdbファイルで出力しました。

構造を表示
リガンドを選択　→　(sele) になる
pdbファイルで出力「save aligned_5nix.pdb, 'sele', state=-1, pdb」

「aligned_5nix.pdb」という名前で、5nixのリガンドのみを含むpdbファイルが出力されました。

出力の方法に関してはPymolWikiの Save *4を参照しました。

「save file [,(selection) [,state [,format]] ]」 という形でつかう
formatとしてpdb以外にもmol、sdfなど多数選択可能
stateは「Default = -1」（現在の状態のみを保存）

stateに関しては、formatを指定する場合に省略すると、引数がおかしいと怒られたので、Defaultで良くても書き加える必要がありそうです。

出力するリガンドを選択する際に、pymolのwindow上でクリックすると、意図しない残基を選択してしまうことがあったので、「Sequence」でリガンドIDを選択した方がやりやすいと思います。

f:id:magattaca:20190323123812p:plain

出力したpdbファイルを再度読み込んで確認します。

f:id:magattaca:20190323123853p:plain

アラインメントの情報（変換された座標）を保ったまま、構造が出力されていそうです。

今度はまとめてsdfで出力して見ます。拡張子をファイル名に記載していれば自動でformatを認識して処理してくれる、ということだったので
「save aligned_ligands.sdf」
とだけ、打ち込みました。

ちなみに、pymol上で座標を確認するには目的の構造を選択した後、

xyz = cmd.get_coords('sele', 1)
print xyz

とすることで、各原子のxyz座標の３要素のリストを要素とする、リストのリストとして座標が出力されました。PyMolWikiの Get Coordinates I *5によると、他にも色々と方法があるそうです。

4. 重なりの比較

準備ができたのでリガンドの重なりの比較を行なっていきたいと思います。

4-1. 5J89と5N2D (グループ１)

5J89(リガンドID: 6GX)と 5N2D(リガンドID: 8J8)を表示させました。前回もちいたファーマコフォアの３点をメインとする構造のように見えます。

f:id:magattaca:20190323124033p:plain

芳香族の重なりはかなり良いようですが、末端の親水性領域についてはズレが大きくなっています。この辺りは溶媒露出面にもなっていそうなので、許容される動きの幅が大きいのかもしれません。

・・・ということは、この親水性部位の座標を狭い範囲で厳格に決めすぎてファーマコフォアを作ることはよくない・・・、ということでしょうか？フィーチャー３点の一つとして３角形をつくる前にもう少し慎重になった方が良かったかもしれません。

4-2. 5N2D (グループ１) と5N2F (グループ2)

同一の論文で報告されている5N2D と 5N2F(リガンドID: 8HW)を表示させました。

f:id:magattaca:20190323124114p:plain

5N2Dと5N2Fの大きな違いは、後者において図右端のフェニル環にさらに環構造が付加していることです。これにともない、タンパク側の chain A 残基56 Tyrの位置が大きく動いています。押しのけられた、、、という感じでしょうか？（Induced Fit？）

残念ながら構造が報告されている論文（J. Med. Chem. 2017(60)5857）*6は読めなかったのですが、オープンアクセスの Oncotarget 2017(8)72167で、構造の伸長に伴うTyrの動きが議論されていました。

4-3. 5J89(グループ１) と5J8O (グループ2)

もう一つのペア、5j89 と 5j8O(リガンドID: 6GZ)を表示させました。

f:id:magattaca:20190323124149p:plain

興味深いことに、このペアではビフェニル骨格で重なっているものの対称的な位置に広がっており、リガンド分子全体としての重なりはとても悪くなっています。リガンドベースではなく、タンパク質を基準としたアラインメントになっているため、このような結果になっていると思われます。共結晶構造がPD-L1の対称的な二量体となっていることを踏まえると、確かに対照的な位置に結合してもおかしくないという気もします。（結晶の群とか対称とかわかりません。すみません）

予期せぬ重ね合わせ構造でびっくりしましたが、これはこれで示唆に富む結果だと思います。なぜかというと、共闘用シェアデータに抜き出されている特許構造をながめても、最近のものに移るにつれて、中心に疎水性、両末端に親水性基をもつ対称性のある構造へと変化しているように見えるからです。（構造式は MarvinSketch 18.24.0で作成）

f:id:magattaca:20190323124305p:plain

Twitterで拝見したDeNA 望月氏の講演資料「コンペティションから見るAI創薬」において、創薬レイドバトルに関して化合物の線対称に近い構造に着目しているとのコメントがありました（スライド p27）。*7 ご自身で記述子を開発されているとのこと、一体どんな工夫をされているのか？とても結果が楽しみです。

"コンペティションから見るAI創薬" というタイトルで DeNA の AIエンジニアが日本オミックス医学会シンポジウムに登壇しました。

コンペという軸で、AI創薬の近年の動向について紹介しています！https://t.co/EQx48elkL6
— DeNA Tech (@DeNAxTech) 2019年3月20日

対称性のある構造の医薬品

ちょっと脱線しますが、上述の対称性のある構造を見た際にC型肝炎治療薬を思い出しました。C型肝炎治療薬といえば、ハーボニー配合錠が薬価や偽造品流通などで話題となりました。こちらはNS5A阻害剤レジパスビル　と　NS5B阻害剤ソホスブビルを含む合剤となっています。 NS5A阻害剤は複数ありますが、例えば以下のように、分子量が大きく対称的な構造をとっている、という印象です。

f:id:magattaca:20190323124338p:plain

経口薬というのが驚きの巨大分子ばかりです。

Wikipedia(英語)の記事「Discovery and development of NS5A inhibitors」に、BMSの研究者によるNS5A 結晶構造の論文（PDB id: 4CL1）（Protein Sci. 2014(23)723）が引用されていました。共結晶構造ではなかったのですが、

２量体の構造解析
ドッキングシミュレーション

等の結果を示すFigureがありました（・・・詳しく読んでません。すみません）。
２量体の形成する疎水性の溝に、対称性の高いリガンド（Daclatasvir : DCV）が結合している図は、これまで見てきたPD-L1の共結晶構造を彷彿とさせるもので、ひょっとしてNS5B阻害剤、PD-L1にもヒットするのでは？？？という気もしてきます。・・・まあ、無理か。。。

当てずっぽうはさておき、LipinskiのRule of 5からは真っ先に弾かれてしまいそうなこれらの分子が経口薬として上市に至っている、という成功例がPD-1/PD-L1 結合阻害剤においても巨大化、対称化の流れを進める背景にもあったりするのかな、と想像してしまいます。創薬の現場からみたらどうなんでしょう？？？

4-4. 5J89(グループ１) と5NIU、5NIX (グループ3)

閑話休題。

残りの5NIU(リガンドID: 8YZ)、5NIX(リガンドID: 8YQ)を5J89と重ね合わせます。

f:id:magattaca:20190323124409p:plain

全体として重なりがよく、グループ3の大きな違いはグループ1、2と異なりベンジル基が左下方向にむかって付け足されていることです。文献（Oncotarget 2017(8)72167）*8中では、この置換基が新しい相互作用の獲得に寄与しているとの指摘もなされていました。

重ね合わせ結果の考察

以上、重ねあわせ構造を順番にながめてみました。６つの構造で、前回ファーマコフォアとして選択した３つのフィーチャーを満たしているようにみえ、とりあえずの選択としては十分要件を満たしていたのではないかな？という印象です。

・・・つまり、選択に根拠はあった！（・・・後付け）

共結晶構造中のリガンドをベースに考える場合、更なるポイントとして前回のファーマコフォアに加えて、

右末端フェニルからの置換基伸長（疎水性ポケット）
左下方向への芳香環伸長
左端、親水性基は座標の自由度がある

といった点を考慮していくと良いかもしれません。

5. ファーマコフォアを再作成

5-1. 今回の目標

観察をもとに、更なる絞り込みのためのファーマコフォアを作りたいと思います。
2.左下への芳香環伸長は新しいポケットを埋めるという意味で魅力的ですが、「共闘用シェアデータ」の活性化合物群の流れをみる限り、こちらのポケットは活用しない方向に流れがシフトしているように見えます。そこで、今回は1.右端疎水性ポケットと3. 親水性基の自由度を考慮したファーマコフォアを作成したいと思います。

5-2. Pymolで出力したSDFをRDKitで読み込む

まずは、Pymolで出力したファイルをRDkitで読み込めるか確認して見ます。

from rdkit import Chem
aligned_lig_sup = Chem.SDMolSupplier('./aligned_ligands.sdf')
aligned_ligs = [x for x in aligned_lig_sup if x is not None]

以下のようなエラーがたくさん出てきました。

RDKit ERROR: [20:29:00] ERROR: non-ring atom 14 marked aromatic
RDKit ERROR: [21:10:06] non-ring atom 4 marked aromatic
RDKit ERROR: [21:10:06] ERROR: Could not sanitize molecule ending on line 72

MarvinViewでは開き、構造を確認することができたので、どうもRDKitとpymolの相性が悪そうです。
RDKitのエラーメッセージによると、環に含まれていない原子（non-ring atom）に芳香族性が紐づけられていることが問題のようです。SDFの中身を色々とみた結果、エラーがでた構造ではすべてカルボキシ基を有しており、結合表(bond block)においてこの官能基を構成する原子間の結合タイプが「4 = aromatic」となっていることが原因（？）、という結論に達しました。*9

aromatic bondと紐づけられているためaromatic atomと認識された原子が、ring atomではないためエラーとなり、sanitizeもできなかったと思われます。（推定）

5-3. Pymolで出力したpdbをRDKitで読み込む

SDFの中身を修正するのは私には難しいので、PDBファイルとして出力した個々の構造の検証に移りたいと思います。

from rdkit.Chem import Draw
pdb_5j8o = Chem.MolFromPDBFile('./aligned_ligands/aligned_5j8o.pdb')
Draw.MolToImage(pdb_5j8o)

f:id:magattaca:20190323124646p:plain

こちらは読み込み自体はうまくいきましたが、PDBフォーマットに結合次数が記載されていないため、全て単結合、カルボキシル基に至ってはジェミナルジオール構造のように認識されています。

RDKitメーリングリストのディスカッション *10を参考に結合次数を割り当てたいと思います。

手順としては

smilesからMolオブジェクトを作成（テンプレート）
AssignBondOrdersFromTemplate(refmol, mol)でテンプレートから次数を割り当て

となります。

from rdkit.Chem import AllChem
smi_5j8o = 'Cc1c(COc2ccc(CN3CCCC[C@@H]3C(O)=O)cc2Br)cccc1-c1ccccc1'
template_5j8o = Chem.MolFromSmiles(smi_5j8o)
mol_5j8o = AllChem.AssignBondOrdersFromTemplate(template_5j8o, pdb_5j8o)
Draw.MolToImage(mol_5j8o)

f:id:magattaca:20190323124742p:plain

うまくいったようです。残りの５つも同様に処理します。

#pdb idをkey、smilesをvalueとする辞書を作成
smi_dict = {}
smi_dict['5j89'] = 'COc1nc(OCc2cccc(c2C)-c2ccccc2)ccc1CNCCNC(C)=O'
smi_dict['5n2d'] = 'COc1cc(OCc2cccc(c2C)-c2ccccc2)cc(OC)c1CNCCNC(C)=O'
smi_dict['5n2f'] = 'Cc1c(COc2cc(F)c(CNCC(=O)CC(O)=O)cc2F)cccc1-c1ccc2OCCOc2c1'
smi_dict['5niu'] = 'Cc1cc(CN[C@H](CO)C(O)=O)c(OCc2cccc(c2)C#N)cc1OCc1cccc(c1C)-c1ccc2OCCOc2c1'
smi_dict['5nix'] = 'Cc1c(COc2cc(OCc3cccc(CN)c3)c(C[C@H]3N[C@H](O)C[C@H]3C(O)=O)cc2Cl)cccc1-c1ccc2OC=COc2c1'

mol_5j8o.SetProp('PDB_id', '5j8o')
#結合次数を割り当てたMolオブジェクトを格納するリストを作成
mols = [mol_5j8o]

for k, v in smi_dict.items():
    # aligned_xxxx.pdb という名前（xxxxはPDB id）のファイルから作成
    path = './aligned_ligands/aligned_' + k + '.pdb'
    tmp = Chem.MolFromPDBFile(path)
    template = Chem.MolFromSmiles(v)
    mol = AllChem.AssignBondOrdersFromTemplate(template, tmp)
    # PDB idをプロパティとして持たせる
    mol.SetProp('PDB_id', k)
    mols.append(mol)

#確認
Draw.MolsToGridImage(mols)

f:id:magattaca:20190323124951p:plain

#SDFとして書き出しておきます
writer=Chem.SDWriter('aligned_pdb_molecules.sdf')
writer.SetProps(['PDB_id'])
for mol in mols:
    writer.write(mol)
writer.close()

5-4. 重ね合わせでフィーチャーを眺める

出力したSDFをつかって、pymol + RDKitでフィーチャーを描画して見ます。

import os
from rdkit import RDConfig
fdef = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'))
from rdkit.Chem.Features import ShowFeats
from rdkit.Chem import PyMol

import subprocess
from IPython import display
showfeatpath = os.path.join(RDConfig.RDCodeDir, 'Chem/Features/ShowFeats.py')

v = PyMol.MolViewer()
v.DeleteAll()
process = subprocess.Popen(['python', showfeatpath, '--writeFeats', 'aligned_pdb_molecules.sdf'], stdout=subprocess.PIPE)
stdout = process.communicate()[0]

png=v.GetPNG()
display.display(png)

f:id:magattaca:20190323125155p:plain

この構造をもとに、5-1. 今回のファーマコフォアの狙いと合わせると、下図のようなフィーチャーの選択が考えられそうです。

f:id:magattaca:20190323125259p:plain

"Aromatic" ２つと "Acceptor" 1つの３点という点は以前と同じですが、選択する場所を変えています。

右側疎水性ポケットを意識した"Aromatic"の選択
複数のリガンドの親水性基が集まる場所を意識した"Acceptor"の選択

としています。

Acceptorの位置を決めるにあたっては、複数のリガンドの親水性基の中心にPseudoatomを作成することとしました。具体的には各リガンドで、該当部位の近辺の親水性基をなす原子を選択し、すべて選んだあとで「pseudoatom hyd, sele」とコマンドを打ち込みました。

Aromatic 2点についてもそれぞれ中心にPseudoatomを作成し、meshで表示すると以下の図のようになりました。

f:id:magattaca:20190323125335p:plain

「Wizard → Measurement」で測定した距離と、別に取得した座標情報「例) xyz = cmd.get_coords('hyd', 1)」とを合わせてまとめると以下のようになりました。

f:id:magattaca:20190323125400p:plain

なんだか格好良くなってきました！それっぽいぞ！（意味不明）

この情報を使って、再度ファーマコフォアサーチを行います。

5-5. ファーマコフォアを作成

前回と設定値以外変わってませんが念のため・・・

from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.Pharm3D import Pharmacophore
from rdkit import Geometry

feat_1=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-11.796, 10.999, -22.536))
feat_2=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-7.625, 10.407, -21.653))
feat_3=ChemicalFeatures.FreeChemicalFeature('Acceptor',Geometry.Point3D(1.222, 9.145, -28.891)) 
feats = [feat_1,feat_2,feat_3]
pcophore = Pharmacophore.Pharmacophore(feats) # ファーマコフォアの設定
d_upper = 1.5  # 特性基間の距離の許容範囲(上限値)
d_lower = 0.5 # 特性基間の距離の許容範囲(下限値)

# feat_1とfeat_2の距離 4.3Åを基準に、下限と上限を設定
pcophore.setLowerBound(0,1, 4.3-d_lower)
pcophore.setUpperBound(0,1, 4.3+d_upper)

# feat_2とfeat_3の距離 11.5Åを基準に、下限と上限を設定
pcophore.setLowerBound(1,2, 11.5-d_lower)
pcophore.setUpperBound(1,2, 11.5+d_upper)

# feat_1とfeat_3の距離 14.6Åを基準に、下限と上限を設定
pcophore.setLowerBound(0,2, 14.6-d_lower)
pcophore.setUpperBound(0,2, 14.6+d_upper)

print(pcophore)
"""
                   Aromatic      Aromatic      Acceptor 
     Aromatic         0.000         5.800        16.100 
     Aromatic         3.800         0.000        13.000 
     Acceptor        14.100        11.000         0.000 
"""

新しいファーマコフォアの準備ができました。関数は前回作成したPSearch3というものを用います。

6. ファーマコフォアサーチ

6-1. サーチ

前回、ファーマコフォアサーチを実施したSDFファイルを読み込みます。

PSerched_suppl = Chem.SDMolSupplier('./PSearched.sdf')
PSerched_mols = [x for x in PSerched_suppl if x is not None]
print(len(PSerched_mols ))
#4220

Embedできた数をPScoreというプロパティに格納しているので、こちらの値が1以上のものを、ファーマコフォアにマッチしたもの(P_Matched)、0のものをマッチしなかったもの(P_MissMatched)として別々のリストに分けます。

P_Matched = []
P_MissMatched = []

for x in PSerched_mols :
    PScore = int(x.GetProp('PScore'))
    if PScore >= 1:
        P_Matched.append(x)
    else:
        P_MissMatched.append(x)   
print('P_Matched:', len(P_Matched))
#2463
print('P_MissMatched:', len(P_MissMatched))
#1757

前回のファーマコフォアにマッチした約2500個について、今回新しく作成したファーマコフォアスクリーニングを適用します。

for x in P_Matched:
    num_embeds = PSearch3(x)
    x.SetIntProp('PScore2', num_embeds)

マッチしたものを取り出します。（上と同じなので省略）

#取り出したリストをP_Matched_2、P_MissMatched_2としている。
print('P_Matched:', len(P_Matched_2))
#997
print('P_MissMatched:', len(P_MissMatched_2))
#1466

前回と今回のファーマコフォア両者を満たすものは約1000個でした。「P_Matched2.sdf」という名前でSDFを出力しました。

6-2. 確認

PandasのDataFrameで読み込んで見ます。

from rdkit.Chem import PandasTools
import pandas as pd
PSearched2_df = PandasTools.LoadSDF('./P_Matched2.sdf')

今回のスクリーニング(PScore2)の値の分布を取得して見ます。

print(PSearched2_df['PScore2'].describe())
"""
count     997
unique     16
top         1
freq      254
Name: PScore2, dtype: object
"""

最大値や最小値が返ってくるのを期待していたのですが、どうも様子がおかしいです。データ型を確認して見ます。

PSearched2_df.dtypes
"""
ID                   object
MW                   object
MolLogP              object
NumHAcceptors        object
NumHDonors           object
NumRotatableBonds    object
PScore               object
PScore2              object
ROMol                object
TPSA                 object
original_id          object
dtype: object
"""

PandasToolsで読み込んだものは全てobject型になっています。こちらの記事「pandasのデータ型dtype一覧とastypeによる変換 *11」によると、要素に文字列を含むDataFrameの列はobject型となるらしく、また、化学の新しいカタチさんの記事RDKitのPandasToolsでデータ分析を加速する *12にもPandasToolsで読み込んだ場合は「プロパティが全て「文字列」として認識」されると記載されていました。

astypeで整数型(int)に変換します。

PSearched2_df = PSearched2_df.astype({'PScore': int, 'PScore2': int})
PSearched2_df.describe()
"""
      PScore  PScore2
count 997.000000  997.000000
mean  14.343029   3.289870
std   12.517565   2.394612
min   1.000000    1.000000
25%   5.000000    1.000000
50%   10.000000   3.000000
75%   20.000000   4.000000
max   99.000000   31.000000
"""

無事出力できました。それぞれの最大値の構造を確認して見ます。

#最大値のindexを取得
PScore_max = PSearched2_df['PScore'].idxmax()
PScore2_max = PSearched2_df['PScore2'].idxmax()
#Molオブジェクトを取り出し
max_mol = PSearched2_df.loc[PScore_max, 'ROMol']
max_mol2 = PSearched2_df.loc[PScore2_max, 'ROMol']

#スコアを構造とともにlegendとして描画する準備
legends = []

for x in [PScore_max, PScore2_max]:
    score1 =  str(PSearched2_df.loc[x, 'PScore'])
    score2 =  str(PSearched2_df.loc[x, 'PScore2'])
    legend =  'PScore:' + score1 + '/ PScore2:' + score2
    legends.append(legend)

Draw.MolsToGridImage([max_mol, max_mol2], legends = legends)

f:id:magattaca:20190323131748p:plain

前回と今回で同じものがEmbed数最大となっていたようです。

7. まとめ

今回は、より良いファーマコフォアの作成を目指して複数の共結晶構造を重ねあわせる、というステップを加えて見ました。タンパク質をあつかうにはPyMolが非常に使いやすく、Wikiの充実しているので少しずつ機能がわかってきました。

今回のファーマコフォアマッチングの課題としてはEmbedの数しか考慮しておらず、リガンドの３次元構造を実際に発生、最適化させることまでできていないことです。３次元構造発生　→　ドッキング　みたいな感じに持っていけると良いのですが、、、進捗が悪い。。。

行ったり来たりあちらこちらで躓いていたらとても冗長になってしまいました。おかしな点が多数あると思いますのでご指摘いただければ幸いです。

*1:AI創薬のためのケモインフォマティクス入門　https://github.com/Mishima-syk/py4chemoinformatics/blob/master/pdf/py4chemoinformatics.pdf

*2:ファーマコフォアについてのよもやま話 http://blog.kzfmix.com/entry/1552452149

*3: 話題のXofluzaの変異をPymolで見てみる https://keetaneblog.hatenablog.com/entry/2019/02/08/210016

*4:PymolWiki Save https://pymolwiki.org/index.php/Save

*5:PyMolWiki Get Coordinates I https://pymolwiki.org/index.php/Get_Coordinates_I

*6:https://pubs.acs.org/doi/abs/10.1021/acs.jmedchem.7b00293

*7:コンペティションから見るAI創薬 https://speakerdeck.com/m_mochizuki/ai-drug-discovery-in-the-view-of-competitions

*8:Oncotarget 2017(8)72167 http://www.oncotarget.com/index.php?journal=oncotarget&page=article&op=view&path%5B%5D=20050

*9:化学構造の表記 http://www.metabolomics.jp/wiki/Aritalab:Lecture/Compchem/MOL

*10: https://sourceforge.net/p/rdkit/mailman/message/31844451/

*11:pandasのデータ型dtype一覧とastypeによる変換 https://note.nkmk.me/python-pandas-dtype-astype/

*12: RDKitのPandasToolsでデータ分析を加速する https://future-chem.com/rdkit-pandastools/

2019-03-17

ファーマコフォアスクリーニングをする話

前回の記事で共結晶構造を元にしたRDKitのファーマコフォアの作成を実施しました。いよいよこちらを用いてスクリーニングを実施したいと思います。

ファーマコフォアの準備
ファーマコフォアサーチの関数作成（関数1）
エラーを回避した関数作成（関数2）
不斉点を考慮しない関数の作成（関数3）
- 関数3の検証
スクリーニングの実施
- スクリーニング対象の再確認
- スクリーニング本番
まとめ

ファーマコフォアの準備

前回の再掲ですが、SAR News No.19の吉森氏による記事「Chemoinformatics Toolkits を用いた創薬システム開発におけるラピッドプロトタイピング」での流れを確認します。

座標を使って特性基を定義
距離情報を使ってファーマコフォアを設定
ファーマコフォアサーチ
1. 対象のフィーチャーをそもそも含んでいるか？
2. フィーチャー間の距離が条件を満たすか？
3. 距離情報を拘束条件にして3D構造を発生させる。

このうちファーマコフォアの設定までは、前回のものをそのまま利用します。

f:id:magattaca:20190317081633p:plain

import os
from rdkit import RDConfig
from rdkit.Chem import AllChem
from rdkit import Chem
from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.Pharm3D import Pharmacophore
from rdkit import Geometry
from rdkit.Chem.Pharm3D import EmbedLib
from rdkit.Chem import rdDistGeom

Ffactory = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'))

feat_1=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-7.625, 10.407, -21.653))
feat_2=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-1.851, 11.051, -24.019))
feat_3=ChemicalFeatures.FreeChemicalFeature('Acceptor',Geometry.Point3D(-0.100, 13.586, -28.137)) 
feats = [feat_1,feat_2,feat_3]

pcophore = Pharmacophore.Pharmacophore(feats) # ファーマコフォアの設定
d_upper = 1.5  # 特性基間の距離の許容範囲(上限値)
d_lower = 0.5 # 特性基間の距離の許容範囲(下限値)

# feat_1とfeat_2の距離 6.3Åを基準に、下限と上限を設定
pcophore.setLowerBound(0,1, 6.3-d_lower)
pcophore.setUpperBound(0,1, 6.3+d_upper)

# Acceptor(feat_3)の代表点の選び方によってfeat_2との距離は[5.0~5.6]の値をとる
pcophore.setLowerBound(1,2, 5.0-d_lower)
pcophore.setUpperBound(1,2, 5.6+d_upper)

# 同様にfeat1とfeat_3は[9.8~11.2]の値をとる
pcophore.setLowerBound(0,2, 9.8-d_lower)
pcophore.setUpperBound(0,2, 11.2+d_upper)

ファーマコフォアをpcophore、フィーチャーファクトリーをFfactoryとして用意しました。molオブジェクトが与えられた場合に、ファーマコフォアサーチを行う関数を作成します。

ファーマコフォアサーチの関数作成（関数1）

距離情報を拘束条件にした3D構造がうまく発生できた場合、その構造をファーマコフォア合致構造とすることとします。また、合致構造間での順位づけとして、構造がうまく生成できた数（条件を満たすコンフォマーの数）をスコアとしてもちいることとしたいと思います。生成させる構造の最大数(count)は10としました。

具体的な関数を以下のように作成しました。

def PSearch(mol):
    #水素原子を付加（３次元構造利用のため）
    molH = AllChem.AddHs(mol)
    #水素原子を綺麗に整える座標計算
    AllChem.Compute2DCoords(molH)
    
    #フィーチャーをそもそも持っているか？
    match, mList = EmbedLib.MatchPharmacophoreToMol(molH, Ffactory, pcophore)
    print(match)
    
    if match == True:
        #距離の検証
        #距離行列の取得
        bounds = rdDistGeom.GetMoleculeBoundsMatrix(molH)

        #ファーマコフォアのマッチング 
        pList = EmbedLib.GetAllPharmacophoreMatches(mList,bounds,pcophore) 
        
        #pListのそれぞれについてFeatureとマッチする原子のidを取得する
        num_match = len(pList)
        phMatches = []
        for i in range(num_match): 
            num_feature = len(pList[i])
            
            phMatch = []
            
            for j in range(num_feature):
                phMatch.append(pList[i][j].GetAtomIds())
                
            phMatches.append(phMatch)
        
        #原子のidとファーマコフォアをもとに3次元構造を埋め込む    
        Generated_embeds = []

        for phMatch in phMatches:
            bm,embeds,nFail =EmbedLib.EmbedPharmacophore(molH, phMatch, pcophore,count=5, silent=1)
            Generated_embeds.append(len(embeds))        
        return sum(Generated_embeds)           

    else:
        return 0

関数1の検証とエラー

この関数をもちいて動作確認のため、共闘用シェアデータ）から取り出した分子のうち、前回用いなかったものを含めて検証します*1。

chain_suppl = Chem.SDMolSupplier('./chain_compounds.sdf', removeHs=False)
chain_mols = [x for x in chain_suppl if x is not None]

この化合物群に適用したところ早速以下のエラーで止まりました。

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-40-74609285009c> in <module>()
----> 1 PSearch(test_m3)

<ipython-input-38-904fad3ab3c2> in PSearch(mol)
     34 
     35         for phMatch in phMatches:
---> 36             bm,embeds,nFail =EmbedLib.EmbedPharmacophore(molH, phMatch, pcophore,count=20, silent=1)
     37 
     38             Generated_embeds.append(embeds)

~/.pyenv/versions/anaconda3-4.4.0/lib/python3.6/site-packages/rdkit/Chem/Pharm3D/EmbedLib.py in EmbedPharmacophore(mol, atomMatch, pcophore, randomSeed, count, smoothFirst, silent, bounds, excludedVolumes, targetNumber, useDirs)
    427       for idx, stereo in mol._chiralCenters:
    428         if stereo in ('R', 'S'):
--> 429           vol = ComputeChiralVolume(m2, idx)
    430           if (stereo == 'R' and vol >= 0) or (stereo == 'S' and vol <= 0):
    431             keepIt = False

~/.pyenv/versions/anaconda3-4.4.0/lib/python3.6/site-packages/rdkit/Chem/Pharm3D/EmbedLib.py in ComputeChiralVolume(mol, centerIdx, confId)
   1241   nbrRanks = []
   1242   for nbr in nbrs:
-> 1243     rank = int(nbr.GetProp('_CIPRank'))
   1244     pos = conf.GetAtomPosition(nbr.GetIdx())
   1245     nbrRanks.append((rank, pos))

KeyError: '_CIPRank'

EmbedPharmacophoreの部分で止まっているようです。ComputeChiralVolumeという箇所のようです。

エラーの出た構造はこちら、

f:id:magattaca:20190317083045p:plain

不斉点をもつ構造でエラーが出たようです。

エラーの原因は不斉点？

pythonファイルの該当部分の周辺をみると、nbrは不斉中心の近隣（neigbhbor）の原子を指定しているように見えました。

エラー箇所、

1243     rank = int(nbr.GetProp('_CIPRank'))

は、近隣原子の特性(atom property)を取得した際に、マジックプロパティの一つ_CIPRankの情報が見当たらないと言う内容のようです。

不斉に関する設定なので、CIPはカーン・インゴルド・プレローグ順位則（wikipedia: Cahn–Ingold–Prelog priority rulesの略と推測できます。不斉中心の近隣の原子のCIPRank、ということですから、CIP順位則に従って各原子に割り当てられた番号（rank）、を意味していそうです。

描画を見る限り光学活性体として認識されているのに、この情報が見当たらない、設定されていないというのはなぜ？？？

おそらくですが、一番最初に水素原子を付加（AllChem.AddHs()）した際に、不斉中心（３級炭素）の周りの原子のプロパティ情報（_CIPRank）失われているように見えます。

検証してみたいと思います。

水素原子の取り扱いが原因？

まずは水素原子を付与する前のMolオブジェクトでプロパティを確認します。
不斉中心の情報を取得します。

Chem.FindMolChiralCenters(test_m3)
#[(31, 'S')]

AtomID: 31の原子に不斉 S が割り当てられています。GetNeighbors()を使って近隣の原子のidを取得します。

chiral_center = test_m3.GetAtomWithIdx(31)
nbrs = chiral_center.GetNeighbors()

取得した近隣の原子について、AtomIDと元素記号、問題の_CIPrankを出力します。

for x in nbrs:
    print('AtomID: ', x.GetIdx(), 'Symbol: ', x.GetSymbol(), 'CIPRank: ', x.GetProp('_CIPRank'))

#AtomID:  30 Symbol:  C CIPRank:  4
#AtomID:  32 Symbol:  C CIPRank:  39
#AtomID:  26 Symbol:  N CIPRank:  42

水素原子を付与する前のMolオブジェクトでは問題なくCIPRankが割り当てられています。値の意味はわかりません・・・・（1,2,3,4の順位かと思っていましたが、なんらかの基準でもっと広い範囲の値が割り当てられ大小を比較していそう・・・）

続いて水素を付与した場合の確認をします。

test_m3H = AllChem.AddHs(test_m3)
AllChem.Compute2DCoords(test_m3H)
Chem.FindMolChiralCenters(test_m3H)
#[(31, 'S')]

不斉中心を持つと言うことは認識されていそうです。

chiral_centerH = test_m3H.GetAtomWithIdx(31)
nbrsH = chiral_centerH.GetNeighbors()

for i in nbrsH:
    print('AtomID: ', i.GetIdx(), 'Symbol: ', i.GetSymbol())
#AtomID:  30 Symbol:  C
#AtomID:  32 Symbol:  C
#AtomID:  26 Symbol:  N
#AtomID:  79 Symbol:  H

CIPRankを確認します。

nbrsH[0].GetProp('_CIPRank')

"""
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-99-89378a1786f7> in <module>()
----> 1 nbrsH[0].GetProp('_CIPRank')

KeyError: '_CIPRank'
"""

付与した水素原子だけでなく、もともとあった炭素原子からも_CIPRankの情報が失われていました。

立体化学を再度割り当てれば良いのかと思い、Chem.AssignStereochemistry()を使ってみましたが、_CIPRankの情報は付与されていませんでした。

そもそも水素原子を付加していたのは、３次元構造を扱うためでした。今回、embedできたか否かを判定するだけ（３次元構造の最適化といった処理をしない）ですので、水素原子付加をしなくても良さそうです。

したがって、最初の2つの処理

molH = AllChem.AddHs(mol)
AllChem.Compute2DCoords(molH)

を除いてしまえば良さそうです。

エラーを回避した関数作成（関数2）

上記2つを除いた関数をPSearch2として作成しました。（コード省略）

この関数を用いたところエラーで止まることなく最後まで計算できました。

PScore2_list = []

for x in chain_mols:
    num_embeds = PSearch2(x)
    PScore2_list.append(num_embeds)
    print(num_embeds)
    x.SetIntProp('PScore', num_embeds)

print(PScore2_list)
#[4, 62, 0, 7, 7, 2, 11, 78, 0, 0, 0, 0, 0, 60, 0, 2, 0, 0, 51, 0, 75, 0, 0, 29, 157, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14, 0, 0, 20]

0の数が多いようなので確認してみます。

#リストから要素が0となっている数を求める
print( len([i for i, x in enumerate(PScore2_list) if x == 0]) )
#24

活性化合物群のなかでもembedできた数が0となる構造が多数ありました（24個）。

実行結果で気になる点として、

RDKit DEBUG: [22:55:56] DEBUG: Removed embedding due to chiral constraints.

という表示が多く出力されました。不斉の制限のせいで弾かれている構造が多そうです。かなり荒い３次元構造の取り扱いをしているので、不斉点にそこまで強い制約を加えない方が良い気がします。

不斉点を考慮しない関数の作成（関数3）

今回の計算にあたって、不斉点を無視して処理したいと思います。具体的にはSMILESに変換（オプションでisomericSmiles=Falseとする）して不斉の情報を除いてからMolオブジェクトに戻すという処理を加えます。

#そのままSMILESに変換した場合
test_iso_smi = Chem.MolToSmiles(test_m3)
print('chiral:', test_iso_smi)
#出力: chiral: Cc1c(COc2cc(OCc3cncc(C#N)c3)c(CN3CCCC[C@H]3C(=O)O)cc2Cl)cccc1-c1cccc(OCCCN2CCC(O)CC2)c1C

#立体の情報を除いてSMILESに変換した場合
test_smi = Chem.MolToSmiles(test_m3, isomericSmiles=False)
print('achiral: ', test_smi)
#出力: achiral:  Cc1c(COc2cc(OCc3cncc(C#N)c3)c(CN3CCCCC3C(=O)O)cc2Cl)cccc1-c1cccc(OCCCN2CCC(O)CC2)c1C

#もう一度Molオブジェクトに戻す
test_mol_re = Chem.MolFromSmiles(test_smi)
AllChem.Compute2DCoords(test_mol_re)
#立体の確認
print(Chem.FindMolChiralCenters(test_mol_re))
#空のリストが出力された

isomericSmiles=Falseとしたものでは@Hがなくなっており、またChem.FindMolChiralCenters()が空のリストを返すことから、不斉の情報がなくなっていることが確認できました。

上記の処理を加えた関数を再定義します。

def PSearch3(mol):
    #水素原子を付加しない
    
    smi = Chem.MolToSmiles(mol, isomericSmiles=False)
    mol_re = Chem.MolFromSmiles(smi)
    AllChem.Compute2DCoords(mol_re)
    
    #フィーチャーをそもそも持っているか？
    match, mList = EmbedLib.MatchPharmacophoreToMol(mol_re, Ffactory, pcophore)
    
    if match == True:
        #距離の検証
        #距離行列の取得
        bounds = rdDistGeom.GetMoleculeBoundsMatrix(mol_re)

        #ファーマコフォアのマッチング 
        pList = EmbedLib.GetAllPharmacophoreMatches(mList,bounds,pcophore) 
        
        #pListのそれぞれについてFeatureとマッチする原子のidを取得する
        num_match = len(pList)
        phMatches = []
        for i in range(num_match): 
            num_feature = len(pList[i])
            
            phMatch = []
            
            for j in range(num_feature):
                phMatch.append(pList[i][j].GetAtomIds())
                
            phMatches.append(phMatch)
        
        #原子のidとファーマコフォアをもとに3次元構造を埋め込む    
        Generated_embeds = []

        for phMatch in phMatches:
            bm,embeds,nFail =EmbedLib.EmbedPharmacophore(mol_re, phMatch, pcophore,count=5, silent=1)
            Generated_embeds.append(len(embeds))    
        return sum(Generated_embeds)           
    else:
        return 0

関数3の検証

活性化合物群に適用します。

PScore3_list = []

for x in chain_mols:
    num_embeds = PSearch3(x)
    PScore3_list.append(num_embeds)
    #Molオブジェクトのプロパティに結果を保持させる
    x.SetIntProp('PScore', num_embeds)

今回はRDKit DEBUGは出ませんでした。期待した挙動を示していそうです。 embedできなかった数を数えます。

False_list = [i for i, x in enumerate(PScore3_list) if x == 0]

print(False_list)
#出力: [8, 14, 21, 25, 27, 28, 29, 30, 31, 32, 33, 34]
print(len(False_list))
#出力: 12

今回はembedできた数が0となった構造は12個でした。キラルの制約を無くしたことで、活性化合物群のファーマコフォアベースのヒット率が上がったようです。

embedできなかった化合物について具体的な構造を一つ見てみます。

test_m8 = chain_mols[8]
print(test_m8.GetProp('PScore'))
# 0

Draw.MolToImage(test_m8)

f:id:magattaca:20190317085029p:plain

Property「'PScore'」に値が保持されていることが確認できました。

ipywidgetsをつかって順番に眺めてみます。

from ipywidgets import interact,fixed,IntSlider
import ipywidgets

def False_mols(idx):
    mol_idx = False_list[idx]
    mol = chain_mols[mol_idx]
    return(display(Draw.MolToImage(mol)))
    
interact(False_mols, idx=ipywidgets.IntSlider(min=0,max=len(False_list)-1, step=1));

f:id:magattaca:20190317085343g:plain

最初のいくつかはファーマコフォアを満たすとして認識されても良さそうですが、後半は芳香環を一つしかもたないなど、確かに要件を満たしていない、ということがわかります。

スクリーニングの実施

スクリーニング対象の再確認

関数の動作確認ができたので、いよいよ本題のスクリーニングに進みたいと思います。

スクリーニング対象の化合物群に関しては、以下の処理をこれまでに行いました。（だいぶ日があいたので再確認）

	処理	内容
①	指標による絞り込み	Rule of Five & フラグメントライクな化合物の除去
②	部分構造で絞り込み-1	ビフェニル構造をもつもの
③	部分構造で絞り込み-2	オルト位に置換基を持つもの

指標としては以下を用いています。

指標	分子量	LogP	水素結合供与体数	水素結合受容体数	回転可能結合数	極性表面積
Rule of 5	≦500	≦5	≦5	≦10
基準	>300	>3		>3	>3	>60

上記の結果として4220個の化合物まで絞り込み、biphenyl_library.sdfというSDFで保存しています。

スクリーニング本番

この化合物群にファーマコフォアサーチを実施します。

biphenyl_suppl = Chem.SDMolSupplier('./biphenyl_library.sdf')
biphenyl_mols = [x for x in biphenyl_suppl if x is not None]
biphenyl_PScores = []

for x in biphenyl_mols:
    num_embeds = PSearch3(x)
    biphenyl_PScores.append(num_embeds)
    x.SetIntProp('PScore', num_embeds)

ヒット（一つ以上ファーマコフォアを満たす構造がembedできた化合物）の数を確認します。

Hit_list = [i for i, x in enumerate(biphenyl_PScores) if x != 0]

print(len(Hit_list))
# 2463

合計2463個の化合物が残りました。約半数に絞り込めたようです。

上記の計算過程で、

RDKit ERROR: [23:37:56] DEBUG: Removed embedding due to c[00:09:28] UFFTYPER: Unrecognized charge state for atom: 1

といったエラーが出ましたが、4つだったので今回は深追いするのをやめようと思います。（困難から逃げるダメな大人）

embedの数が最も多かった化合物を確認してみます。

print([i for i, x in enumerate(biphenyl_PScores) if x == max(biphenyl_PScores)])
# [1351]

max_mol = biphenyl_mols[1351]
max_mol.GetProp('PScore')
#'99'
Draw.MolToImage(max_mol)

数は99で、構造は以下でした。

f:id:magattaca:20190317085803p:plain

・・・ほとんどsp2です。溶解度悪そう、、、

とりあえず結果をSDFで出力しておきます。

#保持するためのpropertyを確認
check_mol = biphenyl_mols[0]
prop_names = [p for p in check_mol.GetPropNames()]
print(prop_names)
# ['original_id', 'MW', 'MolLogP', 'NumHDonors', 'NumHAcceptors', 'NumRotatableBonds', 'TPSA', 'PScore']

writer = Chem.SDWriter('PSearched.sdf')

writer.SetProps(prop_names)
for mol in biphenyl_mols:
    writer.write(mol)
writer.close()

おしまい。

まとめ

今回は、前回作成したファーマコフォアをもちいて、スクリーニングを実施してみました。不斉点を持つ化合物に対して期待した挙動をしめさず、いくつか修正が必要となりましたが、最終的に約4000個の化合物から、2500個程度に絞ることができました。思っていたほど減りませんでしたが、ファーマコフォアの設定ポイント（フィーチャー）として、ビフェニル骨格の芳香環２つを選んでいることに由来しているかもしれません。すでに「ビフェニル構造を部分構造として持つ」という絞り込みを行なったリストに対してスクリーニングを実施しているので、「より数を絞り込む」と言う観点からはビフェニル以外のフィーチャーを選択すべきであったかもしれません。

*1:取り出した分子については以前の記事をご参照ください。

2019-03-11

ファーマコフォアを作ろうとする話

前回の記事でRDKitのファーマコフォアがどのような化学的構造を認識しているのか、その定義を眺めました。実際の化合物に適応するとどうなるのか、活性化合物のデータセット（共闘用シェアデータ）から取り出した分子（マクロサイクル型を除いたもの）で試したいと思います。

1. 活性化合物で検証
2. 共結晶構造をつかって検証
3. 共結晶構造をもとにファーマコフォアを作成
- 3-1. フィーチャーの選択と相対配置
- 3-2. 共結晶構造におけるフィーチャーの座標の取得
4. RDKitでファーマコフォアを設定
5. ファーマコフォアを基準に共結晶構造に重ね合わせ
- 5-1. PDB座標にアラインメント
まとめ

1. 活性化合物で検証

1-1. フィーチャーの探索

まずは活性化合物群を読み込みます。

import os
from rdkit import RDConfig
from rdkit.Chem import AllChem
from rdkit import Chem

chain_suppl = Chem.SDMolSupplier('./chain_compounds.sdf', removeHs=False)
chain_mols = [x for x in chain_suppl if x is not None]

FDefファイルを読み込んでフィーチャーファクトリーを作成します。

fdef = AllChem.BuildFeatureFactory(os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'))
print(fdef.GetFeatureFamilies())
#('Donor', 'Acceptor', 'NegIonizable', 'PosIonizable', 'ZnBinder', 'Aromatic', 'Hydrophobe', 'LumpedHydrophobe')

一つ目の化合物を取り出し、特性基（フィーチャー）の探索を行います。

test_m = chain_mols[0]
# そのままfeatureを検索
rawFeats = fdef.GetFeaturesForMol(test_m)
len(rawFeats)
#17

17個のフィーチャーが認識されました。各フィーチャーがどのフィーチャーファミリーのものか確認します。

feat_fam = []
for feat in rawFeats:
    feat_fam.append(feat.GetFamily())

Python標準ライブラリcollectionsを使って、上記のリストの各要素の個数を確認します。*1

import collections
c = collections.Counter(feat_fam)
print(c)
# Counter({'Hydrophobe': 5, 'Acceptor': 4, 'Aromatic': 3, 'Donor': 2, 'LumpedHydrophobe': 2, 'PosIonizable': 1})

用意されている８種類ファミリーのうち「'Hydrophobe', 'LumpedHydrophobe', 'ZnBinder'」の３つは特殊なように思えたので、より基本的な残りの５種類「'Donor', 'Acceptor', 'NegIonizable', 'PosIonizable', 'Aromatic'」を使うことにします。*2

keep = ('Donor','Acceptor','NegIonizable','PosIonizable','Aromatic')
featLists = []
for feat in rawFeats:
    # ファミリーを取得し、keepに含まれている時のみリストに追加
    if feat.GetFamily() in keep:
        featLists.append(feat)
len(featLists)
# 10

Hydrophobe（5個）とLumpedHydrophobe（2個）を除くと、基本的なフィーチャーは10個となりました。

1-2. フィーチャーの可視化

@iwatobipen先生の記事visualize-pharmacophore-in-rdkitを参考に、フィーチャーとして認識された構造を眺めます。

atom_ids_list = [] 
legend_list = []
for feat in featLists:
    # feat_family_list.append(feat.GetFamily())
    atom_ids_list.append(feat.GetAtomIds())
    # feat_type_list.append(feat.GetType())
    legend_list.append(feat.GetFamily() + ':' + feat.GetType())

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

Draw.MolsToGridImage([test_m]*10, molsPerRow=5, subImgSize=(200, 200), legends=legend_list, highlightAtomLists=atom_ids_list)

f:id:magattaca:20190310235324p:plain

水素結合ドナー(Donor)としてはアミンやアミドのNH、アクセプターとしてはアミド結合のカルボニル酸素、エーテル結合の酸素原子が認識されています。また、ピリジンのNは塩基性があるので、プラス電荷を帯びる塩基性グループ（PosIonizable:BasicGroup）としても良さそうですが、デフォルトのフィーチャーファクトリーでは、水素結合アクセプターとして認識されるようです（下段左端）。

1-3. SMARTSの表現を再確認

前回フィーチャーのSMARTSを眺めましたが、ピリジンのNを認識していると思われる部分を再度確認しましょう。 AcceptorファミリーのフィーチャーAcceptor.SingleAtomAcceptor、「[n;+0;!X3;!\$([n;H1](cc)cc)]」あたりでしょうか?

SMARTSviewerを使ってみます。*3

f:id:magattaca:20190310235402p:plain

「1」で囲まれた芳香族NH（ピロールのような構造？）ではない（「!:否定」）、芳香族Nと解釈できそうなので、ピリジンNを認識しそうです。

もっとかっこよく眺めたい！ということで、RDKitとPyMolの組み合わせを使ってみたいと思います。

1-4. PyMol x RDKit

@iwatobipen先生の記事Visualize pharmacophore with RDKitとGreg先生のノートShow_Ph4_Features_in_PyMOL.ipynbを参考にさせていただきました。

まず、ターミナルでpymolを xml-rpc（remote procedure call）サーバーモードで起動します。

# ターミナル
pymol -R

jupyter notebookに戻って以下を実行します。

from IPython import display
from rdkit.Chem import PyMol
IPythonConsole.ipython_useSVG=True

# viewerを作成
v = PyMol.MolViewer()

# viewerを一旦綺麗に掃除
v.DeleteAll()

# 先の構造を表示
v.ShowMol(test_m)

先に起動していたpymolを眺めると構造式が出ていました。グリグリ回してポーズを決めてから次のコマンドをjupyter notebookで実行するとノート上に静止画（PNG）が表示されました。

png=v.GetPNG()
display.display(png)

f:id:magattaca:20190310235618p:plain

PyMol綺麗・・・

PyMolが無事動いたので、フィーチャーの可視化に進みます。

1-5. PyMolでフィーチャーを可視化

フィーチャーを表示させるための関数を作成します。さっぱりわからないのでGreg先生のipynbからコピペさせていただきます。

RDKitのFeatures.ShowFeatsモジュールを使っているようですが、CGOなど全くわからない・・・

# コピペ
from rdkit.Chem.Features import ShowFeats

def showMolFeats(mol,factory,viewer,name="mol",showOnly=True):
    featLabel = f'{name}-feats'
    dirLabel = f'{name}-feats-dirs'
    if(showOnly):
        ShowFeats._globalSphereCGO = []
        ShowFeats._globalArrowCGO = []
    # workaround for what is either a bug in the way pymol handles CGOs
    # or a gap in my understanding of how it should work
    viewer.server.resetCGO(featLabel)
    viewer.server.resetCGO(dirLabel)
    viewer.server.sphere((0, 0, 0), .01, (1, 0, 1), featLabel)
    viewer.server.cylinder((0, 0, 0), (.01, .01, .01), .01, (1, 0, 1), dirLabel)
    ShowFeats.ShowMolFeats(mol,factory,v,name=name,featLabel=featLabel,dirLabel=dirLabel,
                           useFeatDirs=False,showOnly=showOnly)
    viewer.server.renderCGO(ShowFeats._globalSphereCGO,featLabel,1)
    viewer.server.renderCGO(ShowFeats._globalArrowCGO,dirLabel,1)

フィーチャーを眺めたい分子（molオブジェクト）はtest_m、フィーチャーファクトリーはfdef、ビューワーはv、という名前にしているので、これらを引数にして上の関数を使います。

showMolFeats(test_m,fdef,v)

pymolに移動して回転してから、jupyter notebookに表示します。

png=v.GetPNG()
display.display(png)

f:id:magattaca:20190310235812p:plain

格好いい！！

認識されたフィーチャーが球として表示されています。

先にMolsToGridImageで並べて見ていた時との違いは、黄土色の球の部位です。

表示されないようにFeatureFamilyから除いた「Hydrophobe」や「LumpedHydrophobe」のようです。

2. 共結晶構造をつかって検証

2-1. 共結晶構造のリガンドでフィーチャーを確認

かなりいい感じでフィーチャーを眺められるようになってきたので、PLIPを検証した記事で眺めた共結晶構造（PDB id: 5NIX）のリガンド（ID: 8YQ）でも、同じことをしてみようと思います。

まず、PDBからリガンドの構造をSDFでとってきて、Chain Aに含まれるリガンドの座標のみ残して保存しました（8YQ_noH.sdf）。

suppl_8YQ = Chem.SDMolSupplier('./8YQ_noH.sdf')
mol_8YQ = [x for x in suppl_8YQ if x is not None][0]

showMolFeats(mol_8YQ,fdef,v)
png=v.GetPNG()
display.display(png)

f:id:magattaca:20190310235941p:plain

なかなかごつい感じになりました。ほとんどの原子がフィーチャーに認識されていて、単なるBall&Stick表示のようにも見えてしまいます。

このフィーチャーのなかから、PLIPで探した「複数の共結晶構造に共通する相互作用」に該当する「フィーチャー」を取り出して組み合わせればファーマコフォアのモデルができるはず！

2-2. 共結晶構造の相互作用をPyMolで眺める

まずは、取り出した相互作用形式を再度確認します。

「PDB id:5NIX(Ligand:8YQ)」の場合・・・

	残基番号	アミノ酸	距離	相互作用		残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)					8YQ(chain C/D)	54C	ILE	3.89	Hydrophobic
	56B	TYR	3.82	Hydrophobic		56C	TYR	3.42/3.62	Hydrophobic
	66B	GLN	3.72	Hydrophobic
	115A/B	MET	3.70/3.80	Hydrophobic		115C/D	MET	3.85/3.82	Hydrophobic
	121A/B	ALA	3.69/3.75	Hydrophobic		121C/D	ALA	3.96/3.57	Hydrophobic
	123A	TYR	3.80/3.52/3.85	Hydrophobic		123D	TYR	3.72	Hydrophobic
	124A	LYS	5.03/3.18	π-cation/Salt Bridges		124D	LYS	3.70	Water Bridges/Salt Bridges

このうちChain A/B の残基（左半分）についてpymol上でハイライトして眺めて見ます。

pymolの設定などに関しては或る化みす途さんのブログ話題のXofluzaの変異をPymolで見てみるを参考にさせてきただきました。

「5nix.pdb」を表示させた後、

set cartoon_transparency, 0.8
リガンドを中心に持ってくる（Command+左クリック）
水を消す（object panelのallでH(hide)からwatersを選択）
見たい残基を選択（sequenceを表示「右下Sをクリック」して順番に選択）
残基をLine表示（object panelの(sele)でS(show)からas→wire→lines）

とします。スナップショットを保存するには、右上の「Draw/ray」というボタンから出力したいサイズを選んでDraw→save image to fileとしてやれば良いそうです。以下のようになりました。

f:id:magattaca:20190311000053p:plain

2-3. ファーマコフォア作成に向けた偽アトムの設定

ファーマコフォアを作るには、フィーチャー間の距離や角度といった情報をつかって定義するようです。フィーチャーが原子一つに帰着できるなら良いですが、芳香環みたいな複数の原子からなるグループ（官能基）の場合は、どの点を基準にして距離計測を行うかで結果が変わりそうです。

先に描画した図では芳香環（Aromatic フィーチャー）の真ん中に球が表示されていました。同じような点（偽原子）を作って見たいと思います。

PyMolWikiを参考にPseudoatomというのを作って見ました。

右下緑文字のSelectingをクリックしてAtomsに変更
芳香環の各原子を選択
pseudoatom arom1, sele

と打ち込みました。これで真ん中の芳香環の中に、新しい点(arom1という名前のobjectにしました)ができました

f:id:magattaca:20190311000159p:plain

同様にして他の環にもpseudoatomを作っていきました。

2-4. 相互作用を順番に

~ Tyr / Gln ~

だいたい準備ができた感じがするので一つ一つの相互作用部位を見ていきます。

	残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)	56B	TYR	3.82	Hydrophobic
	66B	GLN	3.72	Hydrophobic

f:id:magattaca:20190311000228p:plain

距離（点線）はWizardのMeasurementを選択してから、atomをクリックしていくと表示されました。 Chain Bの残基56 Tyrと66 GLNは芳香環と疎水性相互作用をしているようです。

~ Ala / Met ~

	残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)	115A/B	MET	3.70/3.80	Hydrophobic
	121A/B	ALA	3.69/3.75	Hydrophobic

f:id:magattaca:20190311000248p:plain

Chain Aの残基115 METとChain Bの残基121 ALA、Chain Bの残基115 METとChain Aの残基121 ALAがそれぞれ組み合わせとなって別々の芳香環と疎水性相互作用をしているようです。

~ Tyr / Lys ~

	残基番号	アミノ酸	距離	相互作用
8YQ(chain A/B)	123A	TYR	3.80/3.52/3.85	Hydrophobic
	124A	LYS	5.03/3.18	π-cation/Salt Bridges

f:id:magattaca:20190311000259p:plain

Chain Aの残基123 TYR（黄色）は芳香環と疎水性相互作用、Chain Aの残基124 LYS（オレンジ色）は芳香環との相互作用に加えて、末端カルボキシ基とも相互作用が見られるようです。

3. 共結晶構造をもとにファーマコフォアを作成

3-1. フィーチャーの選択と相対配置

タンパク質残基との相互作用を眺めることで、リガンド側で大事そうな部位がわかってきました。

SAR news No.19の記事を参考に、aromatic ２つとacceptor 1つの３点のフィーチャーをファーマコフォアとして選んで見ました。

（大事な点が３つあります！っていうと格好いいから３にした）

f:id:magattaca:20190311000349p:plain

距離をPymol上で確認して見ます。

Mouse Modeを3-Button Editingに変更
２点を選ぶと距離、３点では角度、４点で２面角が表示される

f:id:magattaca:20190311000432p:plain

結果をまとめると以下のようです。

f:id:magattaca:20190311000501p:plain

便宜的にAcceptorはカルボキシ基の炭素原子との距離を測りました。酸素原子を使った場合の距離はそれぞれ、5.1Åが[5.0, 5.6] 、10.4Åが[9.8, 11.2]となります。

着目するフィーチャー３点と、その位置関係が決まりました。このファーマコフォアを定義するため、フィーチャーの座標（位置情報）を確認します。

3-2. 共結晶構造におけるフィーチャーの座標の取得

pymolでshowfeatを使う際に、--writeFeatsのオプションをつかうことで、フィーチャーの座標を出力できるようになるそうです。@iwatobipen先生の記事Visualize pharmacophore with RDKitのコードをまたしてもそのままコピペさせていただきました。

(Greg先生の先の関数「showMolFeats」のどこかに--writeFeatsを書き込めばよかったのかもしれませんが私の能力ではわかりませんでした・・・)

import subprocess
import pprint
showfeatpath = os.path.join(RDConfig.RDCodeDir, 'Chem/Features/ShowFeats.py')

v = PyMol.MolViewer()
v.DeleteAll()
process = subprocess.Popen(['python', showfeatpath, '--writeFeats','./8YQ_noH.sdf'], stdout=subprocess.PIPE)
stdout = process.communicate()[0]

res = stdout.decode('utf-8').replace('\t', ' ').split('\n')
pprint.pprint(res)

下の図のように、座標が出力されます。（一部抜粋。全フィーチャーが並んでいる。）

f:id:magattaca:20190311000940p:plain

ここから必要なフィーチャーの座標を抜き出します。フィーチャーを間違えないようにAtom_Idを確認しておきます。

Draw.MolToImage(mol_8YQ, includeAtomNumbers = True)

f:id:magattaca:20190311001255p:plain

Atomid と Familyから該当の座標を取り出してきます。Draw.MolToImageで表示させたAtomNumberと、ShowFeatのAtom Idは何故か１ずつ番号がずれていました。理由はよくわかりません・・・とりあえず先に進みます。

芳香環２つの座標は以下・・・

'Aromatic -7.625 10.407 -21.653 1.0 # 1 10 11 36 35 34'
'Aromatic -1.851 11.051 -24.019 1.0 # 6 8 37 38 39 19'

アクセプター２点（カルボキシル基酸素原子）の間が、ちょうど'NegIonizable'として認識されているのでこちらの座標を使いたいと思います。

'NegIonizable -0.100 13.586 -28.137 1.0 # 45 28 20'
'Acceptor 0.862 14.067 -27.897 1.0 # 20', 'Acceptor -1.171 13.517 -28.363 1.0 # 28'

フィーチャーの座標位置とその位置関係が決まりました。いよいよRDKitを使ってファーマコフォアを設定したいと思います。

4. RDKitでファーマコフォアを設定

4-1. ファーマコフォア設定方法の流れ

SAR News No.19の吉森氏による記事「Chemoinformatics Toolkits を用いた創薬システム開発におけるラピッドプロトタイピング」での流れを確認します。

座標を使って特性基を定義
距離情報を使ってファーマコフォアを設定
ファーマコフォアサーチ
1. 対象のフィーチャーをそもそも含んでいるか？
2. フィーチャー間の距離が条件を満たすか？
3. 距離情報を拘束条件にして3D構造を発生させる。

以上の流れを辿って見ます。テストとして、記事冒頭に用共闘用シェアデータから取り出した分子を用います。

4-2. ファーマコフォアの設定を実践

フィーチャーの定義にはChemicalFeaturesモジュールを使うようです。

ChemicalFeatures.FreeChemicalFeature(Feature Family, Feature Type, 位置) とすることで、フィーチャーを定義できます（Typeは省略可能）。位置の指定にはGeometryのPoint3Dをもちいます。

from rdkit.Chem import ChemicalFeatures
from rdkit.Chem.Pharm3D import Pharmacophore
from rdkit import Geometry

feat_1=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-7.625, 10.407, -21.653))
feat_2=ChemicalFeatures.FreeChemicalFeature('Aromatic',Geometry.Point3D(-1.851, 11.051, -24.019))
feat_3=ChemicalFeatures.FreeChemicalFeature('Acceptor',Geometry.Point3D(-0.100, 13.586, -28.137)) 
feats = [feat_1,feat_2,feat_3]

Aromatic ２つとAcceptor１つが定義できました。

続いてPharm3D.Pharmacophoreモジュールを使って距離情報を設定します。先に定義したフィーチャーのリストでファーマコフォアを設定した後、フィーチャー間の距離のとりうる範囲を下限値（setLowerBound）、上限値（setUpperBound）という形で与えます。

SAR Newsの記事では距離の許容範囲の値として、下限値に0.5（d_lower）、上限値に1.5（d_upper）が用いられていました。

pcophore = Pharmacophore.Pharmacophore(feats) # ファーマコフォアの設定
d_upper = 1.5  # 特性基間の距離の許容範囲(上限値)
d_lower = 0.5 # 特性基間の距離の許容範囲(下限値)

# feat_1とfeat_2の距離 6.3Åを基準に、下限と上限を設定
pcophore.setLowerBound(0,1, 6.3-d_lower)
pcophore.setUpperBound(0,1, 6.3+d_upper)

# Acceptor(feat_3)の代表点の選び方によってfeat_2との距離は[5.0~5.6]の値をとる
pcophore.setLowerBound(1,2, 5.0-d_lower)
pcophore.setUpperBound(1,2, 5.6+d_upper)

# 同様にfeat1とfeat_3は[9.8~11.2]の値をとる
pcophore.setLowerBound(0,2, 9.8-d_lower)
pcophore.setUpperBound(0,2, 11.2+d_upper)

print(pcophore)

f:id:magattaca:20190311001409p:plain

以上でファーマコフォアの設定が終わりました。printで出力したPharmacophoreを見てみると、距離行列となっており、対角線右上三角に上限値、左下三角に下限値の情報を持つようです。

4-3. ファーマコフォアサーチ ~step 1. フィーチャーの有無~

それではテスト分子を使ってファーマコフォアサーチを行います。３次元構造を扱うので水素原子を付加しておきます。

test_mH = AllChem.AddHs(test_m)
AllChem.Compute2DCoords(test_mH)

まずはフィーチャーをそもそも持っているか？

Parm3D.EmbedLibモジュールのMatchPharmacophoreToMolを使います。

EmbedLib.MatchPharmacophoreToMol(molオブジェクト,フィーチャーファクトリー, ファーマコフォア)と、することで対象のmolオブジェクトにファーマコフォアのマッピングを行い、可能なマッピングのリストを作成します。

戻り値は2要素のタプルで、

全てのフィーチャーが見つけられたか否かのブール値
フィーチャーの並びのリスト

となっています。

テスト分子をtest_mH、フィーチャーファクトリーをfdef、ファーマコフォアをpcophoreとして実行します。

from rdkit.Chem.Pharm3D import EmbedLib

match, mList = EmbedLib.MatchPharmacophoreToMol(test_mH, fdef, pcophore)
print(match)
# True

どうやら全てのフィーチャーのマッチングはOKだったみたいです。確認のためリストからフィーチャーの情報を取り出して見ます。

print(len(mList))
#3
print(len(mList[2]))
#4
print(mList[2][0].GetFamily())
#Acceptor
print(mList[2][0].GetType())
#SingleAtomAcceptor
print(mList[2][2].GetAtomIds())
#(15,)

心配なのでAcceptorについて描画して見ます。

highlight_list = [mList[2][x].GetAtomIds() for x in range(len(mList[2]))]
Draw.MolsToGridImage([test_m]*4, molsPerRow=4, subImgSize=(200, 200), highlightAtomLists=highlight_list)

f:id:magattaca:20190311001756p:plain

うまく機能していそうです。

4-4. ファーマコフォアサーチ ~step 2. 距離条件~

ついで距離条件の検証を行います。

まず、rdDistGeomモジュールのGetMoleculeBoundsMatrixを使って分子の距離行列を取得します。

Parm3D.EmbedLibモジュールのGetAllPharmacophoreMatchesを使って、先に取得した分子の距離行列が、ファーマコフォアで定義されている距離の条件を満たすか否かを判定します。

from rdkit.Chem import rdDistGeom
# 距離行列の取得
bounds = rdDistGeom.GetMoleculeBoundsMatrix(test_mH)
#ファーマコフォアのマッチング 
pList =EmbedLib.GetAllPharmacophoreMatches(mList,bounds,pcophore) 
print(len(pList))
# 4

ファーマコフォア（３つのフィーチャー）の条件を満たす組み合わせは４つあったようです。一つ目の組み合わせのフィーチャーを構成する原子のIDを取得して見ます。

AtomIds_list = []

for i in range(len(pList[0])):
    p = pList[0][i]
    print(p.GetFamily(), ':', p.GetAtomIds())
    AtomIds_list.append(p.GetAtomIds())

f:id:magattaca:20190311002002p:plain

Draw.MolsToGridImage([test_m]*3, subImgSize=(200, 200), highlightAtomLists=AtomIds_list)

f:id:magattaca:20190311002048p:plain

他の組み合わせでは以下になりました。

f:id:magattaca:20190311002120p:plain

想定していたのと異なり、一番遠くの芳香環を使った組み合わせとなっていました。
ところで、このマッチングは距離を条件につかっているのですが、検索対象としているmolオブジェクトは３次元構造を生成させていません（pymolで眺めた時もぺちゃんこだった）。

SAR Newsの記事の該当部分を引用すると

一般的なファーマコフォアサーチにおいては、分子の3D構造生成後、ファーマコフォアのサーチを行うが、RDKitにおいては、事前に分子の3D構造を生成させるのではなく、ファーマコフォアを満たす3D構造が生成できるか否かを判定基準として、ファーマコフォアのサーチを行うことができる。

とのことだそうです。それでは距離条件をみたすような３次元構造が本当に作れるのかどうか、検証したいと思います。

4-5. ファーマコフォアサーチ ~step 3. 3D構造の発生~

距離情報を拘束条件にして3D構造を発生させます。まずはマッチしたフィーチャーの原子IDのリストを作成します。コピペ・・・

num_match = len(pList)
phMatches = []
for i in range(num_match): 
    num_feature = len(pList[i])
    phMatch = []
    for j in range(num_feature):
        phMatch.append(pList[i][j].GetAtomIds())
    phMatches.append(phMatch)

距離情報を拘束条件にした3D構造の発生にはPharm3D.EmbedLibのEmbedPharmacophoreを使います。

引数がたくさんありますが、対象のmolオブジェクト(mol)に対して、ファーマコフォアのフィーチャーを構成する原子のid(atomMatch)をつかって、ファーマコフォア(pcophore)の条件を満たすように3次元構造を生成(embedding)することができるようです。生成させる構造の最大数をcountで設定します。（silentはよくわかりません・・・）

戻り値は３要素のタプルで、

ファーマコフォアに合うように調整された拘束条件の距離行列
３次元構造（コンフォマー１つを有する分子）のリスト
３次元構造の生成に失敗した回数

となっています。

生成させたのち、ETKDG法を使って３次元構造の最適化を実施します。

confs = []

for phMatch in phMatches:
    bm,embeds,nFail =EmbedLib.EmbedPharmacophore(test_mH, phMatch,
                                                 pcophore,count=20,
                                                 silent=1)
    print("Generated embeds: ",len(embeds))
    print("Failed attempts: ",nFail)
    for embed in embeds:
        AllChem.EmbedMolecule(embed, AllChem.ETKDG())
        confs.append(embed)

f:id:magattaca:20190311002232p:plain

フィーチャーの組み合わせによって、構造生成の成否回数が異なっているようです。

SAR newsの記事ではUniversal Force Field法を使って構造最適化を行うと記載してあったのですが、以下のようなエラーメッセージが出てしまいました。

RuntimeError                              Traceback (most recent call last)
<ipython-input-254-6895836e8c95> in <module>()
      8     print("Failed attempts: ",nFail)
      9     for embed in embeds:
---> 10         AllChem.UFFOptimizeMolecule(embed)
     11         confs.append(embed)

RuntimeError: Invariant Violation
    bad direction in linearSearch
    Violation occurred on line 234 in file Code/Numerics/Optimizer/BFGSOpt.h
    Failed Expression: status >= 0
    RDKIT: 2018.09.1
    BOOST: 1_65_1

よくわからない・・・。飛ばそう・・・

念のためETKDGで最適化した３次元構造を一つ眺めて見ます。

from rdkit.Chem.Draw import IPythonConsole
import py3Dmol
IPythonConsole.drawMol3D(confs[0])

f:id:magattaca:20190311002429p:plain

いい感じにできていそうです。

5. ファーマコフォアを基準に共結晶構造に重ね合わせ

活性化合物に関して、ファーマコフォアの基準を満たす３次元構造の取得までが確認できました。この3D構造は、共結晶構造のリガンド構造と類似しているはず・・・なので、共結晶構造に重ね合わせてタンパク質と一緒に描いて見たいと思います。

5-1. PDB座標にアラインメント

ファーマコフォア作成に用いたリガンドの座標はPDBの共結晶構造を元にしているので、この座標をテンプレートとしてアライメントをとります。

Numerics.rdAlignmentモジュールのGetAlignmentTransformを使って、参照分子の座標(refPoints)と目的の分子の座標(probePoints)とで、RMSDが最小となる最適なアラインメントを計算します。

rdkit.Numerics.rdAlignment.GetAlignmentTransform(refPoints, probePoints)の戻り値として、SSD値（Sum of Squared Difference）と、4x4変換行列（transform matrix）のアレイを要素として持つ、2-タプルが得られます。

得た変換をAllChem.TransformMolを使って、適応することで、アラインメントをとった座標を持つmolオブジェクトを得ます。*4

# 生成した構造のフィーチャーを取得
match, mList_0 = EmbedLib.MatchPharmacophoreToMol(confs[0], fdef, pcophore)
bounds_0 = rdDistGeom.GetMoleculeBoundsMatrix(confs[0])
pList_0 =EmbedLib.GetAllPharmacophoreMatches(mList_0,bounds_0,pcophore) 

# ファーマコフォアで設定したフィーチャーのリストを参照(ref)として使用する
refFeats = feats

# pList[0]で認識されているフィーチャーを
#ファーマコフォアのフィーチャーの定義の順番に並べ替えたリストにする
probFeats = [pList_0[0][1], pList_0[0][0], pList_0[0][2]]

# ref、prob、それぞれのフィーチャーの座標を取得し、リストにする
probPts = [list(x.GetPos()) for x in probFeats]
refPts = [list(x.GetPos()) for x in refFeats]

# ２つの座標をつかって変換行列を取得
ssd,tform = rdAlignment.GetAlignmentTransform(refPts, probPts)
# 変換行列を使って、生成した構造の座標を変換
AllChem.TransformMol(confs[0], tform)

#変換したmolオブジェクトをmolファイルとして出力
Chem.MolToMolFile(confs[0], 'conf.mol')

出力したmolファイルと、元にした共結晶構造 5NIX.pdb を同時にpymolで表示させました。

f:id:magattaca:20190311002959p:plain

もともとのリガンドの色がシアンで、重ね合わせたリガンドが緑色です。なかなかいい感じで重なっているようにみえます。今回ファーマコフォアとして指定しなかった末端の芳香環についてはズレが大きいように見える一方で、ファーマコフォアマッチングした部分は重なりがよく見えます。うまくいった！！！　・・・のか？？？

まとめ

以上、今回はファーマコフォアの作成とテスト分子を使った動作検証を行って見ました。ファーマコフォアの基準を満たす３次元構造の取得までが確認できましたので、この取得ができるか否かまでを、ファーマコフォアを満たすことができるか否かという判定基準として活用していけば多分いい感じでスクリーニングになるはず。。。

あまりにもわからないことが多すぎて先生方の記事からひたすらコピペを繰り返しましたが、それでも正しく使えているのか自信がまったくありません。特にBounds Matrixが理解できなかった・・・（3次元に関する条件にしては行列の各要素が3個でもないし・・・）。

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

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

*1:参考: PythonのCounterでリストの各要素の出現個数をカウント

*2:参考: RDKitブログ記事Using Feature Maps

*3:citatition:

SMARTS viewer https://smartsview.zbh.uni-hamburg.de
1. Schomburg, H.-C. Ehrlich, K. Stierand, M.Rarey; From Structure Diagrams to Visual Chemical Patterns, J. Chem. Inf. Model., 2010, 50 (9), pp 1529-1535
  
  *4:参考: RDKitのメーリングリスト