今回の話題は前回から引き続き類似構造検索（Similarity Search）についてです。タイトル通り、RDKitでC++プログラミングするという内容です。RDKitを取り上げるのは初めてで、ブログって何を書くか毎回悩むのですが、色々と忘れないうちに類似構造検索するというネタにしました。
はじめに
ケモインフォマティクスのツールキットとして良く利用されるのはRDKitであり、OpenBabelは少数派な気がします。RDKitはPythonと非常に相性が良いので、使いやすいのことが理由の一つに挙げられるでしょう。それは恐らく、RDKitはPythonを通して利用することを前提とした設計になっており、Pythonが昨今のAIブームで人気プログラミング言語になったこととも無関係では無いでしょう。ではなぜOpenBabelばかり取り上げてきたかというと、C++でRDKitを取り扱うのは少々面倒だからです。先ず、C++に関するドキュメントが少ないです。仕方ないからソースコードを読むのですが、RDKitの場合はソースコードが色々なディレクトリに散らばっていて、OpenBabelよりもクラスの使い方を調べるのに苦労します。そもそもどのようなヘッダーファイルをインクルードすれば良いか？そしてビルドするときにリンクすべきライブラリは何なのか？が分かりずらいのです。Pythonの場合はこのような苦労は一切なく、しかも日本語の情報も結構ネット上にあり、C++との差は歴然です。しかしながら、敢えてC++プログラミングに挑戦するのがこのブログなのです。
インストール
OpenBabel 3.1のインストールにAnacondaを利用しましたが、RDKitのインストールもAnacondaで簡単にできます。コマンドは以下のようになります。
$ conda install rdkit -c rdkit通常は仮想環境を作成してからインストールします。例えばOpenBabel 3.1のとき作成した”ob3″という仮想環境にインストールする場合は、以下のように一旦仮想環境を有効化してからインストールします。
$ conda activate ob3
$ conda install rdkit -c rdkitヘッダーファイルとライブラリはそれぞれ$HOME/anaconda3/envs/ob3/includeの下と$HOME/anaconda3/envs/ob3/libの下に格納されます。
今回のプログラム
今回のC++プログラムは前回と同じ類似構造検索計算をRDKitを利用して作成したものです。プログラム名は”SimSearch_RD.cpp“としています。入力データは”Fexofenadine.smi“と”test19.smi“です。それぞれのリンクをクリックしてダウンロードして下さい。
C++プログラムの内容は以下の通りです。
1 // C++
2 #include <iostream>
3 #include <fstream>
4 #include <sstream>
5 #include <string>
6 #include <vector>
7 #include <algorithm>
8
9 // C
10 #include <cstdio>
11 #include <cstdlib>
12 #include <ctime>
13 #include <cstring>
14 #include <cmath>
15 #include <unistd.h>
16
17 // RDKit
18 #include <GraphMol/RDKitBase.h>
19 #include <GraphMol/FileParsers/MolSupplier.h>
20 // RDKit FingerPrint
21 #include <GraphMol/Fingerprints/MorganFingerprints.h>
22
23 using namespace std;
24 using namespace RDKit;
25
26 int main(int argc, char* argv[]) {
27 ifstream isd1("Fexofenadine.smi");
28 MolSupplier* ms1 = new SmilesMolSupplier(&isd1, true, "\t", 0, 1, false, true);
29 ROMOL_SPTR ref = ROMOL_SPTR(ms1->next());
30
31 // calculate fingerprint
32 SparseIntVect<boost::uint32_t>* fp_ref = MorganFingerprints::getFingerprint(*ref, 2);
33
34 ifstream isd2("test19.smi");
35 MolSupplier* ms2 = new SmilesMolSupplier(&isd2, true, "\t", 0, 1, false, true);
36
37 while (true) {
38 if (ms2->atEnd()) break;
39 ROMOL_SPTR mol = ROMOL_SPTR(ms2->next());
40 if (!mol) continue;
41
42 SparseIntVect<boost::uint32_t>* fp_mol = MorganFingerprints::getFingerprint(*mol, 2);
43 double Tc = TanimotoSimilarity(*fp_ref, *fp_mol);
44 cout << mol->getProp<string>("_Name") << "\t" << Tc << endl;
45 }
46
47 return 1;
48 } いつものように簡単にプログラムの解説をします。先ずヘッダーファイルは17–21行目でインクルードしています。このプログラムでは最低限のヘッダーファイルしかありませんが、RDKitの利用したい機能のよって適切なヘッダーファイルを追加しないといけません。しかしながら、どこかのドキュメントに分かりやすくまとまっていないので、RDKitのソースコードを読みながら探すことになります。
SMILESファイルの読み込みについては27–29行目に記述しています。先ず27行目のように、普通にファイルをオープンします。次の28行目では、分子を読み込むためのクラスである”MolSupplier”のインスタンスを作成します。MolSupplierクラスを継承する形でSMILESを読み込むクラスであるSmilesMolSupplierが用意されています。SDファイルの場合はSDMolSupplierクラスを利用して、PDBの場合はPDBMolSupplierクラスを利用します。最後に29行目でrefという変数に分子（Fexofenadine）情報を格納しています。ここでROMOL_SPTRという型が出てきています。RDKitにはROMolという分子用のクラス（OpenBabelでいうところのOBMolのことです）があり、このクラスのスマートポインタを示す型としてROMOL_SPTRが定義されています。自分でメモリ管理する自信があるならROMol* ref=ms1->next()としても良いかと思いますが、なるべくスマートポインタを利用しましょう。
フィンガープリントの計算は32行目で記述しています。フィンガープリントを格納する変数型は”SparseIntVect”です。”MorganFingerprints::getFingerprint”という関数でフィンガープリントを計算して、結果をポインタとして返しています。そして43行目で、フィンガープリントから谷本係数を計算しています。Similarity metricについては、他にもDiceSimilarityとTverskySimilarityが利用できます。
コンパイルと実行
このプログラムのコンパイルは以下のようなコマンドで実行できます。
g++ -std=c++11 -w SimSearch_RD.cpp -I ${HOME}/anaconda3/envs/ob3/include/rdkit \
-L ${HOME}/anaconda3/envs/ob3/lib \
-Wl,-rpath=${HOME}/anaconda3/envs/ob3/lib \
-lRDKitFileParsers -lRDKitRDGeneral -lRDKitFingerprints注意してほしいのはg++(gcc)のバージョンです。CentOS7で標準的にインストールされるgccのバージョンではコンパイル（リンク）できません。もっと新しいバージョンを利用してください。ちなみに私が使用しているgccのバージョンは9.3.0です。Ubuntu 20.04 LTSで普通にインストールされるものです。
コンパイルして出来たロードモジュールを実行すると次のようになります。
$ ./a.out
101845628 0.267974
123506854 0.735537
139808122 0.486301
141260304 0.827586
14677762 0.562044
150332786 0.740458
150598871 0.692857
150622200 0.773109
150794865 0.736842
22667272 0.47619
3348 1
57619841 0.177033
57619845 0.21
57619846 0.164319
67722853 0.653543
67767321 0.620155
68320883 0.521739
71141694 0.934579
9936344 0.884956ちなみに
上記のC++プログラムをPythonで作成したら以下のようになります。
1 from rdkit import Chem
2 from rdkit.Chem import AllChem
3 from rdkit import DataStructs
4
5 suppl1 = Chem.SmilesMolSupplier("Fexofenadine.smi", delimiter="\t", titleLine=False)
6 ref = next(suppl1)
7 fp_ref = AllChem.GetMorganFingerprint(ref, 2)
8
9 suppl2 = Chem.SmilesMolSupplier("test19.smi", delimiter="\t", titleLine=False)
10 for mol in suppl2:
11 fp_mol = AllChem.GetMorganFingerprint(mol, 2)
12 Tc = DataStructs.TanimotoSimilarity(fp_ref, fp_mol)
13 print("%s\t%f" % (mol.GetProp("_Name"), Tc))単純に行数だけを比較しても3倍以上短いです。しかも、コンパイルの必要もありません。適当に”SimSearch.py”というファイル名で上記のプログラムを保存したとして、それを実行するには以下のコマンドだけです。
$ python ./SimSearch.py最後に
やっぱりRDKitはPythonで利用すべきです。私ももっとPythonを勉強したいと思います。
Category: Linux関連, RDKit, プログラミング関連