in silico 創薬のためのC/C++プログラミング環境

ブログ開始

 こんにちは、私はインテージヘルスケアの金井といいます。in silico創薬における受託解析計算などを仕事にしています。今回、私もブログなるものを書くことになりました。何分初めてのことで、in silico創薬関連のどの様な事を書けば良いのか思いつきません。創薬からは少し遠いですが、プログラミング関連の話題を扱っていこうと思います。私が良く使うプログラミング言語は C++で、良く使うOSはLinuxです。というか、プログラミングに関して言えばLinuxしか使っていません。C++は良い言語だと思っていますが、最近は Pythonも意識して使うようにしています。実はPythonが結構苦手です。もしかしたら、Pythonの勉強の成果みたいなこともブログにするかもしれません。

 ブログに慣れてきたら有名なChemoinformaticsのツールキットであるOpenBabelやRDKitを使ったプログラミングの話とかも書いてみようと思います。もちろん私が苦手なPythonは使いませんよ。C++ですよ。Pythonを使ったRDKitの話題はウェブ上に色々あるけど、C++ならあまり見あたらないし、ネタ被りしなそうだから丁度良いかもしれませんね。

プログラミングって

 色々なプログラミング言語についての入門書やウェブサイトが世の中にはたくさんあります。新しい言語を習得しようとするとき、このような書籍などを読んで、サンプルプログラムなど実行してみることでしょう。しかしながら、そのサンプルが正常動作するまでの環境構築が非常に面倒だったりします。特にプログラミング超初心者の方などは環境構築の時点で挫折してしまうのではないでしょうか。というかエディタの使い方の時点で挫折しそうですよね。私もPythonを勉強しているのですが、これがまた環境構築が厄介です。Anacondaを利用すると非常に楽に構築出来るらしいのですが、未だ手を出していません。Pythonは是非とも習得したいので、少しづつでも努力しています。

 ところで、私が普段使っているC++も環境依存の激しい言語です。だってLinux上で作ったプログラムをWindows上でコンパイルすることすら簡単に出来ないのですから。プログラミング環境は人によって千差万別です。そんなわけで、今回は自分がC/C++のプログラミングで使用しているツールやライブラリなど少し紹介します。そもそも何故C/C++ なのかというと、数値計算や科学技術計算で広く利用される言語で、私が大学生時代から良く使っていたからです。(あと、数値計算ではFortran もよく使われます)科学技術計算で良く使われるので、何処のスーパーコンピュータにも専用コンパイラ(C/C++, Fortran)がインストールされています。あと、in silico創薬関連ツールの核心部分は大体C/C++で書かれているのではないかと思います。前述したOpenBabelやRDKitもC++です。このへんが私がC++を愛用する理由ですね。

C/C++プログラミングツール

 私がC/C++プログラミングで使っているツールを思いつくままに列挙してみました。

ツール名説明
Emacsプログラムを書くためのエディタ
Vimこれもエディタだが自分はviewerとして使っている
GCCコンパイラ
Clangコンパイラ
CMakeビルドの自動化ツール
DoxygenMakefileの自動生成ツール
Git, Subversionコードのリビジョン管理ツール
GDBデバッガー

 パッと思いつくだけあって、このリストにあるものはC++プログラミングには全部重要だし、全部使います。むしろ未だ何か見落としていると思います。恐らく初心者が「先ずこれらのツールを全部覚えてくださーい」とか言われたら、ちょっと引くんじゃないでしょうか。もちろん入門書にあるサンプルプログラムを動かす程度であれば、CMakeなど全く必要ありません。しかしながら、結局は勉強することになるのではないかと思います。各ツールに関しての話題は今後のブログのネタに出来ればと思います。

 Linux上でプログラミングする利点の一つとして上記のようなツール初めからインストールされているか、ちょっとしたコマンド一発で簡単にインストールできてしまう点にあります。私は普段の仕事でLinuxばかり利用しますが、プログラミングを勉強してみたいという方にもLinuxはおすすめです。

Chemoinformatics関連でよく使うライブラリ

 プログラミングでは必ず何かしらのライブラリを利用します。Hello Worldのような極めて単純なプログラムでも標準ライブラリを使用しています。C/C++でのプログラミングが面倒な理由の一つは必要なライブラリを揃えるところです。しかもただ揃っていれば良いわけでなく、バージョンも必要とされるものでなければなりません。面倒くさすぎですよね。多分Pythonの人気がある理由として、C/C++よりはライブラリに関する面倒くささが軽減されているからかもしれません。

 私がChemoinformatics関連のC/C++プログラミングで良く使うライブラリは以下の通りです

ライブラリ名説明
Boost便利なC++ライブラリ群
RDKitChemoinformatics用のツールキット
OpenBabelChemoinformatics用のツールキット
Eigen線形代数ライブラリ
BLAS, LAPACK線形代数ライブラリ
FFTW高速フーリエ変換
SQLite軽量データベース
TinyXML2XMLパーサー

 RDKitとOpenBabelの2つのライブラリは我々にとって必須です。これらが無いと全く仕事になりません。どちらも同じようなことが出来るライブラリですが、プログラミングという観点からするとOpenBabelの方が多機能でAPIも使い易いです。ただし、大きな弱点があります。それはスレッドセーフでは無いという点です。これはマルチスレッド化が当たり前になってきている現在では致命的です。一方RDKitはスレッドセーフです。しかも分子を扱うクラスのシリアライズが標準的に出来るのでMPIでの並列化もやりやすい気がします。でもちょっとAPIが使いづらいです。どちらも一長一短ですよね。

 他にも色々ライブラリを利用しています。特に計算屋であればBLASやLAPACKのような線型代数ライブラリは基本中の基本です。C++ではLAPACKよりもEigenが使い易いので、現在の私はこちらを良く使っています。大変オススメです。まあ今後このような話題についてもちょくちょく触れていければ、ブログのネタ切れを回避できるのではないでしょうか。

