ケモインフォマティクス(Chemoinformatics)とは、主に化学に関連するデータをコンピュータを用いて解析し、様々な課題を解くために用いる方法論のことです。ケモインフォマティクスという言葉は1990年代終わりから2000年代はじめに定義され、現在では製薬業界や薬学系のアカデミアにおいて、薬剤の効果と化合物特性の関連性を解析したり、大量の化合物情報の可視化、化合物の類似性に基づいたクラスタリングなど、多岐にわたるプロセスで利用されています。
近年、ディープラーニング(Deep Learning)の創薬応用が模索されていますが、活性や物性を予測するQSAR（構造活性相関）だけではなく、新規デザイン提案や合成経路提案といった従来のケモインフォマティクスではあまり扱わなかった領域への応用研究も盛んに行われています。
そもそもどのような化合物を作るべきか？、またそれをどのように合成するか？を考えるプロセスは背景知識と想像力が求められる領域であり、従来は人以外が担うのは難しい領域であると認識されていましたが、このような領域に対してもAIと呼ばれるものの進出がここ数年(2017-2019)で急速に進みました。
ケモインフォマティクスはすでに様々な場面で利用されていますが関連情報はあまりありませんでした。この理由として考えられることはいくつかありますが、オープンソースのツールキットが存在しなかった、公開データベースが存在しなかった、という二点が最も大きな原因なのは間違いありません。しかし、RDKitというオープンソースのケモインフォマティクスツールキットとChEMBLという公開データベースの登場により、この点は解消しました。
近年では、ケモインフォマティクスもバイオインフォマティクス同様、ウェブで検索すれば多くの情報がすぐに得られ、自己学習することは十分可能となっていますが、最初の一歩を踏み出すためのまとまった情報として、「ケモインフォマティクスに関する基礎を学び、それらを応用できるようになるコンテンツ」を用意することにしました。近年のAI創薬ブームを考慮して後半の章では「AI創薬」の文脈で用いられるディープラーニングを利用した化合物の活性予測と化合物提案の章を収録していますので、一通り学習することで、最近のトレンドについていけるようになっているはずです。
- warning
-
ここは@iwatobipenがRDKitについて熱く語るサブセクションです。下書き段階での「申し上げる」とか「踏まえて」等の言い回しをそのまま生かし、自称が「拙者」で「ござる」調の＠iwatobipen風文体です。
拙者、本書の一部を執筆する＠iwatobipenと申す。ここではRDKitに関して熱く語ってみようと思うでござる。
RDKitのRDはなんなのか？実はRational Discoveryの略称であり現在のオープンソースの前身となるフレームワークは2000年に開発されたでござる。随分と古いでござるね。その後、2006年にコードがオープンソースになりsourceforgeから公開されたでござる。PythonのケモインフォマティクスツールキットはRDKit以外にOpenBabelもあるぞと思われる読者もござろう。OpenBabelは2005年に最初のリリースがなされておる。いずれも、もう10年以上の歴史があるツールキットでござる。拙者がこの辺りに興味を持ち始めた2012年ころはどちらかというとOpenBabelの方がメジャーだったように記憶しておる。当時、日本語の記事はほぼ皆無であり、拙者は本書共著者であり業界のパイオニアでもある@fmkz___殿のケモインフォクックブックなどを参考にRDKitのコードを書いて試行錯誤していたでござるよ。なお、ケモインフォ関連のヒストリを追いたい御仁はこちらの記事を一読されるとよかろう。
おっと話が横道に逸れてしまった。本題に戻ろう。
開発者のGreg Landrum氏いわく
RDKitはケモインフォマティクスにおけるSwiss Army Knifeであり、様々な機能ピースの集合体である
これはまさに的を得た表現でござる。公式ドキュメントを見ればわかるでござろうが、既に色々な機能が用意されておるのだ。 化合物情報の読み込み、書き込みに始まり、構造の描画、３次元構造配座発生、Rグループ分解、記述子、フィンガープリント計算、ファーマコフォア算出などなど、挙げればきりがないほどの機能が実装されておる。解析から可視化まで幅広い範囲をカバーできるのだ。 さらにContributerらがRDKitを利用して開発したツール群がその熱い想いとともにContribフォルダーに詰められておるのだ。どうじゃ使ってみたくならんか？。拙者はもう書きながらも早くRDKitに触りたくなってきたでござる。
|
Note
|
@iwatobipenももちろんContributerの1人で、Fastclusterという大量の化合物ライブラリを高速にクラスタリングするコードを提供しています。(by @fmkz___) |
RDKitは開発やユーザーコミュニティの活動も活発で、どんどん機能追加がされておる。世界中の有能な研究者が全体で盛り上げ開発していくスタイルはオープンソースの強みであり、魅力であろう。もしチャンスがあれば毎年開催されるRDKit User Group Meetingへの参加を検討するのもよかろう。Face2Faceでユーザー同士議論ができるのは何事にも代え難いものがあるでござる。 また、先ほど拙者が使い始めた当時は日本語の情報ほぼ皆無であったと申したが、近年は非常に良質な日本語記事もたくさん増えておる。下記に何個か例を挙げたでござる。Qiitaにも多くの記事が掲載されているでござるよ。
また、有志によるRDKit-users-jpも立ち上がっておる。英語での質問がちょっと、、、と思われる御仁はこちらに質問を投げかけるとよかろう。また、最新版のRDKitのリポジトリには日本語のドキュメントもマージされておる。こちらも参考になるであろう。 本書ではRDKitの一部の機能しか使わん。それでも非常に多くのことができると感じていただけるはずじゃ。興味持ちはじめの一歩を踏み出したら後はどんどん自分の興味、意欲のままに足を進めていけばよかろう。何かわからないことがあれば上記のコミュニティに問いかけ、本書のリポジトリへIssueとして投稿してみるのもよかろう。 さあそれでは始めよう！
次のような方々を読者として想定しています。
-
医学薬学系の大学院生及び薬学系のデータ解析を行いたいポスドク
-
製薬企業の薬理研究者で自分のデータを自分で解析したい人
-
創薬化学者でケモインフォマティクスの必要性を感じている方や謎の力により突然アサインされた方
-
ケモインフォマティクスを学んでみようと考えているバイオインフォマティシャン
-
AI創薬に興味があるがなにからはじめたらいいかわからない人
本書で使用したプログラミングコードは全てMishima.sykのpy4chemoinformaticsリポジトリのnotebooksディレクトリに置いてありますので利用してください。またそれぞれの章の最初の
にその章のJupyter notebookへのリンクを張っていますので適宜参照してください。
2章のインストールを行うとgitコマンドが使えるようになりますので、以下のコマンドでpdfを含む本書の全てのデータがダウンロードできます
$ git clone https://github.com/Mishima-syk/py4chemoinformatics.gitもともとはBioに対してChemoと語感を合わせて登場してきたように記憶しているが、Journal of Cheminformaticsの創刊により一時期Chemに大きく離されていました。
最近のGoogle trendによるとどちらでもいいようですが個人的にはRhymeを重視したほうが良いと思うので本書ではChemoの方を使うことにします。
本書を執筆するにあたり、バグフィックスや改善のための助言をしてくれた以下の方々に感謝いたします。
@antiplastics, @bonohu, @ReLuTropy, @ski_nanko, @torusengoku, @yamasaKit_, keisuke-yanagisawa, Masaru Nagaku
ここから先は(Nujabes - reflection eternalを聴きながら書きました by @fmkz___ 2019/03/03)
まず、本書を書くきっかけとなった@bonohuに感謝したいと思います。@bonohuのDr. Bonoの生命科学データ解析の出版後のMishima.sykのミーティングで「Bono本のChemoinformatics版あったらいいよね」という話がどこからともなくでた際に、「書けばええんちゃう、むしろなんで書かんの？」と言ってくれたことが本書を執筆するきっかけであることは間違いありません。また@souyakuchanの創薬 Advent Calendar 2018も執筆のいい刺激になりました。というより、ここで章立てしなかったら具体的に動き出さなかったと思います。
また、忘れてはいけないのはy-samaの存在です。Mishima.sykを初期から盛り上げてきたy-samaは2019/01/06に永眠しました。彼はデータサイエンティストを目指す人のpython環境構築 2016やDruglikenessについてのよもやま話といった素晴らしいエントリを残しました。彼が存命であればきっと3人で執筆していたし、内容ももっと充実していたことでしょう。この出来事も我々に執筆しようという強い動機を与えました。
最後にMishima.sykに参加して美味しいワインやビールを飲みなから毎度熱い議論を交わしていただいた参加者の方々にも感謝します。いくつかのコンテンツはMishima.sykでの発表をもとにしており、みなさんのフィードバックをもとに加筆訂正してあります。
もし、本書を読んで、ケモインフォマティクスって面白いなと感じたり、創薬やってみたいなと感じる方がいたら、是非Mishima.sykに参加してみてください。きっと楽しいと思います。今後の創薬研究では所属を超えてお互いにプッシュしあって自身のスキルを高めていくことが重要になるでしょう。というより、既にそういう社会になっているのだと思います。本書が皆さんの楽しい研究生活を送る役に立てば幸いです。
やりたいことをやって生きてきて 私自身は自分の人生に後悔はありません 人生は楽しんだもの勝ち 皆さんも嫌なことは嫌だと言って自分の喜びを最大限に追い求めて人生を満喫した方が楽しいと思いますよ 皆様の人生に幸多い事を願っています
This document is copyright © 2019 by @fmkz___ and @iwatobipen
This document is Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International Public License.
本書に必要な環境構築を行います。
Anacondaは機械学習を行うための準備を楽にするためのパッケージです。また、後ほど説明するRDKitを簡単にインストールできます。
- なぜAnacondaを利用するのか？
-
プログラミング言語Pythonは比較的多くの標準ライブラリが用意されていますが、ケモインフォマティクス用のライブラリは自分でインストールする必要があります。この作業は慣れれば大した問題ではないですが、初学者にとっては面倒くさいでしょう。この手間を軽減するのが理由です。
- Pythonには大きく2.x系のバージョンと3.x系のバージョンとがありますが？
-
2.x系のサポートは2020年1月1日をもって終了するため、新しく学ぶ方は2.x系を使う必要はありません。
では早速Anacondaをインストールしましょう。公式サイトにアクセスして、ご自身の環境にあったPython3系インストーラーをダウンロードします。OSがLinux/MacであればGUI/CUIのインストーラーを選べるのでPython3.7の64-Bit Command Line Installerをダウンロードします。
$ bash ~/Downloads/Anaconda3-2019.03-MacOSX-x86_64.sh #インストーラー名は適宜変更してくださいEnterを押します
Welcome to Anaconda3 2018.12
In order to continue the installation process, please review the license
agreement.
Please, press ENTER to continue
>>>Enterを押し続けてYes,Noでyesを入力します
Do you accept the license terms? [yes|no]
[no] >>>インストール先を聞かれるのですがデフォルトで問題ないのでリターンを押します
Anaconda3 will now be installed into this location:
/Users/kzfm/anaconda3
- Press ENTER to confirm the location
- Press CTRL-C to abort the installation
- Or specify a different location belowインストール後にVSCodeをさらにインストールするかを聞かれるのでNoを押します。
Thank you for installing Anaconda3!
===========================================================================
Anaconda is partnered with Microsoft! Microsoft VSCode is a streamlined
code editor with support for development operations like debugging, task
running and version control.
To install Visual Studio Code, you will need:
- Internet connectivity
Visual Studio Code License: https://code.visualstudio.com/license
Do you wish to proceed with the installation of Microsoft VSCode? [yes|no]
>>> Please answer 'yes' or 'no':
>>>Anacondaのインストールが完了すると、コマンドプロンプトまたはターミナルから’conda’コマンドが使えるようなにります。
コマンドの-nの後ろは"py4chemoinformatics"としていますが皆さんの好きな名前でも構いません。 仮想環境構築後、本章以降で利用するパッケージをインストールします。
$ conda create -n py4chemoinformatics
$ source activate py4chemoinformatics # Mac/Linux
$ activate py4chemoinformatics # Windows
# install packages
$ conda install -c conda-forge rdkit
$ conda install -c conda-forge seaborn
$ conda install -c conda-forge ggplot
$ conda install -c conda-forge gitもし仮想環境でjupyter notebookが動かない場合はもう一度インストールする必要があります。
$ conda install notebook ipykernelRDKitはケモインフォマティクスの分野で最近よく用いられるツールキットの一つです。オープンソースソフトウェア(OSS)と呼ばれるものの一つで、無償で利用することが可能です。詳しくははじめにを参照してください。
統計データの視覚化のためのパッケージの一つです。
グラフ描画パッケージの一つで一貫性のある文法で合理的に描けることが特徴です。もともとはRという統計解析言語のために開発されましたが、yhatという会社によりPythonに移植されました。
バージョン管理システムです。本書ではGitについては説明しませんのでもしGitについて全然知らないという方はサルでもわかるGit入門でも読みましょう。
「はじめに」でも説明しましたが、以下のコマンドでpdfを含む全てのデータがダウンロードされますので必要に応じてダウンロードしてください。
$ git clone https://github.com/Mishima-syk/py4chemoinformatics.git- なぜ仮想環境を作るのでしょうか
-
いくつかのシステムでは様々な機能を提供するために内部的にPythonを利用しているため、特定のパッケージのためにPythonのバージョンを変更してしまうと問題が起こることがあります。仮想環境はこのような問題を解決します。もし、パッケージが異なるライブラリのバージョンを要求しても仮想的なPython環境を準備して試行錯誤できます。不要になれば仮想環境を簡単に削除でき、もとの環境にトラブルを持ちこむこともありません。このように、ひとつのシステム内にそれぞれ個別の開発環境を作成できるようにすることで開発時によく起こるライブラリの依存問題やPythonのバージョンの違いに悩まされることがなくなります。
本書では本書用に一つだけ仮想環境を用意しますが、実際はいくつもの仮想環境をつくって開発することが多いです。そのため、よく利用するcondaのサブコマンドを挙げておきます。
$ conda install <package name> # install package
$ conda create -n 仮想環境の名前 python=バージョン # 仮想環境の作成。
$ conda info -e # 作った仮想環境一覧の表示
$ conda remove -n 仮想環境の名前 # 仮想環境の削除
$ conda activate 仮想環境の名前 # 仮想環境を使う(mac/linux)
$ activate 仮想環境の名前 # 仮想環境を使う(Windows)
$ conda deactivate # 仮想環境から出る(source deactivateは廃止されました)
$ conda list # 今使っている仮想環境にインストールされているライブラリの一覧を表示この章ではPythonに触れたことのない読者のために効率的に勉強するためのサイトや本などを紹介します。 もしこれ以降の章でわからないことなどがあったら、この章のサイトや本を参考に学んでみてください。
- Pythonスタートブック増補改訂版
-
プログラミング自体が初心者であればこの本が良いでしょう。
- みんなのPython 第4版
-
JavascriptやJavaなどのなにかプログラミングを少しかじっていて、これからPythonを覚えたいのであればこちらの本をおすすめします。
- Python Boot Camp(初心者向けPythonチュートリアル)
-
一般社団法人PyCon JPが開催している初心者向けPythonチュートリアルイベントです。全国各地で行われているので近くで開催される場合には参加するとよいでしょう
- その他ローカルコミュニティなど
-
あちこちで入門者向けからガチのヒト向けまでの勉強会やコミュニティなどもあるので、そういうのに参加してモチベーションを高めるのもよい方法です。
- udemy/python
-
オンライン学習サービスを利用するのも効果的な手段のひとつですが、筆者は試したことがないのでわかりません。 周りの評判を聞いてみても良いでしょう。YouTubeを探すのもありです。
- py4chemoinformaticsのissues
-
py4chemoinformaticsのissuesに質問していただければお答えします。わかりにくい場合だったら修正しますので、よりよくなってみんなハッピー。
- Qiita
-
Qiitaで探せば大抵答えが見つかるはずです。
- stackoverflow
-
それでも答えが見つからなかったらsofで探すか質問しましょう
- Mishima.syk
-
本書を書いている人たちが集まるコミュニティです。特に話題をPythonに限定していませんが、Pythonを使ったネタが多めです。かなりガチですが、初心者対応も万全でハンズオンに定評があります。質問されれば大体答えられます。
Jupyter notebookを利用すると、コードを書いて結果を確認するということがとても簡単にできるようになります。
Jupyter notebookはWebブラウザーベースのツールで、コードだけではなくリッチテキスト、数式、なども同時にノートブックに埋め込めます。また結果を非常に綺麗な図として可視化することも容易にできます。つまり、化学構造やグラフも描画できるため、ケモインフォマティクスのためのプラットフォームとして使いやすいです。さらに、プログラミングの生産性を上げるような、ブラウザ上でコードを書くとシンタックスハイライトや、インデント挿入を自動で行ってくれたりという便利な機能もついているので、特に初学者は積極的に使うべきでしょう。
ケモインフォマティクスに限らず、インフォマティクスを学ぶにあたり、機械学習は外せません。本書でもある程度の機械学習の知識があることを前提に進めていきます。Pythonで機械学習をするにはScikit-learnというライブラリを利用するのが定番であり、本書でも特に説明せずに利用していきますが、初学者のために参考となる書籍などをすすめておきます。
- Pythonではじめる機械学習 ―scikit-learnで学ぶ特徴量エンジニアリングと機械学習の基礎
-
Pythonで機械学習をやるための基礎を学べます。数学的な表現があまりないので読みやすいです。
- sklearn-tutorial
-
y-samaによるsklearnのチュートリアルハンズオンのjupyter notebookです。
この章ではケモインフォマティクスでよく使うデータベースを紹介します。
ChEMBLはEBIのChEMBLチームにより維持管理されている医薬品及び開発化合物の結合データ、薬物動態、薬理活性を収録したデータベースです。データは主にメディシナルケミストリ関連のジャーナルから手動で抽出されており、大体3,4ヶ月に一度データの更新があります。
メディシナルケミストリ関連のジャーナルからデータを収集しているため、QSARに関連する情報や背景知識を論文そのものに求めることが可能であり、創薬研究をする際には有用です。
|
Note
|
ChEMBLはもともとはStARliteという商用データベースでした。詳しくは慶応大学池田先生のChEMBLに関する資料を参照してください。 |
PubChemはNCBIにより維持管理されている低分子化合物とその生物学的活性データを収録している公開リポジトリです。5000万件以上の化合物情報と、100万件を超えるアッセイデータを含みそのデータ量の多さが特徴とも言えます。もうひとつの特徴はデータをアカデミアからの化合物登録やアッセイ結果の登録により成長することであり、ここが先のChEMBLとの大きな違いです。
特にPubChemは初期スクリーニングのデータが多いため、そのようなデータに対しなんらかのマイニングや分析を行いたい場合は有用だと考えられます。
- どちらを使うべき?
-
QSARをやりたい場合にはやはりChEMBLのデータを利用することが多いです。IC50のようなデータが得られていることが多いですし、モデルの解釈に元論文をあたることができるというのが大きな理由です。
|
Note
|
ChEMBLはユーザーインターフェースを刷新中で現在beta版のテストを行っていますが、いずれこちらに置き換わると思うので新バージョンのインターフェースでの検索方法を紹介します。 |
まずはChEMBLにアクセスし、画面上部のCheck out our New Interface (Beta). というリンクをクリックして新しいインターフェース画面に移行します。
ChEMBLのデータは主に4つのカテゴリに分かれていて、一意なIDが振られており相互に関連付けされています。それぞれのカテゴリについて簡単に説明すると
- Targets
-
ターゲット分子についてその分子を対象としてアッセイされた論文に関してまとめられており、どういったジャーナルに投稿されているかや、どの年に投稿されたのかといった情報がまとめられています。また、アッセイに関しても同様にまとめられています。
- Compounds
-
化合物に関する基本的な物理量（分子量など）のほか、Rule of 5を満たしているかといった分子の特性情報や、臨床情報などの創薬関連情報のほか、ChEMBLでの関連アッセイ、関連論文のサマリがまとめられています。
- Assays
-
アッセイに関する情報と元論文との関連付けがされているほか、アッセイに供された化合物データへのリンクが貼られています。
- Documents
-
論文のタイトル、ジャーナル名、アブストラクトの他に関連論文データへのリンクと、その論文中で行われたアッセイへのリンクと使われた化合物データへのリンクが貼られています。
ある創薬ターゲット分子がどのくらい研究開発されているかを知るために、それをターゲットとしてどのくらいの化合物が合成されたのか？さらに骨格のバリエーションはどのくらい存在するのかを調べたい場合がよくあります。ChEMBLを利用するとターゲット名で探索して関連化合物をダウンロードすることができます。
ここでは抗がん剤のターゲットとして知られているTopoisomerase2を検索します。画面上部のフォームにtopoisomeraseと入力して検索するとスクリーンショットのように表示されるはずです。
サジェスト機能による絞り込みでいくつか候補をリスト表示してくるのでTOP2Bを選んでください。画面をスクロールするとAssociated Compoundsセクションがありますのでグラフのタイトル(Associated Compounds for Target CHEMBL3396)をクリックすると関連化合物一覧画面が開きます。
259化合物存在することがわかります。スクロールすると全体をみることができます。画面右のアイコンをクリックするとそれぞれCSV(カンマ区切りテキスト),TSV(タブ区切りテキスト),SDF(5章で説明しています)の形式でダウンロードできます。
QSARモデルを作る場合、アッセイの活性値と対応する化合物の構造情報が必要です。ChEMBLの場合アッセイのページからダウンロードすることでQSARモデル作成のためのデータを得ることができます。
大体次のような手順を辿ることがおおいです。
-
論文データを検索してからそれに関連付けられているアッセイデータを辿る
-
ターゲットを検索してそれに紐付いているアッセイデータからQSARに使えそうなものを選ぶ
ここでは後者のターゲットから検索してQSARモデルに使えそうなアッセイデータを探します。心毒性関連ターゲットとしてよく知られているhERGのQSARモデルを作りたいという状況を想定しています。
検索フォームにhERGと入力して、Search hERG for all in Assaysを選びます。361件ヒットしました。
モデル構築のためのデータが欲しいのでデータ数が多い順に並べ替えます。ヘッダーのCompoundsをクリックして降順に並べ替えます。
論文由来で最もアッセイ数の多いCHEMBL829152を選んでクリックしてアッセイページを開きます。Activity chartの円グラフをクリックすると詳細画面が開くのでSelect allで全選択してTSV形式でダウンロードします。
- NOTE
ダウンロードしたファイルをエディタで開くと^@C^@h^@E^@M^@B^@L^@と文字化けすることがあります。これはutf-16-leでエンコードしているためです(こうしないとExcelで問題が発生するようです)。
viの場合':e ++enc=utf16le’と打てばきちんと表示されるようになります。
ZINCは購入可能な試薬をコレクションしたデータベースです。現在のバージョンは15で約7億5000万の構造が収載されています。 もともとがドッキングシミュレーションでの利用を想定して開発されているため、三次元化したデータをダウンロードすることも可能です。ZINCのデータでバーチャルスクリーニング(6章で説明します)を行い、ヒットした化合物を購入し実際のアッセイに供するというのが主な使い方だと思います。
データのダウンロード方法は上部のTranchesタブをクリックすると次の画面に縦軸にLogP横軸に分子量の大きさで分類されそれぞれの区画にいくつの化合物が収載されているかの表が表示されます。
ここから必要なデータセットを選んでダウンロードボタンを押すと、実際にデータセットのURLが列挙されたテキストファイルが得られますのでそれぞれにアクセスしてデータをダウンロードします。
統合TVは生命科学分野の有用なデータベースやツールの使い方を動画で紹介するサイトで、ライフサイエンス統合データベースセンター(DBCLS)により管理、運用されています。その名の通りバイオインフォマティクス関連の動画が多いですが、ケモインフォマティクスを紹介した動画もいくつかありますので参考にしてみてください。文献・辞書・プログラミングのカテゴリも役に立つはずです。
|
Note
|
生命科学データベース・ウェブツール 図解と動画で使い方がわかる! 研究がはかどる定番18選という書籍も出版されています。 |
- NOTE
-
これ以外にもケモインフォマティクスに有用なデータベースがあればお知らせください。IssueやPRでも受け付けてます。
この章ではRDKitを使って分子の読み込みの基本を覚えます。
Simplified molecular input line entry system(SMILES)とは化学構造を文字列で表現するための表記方法です。 詳しくはSMILES Tutorialで説明されていますが、例えばc1ccccc1は6つの芳香族炭素が最初と最後をつないでループになっている構造、つまりベンゼンを表現していることになります。
SMILESで分子を表現することがわかったので、SMILESを読み込んで分子を描画させてみましょう。まずはRDKitのライブラリからChemクラスを読み込みます。二行目はJupyter Notebook上で構造を描画するための設定です。
from rdkit import Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import DrawRDKitにはSMILES文字列を読み込むためにMolFromSmilesというメソッドが用意されていますので、これを使い分子を読み込みます。
mol = Chem.MolFromSmiles("c1ccccc1")続いて構造を描画しますが、単純にmolを評価するだけで構造が表示されます。
mol図のように構造が表示されているはずです。
上のように原子を線でつなぎ構造を表現する方法（構造式）と、SMILES表記はどちらも同じものを表現しています。構造式は人が見てわかりやすいですが、SMILESはASCII文字列で表現されるのでより少ないデータ量で表現できるというメリットがあります。
|
Note
|
文字列で表現できるということは、文字列生成アルゴリズムを応用することで新規な化学構造を生成することも可能ということです。この内容に関しては12章で詳しく説明します。 |
複数の化合物を一つのファイルに格納する方法にはいくつかありますが、sdというファイル形式を利用するのが一般的です。sdファイルということで、ファイルの拡張子は.sdfとなることが多いです。
MDL社で開発された分子表現のためのフォーマットにMOL形式というものがあります。このMOL形式を拡張したものがSDファイルです。具体的にはMOL形式で表現されたものを$$$$という行で区切ることにより、複数の分子を取り扱えるようにしてあります。
MOL形式は分子の三次元座標を格納することができ二次元だけでなく立体構造を表現できる点はSMILESとの大きな違いです。
4章を参考にChEMBLのトポイソメラーゼII阻害試験(CHEMBL669726)の構造データをsdfファイル形式でダウンロードします。
- NOTE
具体的な手順はリンクのページを開いて、検索フォームにCHEMBL669726を入力すると検索結果が表示されるので、Compoundsタブをクリックします。その後、全選択してSDFでダウンロードするとgzip圧縮されたsdfがダウンロードされるので、gunzipコマンドまたは適当な解凍ソフトで解凍してください。それをch05_compounds.sdfという名前で保存します。
RDKitでsdfファイルを読み込むにはSDMolSupplierというメソッドを利用します。複数の化合物を取り扱うことになるのでmolではなくmolsという変数に格納していることに注意してください。どういう変数を使うかの決まりはありませんが、見てわかりやすい変数名をつけることで余計なミスを減らすことは心がけるとよいでしょう。
mols = Chem.SDMolSupplier("ch05_compounds.sdf")何件の分子が読み込まれたのか確認します。数を数えるにはlenを使います。
len(mols)34件でした。
創薬の化合物最適化ブロジェクトで、分子の形を変更しないで化合物の特性を変えたいということがあります。このような場合、芳香環を形成する炭素、窒素、硫黄、酸素などの原子種を入れ替えることでより良い特性の化合物が得られることがありますがこのようにヘテロ原子(水素以外の原子)を入れ替えるアプローチをヘテロシャッフリングといいます。
ヘテロシャッフリングを行うことで、活性を維持したまま物性を変化させて動態を良くする、活性そのものを向上させる、特許クレームの回避といった効果が期待できます。
少しの構造の違いが選択性や薬物動態が影響を与える有名な例として、Pfizer社のSildenafilとGSK社のVardinafilが挙げられます。
二つの構造を比較すると中心の環構造部分の窒素原子の並びが異なっているだけで極めて似ています。両分子は同じ標的蛋白質を阻害しますが、その活性や薬物動態は異なります。
上記の画像を生成するコードを示します。単にDraw.MolsToGridImageを適用するのではなく Core構造をベースにアライメントしていることとDraw.MolToGridImageのオプションにlegendsを与え、分子名を表示していることに注意してください。
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import rdDepictor
from rdkit.Chem import rdFMCS
from rdkit.Chem import TemplateAlign
IPythonConsole.ipython_useSVG = True
rdDepictor.SetPreferCoordGen(True)
sildenafil = Chem.MolFromSmiles('CCCC1=NN(C)C2=C1NC(=NC2=O)C1=C(OCC)C=CC(=C1)S(=O)(=O)N1CCN(C)CC1')
vardenafil = Chem.MolFromSmiles('CCCC1=NC(C)=C2N1NC(=NC2=O)C1=C(OCC)C=CC(=C1)S(=O)(=O)N1CCN(CC)CC1')
rdDepictor.Compute2DCoords(sildenafil)
rdDepictor.Compute2DCoords(vardenafil)
res = rdFMCS.FindMCS([sildenafil, vardenafil], completeRingsOnly=True, atomCompare=rdFMCS.AtomCompare.CompareAny)
MCS = Chem.MolFromSmarts(res.smartsString)
rdDepictor.Compute2DCoords(MCS)
TemplateAlign.AlignMolToTemplate2D(sildenafil, MCS)
TemplateAlign.AlignMolToTemplate2D(vardenafil, MCS)
Draw.MolsToGridImage([sildenafil, vardenafil], legends=['sildenafil', 'vardenafil'])ヘテロシャッフルした分子を生成するためにHeteroShuffleというクラスを定義します。オブジェクトの生成にはシャッフルしたい分子と変換したい部分構造（Core）を与えます。クラス内のコードではまず、分子をCoreで切断し、Coreとそれ以外に分けます。CoreのAromatic原子で、置換基がついてない原子のみが置換候補になります。シャッフル後のCoreとCore以外のパーツを再結合するための反応オブジェクトを生成するメソッドがmake_connectorです。このメソッドで作られた反応オブジェクトを利用してre_construct_molで分子を再構築しています。
考えられる原子の組み合わせを構築するために、itertools.productに、候補原子（C, S, N, O）の原子番号と、環を構成する原子数target_atomic_numsを与えます。その後に分子として生成できないものは排除するのでここでは考えられる全部の組み合わせを出します。
import copy
import itertools
from rdkit import Chem
from rdkit.Chem import AllChem
class HeteroShuffle():
def __init__(self, mol, query):
self.mol = mol
self.query = query
self.subs = Chem.ReplaceCore(self.mol, self.query)
self.core = Chem.ReplaceSidechains(self.mol, self.query)
self.target_atomic_nums = [6, 7, 8, 16]
def make_connectors(self):
n = len(Chem.MolToSmiles(self.subs).split('.'))
map_no = n+1
self.rxn_dict = {}
for i in range(n):
self.rxn_dict[i+1] = AllChem.ReactionFromSmarts('[{0}*][*:{1}].[{0}*][*:{2}]>>[*:{1}][*:{2}]'.format(i+1, map_no, map_no+1))
return self.rxn_dict
def re_construct_mol(self, core):
'''
re construct mols from given substructures and core
'''
keys = self.rxn_dict.keys()
ps = [[core]]
for key in keys:
ps = self.rxn_dict[key].RunReactants([ps[0][0], self.subs])
mol = ps[0][0]
try:
smi = Chem.MolToSmiles(mol)
mol = Chem.MolFromSmiles(smi)
Chem.SanitizeMol(mol)
return mol
except:
return None
def get_target_atoms(self):
'''
get target atoms for replace
target atoms means atoms which don't have anyatom(*) in neighbors
'''
atoms = []
for atom in self.core.GetAromaticAtoms():
neighbors = [a.GetSymbol() for a in atom.GetNeighbors()]
if '*' not in neighbors and atom.GetSymbol() != '*':
atoms.append(atom)
print(len(atoms))
return atoms
def generate_mols(self):
atoms = self.get_target_atoms()
idxs = [atom.GetIdx() for atom in atoms]
combinations = itertools.product(self.target_atomic_nums, repeat=len(idxs))
smiles_set = set()
self.make_connectors()
for combination in combinations:
target = copy.deepcopy(self.core)
for i, idx in enumerate(idxs):
target.GetAtomWithIdx(idx).SetAtomicNum(combination[i])
smi = Chem.MolToSmiles(target)
target = Chem.MolFromSmiles(smi)
if target is not None:
n_attachment = len([atom for atom in target.GetAtoms() if atom.GetAtomicNum() == 0])
n_aromatic_atoms = len(list(target.GetAromaticAtoms()))
if target.GetNumAtoms() - n_attachment == n_aromatic_atoms:
try:
mol = self.re_construct_mol(target)
if check_mol(mol):
smiles_set.add(Chem.MolToSmiles(mol))
except:
pass
mols = [Chem.MolFromSmiles(smi) for smi in smiles_set]
return mols上のコードで使われているcheck_molという関数はc1coooo1のような６員環の構造もAromaticだと判定されてしまうのでそれを避けるために使っています。O, Sが許容されるのは５員環のヘテロ芳香環のみにしました。
def check_mol(mol):
arom_atoms = mol.GetAromaticAtoms()
symbols = [atom.GetSymbol() for atom in arom_atoms if not atom.IsInRingSize(5)]
if not symbols:
return True
elif 'O' in symbols or 'S' in symbols:
return False
else:
return True実際に使ってみます。
# Gefitinib
mol1 = Chem.MolFromSmiles('COC1=C(C=C2C(=C1)N=CN=C2NC3=CC(=C(C=C3)F)Cl)OCCCN4CCOCC4')
core1 = Chem.MolFromSmiles('c1ccc2c(c1)cncn2')
# Oxaprozin
mol2 = Chem.MolFromSmiles('OC(=O)CCC1=NC(=C(O1)C1=CC=CC=C1)C1=CC=CC=C1')
core2 = Chem.MolFromSmiles('c1cnco1')元の分子
ht = HeteroShuffle(mol1, core1)
res = ht.generate_mols()
print(len(res))
Draw.MolsToGridImage(res, molsPerRow=5)Gefitinibを入力とした場合の変換結果の一部です。芳香環を形成する原子が元の化合物から変化した分子が出力されています。 また、Coreで指定したキナゾリン部分のみが変換されています。
ht = HeteroShuffle(mol2, core2)
res = ht.generate_mols()
print(len(res))
Draw.MolsToGridImage(res, molsPerRow=5)Oxaprozinを入力とした場合の変換結果です。こちらは中心に、オキサゾールと呼ばれる5員環構造を有してます。５員環を形成する芳香環にはチオフェン、フランなどのように窒素や酸素を含むものもあります。以下の例でもS、Oが5員環の構成原子に含まれている分子が出力されています。
どうでしょうか。二つの分子の例を示しました。一つ目の例、Gefitinibは、分子を構成する芳香環が、キナゾリンとベンゼンでした。キナゾリンは、ベンゼンとピリミジンという二つの６員環が縮環した構造です。６員環をベースに構成される芳香環を形成する原子の候補は炭素と窒素になります。（ピリリウムイオンなど電荷を持つものも考慮すれば酸素や硫黄も候補になりますが、通常このような構造をDrug Designで使うことは少ないので今回の説明からは外しています。複素環式化合物の説明） Oxaprozinはオキサゾールを有しています。５員環の芳香環を形成する原子の候補は炭素、窒素、硫黄、酸素が挙げられます。このような分子の場合の例として紹介しました。 いずれのケースでも上記のコードでヘテロ原子がシャッフルされたものが生成されています
|
Tip
|
2018.03.1 以降のRDKitではEnumerateHeterocyclesというメソッドが実装されているので上記のコードを書かずともヘテロシャッフルを容易に行うことができます。また、全部の芳香族原子を変換の対象にしたくないという場合があるかと思います。その時は、除外したい原子に_protectedという属性持たせておくと良いです。実際の例をコードで見てみましょう。 |
# EnumerateHeterocyclesはジェネレーターを返します。これは組み合わせが多くリストで返すとメモリを多く消費するためです。以下の例ではわかりやすくするためリストに変換しています。
from rdkit import Chem
from rdkit.Chem import EnumerateHeterocycles
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import AllChem
import copy
omeprazole = Chem.MolFromSmiles('CC1=CN=C(C(=C1OC)C)CS(=O)C2=NC3=C(N2)C=C(C=C3)OC')
enumerated_mols = EnumerateHeterocycles.EnumerateHeterocycles(omeprazole)
enumerated_mols = [m for m in enumerated_mols]
print(len(enumerated_mols))
# 384
Draw.MolsToGridImage(enumerated_mols[:6], molsPerRow=3)
# 今度は56縮環部分は保護してみます。
ringinfo = omeprazole.GetRingInfo()
ringinfo.AtomRings()
# ((1, 6, 5, 4, 3, 2), (14, 13, 17, 16, 15), (18, 19, 20, 21, 15, 16))
protected_omeprazole = copy.deepcopy(omeprazole)
for idx in ringinfo.AtomRings()[1]:
atom = protected_omeprazole.GetAtomWithIdx(idx)
atom.SetProp("_protected", "1")
for idx in ringinfo.AtomRings()[2]:
atom = protected_omeprazole.GetAtomWithIdx(idx)
atom.SetProp("_protected", "1")
enumerated_mols2 = EnumerateHeterocycles.EnumerateHeterocycles(protected_omeprazole)
enumerated_mols2 = [m for m in enumerated_mols2]
print(len(enumerated_mols2))
# 4
Draw.MolsToGridImage(enumerated_mols2)
J. Med. Chem. 2012, 55, 11, 5151-5164ではPIM-1キナーゼ阻害剤におけるNシャッフリングの効果をFragment Molecular Orbital法という量子化学的なアプローチを使って検証しています。さらにJ. Chem. Inf. Model. 2019, 59, 1, 149-158ではAsp–Arg塩橋とヘテロ環のスタッキングのメカニズムを量子化学計算により探っており、置換デザインの指標になりそうです。
また、バイオアベイラビリティ改善のためにヘテロシャッフリングを行った例としてはJ. Med. Chem. 2011, 54, 8, 3076-3080があります。
２つの化合物が似ているとはどういうことでしょうか？なんとなく形が似ている？という表現は科学的ではありません。ケモインフォマティクスでは類似度(一般的に0-100%値を取ります)や非類似度(距離)といった定量的な尺度により似ているかどうかを評価します。
ここでは主に２つの代表的な尺度を紹介します。
分子の全体的な特徴を数値で表現するものを記述子と呼びます。分子量や極性表面性(PSA)、分配係数(logP)などがあり、現在までに多くの記述子が提案されています。これらの記述子の類似性を評価することで２つの分子がどのくらい似ているかを表現することが可能です。また分子全体の特徴を1つの数字で表現しており局所的な特徴ではないということに注意してください。
|
Note
|
いくつかの記述子に関しては市販ソフトでないと計算できない場合があります。 |
もう一つがフィンガープリントです。フィンガープリントとは分子の部分構造を0,1のバイナリーで表現したもので部分構造の有無とビットのon(1),off(0)を対応させたものになり、部分構造の集合を表現することで分子の特徴を表現しています。フィンガープリントには固定長FPと可変長FPの二種類が存在し、古くはMACCS keyという固定長FP(予め部分構造とインデックスが決められているFP)が使われていましたが、現在ではECFP4(Morgan2)という可変長FPが利用されるのが普通です。
RDKitのフィンガープリントに関しては開発者のGregさんのスライドが詳しいので熟読してください。
今回はこのECFP4(Morgan2)を利用した類似性評価をしてみましょう。
SMILESは構造をASCII文字列で表現したものでフィンガープリントは部分構造の有無をバイナリで表現したものです。違いは前者は構造表現の一つであるのに対し、後者は特徴表現の一つだということです。 部分構造の有無だけしか表現していないため、部分構造間の関係性(どう位置関係でつながっているのか)といった情報が失われ、もとの構造に戻ることはありません。
テキストマイニングでよく用いられるBag-of-Wordsに対応するのでBag-of-Fragmentsと呼ぶ人もいます。
簡単な分子としてトルエンとクロロベンゼンの類似性を評価してみましょう。
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsoleSMILESで分子を読み込みます。
mol1 = Chem.MolFromSmiles("Cc1ccccc1")
mol2 = Chem.MolFromSmiles("Clc1ccccc1")一応目視で確認しておきます。
Draw.MolsToGridImage([mol1, mol2])ECFP4に相当する半径2のモルガンフィンガープリントを生成します。
fp1 = AllChem.GetMorganFingerprint(mol1, 2)
fp2 = AllChem.GetMorganFingerprint(mol2, 2)類似度の評価にはタニモト係数を使います。
DataStructs.TanimotoSimilarity(fp1, fp2)
# 0.5384615384615384類似度の補数(すなわち1 - 類似度)である距離を求めたい場合はreturnDistance=Trueとします。 後述のクラスタリングなどでは類似度よりも距離を用いることの方が多いです。
DataStructs.TanimotoSimilarity(fp1, fp2, returnDistance=True)
# 0.46153846153846156ここまでで化合物の類似性の評価方法について説明しました。この類似性の指標を用い多くの化合物の中から特定の化合物群を選び出すことをバーチャルスクリーニングと呼びます。
例えば薬になりそうな化合物が特許や論文で発表されたり、自社のアッセイ系で有望そうな化合物が見つかった場合、自社の化合物ライブラリデータベースや市販化合物のデータベースの中に類似の化合物で、より有望そうなものがあるかどうかを探したいことがあります。ここではノイラミニダーゼ阻害薬として知られるインフルエンザ治療薬であるイナビルの類似体が購入可能であるかをZINCを利用して調べます。
イナビルの分子量が約350,ALogPが-2.92だったので、ZINCの分子量350-375,LogPが-1未満の340万化合物の区画を選択しました。この区画は16のファイルに分かれていますが、最初の1セットだけダウンロードして使ってみます。
|
Note
|
データのダウンロード方法は4章で説明しています。 |
jupyter notebookでは!で始めるとShellコマンドを実行できます。以下はjupyter notebook上でwgetコマンドでZINCのデータセットをダウンロードする例です。
!wget http://files.docking.org/2D/EA/EAED.smi- NOTE
-
もし、Windowsなどでwgetコマンドがない場合はブラウザにURLを直接入力するとファイルのダウンロードができるので、ダウンロードしたファイルをjupyterを起動しているディレクトリに移動してください。
SmilesMolSupplier関数を使い、ファイルからまとめてSMILESを読み込みmolオブジェクトにします。5章で紹介したSDMolSupplierのSMILESバージョンと考えてください。
spl = Chem.rdmolfiles.SmilesMolSupplier("EAED.smi")
len(spl)
# 195493続いてイナビル(Laninamivir)との類似度を調べるための関数を用意します。
laninamivir = Chem.MolFromSmiles("CO[C@H]([C@H](O)CO)[C@@H]1OC(=C[C@H](NC(=N)N)[C@H]1NC(=O)C)C(=O)O")
laninamivir_fp = AllChem.GetMorganFingerprint(laninamivir, 2)
def calc_laninamivir_similarity(mol):
fp = AllChem.GetMorganFingerprint(mol, 2)
sim = DataStructs.TanimotoSimilarity(laninamivir_fp, fp)
return sim調べてみます。
similar_mols = []
for mol in spl:
sim = calc_laninamivir_similarity(mol)
if sim > 0.2:
similar_mols.append((mol, sim))結果を類似度の高い順に並べ替えて最初の10件だけ取り出します。
similar_mols.sort(key=lambda x: x[1], reverse=True)
mols = [l[0] for l in similar_mols[:10]]描画してみます。
Draw.MolsToGridImage(mols, molsPerRow=5)
類似度を確認すればわかりますが、今回調べた約20万件の化合物は最高でも23%の類似度の化合物しか見いだせませんでした。しかしZINCは7億5000万件のデータを収録してあるのでその中にはもっと似ている化合物はたくさんあるはずです。
例えば市販化合物を購入してライブラリを作る場合にはできるだけ多様性をもたせたいので、似ている化合物ばかりが偏らないように類似化合物同士をまとめ、その中の代表を選びます。このように化合物を構造の類似性でまとめたい場合、クラスタリングという手法を使います。
Novrtisの抗マラリアアッセイの5614件のヒット化合物をクラスタリングします。
クラスタリング用のライブラリをインポートし、データを読み込みます。
from rdkit.ML.Cluster import Butina
mols = Chem.SDMolSupplier("ch06_nov_hts.sdf")何らかの理由でRDKitで分子の読み込みができない場合、molオブジェクトではなくNoneが生成されます。このNoneをGetMorganFingerprintAsBitVect関数にわたすとエラーになるので、Noneを除きながらフィンガープリントを生成します。
fps = []
valid_mols = []
for mol in mols:
if mol is not None:
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
fps.append(fp)
valid_mols.append(mol)フィンガープリントから距離行列(下三角の距離行列)を生成します。ある化合物を参照化合物とし、複数の化合物の間の類似度や距離を計算する時はDataStructs.BulkTanimotoSimilarityを用いるとまとめて計算ができます。
distance_matrix = []
for i, fp in enumerate(fps):
distance = DataStructs.BulkTanimotoSimilarity(fps[i], fps[:i+1],
returnDistance=True)
distance_matrix.extend(distance)距離行列を用いて化合物をクラスタリングします。3番目の引数は距離の閾値です。この例では距離0.2つまり80%以上の類似度の化合物でクラスタリングしています。
clusters = Butina.ClusterData(distance_matrix, len(fps), 0.2, isDistData=True)クラスタ数を確認します。
len(clusters)
#2492最初のクラスタの構造を表示してみます
mols_ = [valid_mols[i] for i in clusters[0]]
Draw.MolsToGridImage(mols_, molsPerRow=5)
今回はRDKitに用意されているライブラリでクラスタリングを行いましたが、Scikit-learnでも幾つかの手法が利用できますし、実際にはこちらの方を使うことが多いです。
ここでは抗凝固薬として上市されているapixaban, rivaroxabanの類似性を評価します。
apx = Chem.MolFromSmiles("COc1ccc(cc1)n2nc(C(=O)N)c3CCN(C(=O)c23)c4ccc(cc4)N5CCCCC5=O")
rvx = Chem.MolFromSmiles("Clc1ccc(s1)C(=O)NC[C@H]2CN(C(=O)O2)c3ccc(cc3)N4CCOCC4=O")Draw.MolsToGridImage([apx, rvx], legends=["apixaban", "rivaroxaban"])
構造を見るとわかる通りなんとなく似ていますが、この２つの化合物は両方共FXaというセリンプロテアーゼの同じポケットに同じように結合し、蛋白質の働きを阻害することが知られています。
apx_fp = AllChem.GetMorganFingerprint(apx, 2, useFeatures=True)
rvx_fp = AllChem.GetMorganFingerprint(rvx, 2, useFeatures=True)
DataStructs.TanimotoSimilarity(apx_fp, rvx_fp)
# 0.4062540%くらいの類似度ということになりました。実はapixabanもrivaroxabanも複合体結晶構造が解かれており、PyMOLを利用して重ねて描画しました。
- NOTE
-
本書の内容を超えるので、PyMOLの使い方などについては説明しませんが、興味があればこの辺りを参考にしてください。
図からわかるように、apixaban,rivaroxabanは立体的に綺麗に重なっています。特にメトキシフェニルとクロロチオールはS1ポケットと呼ばれる部位に位置しており何らかの強い相互作用をしていると言われています。蛋白質のリガンド結合部位(ポケット)が明確になれば、メディシナルケミストが次にどの辺りを修飾するかの戦略が立てやすくなり、プロジェクトの成功確率や進捗速度も上がります。
このように、X線やクライオ電験などにより決められた蛋白質の形状に基づいて構造最適化を行うアプローチをStructure Based Drug Design(SBDD)と呼びます。また、ポケットがわかれば、物理的にポケットに結合するような化合物をスクリーニングすることができ、これを構造ベースのバーチャルスクリーニング(SBVS)と呼び、前の章でやったようなリガンドベースのバーチャルスクリーニング(LBVS)と区別することもあります。
本書のケモインフォマティクスの内容とは大きくかけ離れるのですが、FXa阻害剤の歴史を辿り、世代を経てどのような改良がなされていったのかを理解するのは分子設計をする上でかなり役に立ちます。また、S1ポケットの相互作用の解釈は目視や古典力学では非常に難しく、Fragment Molecular Orbital Methodのような量子化学計算によってのみ解釈できるため、今後の分子設計において量子化学計算は必須になることは間違いないと思います。
グラフとはノード（頂点）群とノード間の連結関係を示すエッジ（枝）群で構成されるデータのことを指します。化学構造はこのグラフで表現できます。つまり原子をノード、結合をエッジとしたグラフ構造で表せます。
通常、6章で紹介したようなフィンガープリントを使い分子同士の類似性を評価することが多いですが、グラフ構造を利用して類似性を評価する手法もあります。次に紹介するMCS（Maximum Common Substructure）は対象となる分子集合の共通部分構造のことを指します。共通部分構造が多いほとそれらの分子はより似ていると考えます。
最大共通部分構造Maximum Common Substructure(MCS)とは与えられた化学構造群において共通する最大の部分構造のことです。RDKitではMCS探索のためにrdFMCSというモジュールが用意されています。
今回はMCS探索のサンプルデータとしてrdkitに用意されているcdk2.sdfというファイルを利用します。RDConfig.RDDocsDirが、サンプルデータのディレクトリを表す変数で、そのディレクトリ以下のBooks/data/にcdk2.sdfというファイルが存在するので、os.path.joinメソッドでファイルパスを設定します。尚、os.path.joinはosのパスの違いを吸収するためのpythonの組み込みモジュールです。
import os
from rdkit import Chem
from rdkit.Chem import RDConfig
from rdkit.Chem import rdFMCS
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
filepath = os.path.join(RDConfig.RDDocsDir, 'Book', 'data', 'cdk2.sdf')
mols = [mol for mol in Chem.SDMolSupplier(filepath)]
# 構造を確認します
Draw.MolsToGridImage(mols[:7], molsPerRow=5)
読み込んだ分子を使ってMCSを取得します。RDKitではMCSの取得方法に複数のオプションが指定できます。 以下にそれぞれのオプションでの例を示します。
-
デフォルト
-
原子がなんであっても良い（構造とボンドの次数があっていれば良い）
-
結合次数がなんでも良い（例えば、ベンゼンとシクロヘキサンは同じMCSとなる）
result1 = rdFMCS.FindMCS(mols[:7])
mcs1 = Chem.MolFromSmarts(result1.smartsString)
mcs1
print(result1.smartsString)
#[#6]1:[#7]:[#6](:[#7]:[#6]2:[#6]:1:[#7]:[#6]:[#7]:2)-[#7]
result2 = rdFMCS.FindMCS(mols[:7], atomCompare=rdFMCS.AtomCompare.CompareAny)
mcs2 = Chem.MolFromSmarts(result2.smartsString)
mcs2
print(result2.smartsString)
#[#6]-,:[#6]-,:[#6]-[#6]-[#8,#7]-[#6]1:[#7]:[#6](:[#7]:[#6]2:[#6]:1:[#7]:[#6]:[#7]:2)-[#7]
result3 = rdFMCS.FindMCS(mols[:7], bondCompare=rdFMCS.BondCompare.CompareAny)
mcs3 = Chem.MolFromSmarts(result3.smartsString)
mcs3
print(result3.smartsString)
#[#6]1:[#7]:[#6](:[#7]:[#6]2:[#6]:1:[#7]:[#6]:[#7]:2)-[#7]
RDKitではMCSに基づく類似性を数値化するアルゴリズムのひとつにFraggle Similarityが実装されています。これを利用することでクラスタリングや、類似性に基づいた解析が行なえます。
from rdkit.Chem.Fraggle import FraggleSim
sim, match = FraggleSim.GetFraggleSimilarity(mols[0], mols[1])
print(sim, match)
#0.925764192139738 *C(C)C.*COc1nc(N)nc2[nH]cnc12
match_st = Chem.MolFromSmiles(match)
match_st
このようにFraggleSimilarityは類似性及びマッチした部分構造を返します。ECFPを利用した類似性よりもケミストの感覚に近いことが多いです。詳しくは参考リンクを参照してください。
参考リンク
創薬研究の構造最適化ステージにおいて、起点となる化合物（リード化合物）をどのように構造変換していくかは非常に重要ですが、ステージが進んだ場合どの構造変換が活性や物性に影響を及ぼしたかというレトロスペクティブな解析することもまた非常に大切です。
|
Tip
|
興味があればhttps://sar.pharm.or.jp/wp-content/uploads/2018/09/SARNews_19.pdfを読むとよいです。 |
Matched Molecular Pair(MMP)は、二つの分子のうち一部の部分構造だけが異なりそれ以外は同一な分子のペアのことです。例として、クロロベンゼンとフルオロベンゼンはCl基とF基のみが異なるのでMMPです。このようなペアの特性の変化を大量に解析することで、置換基変換のトレンドを掴むことができます。これをMatched Molecular Pair Analyisis（MMPA)と呼びます。大規模なデータでMMPAを行うことにより、置換基の変化がもたらす特性変化の普遍的なルールを抽出できます。このようなルールを理解していれば構造最適化を効率的に進められることになります。
ここではRDKitのContribに提供されているRDKit/Contrib/MMPA[mmpa]を使ってMMP解析を行います。
RDKitインストール先の下にあるContrib/mmpaに移動し、pythonスクリプトを順次実行します。
python rfrag.py <MMPAを実施したいFileの名前 >フラグメント化したデータの保存ファイル名
# 例えば
# python rfrag.py <data/sample.smi >data/sample_fragmented.txt
python indexing.py <先のコマンドでできたフラグメントのファイル >MMP_アウトプットファイル.CSV
# 例えば
# python indexing.py <data/sample_fragmented.txt >data/mmp.csv以上のコマンドを実行すると分子A, 分子B, 分子AのID, 分子BのID, 変換された構造のSMIRKS, 共通部分構造（context）のcsvファイルが生成されます。このデータに基づき活性や物性などを紐つけることでMMPAが行えます。
|
Note
|
SMIRKSは分子の変換をSMILESのように文字列表記によって表現する手法です。 |
MMPの拡張としてMatched Molecular Series(MMS)という手法も提案されています。MMPは分子のペアですが、このペアを共通構造をもった3つ以上の集団としてリスト化したものがMMSです。
実際にMMSを作ってみましょう。以下の例ではChEMBLよりFactor Xaのデータをダウンロードし、一例として使いました。MMSの実装に関してはNoel O’Boyle氏のRDKit UGMでのプレゼンテーションのコードを利用しています。
まず利用するライブラリの読み込みと、データの読み込みを行い、SaltRemoverを使い脱塩しています。
import sys
import os
import pandas as pd
from rdkit import Chem
from rdkit.Chem import rdMMPA
from rdkit.Chem import RDConfig
from rdkit.Chem import rdBase
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem import Draw
from rdkit.Chem import SaltRemover
# RDKit contrib/mmpaをPythonのSystem Pathに追加して次のパートで分子をフラグメント化するためのrfragをインポートできるようにします
mmpapath = os.path.join(RDConfig.RDContribDir, 'mmpa')
sys.path.append(mmpapath)
# データをPandasのDataFrameとして読み込み、SmilesをMolオブジェクトに変換、CHEMBLのIDをプロパティとして登録します。
df = pd.read_csv('Chembl_FXa.txt', sep='\t')
remover = SaltRemover.SaltRemover()
mols = []
for smi, chembl_id in zip(df.CANONICAL_SMILES, df.CMPD_CHEMBLID):
try:
mol = Chem.MolFromSmiles(smi)
mol.SetProp('CMPD_CHEMBLID', chembl_id)
mol = remover.StripMol(mol)
mols.append(mol)
except:
print(smi)続いてRDKit contribに登録されているmmpaのrfragをインポートして、分子をフラグメントに分割します。 以下のコードで用いているrfrag.fragment_molで呼ばれている関数の実態はrdMMPA.FragmentMolです。この関数は、与えたMolオブジェクトを回転可能ボンドにて切断したフラグメントを生成します。 出力の形式はfragment(SMILES, COMPOUND_ID)とした場合、SMILES, COMPOUND_ID, CORE, CHAINS というカンマ区切りのテキストのリストになります。COREというのは例えばA-B-Cという分子（ーが回転可能ボンドを意味します）を考えると、切断の仕方がA/B-C, A/B/C, A-B/Cという三つのケースがあります。二番目のケースの時にBがCoreにA, CがCHAINSになります。一番目と、三番目はどちらも (A,B-C), (A-B,C) というフラグメントはCHAINSにアサインされ、COREはNONEになります。
import rfrag
remover = SaltRemover.SaltRemover()
rfragdata = []
for smi, chembl_id in zip(df.CANONICAL_SMILES, df.CMPD_CHEMBLID):
try:
mol = Chem.MolFromSmiles(smi)
mol = remover.StripMol(mol)
smi = Chem.MolToSmiles(mol)
out = rfrag.fragment_mol(smi, chembl_id)
rfragdata.append(out)
except:
print(smi, chembl_id)MMSを作成する関数を定義します。コードはUGMの資料に記載されているものをほぼそのまま利用しますが、Jupyter上で全ての処理をおこなうため読み込み先をファイルからリストに変更しました。
以下MMSの作成プロセスの概要です。
-
各分子を一定のルール（回転可能結合で切断など）でカット
-
カットしたフラグメントがkeyの辞書を作成、同じキーを持つ分子のフラグメントを辞書のvalueに格納
上記の作業を繰り返すとことで共通のスキャフォールドを持つ分子をまとめられます。共通のスキャフォールドでまとまった分子は、スキャフォールド以外の置換基が異なる分子となります。
創薬において、前臨床試験の前のステージに構造最適化というステージがあり、そこでは化合物の主要骨格以外の部分をちょこまかと変換して薬にふさわしいバランスの取れたプロパティにします。
この主要骨格のことをスキャフォールドと呼びます。例えばこの特許ではRを除いた部分は固定されておりこの主要骨格をスキャフォールドと呼びます。
以下のコードの簡単な説明をします。
Fragというnamedtupleはid, scaffold, rgroupを属性に持ちます。Namedtupleを利用することで簡単なクラスと同じ働きをするオブジェクトを作ることができます。
getFrags関数は、フラグメント化してあるデータを元にFragクラスにデータを登録します。登録の際に、フラグメントの原子数を比較しどちらをScaffoldにするかを決めています。この関数を呼ぶことで各分子のid, scaffold, rgroupのペアを作ります。
次のgetSeries関数はFragオブジェクト内に格納されたscaffoldと与えられたデータを照合し、同じscaffoldを持つrgroupのデータをまとめます。これによってペアからシリーズにしています。Seriesクラスはscaffoldが文字列（SMILES)、rgroupsはリストで複数のフラグメントを持てるようになっています。
from collections import namedtuple
Frag = namedtuple( 'Frag', ['id', 'scaffold', 'rgroup'] )
class Series():
def __init__( self ):
self.rgroups = []
self.scaffold = ""
def getFrags(rfrags):
frags = []
for lines in rfrags:
for line in lines:
broken = line.rstrip().split(",")
if broken[2]: # single cut
continue
smiles = broken[-1].split(".")
mols = [Chem.MolFromSmiles(smi) for smi in smiles]
numAtoms = [mol.GetNumAtoms() for mol in mols]
if len(numAtoms) < 2:
continue
if numAtoms[0] > 5 and numAtoms[1] < 12:
frags.append(Frag(broken[1], smiles[0], smiles[1]))
if numAtoms[1] > 5 and numAtoms[0] < 12:
frags.append(Frag(broken[1], smiles[1], smiles[0]))
frags.sort(key=lambda x:(x.scaffold, x.rgroup))
return frags
def getSeries(frags):
oldfrag = Frag(None, None, None)
series = Series()
for frag in frags:
if frag.scaffold != oldfrag.scaffold:
if len(series.rgroups) >= 2:
series.scaffold = oldfrag.scaffold
yield series
series = Series()
series.rgroups.append((frag.rgroup, frag.id))
oldfrag = frag
if len(series.rgroups) >= 2:
series.scaffold = oldfrag.scaffold
yield seriesMMSを作る準備ができたので実行します。同じスキャフォールドに対して４つ以上置換基の変換があったデータのみを可視化します。
frags = getFrags(rfragdata)
series = getSeries(frags)
series =[i for i in series]
from IPython.display import display
for s in series[:50]:
mols = [Chem.MolFromSmiles(s.scaffold)]
ids = ['scaffold']
for r in s.rgroups:
rg = Chem.MolFromSmiles(r[0])
mols.append(rg)
ids.append(r[1])
if len(mols) > 5:
display(Draw.MolsToGridImage(mols, molsPerRow=5, legends=ids))
print("########")
スキャフォールドに対して5つの置換基のMMSが表示されました。
|
Note
|
このMMSを利用して活性予測を行うこともできます。 |
|
Warning
|
この内容は入門の内容を超えるので興味がなければ飛ばしてください |
MMPは変換前、変換後の情報をノード、変換ルールをエッジとするグラフ構造と考えることができます。Cytoscapeなどのネットワーク可視化ツールを利用するとこのグラフ構造を直感的に把握できます。
RDKitには先に紹介したMMPAの他にmmpdbという別プロジェクトがあります。 こちらはコマンドラインのツール群とデータベースシステムとして提供されているため、長期的な管理がしやすいという特徴があります。本セクションではこのmmpdbとCytoscapeを利用したMMPの可視化を紹介します。
Cytoscapeはオープンソースのネットワーク可視化ソフトで色々なシーンで広く使われています。化合物の構造表示用プラグインを使うことで構造のネットワークを表示することができます。
インストールは簡単でダウンロードサイトから対応するOSのインストーラをダウンロードして指示のとおりにインストールするだけです。
インストールが完了したらCytoscapeを起動して化合物構造描画用のChemviz2プラグインをインストールします。手順は簡単でApps→App Managerからchemviz2を選択してインストールします。
今回利用するデータは<Inhibition of recombinant GSK3-beta> J. Med. Chem. (2008) 51:2062-2077の151化合物です。MMPAを行うにはHTSのような探索データではなくて構造最適化のようにスキャフォールドが決まっているものを使うのが原則です。
コマンドの流れを載せておきます。SMILESのtextと活性や物性値のデータは別々にデータベースに登録する必要があります。
$ mmpdb fragment smiles.txt -o CHEMBL930273.fragments # fragmentation
$ mmpdb index CHEMBL930273.fragments -o CHEMBL930273.db # make db
$ mmpdb loadprops -p act.txt CHEMBL930273.db # load propertiesそのあとCytoscapeで読み込むためのgmlファイルを作成しますが、これは本書の範囲を超えるので割愛します。もし興味があるのであればコードを直接読んでもらうといいのですが 流れは以下のとおりです。
-
gmlファイルをCytoscapeで読み込む
-
Cytoscapeで属性を各パラメータにアサインして視覚的に理解しやすくする
-
ノードの大きさを物性値に対応
-
エッジの色を活性差に対応
-
chemviz2 pluginで構造を描画してノードに貼り付ける
-
MMPネットワークを見てみましょう。あまり活性差のないMMPが左上の方に固まっています。右下の方にはエッジが赤い（活性差が大きい）ものが観測されます。このような小さな置換基変化が大きな活性差を生むもMMPをActivity Cliffと呼びます。一般的にActivity Cliffは創薬プロジェクトにおいてブレークスルーとなることが多いため、このような活性変化を見逃さないことは大切です。
実際にどういう置換が行われたのかを確認すると、OH基がMeO基に置換されることで活性の消失が起こっています。
MMPだけではこのように単純に事実しかわからないので、もう少し深く考察するために類似体の複合体結晶構造を探してみました。するとPDBID:5OY4というGSK3βと類似化合物の複合体が見つかりました。
OH基をMeO基に置換するとポケットの壁にぶつかりそうですね。つまりこのActivity Cliffはリガンドと蛋白質の立体障害 により引き起こされたと考えられます。
沢山のデータがどのように分布しているのかを見るには適当な空間にマッピングするのが一般的です。特にケモインフォマティクスではケミカルスペースという言葉が使われます。
ケミカルスペースとは化合物を何らかの尺度でn次元の空間に配置したものを指します。一般に、2次元または3次元が使われることが多いです（人間の理解のため）。尺度つまり類似性に関しては色々な手法が提案されていますが、うまく化合物の特徴を表すような距離が定義されるように決められることが多いです。
今回は睡眠薬のターゲットとして知られているOrexin Receptorのアンタゴニストについて、どの製薬企業がどういった化合物を開発しているのかを視覚化してみます。データのダウンロード方法は4章を参照してください。今回は表の10個の論文のデータを利用しました。
今回知りたいことは主に以下の２つです。
-
似たような化合物を開発していた会社はあったのか？
-
Merckは似たような骨格ばかり最適化していたのか、それとも複数の骨格を最適化したのか？
Doc ID |
Journal |
Pharma |
CHEMBL3098111 |
Merck |
|
CHEMBL3867477 |
Merck |
|
CHEMBL2380240 |
Rottapharm |
|
CHEMBL3352684 |
Merck |
|
CHEMBL3769367 |
Merck |
|
CHEMBL3526050 |
Actelion |
|
CHEMBL3112474 |
Actelion |
|
CHEMBL3739366 |
Heptares |
|
CHEMBL3739395 |
Actelion |
|
CHEMBL3351489 |
Eisai |
描画ライブラリにはggplotを使います。主成分分析(PCA)を利用して、化合物が似ているものは近くになるように分布させて可視化します。まずは必要なライブラリをインポートします
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
import numpy as np
import pandas as pd
from ggplot import *
from sklearn.decomposition import PCA
import osダウンロードしたsdfを読み込んで、製薬企業とドキュメントIDの対応が取れるようにしてそれぞれの化合物についてフィンガープリントを構築します。もし不明な点があれば6章を確認してください。
oxrs = [("CHEMBL3098111", "Merck" ),("CHEMBL3867477", "Merck" ),
("CHEMBL2380240", "Rottapharm" ),("CHEMBL3352684", "Merck" ),
("CHEMBL3769367", "Merck" ),("CHEMBL3526050", "Actelion" ),
("CHEMBL3112474", "Actelion" ),("CHEMBL3739366", "Heptares" ),
("CHEMBL3739395", "Actelion" ), ("CHEMBL3351489", "Eisai" )]
fps = []
docs = []
companies = []
for cid, company in oxrs:
sdf_file = os.path.join("ch08", cid + ".sdf")
mols = Chem.SDMolSupplier(sdf_file)
for mol in mols:
if mol is not None:
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
docs.append(cid)
companies.append(company)
fps.append(arr)
fps = np.array(fps)
companies = np.array(companies)
docs = np.array(docs)フィンガープリントの情報を確認すると10の論文から293化合物のデータが得られていることがわかります。
fps.shape
# (293, 2048)これで主成分分析の準備完了です。主成分の数はn_componentsで指定できますが今回は二次元散布したいので2にします。
pca = PCA(n_components=2)
x = pca.fit_transform(fps)描画します。colorオプションを変えると、それぞれのラベルに応じた色分けがされるので、OMPANYとDOCIDの2つの属性を選んでみました。
d = pd.DataFrame(x)
d.columns = ["PCA1", "PCA2"]
d["DOCID"] = docs
d["COMPANY"] = companies
g = ggplot(aes(x="PCA1", y="PCA2", color="COMPANY"), data=d) + geom_point() + xlab("X") + ylab("Y")
g各製薬会社がどのような化合物を最適化したのかがわかるようになりました。ケミカルスペースの中心部に各社重なる領域があるので、Merck, Acterion, Eisai, Heptaressは似たような化合物を最適化していたと思われます。Acterionはうまく独自性のある方向(左下)に展開できたのか、展開できなくてレッドオーシャン気味の中心部に進出してきたのかは興味深いです。
またMerckは色々な骨格を最適化していたようです。同時に最適化したのか先行がこけてバックアップに走ったのかわかりませんが、多数の骨格の最適化が動いていたのは間違いないので、ターゲットとしての魅力が高かったということでしょう。実際SUVOREXANTは上市されましたしね。
本章では論文データを利用しましたが、実際の現場でこのような解析をする場合には論文データは使いません。なぜなら企業が論文化するときはそのプロジェクトが終わったこと（成功して臨床に進んだか、失敗して閉じたか）を意味するからです。実際の場面では特許データを利用して解析をします。
このような解析とメディシナルケミストの経験と洞察力をもとに他社状況を推測しながら自分たちの成功を信じてプロジェクトは進んでいきます。
PCAよりもtSNEのほうが分離能がよく、メディシナルケミストの感覚により近いと言われています。sklearnではPCAをTSNEに変更するだけです。
from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, random_state=0)
tx = tsne.fit_transform(fps)描画するとわかりますが、PCAに比べてよく分離されています。
d = pd.DataFrame(tx)
d.columns = ["PCA1", "PCA2"]
d["DOCID"] = docs
d["COMPANY"] = companies
g = ggplot(aes(x="PCA1", y="PCA2", color="COMPANY"), data=d) + geom_point() + xlab("X") + ylab("Y")
g
今回紹介したPCA,tSNEの他にも色々な描画方法があるので調べてみるとよいでしょう。
今回取り上げた事例では、化合物のフィンガープリントをPCAやtSNEを用いて次元圧縮した結果に企業名を載せると綺麗にマッピングができました。ですが、毎回情報があるとは限りません。このような場合、クラスタリングを行い化合物を適当なグループに分割します。
クラスタリングには、ｋ近傍法や、階層別クラスタリングなどがありますが、これらはどれくらいのサイズのクラスタに分けるべきかを自動的に決定することはできません。今回は、付与するクラスタラベルがほぼ企業名に対応するというように、いい感じに分割したいので密度に基づくクラスタリング手法であるHDBSCANを利用します
インストールはcondaまたは、pipコマンドでインストールできます。可視化にはseabornを使いました。データの読み込みは上のパートと同じなので省略しますが、コード実行時にWarningが出るのを避けるため、RDLogger.DisableLog('rdApp.*')を入れました(jupyter notebook参照)。
pip install hdbscan
#or
conda install -c conda-forge hdbscan
conda install -c conda-forge seabornインストールが終わったら早速使ってみます。
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit import RDLogger
from sklearn.manifold import TSNE
from hdbscan import HDBSCAN
## 以下のパッケージは後半で利用します
from sklearn.model_selection import train_test_split
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.ensemble import RandomForestRegressor
from mlinsights.mlmodel import PredictableTSNE
sns.set_context('poster')
sns.set_style('white')
sns.set_color_codes()
plot_kwds = {'alpha' : 0.5, 's' : 80, 'linewidths':0}
RDLogger.DisableLog('rdApp.*')
seed = 794
oxrs = [("CHEMBL3098111", "Merck" ), ("CHEMBL3867477", "Merck" ), ("CHEMBL2380240", "Rottapharm" ),
("CHEMBL3352684", "Merck" ), ("CHEMBL3769367", "Merck" ), ("CHEMBL3526050", "Actelion" ),
("CHEMBL3112474", "Actelion" ), ("CHEMBL3739366", "Heptares" ), ("CHEMBL3739395", "Actelion" ),
("CHEMBL3351489", "Eisai" )]
fps = []
docs = []
companies = []
mol_list = []
for cid, company in oxrs:
sdf_file = os.path.join("ch08", cid + ".sdf")
mols = Chem.SDMolSupplier(sdf_file)
for mol in mols:
if mol is not None:
mol_list.append(mol)
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
docs.append(cid)
companies.append(company)
fps.append(arr)
fps = np.array(fps)
companies = np.array(companies)
docs = np.array(docs)
trainIDX, testIDX = train_test_split(range(len(fps)), random_state=seed)まずtSNEで企業のラベルとともにケミカルスペースを眺めましょう。
tsne = TSNE(random_state=seed)
res = tsne.fit_transform(fps)
plt.clf()
plt.figure(figsize=(12, 6))
sns.scatterplot(res[:,0], res[:,1], hue=companies, **plot_kwds)
きれいに分かれているのが確認できたので、次はHDBSCANを実行し、その結果を可視化してみます。
HDBSCANのメソッドはScikit-learnに準拠しているので、Scikit-learn同様にオブジェクトを作成し、fitを呼び出すだけで完了ですが、クラスタリング実施の際に、各データの特徴量に基づく距離を指定する必要があります。HDBSCANにはデフォルトでユークリッド距離、コサイン距離など多くのメトリクスが用意されています。
今回は化合物の距離なのでタニモト距離を使いますが、これは自分で実装する必要があります(tanimoto_dist)。このようにユーザー定義関数を利用する場合は、metric='pyfunc’として、func引数に定義した関数を渡します。
def tanimoto_dist(ar1, ar2):
a = np.dot(ar1, ar2)
b = ar1 + ar2 - ar1*ar2
return 1 - a/np.sum(b)
clusterer = HDBSCAN(algorithm='best', min_samples=5, metric='pyfunc', func=tanimoto_dist)
clusterer.fit(fps)計算が終わると、clustererオブジェクトのlabels_というプロパティにラベルが付与されるので、それを使ってプロットしてみましょう。今回はmin_sample=5とし、クラスタを形成時に最低限含まれるべき化合物数を5としました。このパラメータを増減すると、異なる出力が得られますので、ある程度試行錯誤してみるのも良いでしょう。 min_cluster_sizeは、クラスターサイズの規定をするパラメータです。この値より小さいクラスターはノイズとして扱われます。min_samplesは、考慮するコア（クラスタ中心）近傍のサンプル数です。このパラメータを大きくするほど、コンサバティブ（より近傍が多いもののみをクラスタにする）なモデルを作ります。クラスタリングの結果に正解不正解といったものはなく、実務に落とし込んだ際に、メディシナルケミストが見てリーズナブルかどうかが大事かと思います。従って、プロジェクトを一緒に勧めてるメディシナルケミストの方と結果を見ながら良い塩梅のところを見つけるのが良いかなと思います。公式ドキュメントに例があるので興味のある方はご覧ください。
*プロットで灰色になっている部分は、どのクラスタにも属さないと判断されたものです。
plt.clf()
plt.figure(figsize=(12, 6))
palette = sns.color_palette()
cluster_colors = [sns.desaturate(palette[col], sat)
if col >= 0 else (0.5, 0.5, 0.5) for col, sat in zip(clusterer.labels_, clusterer.probabilities_)]
plt.scatter(res[:,0], res[:,1], c=cluster_colors, **plot_kwds)
おおよそtSNEで分離されているクラスタとラベルが対応しているのが確認できます。このようにクラスタ数を予め決めずに、密度に基づいてクラスタを分割する手法としてHDBSCANは便利です。
さて、ここまでで化合物の特徴を元にそれを低次元のケミカルスペースにマッピングする手法を学びました。
PCAやtSNEを使ってマップしたときに、次に作る化合物はこの空間のどこに位置するのだろうか？と考えることは多いです。また、HTSライブラリの拡張など際に、今までにないケミカルスペースを埋めたい場合も既存のケミカルスペースとの比較がしたくなります。
Scikit-learnのPCAには、transform()メソッドがあり、作成した主成分空間に新しい化合物をマッピングできます。一方でtSNEはfit_transform()メソッドはありますが、transform()メソッドはありません。したがって、新しい化合物をマップしたい場合、全化合物で計算する必要があります。手法によらず、新しいデータを既存のケミカルスペースへマッピングする方法はないのものでしょうか。予測モデルの作成と同じ要領で、ケミカルスペースの予測モデルがあればいいのではないでしょうか。 つまり、トレーニングデータを使い特徴量＝＞低次元空間への次元圧縮の予測モデルを作ります。こうすることで、既存のケミカルスペースに新しい化合物をマッピングできるようになります。これを実装してるパッケージがmlinsightsです。このパッケージはpipコマンドでインストールできます。
pip install mlinsightsmlinsightsは複数の手法を実装していますが今回は、PredictableTSNEのみの事例紹介に絞ります。 ランダムにトレーニングとテストデータに分割しトレーニングデータでTSNEを実施、その結果をRandomForestとGausianProcessRegressorの2つの手法で学習、このモデルを用いてテストデータをマッピングするというテストをします。
trainFP = [fps[i] for i in trainIDX]
train_mol = [mol_list[i] for i in trainIDX]
testFP = [fps[i] for i in testIDX]
test_mol = [mol_list[i] for i in testIDX]
allFP = trainFP + testFP
tsne_ref = TSNE(random_state=seed)
res = tsne_ref.fit_transform(allFP)
plt.clf()
plt.figure(figsize=(12, 6))
sns.scatterplot(res[:,0], res[:,1], hue=['train' for i in range(len(trainFP))] + ['test' for i in range(len(testFP))])ランダムスプリットなので、まんべんなくサンプリングできています。
PredictableTSNEもscikit-learnと同じなので、オブジェクトを作ってfitで学習、Transformデータの次元削減をします。transformerにはTSNEやPCA次元圧縮用のオブジェクトを、Estimatorにはその結果を元にモデルを作る学習機を渡します。後でトレーニングデータとテストデータを重ねるのでkeep_tsne_outputs=Trueとしています。
一点注意としてfitを呼び出すときに本来不要なyに対応する値を渡す必要があります。下のコードではダミーの値としてただのIndexを与えています。
rfr = RandomForestRegressor(random_state=seed)
tsne1 = TSNE(random_state=seed)
pred_tsne_rfr = PredictableTSNE(transformer=tsne1, estimator=rfr, keep_tsne_outputs=True)
pred_tsne_rfr.fit(trainFP, list(range(len(trainFP))))
pred1 = pred_tsne_rfr.transform(testFP)
plt.clf()
plt.figure(figsize=(12, 6))
plt.scatter(pred_tsne_rfr.tsne_outputs_[:,0], pred_tsne_rfr.tsne_outputs_[:,1], c='blue', alpha=0.5)
plt.scatter(pred1[:,0], pred1[:,1], c='red', alpha=0.5)
赤いポイントがテストデータで意外といい感じです。 次はGaussianProcessで同じことをやってみます。
gbr = GaussianProcessRegressor(random_state=seed)
tsne2 = TSNE(random_state=seed)
pred_tsne_gbr = PredictableTSNE(transformer=tsne2, estimator=gbr, keep_tsne_outputs=True)
pred_tsne_gbr.fit(trainFP, list(range(len(trainFP))))
pred2 = pred_tsne_gbr.transform(testFP)
plt.clf()
plt.figure(figsize=(12, 6))
plt.scatter(pred_tsne_gbr.tsne_outputs_[:,0], pred_tsne_gbr.tsne_outputs_[:,1], c='blue', alpha=0.5)
plt.scatter(pred2[:,0], pred2[:,1], c='red', alpha=0.5)
GPモデルはあまり性能が出ませんでした。PredictableTSNEは新しい化合物の既存スペースへマッピングの一つのツールとなり得ますが、この例のようにモデルにより結果が大きく左右されることがあるようです。 実プロジェクト投入時にはよく検証しないとメディシナルケミストの方をミスリードするリスクがあるかもしれませんが一つのアプローチとして紹介しました。
|
Tip
|
新しいバージョンのpandasを利用しているとggplotを呼ぶ際にエラーが起きることがあるかもしれません。github_issue 解決策はこのリンクにあるようにインストールされたggplotのフォルダ内のコードggplot/utils.pyのpd.tslib.Timestampをpd.Timestampに、ggplot/stats/smoothers.pyのfrom pandas.lib import Timestamp を from pandas import Timestampに変えると動くと思います。 |
化学構造と生物学的活性における相関関係をStructure Activity Relationship(SAR)またはQuantative SAR(QSAR)と呼びます。一般的には似たような化合物は似たような生物学的活性を示すことが知られており、この相関関係を理解しドラッグデザインに活かすことが創薬研究において大変重要です。
また、このような問題には細胞の生死、毒性の有無といった化合物がどのクラスに入るのかを推定する分類問題と阻害率（%inhibition）といった連続値を推定する回帰問題の2つがあります。
ChEMBLからhERG阻害アッセイの73データを用いてIC50が1uM未満のものをhERG阻害あり、それ以外をhERG阻害なしとラベルします。
まずは必要なライブラリをインポートします。
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix, f1_score
from sklearn.ensemble import RandomForestClassifierChEMBLでダウンロードしたタブ区切りテキストの処理は8章とほぼ同じですが、今回は活性データが欲しいのでSTANDARD_VALUEという列を探して数値を取り出します。この値が1000nM未満であればPOSというラベルを、そうでなければNEGというラベルを振ります。最後にラベルをnumpy arrayにしておきます。
mols = []
labels = []
with open("ch09_compounds.txt") as f:
header = f.readline()
smiles_index = -1
for i, title in enumerate(header.split("\t")):
if title == "CANONICAL_SMILES":
smiles_index = i
elif title == "STANDARD_VALUE":
value_index = i
for l in f:
ls = l.split("\t")
mol = Chem.MolFromSmiles(ls[smiles_index])
mols.append(mol)
val = float(ls[value_index])
if val < 1000:
labels.append("POS")
else:
labels.append("NEG")
labels = np.array(labels)続いてmolオブジェクトをフィンガープリントに変換します。このフィンガープリントからhERG阻害の有無を予測するモデルを作成します。
fps = []
for mol in mols:
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
fps.append(arr)
fps = np.array(fps)データセットを訓練セットテストセットの2つに分けます。テストセットは作成した予測モデルの精度を評価するためにあとで使います。
x_train, x_test, y_train, y_test = train_test_split(fps, labels)予測モデルを作成するにはインスタンスを作成してfitメソッドで訓練させるだけです
rf = RandomForestClassifier()
rf.fit(x_train, y_train)先程分割しておいたテストセットを予測します。
y_pred = rf.predict(x_test)Confusion matrixを作成します。
Confusion matrixとはクラス分類の結果をまとめた表です。クラスを正しく分類できているかをわかりやすく視覚化できTP,TNが多くFP,FNが少ないほどよく分類できているということになります。
| Actual class | |||
|---|---|---|---|
Positive |
Negative |
||
Predicted class |
Positive |
True Positive(TP) |
False Positive(FP) |
Negative |
False Negative(FN) |
True Negative(TN) |
|
confusion_matrix(y_test, y_pred)
#array([[11, 1],[ 5, 2]])| 11 | 1 |
|---|---|
5 |
2 |
F1スコアを見てみましょう。
f1_score(y_test, y_pred, pos_label="POS")
#0.4あまりよくないですね。
|
Note
|
train_test_split関数がランダムに訓練セットとテストセット分割するので、Confution matrix,F1スコアの値は実行するたびに変わります。 |
-
正しいと予測されたもののうち本当に正しいものの割合を適合率と呼ぶ precision=TP/(TP+FP)
-
正しいものが正しいと予測された割合を再現率と呼ぶ recall=TP/(TP+FN)
F1スコアは適合率と再現率の調和平均となり
F1 = 2 * (precision * recall) / (precision + recall)
で計算されます。
回帰モデルは最初に説明したとおり、連続値を予測するモデルとなります。今回はRandomForestの回帰モデルを作成して、その精度をR2で評価します。データは分類問題で使ったhERGのアッセイデータを利用することにしましょう。最初に必要なライブラリをインポートします。
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from math import log10分類問題のときにはラベル化しましたが、今度は連続値を予測したいのでpIC50に変換します。(なぜpIC50にすると都合が良いのかはそのうち補足する)
pIC50s = []
with open("ch09_compounds.txt") as f:
header = f.readline()
for i, title in enumerate(header.split("\t")):
if title == "STANDARD_VALUE":
value_index = i
for l in f:
ls = l.split("\t")
val = float(ls[value_index])
pIC50 = 9 - log10(val)
pIC50s.append(pIC50)
pIC50s = np.array(pIC50s)データセットをトレーニングセットとテストセットの2つに分割します。フィンガープリントは分類モデルのときに作成したものを流用します。
x_train, x_test, y_train, y_test = train_test_split(fps, pIC50s)訓練します。Scikit-learnの場合はこの手順はどの手法でもほぼ同じメソッドでfitしてpredictです。
rf = RandomForestRegressor()
rf.fit(x_train, y_train)予測しましょう。
y_pred = rf.predict(x_test)予測精度をR2で出してみます。
r2_score(y_test, y_pred)
#0.52まずまずといったところでしょうか。
回帰の当てはまりの良さへの1つの評価指標としてよく利用されるもので決定係数とも呼ばれます。
今回紹介した手法は似たような化合物は似たような生物学的活性を示すという仮設に基づいて生成されるモデルです。もしトレーニングセットに似ている化合物が含まれなかった場合の予測精度はどうなるのでしょうか？
当然その場合は予測された値は信頼できませんよね。つまり、予測値にはその予測が確からしいか？という信頼度が常についてまわります。そのようなモデルが信頼できる、または適用できる範囲をapplicability domainと呼びます。これに関しては明治大学金子先生のモデルの適用範囲・モデルの適用領域が詳しいです。
Hugo Kubinyi先生は似ている化合物は果たして似た活性を示すのか？という疑問を、estradiolのOH基をMethoxy基に変換すると活性が消失する例を上げて説明されていました。その当時のジャーナルが見つからなかったのですが大体次のようなものでした。
Estrogen Receptorのリガンドであるestradiolの水酸基をMeO基に置換した右側の化合物は活性を示しません。両者の類似性は77%と非常に高いにも関わらずです。
この活性消失の理由は複合体結晶構造から綺麗に説明されます。水酸基はGLU353と水素結合を形成していますが、これをMeOに置換すると水素結合能が失われアゴニスト活性がなくなります。
以上が我々が定義した類似性の尺度が蛋白質の分子認識の尺度と一致しないという例です。上記のように必ずしも化合物の類似度から活性が予測できるわけではなく、極めて類似度が高いのに活性が消失してしまうことが多々あります。特にMMPの文脈で説明されたActivity Cliffはこのような事象にそれっぽい名前をつけたものですし、マジックメチルと呼ばれるメチル基導入の効果もコンフォメーションが大きく変化した結果立体的な類似性が減少したという事象に名をつけたような気がします。
applicability domainはそもそも、類似度（距離）という尺度が正しいと仮定をおいて、トレーニングセットの類似性からその予測が信頼できるかという確度を測る手法なのでその辺りは注意したほうがよいでしょう。リード最適化だと、似たような化合物を作るのが普通なのでapplicability domainの問題よりも上記の類似性の尺度が細かいところで狂ってくるほうが問題視されることが多いです。
本章からディープラーニングを利用して、QSARモデルや生成モデルを作成します。
生物の脳には神経細胞が存在し、それらがネットワークを形成することで情報を伝達したり、記憶や学習しています。このネットワーク構造を数理モデル化したものがArtificial Neural Network(ANN)です。
一般的なANNは、学習のための情報を入れる入力層、入力情報のパターンを元に反応（神経シナプスの発火に対応）を学習する中間層（または隠れ層）、最後の出力層の三層から構成されていますが、ディープラーニングは隠れ層を多層にすることで高精度な予測を可能とします。
本書では特にこれ以上の説明はしませんが、自分でゼロからコードを書いて理解したい場合はゼロから作るDeep Learningが助けになります。また、理論についてきちんと学びたい場合は深層学習をお薦めします。
Tensorflowは Googleが開発してOSSとして公開している機械学習のためのフレームワークです。主にディープラーニングのフレームワークとして使われることが多いです。
NOTE: Tensorflowは最近1.xから2.xにメジャーアップデートをしましたが、2.x版はまだ登場したばかりで参考情報が少ないので、1.x系を利用します。また同じ1.xでもバージョンによってAPIが異なるので、動かしたいコードがあった場合、どのバージョンで書かれているかに気をつけてください。
KerasはTensoflowなどの低レベルフレームワークをバックエンドにした高レベルAPIで、より簡単にコードが書けます。KerasはTensorflowとは独立して開発されてきましたが、最近、TensorflowはKerasを同梱するようになりました。そのため別にイントールせずにKerasを利用できます。Tensorflow同梱バージョンのKerasが本家の最新バージョンになっていないことがあります。
どちらのKerasを使うのが良いかは悩ましいところですが、本書では利便性のためにTensorflow統合版Kerasを利用します。
KerasとTensoflowをオフィシャルブログを参照しつつ少し整理しておきます。元々KerasはTensoflowとは別のプロジェクトとして開発されており（もちろん今もです）Kerasを使うにはTensorflowとはインストールする必要がありました。しかし2017年Keras2.xのメジャーバージョンアップのタイミングのあたりでTensoflowプロジェクトがKerasを統合するようになりました。以下の英文は上記リンクの記事の抜粋です。現在はTensorflowからKerasを呼び出し利用することが可能になっています。
TensorFlow integration Although Keras has supported TensorFlow as a runtime backend since December 2015, the Keras API had so far been kept separate from the TensorFlow codebase. This is changing: the Keras API will now become available directly as part of TensorFlow, starting with TensorFlow 1.2. This is a big step towards making TensorFlow accessible to its next million users.
Tensorflow とKerasをインストールしてみましょう。 anacondaでインストールする場合、GPU対応バージョンを使うか、CPUバージョンを使うかでインストールするパッケージが少し異なります。
# CPU版
$ conda install -c conda-forge tensorflow
# GPU版
$ conda install -c anaconda tensorflow-gpu|
Note
|
pipコマンドを利用してTensorFlowをインストールすることもできます。その場合は公式ドキュメントを参照してください。しかし基本的にはCondaで環境を作ったらCondaでパッケージを入れることが望ましいでしょう。 |
参考リンク
Google colaboratoryはクラウド上で実行できるJupyter notebook環境です。Theano, Thensorflow, Keras, Pytorchなどのディープラーニング用のフレームワークがインストール済みなのと時間の制限はありますがGPUが使えるため、手元にGPUマシンがなくてもディープラーニングを利用できる点が非常に魅力的です。
利用にはGoogleのアカウントが必要なので、もしGoogleアカウントを持っていなければこの機会にアカウントを取得し利用してみると良いでしょう。
Google のアカウントがある方はGitHub状のノートブックをそのままColab上で実行することもできます。以前Mishima.sykで使ったScikit-learnハンズオンのノートブックを開いてみます。
- NOTE
-
@y__samaが作成したノートブックですが、データの準備から、AutoSklearn まで一通り学習できるようになっています。
まずGoogle colaboratoryにアクセスします。もし下のような画面にならなければ左上の「ファイル」から「ノートブックを開く」を実行してください
次にGitHubというタブをクリックし、以下のURLをコピーペーストするとJupyter Notebookからコードを動かせます。
ノートブックを開くとJupyter Notebookと同じような画面になります。Shift＋リターンキーでセルのコードを実行できます。
GoogleColabデフォルトで利用できるライブラリを確認するにはセルの中で '!pip freeze’と打つと列記されます。
-
absl-py==0.7.0
-
alabaster==0.7.12
-
たくさん出てくる
-
yellowbrick==0.9.1
-
zict==0.1.3
-
zmq==0.0.0
本章ではDNNを利用して構造活性相関解析をします。
はじめにDNNを利用したシンプルな予測モデルを構築してみます。ここでは9章と同じデータを使います。最初に分類モデルを作成し、Positiveのラベルを[0, 1], Negativeのラベルを[1, 0]の二次元のOneHotベクトルで表します。KerasのModelオブジェクトを利用してモデルを作成した場合、上記の二次元のそれぞれの期待値が得られます。どちらのクラスに属する可能性が高いかを知るにはNumpyのArgmax関数を使えばよいです。
|
Note
|
OneHotベクトルとはある一つの値が1でそれ以外が0になるようなベクトルです。10クラスの分類問題を考えた場合[1, 0, 0, 0, 0, 0, 0, 0, 0, 0]のようにどこかが1で残りの9個が0になるようなベクトルでクラスを表現できます。上の例ではPositive/Negativeの2クラスですので、OneHotベクトルは2次元になっています。 |
必要なライブラリをインポートします。
from rdkit import Chem, DataStructs
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import IPythonConsole
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix, f1_score
from tensorflow.python.keras.layers import Input
from tensorflow.python.keras.layers import Dense
from tensorflow.python.keras.layers import Dropout
from tensorflow.python.keras.layers import Activation
from tensorflow.python.keras.model import Model次にデータを読み込みます。9章では”POS”/”NEG”をlabelsというリストに入れたので一次元表現でしたが、今回はここが二次元になっています。
mols = []
labels = []
with open("ch09_compounds.txt") as f:
header = f.readline()
smiles_index = -1
for i, title in enumerate(header.split("\t")):
if title == "CANONICAL_SMILES":
smiles_index = i
elif title == "STANDARD_VALUE":
value_index = i
for l in f:
ls = l.split("\t")
mol = Chem.MolFromSmiles(ls[smiles_index])
mols.append(mol)
val = float(ls[value_index])
if val < 1000:
labels.append([0,1]) # Positive
else:
labels.append([1,0]) # Negative
labels = np.array(labels)続いて分類モデルと回帰モデルを順次作成します。
まずは回帰モデルで、入力は9章と同じECFPを利用しています。DNNの構築には入力データの次元を明示的に指定する必要があるためnBitsという変数を定義しています。
|
Tip
|
train_test_splitにrandom_stateで適当な整数を指定すると毎回同じデータが得られるので検証の際に有用です。 |
nBits = 2048
fps = []
for mol in mols:
fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=nBits)
arr = np.zeros((1,))
DataStructs.ConvertToNumpyArray(fp, arr)
fps.append(arr)
fps = np.array(fps)
x_train1, x_test1, y_train1, y_test1 = train_test_split(fps, labels, random_state=794)入力が2048次元、300ニューロンの全結合層が三層、最後の出力層が2となるニューラルネットワークを作成します。活性化関数にはReLU, 出力層には二次元の多クラス分類のためにSoftmaxを用いました。
Dropout層はランダムにニューロンを欠損させることにより過学習を防ぐ役割を果たします。
Kerasではモデルを定義した後compile関数を呼ぶことでモデルを構築します。optimizer, lossは目的に応じて変更する必要があります。今回はlossとして’categorical_crossentropy', optimizerとして’adam’を使いましたがoptimzerはadam以外にも多くあるのでどれが適切かは実際は試行錯誤が必要となるでしょう。
# Define DNN classifier model
epochs = 10
inputlayer1 = Input(shape=(nBits, ))
x1 = Dense(300, activation='relu')(inputlayer1)
x1 = Dropout(0.2)(x1)
x1 = Dense(300, activation='relu')(x1)
x1 = Dropout(0.2)(x1)
x1 = Dense(300, activation='relu')(x1)
output1 = Dense(2, activation='softmax')(x1)
model1 = Model(inputs=[inputlayer1], outputs=[output1])
model1.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'])|
Note
|
KerasにはSequentialモデルが用意されており、これを使うことで上記の例（Functional API）よりもシンプルにネットワークを記述できます。今回Functional APIでモデルを定義したのは、こちらに慣れておくと入力が複数の場合やより複雑なモデルの構築にも対応しやすいからです。もしSequentialの書き方に興味がある方は公式サイトやQiitaを調べてください。 |
|
Note
|
DNNは初期のランダムに発生させた重みに基づいて予測した予測値と実際の値を比較し、その差（LOSS)を最小化するように重みを更新するBackpropagationという手順を繰り返しながらモデルを最適化します。この繰り返しの回数を指定するのがEpochsです。Epochsを増やすとどんどん賢くなるように思われるかもしれませんが、計算コストがかかることと過学習のリスクもあるので長ければ良いわけではないのでLoss/Accuracyなどを観測して適切なEpoch数を考えましょう。 |
トレーニングデータを用いてEpoch毎に正解の値と予測値の誤差を小さくするように重みを調整していきます。十分な量のトレーニングデータを用い学習したとしても繰り返しすぎると、同じトレーニングデータで何度も何度も学習することになるのでモデルの汎化性能が落ちてしまいます。
過学習の判断にはEpoch毎、Training set/ Validation setの精度を評価しプロットするとTraining setの精度が向上する一方でValidation setの精度が変化しないまたは悪くなっているかどうかを調べればよいです。KerasにはEarly stoppingという機能があり、ある一定回数以上学習を進めてもモデルの性能が変化しなかったらそこで学習を打ち切るといったこともできます。
詳しくはEarly stoppingの紹介と参考文献をご覧ください。
モデルを構築したら後はScikit-learnと同じ感覚でfit/predictが行えます。
hist1 = model1.fit(x_train1, y_train1, epochs=epochs)最後に結果を可視化してみます。
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(range(epochs), hist1.history['acc'], label='acc')
plt.legend()
plt.plot(range(epochs), hist1.history['loss'], label='loss')
plt.legend()今回の例ではだいたい6Epochくらいでモデルが良い精度になりました。
次にテストデータで検証します。
y_pred1 = model1.predict(x_test1)
y_pred_cls1 = np.argmax(y_pred1, axis=1)
y_test_cls1 =np.argmax(y_test1, axis=1)
confusion_matrix(y_test_cls1, y_pred_cls1)ちょっと微妙でしょうか、、、
回帰モデルも基本的には先ほどの分類問題と同じです。今度は回帰なので最後の出力層は値そのもの、つまり一次元になります。また活性化関数はSigmoidなどでは0-1になってしまうのでLinearとしています。学習データは9章のコードを流用しています。
from math import log10
from sklearn.metrics import r2_score
pIC50s = []
with open("ch09_compounds.txt") as f:
header = f.readline()
for i, title in enumerate(header.split("\t")):
if title == "STANDARD_VALUE":
value_index = i
for l in f:
ls = l.split("\t")
val = float(ls[value_index])
pIC50 = 9 - log10(val)
pIC50s.append(pIC50)
pIC50s = np.array(pIC50s)
x_train2, x_test2, y_train2, y_test2 = train_test_split(fps, pIC50s, random_state=794)次にモデルを定義します。Lossの部分が先ほどの分類モデルとは異なり、MSEになっていることに注意して下さい。
epochs = 50
inputlayer2 = Input(shape=(nBits, ))
x2 = Dense(300, activation='relu')(inputlayer2)
x2 = Dropout(0.2)(x2)
x2 = Dense(300, activation='relu')(x2)
x2 = Dropout(0.2)(x2)
x2 = Dense(300, activation='relu')(x2)
output2 = Dense(1, activation='linear')(x2)
model2 = Model(inputs=[inputlayer2], outputs=[output2])
model2.compile(optimizer='adam', loss='mean_squared_error')ここまでできたら後は同じです。
hist = model2.fit(x_train2, y_train2, epochs=epochs)
y_pred2 = model2.predict(x_test2)
r2_score(y_test2, y_pred2)
plt.scatter(y_test2, y_pred2)
plt.xlabel('exp')
plt.ylabel('pred')
plt.plot(np.arange(np.min(y_test2)-0.5, np.max(y_test2)+0.5), np.arange(np.min(y_test2)-0.5, np.max(y_test2)+0.5))いかがでしょうか。予測モデルはちょっとUnderEstimate気味ですかね。DNNは重ねるレイヤーの数、ドロップアウトの割合、隠れ層のニューロンの数、活性化関数の種類など数多くのパラメータをチューニングする必要があります。今回の例は決め打ちでしたが、色々パラメータを変えてモデルの性能を比較してみるのも面白いです。
さて、ここまで分子のフィンガープリントを入力としてRandomForestやDNNのモデルを作成してきました。DNNが大きく注目を浴びた理由の一つに人が特徴量を抽出しなくてもモデルが特徴量を認識してくれるということが挙げられます。
例えば画像の分類においては、からSIFTという特徴量を人が定義し、これを入力としたモデルが作られていましたが、現在のDNNにおいては基本的に画像のピクセル情報そのものを利用しています。
ケモインフォマティクスに置き換えてみると、SIFTは分子のフィンガープリントに相当します。ですのでここ(入力)をもっとPrimitiveな表現に変えることでDNNの性能が上がるのではないか？と考えるのは至極当然の流れです。2015年、Harvard大学の, Alan Aspuru-Guzikらのグループは一つのチャレンジとしてNeural Finger print/NFPというものを提唱しました。
今まで利用してきたECFPとNFPとの違いを、彼らの論文中の図を引用して示します。
ECFP(Circular Fingerprints)は入力の分子それぞれの原子からN近傍（Nは任意）までの原子までの情報をHash関数（この例ではMod）任意の値に変換、で固定長のベクトルに直すといったものでした。ざっくりいうと部分構造の有無を0/1のビット情報に直したものを利用するといったイメージです。一方、今回紹介するNFPはECFPにコンセプトは似ているのですが、Hash関数の部分がSigmoidに、Modで離散化する部分がSoftmaxになっています。従って入力されるデータセットによりECFPよりも柔軟に分子のフィンガープリントを生成することが期待されます。
この論文が発表されて以降、数多くの実装がGitHubに公開されていますが、各実装ごとにKerasでもBackendがTheanoであったり、Keras/Tensorflowであっても、Keras1.xじゃないと動作しなかったりと意外と環境依存のものが多く扱いにくい状況になっています。残念なことに今回構築した環境で動作するものが公開されていませんのでKeras2.x/Python3.6で動作するものをこちらのコードをベースに作成しました。
SIFTが提案されたのは1999年です。原著論文によると、物体（画像）認識においてピクセルそのものを扱う場合の難しさは、オブジェクトの位置、回転、大きさ（スケール）、光度などが異なるものを扱うところにあるようです。これら、変動する値を普遍的な特徴量に変える方法が種々研究されていたようです。ピクセルそのものを使う方法が全くないわけではなく、私が機械学習を勉強する際に購入したpythonではじめる機械学習には人の顔の画像データを学習して分類する例が載っています。ここではピクセルデータを入力に、主成分分析により顔の特徴を抽出し分類をしています。この問いに関して明確に有効だったという文献を探せていませんが、タスクによっては有効であったと思います。もし詳しい方いたらぜひコメントください。
git clone https://github.com/iwatobipen/keras-neural-graph-fingerprint.gitexample.pyというファイルのコードを眺めるとなんとなく雰囲気がつかめると思います。分子の表現は、これまでの例はフィンガープリントをRDKitを使い生成していましたが、今回はこのフィンガープリントそのものをDNNが学習します。
ということで、分子をグラフとして表現したものが入力になります。Atom_matrixとして(max_atoms, num_atom_features)をEdge_matrixとして(max_atoms, max_degree)をbond_tensorとして(max_atoms, max_degree, num_bond_features)という三つの行列を使います。分子はそれぞれ原子数が異なるためmax_atomsで最大原子数を定義しています。こうすることで分子ごとに同一の行列サイズの入力となりバッチ学習が可能となります。
Exampleを実行するのであれば下記のコマンドを入力してください。
python example.py参考リンク
Deep Learningがメディナルケミストリに大きなインパクトをもたらしたものの一つに生成モデルがあげられます。特にこの数年での生成モデルの進化は素晴らしいです。ここではMarcus Olivecronaにより開発されたREINVENTを使って新規な合成案を提案させてみましょう。
11章で構築した予測モデルを一般に識別モデルといいます。一方で入力の分布をモデル化することでモデルからのサンプリングつまり入力データを生成することが出来るようになります。これを生成モデルといいます。
詳しくはPRMLの1.5.4を読むのをおすすめします
pytorchというディープラーニングのライブラリをcondaでインストールします。新しいバージョンでは動かないのでバージョンを指定してインストールします。
keras同様TensorFlowをより便利に使うためのライブラリです。
$ conda install pytorch=0.3.1 -c pytorch続いてREINVENT本体をGitHubからクローンします。
$ cd <path to your working directory>
$ git clone https://github.com/MarcusOlivecrona/REINVENT.git続いて、ChEMBLの110万件くらいのデータセットで予め訓練済みのモデルをダウンロードしてきて元のデータと置き換えます。 このデータはGTX 1080TiGPUマシンを利用して5,6時間かかっていますのでもしトレーニングを自分で行うのであればGPUマシンは必須です。
$ wget https://github.com/Mishima-syk/13/raw/master/generator_handson/data.zip
$ unzip data.zip
$ mv data ./REINVENT/これで準備が整いました。
ここではJanuviaとして知られる抗糖尿病薬sitagliptinの類似体を生成するようなモデルを作成してみます。
まずはtanimoto係数をスコアとして類似度の高い構造を生成するようにモデルを訓練します。今回は3000ステップ訓練しますが、大体ちょっと前のMacbook Airで7,8時間かかるので気長に待ちましょう。待てない場合はここのデータを使ってください。
./main.py --scoring-function tanimoto --scoring-function-kwargs query_structure 'N[C@@H](CC(=O)N1CCn2c(C1)nnc2C(F)(F)F)Cc3cc(F)c(F)cc3F' --num-steps 3000 --sigma 80ここからはjupyter notebookを立ち上げます。
必要なライブラリを読みこみます。sys.path.appendはREINVENTのディレクトリを指定してください。
%matplotlib inline
import sys
sys.path.append("[Your REINVENT DIR]")
from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs, Draw
import torch
from model import RNN
from data_structs import Vocabulary
from utils import seq_to_smiles続いて、トレーニングしたモデルから50化合物サンプリングします。
voc = Vocabulary(init_from_file="/Users/kzfm/mishima_syk/REINVENT/data/Voc")
Agent = RNN(voc)
Agent.rnn.load_state_dict(torch.load("sitagliptin_agent_3000/Agent.ckpt"))
seqs, agent_likelihood, entropy = Agent.sample(50)
smiles = seq_to_smiles(seqs, voc)実際にどんな構造が生成されたのか見てみましょう。
mols = []
for smi in smiles:
mol = Chem.MolFromSmiles(smi)
if mol is not None:
mols.append(mol)
Draw.MolsToGridImage(mols, molsPerRow=3, subImgSize=(500,400))まずまずといったところでしょうか？
- NOTE
-
どのあたりに興味があって、さらに知りたいというようなリクエストをissueに投げるかtwitterででもreplyしてもらえると。またおすすめの提案も助かります。
Pattern Recognition and Machine Learning(PRML)を一通り読めるようになることを目指すといいでしょう。英語のpdfは先のリンクからフリーでダウンロードできますが、日本語はパターン認識と機械学習 上, 下で販売されています。
PRMLが手強いという方は「PRMLを読む前に」などで検索するともう少し易しい本が出てくるので自分に合ったものを選ぶとよいと思います。
本書はAI創薬を念頭においているため、機械学習や解析手法の基本的なところを説明していますが、本来ケモインフォマティクスはバイオインフォマティクスと同じように、分子の表現方法や、データの効率的な格納方法、高速な検索技術をも含んでいます。このような情報学（ITの側面）としてのケモインフォマティクスに興味がある方はChemoinformatics: Basic Concepts and Methodsから読み進んで興味のあるトピックを深掘りすることをお薦めします。
製薬企業やアカデミアの薬物動態や毒性、薬理に属している方で本書の先を知りたい方は是非ともDrug-Like Properties: Concepts, Structure Design and Methods from ADME to Toxicity Optimizationを一読することをお薦めします。本書は大体製薬企業の合成部門に配属された新入社員が輪読するようなテキストなので、本書を読み終えた方であれば楽しく読めると思います。もし理解が追いつかない部分があれば関連書籍をあたってもいいですし、この本を手がかりにさらに学習しても良いと思います。
また、薬物動態に携わっている方が本書のQSAR/QSPRを理解できていればPBPKモデリングの際に強みとして活かせるはずです。動態プロファイルの最適化は薬品の差別化戦略にとても大きな意味を持つのでQSPR+PBPKに強いというのは非常に重宝されるのではないでしょうか？
本書では低分子化合物をもとにしたインフォマティクスの手法を紹介してきましたが、結果の解釈をする場合にはターゲットとなる蛋白質の理解は必須になります。つまり蛋白質の立体構造の理解無しにはドラッグデザインはできないということです。そのためSBDD関連の本を読んで学ぶと良いでしょう。
|
Note
|
残念ながら私はSBDDを本で勉強したことがないのでだれか良書を教えてください |
さらにSBDDは蛋白質を扱うことになるので、必然的にchemoinformatics,bioinformaticsと括る必要がなくなります。創薬の枠組みのなかでどちらも理解しているほうがより深い考察ができるようになるので両方できるようになりましょう。そのほうが絶対楽しいです。DRY解析教本や生命科学データ解析を支える情報技術もきっとあなたのキャリアの助けになるでしょう。
もう一つ、6章で触れましたが、蛋白質とリガンドの相互作用を理解するのに量子化学計算は重要です。特にこれからのSBDDで量子化学に基づいて相互作用を解釈する能力は必須になると断言できます。毛嫌いせずに基礎量子化学―軌道概念で化学を考えるなどを読んでください。もしGamessを使うなら新版 すぐできる 量子化学計算ビギナーズマニュアルも助けになることでしょう。少なくともenergy decomposition analisysは押さえておきましょう、計算の解釈能力とプロジェクトへの貢献度が上がります。さらにFMOは言うまでもないことですが不可欠のツールですので、各成分を理解しておくとよりよりドラッグデザインの助けになります。