Ligand Expoに対するファーマコフォア探索

今回は、PARP-1のファーマコフォアを使って、Ligand Expoに登録されている化合物データに対して探索を行ってみます。

Step 1: Ligand Expoから化合物データの取得

ここから以下の2つのファイルを取得します。

Components-pub.sdf.gz:
PDBに登録されている低分子化合物のIDと構造情報
cc-to-pdb.tdd:
低分子化合物のIDとPDB IDの関係を示した情報

Components-pub.sdf.gzは解凍して前回の記事の方法(step 3)で配座生成を行います。

Step 2: PARP-1のファーマコフォアデータの保存と読み込み

前回の記事の方法(step 2)で作成したファーマコフォアにExclusion Volumes Coatを追加したファーマコフォアをPML形式で保存します。

Exclusion Volumes Coatの追加方法:
Structure-Based Perspecitveにおいて[Pharmacophore]→[Add Exclusion Volumes Coad]

Pharmacophoreの保存方法:
Structure-Based Perspectiveにおいて[File]→[Save as File]でpmz形式で保存する。ここではPARP1.pmzで保存します。

Pharmacophoreの読み込み方法:
Screening Perspectiveにおいて[File]→[Open]でPARP1.pmzを読み込みます。

Step 3: Ligand Expoに対するファーマコフォア探索

Screening PerspectiveにおいてComponents-pub.ldbを読み込んでファーマコフォア探索を実施します。Components-pub.ldbには約28,000の化合物が登録されています。

上図は、Relatieve Pharmacophore-Fit Scoreの最も良い化合物となります。

上位10化合物のうち6化合物がPARP-1の阻害剤でした。
今回利用したファーマコフォアは、なかなか選択性の高いものだと思います。
一方で、PARP-1阻害剤ではない上位の化合物は、擬陽性だと思われますが、PARP-1を阻害する可能性もあると思います。ドラッグリポジショニングの用途でも利用できるかもしれません。

Name列のIDが低分子化合物のIDとなります。PDBとの関係は、
cc-to-pdb.tddに記載されています。

LigandScoutでバーチャルスクリーニング

今回は、LigandScoutを用いたファーマコフォアの構築とバーチャルスクリーニングを行ってみます。

例として、DUD-EからPARP1の阻害剤とそのデコイを取得してバーチャルスクリーニングを実施します。DUD-Eは、ドッキングソフトウェアのベンチマークを行うためのデータベースであり、102の標的タンパク質に対して22,886の活性分子が登録されています。さらに1つの活性分子あたり50個のデコイ分子が登録されています。まずは、ここからactives_final.sdf.gzとdecoys_final.sdf.gzをダウンロードします。

Step 1: PARP1と阻害剤の複合体構造の取得

LigandScoutを起動し、PDB 4-letter codeである”3l3m”を入力して右側のボタンを押すとPDBからデータを取得することができます。

Step 2: ファーマコフォアの構築

阻害剤の周辺にある黄色の四角で囲まれている領域をクリックします。

Water(赤四角)を選択しdeleteキーを押すことで水分子を削除します。その後、ファーマコフォアを構築するためにpharmacophoreアイコン(青四角)を押します。

ファーマコフォアの構築ができました。Copy to other perspectivesの真ん中のボタンを押して、screening perspectiveを選択します。

Step 3: 活性分子とデコイ分子の配座生成

Screeningボタン(青四角)を押しScreening perspectiveに移動します。その後、右上のデータベースアイコン(赤四角)を押します。Create Screening Database windowが開きますので、Input FileをDUD-Eから取得したactives_final.sdfとし、Output Fileをactives_final.ldbとしてNextボタンを押します。本操作により配座生成を行うことができます。デコイ分子も同様の操作を行います。

Step 4: バーチャルスクリーニング

青四角内のアイコンを押して、Step 3で作成したactives_final.ldbとdecoys_final.ldbを入力します。actives_final.ldbの左側のマークを緑(活性分子)、そして、decoys.final.ldbの左側のマークを赤(デコイ分子)にします。

次に、赤四角内のアイコンを押すことでバーチャルスクリーニングがスタートします。スクリーニングが終わるまでしばらく待ちます。

スクリーニング終了後、緑四角内のアイコンを押すことでスクリーニング精度の評価をAUC、EFにより実施することができます。


残念ながらヒット化合物はゼロでした。さすがにファーマコフォアを構成する8個のFeatureを全て満たす分子は無いようです。

そこで、Hide advanced optionsを表示し、Max. number of omitted features (赤四角)を2とします。これで少なくとも6個のFeatureが満たされればヒットすることになります。

AUC(100%)が0.55となり若干改善されました。

最後にファーマコフォアの修正と設定の変更を行います。赤矢印で示す4つのFeatureをオプションとします([Pharmacophore]→[Mark Features as Optinal])。そして青矢印のFeatureをdeleteキーで削除します。

この操作により、少なくとも3つのFeatureは満たす必要があるが、残りのFeatureに対してはマッチすればスコアに加算されることになります。

PARP1は、NADをADPリボースとニコチンアミドに分解する酵素です。この3つのFeatureはニコチンアミドに相当する部分ですので、オプションにはしませんでした。ここで修正したファーマコフォアを用いて、これまでと同じ手順によりバーチャルスクリーニングを実施します。

AUC(100%)が0.84、EF(1%)が33.5と大幅に性能が改善されました。

スコアトップの化合物は活性化合物の1つである上図の化合物でした。基質の情報、無視できるFeature数、そしてオプション等を利用することで、単独の結晶構造から構築したファーマコフォアでも良い性能を出すことができます。