C++ による RDKit 入門

RDKit は Python/C++ ベースのケムインフォマティクスソフトウェア環境です．
このページは，分子を構築したり構築された分子からいくつかの情報を取得するごく基本的な方法を学んだ結果を再構成したものです．
教材は，RDKit のドキュメント "Getting Started with the RDKit in C++" およびサンプルプログラム（17 あるサンプルプログラムの最初の 3 個分）です．
ご存知の方向けに大雑把に書いておくと，SMILES 形式や MDL MOL 形式から RDKit::ROMol クラスや RDKit::RWMol クラスのインスタンスを作成します．

当面は Python を使って勉強し，その上で，当サイトで開発している C++ ソフトウェアに組み込むべきテーマを選択します（2023 年 12 月 21 日）．

インフォメーション

ウェブサイト

The RDKit Documentation：RDKit のサイト直下にあるドキュメントのトップページです．

Getting Started with the RDKit in C++：C++ で始めるためのドキュメントです．

RDKit: Main Page：C++ API のページです．

rdkit/Docs/Book/C++Examples：GitHub にあるサンプルプログラムです．example1〜17 があります．

目次（ページ内リンク）

考慮すべき点

例１ RDKit::ROMol クラスを使った分子の構築

例２ RDKit::MolFileToMol 関数を使った分子の読み込み

例3 RDKit::SDMolSupplier クラスを使った複数の分子の読み込み

例4 表記法のテスト

インストール

使用している主な Debian パッケージを示します．
librdkit1 は，librdkit-dev をインストールしたら自動的にインストールされます．

パッケージ名(バージョン)	注
g++ (12.2.0)	GNU C++ コンパイラ
librdkit-dev (202209.3)	RDKit（ケムインフォマティクスソフトウェアコレクション）の開発ファイル
librdkit1 (202209.3)	RDKit（ケムインフォマティクスソフトウェアコレクション）の共有ライブラリ

考慮すべき点

Debian パッケージのインストール先 (1)
librdkit-dev のヘッダファイルは，/usr/include/rdkit/ 以下にインストールされます．コンパイルの I オプションではここを指定すればよさそうです．

Debian パッケージのインストール先 (2)
librdkit のライブラリは，/usr/lib/libRDKit*** というファイル名でインストールされます．ビルドの l オプションでライブラリの名前がわからなくなったらここを探します．
librdkit1 (202209.3)なら "ls /usr/lib/libRDKit*.2022.09.3" とでもすれば，何というライブラリがインストールされているか判ります．

以下のパラグラフは，"Getting Started with the RDKit in C++" 冒頭の内容の一部を改変して引用しています．

Python を使うか C++ を使うか
独自のプログラミングロジックを多用し，I/O や SMARTS マッチング，2D 画像の準備などの周辺処理にのみ使用するような，より複雑なことをおこなうのであれば，C++ の方がよりパフォーマンスが向上する可能性が高いでしょう．

Python より C++ の方が API は不完全です
RDKit の高度な機能の多くは Python で開発されており，C++ へのバックポートは需要に応じておこなわれます．そのため，コンフォーマー間の RMS 差の計算（GetConformerRMSMatrix）など，C++ では利用できない機能もあります．

コンパイル
2018 年 3 月リリースの時点で，RDKit はモダン C++ を使用しており，これはリリース時点で C++17 までを意味します．Clang と gcc では，フラグ -std=c++17 を使用してください．

メモリリーク
RDKit では，多くの関数が分子へのポインタを返し，そのようなポインタを引数として受け取ります．std::shared_ptr や std::unique_ptr といったスマートポインタがメモリリークを防ぐのに役立ちます．

例１ RDKit::ROMol クラスを使った分子の構築

最初のプログラムでは，プログラム内で分子を構築し，構築された分子から情報を取得します．
RDKit には C++ で分子を表現するためのクラスが 2 個あります．RDKit::ROMol と RDKit::RWMol です。
ROMol は読み取り専用分子です（Read-Only）．分子を編集する必要がある場合は RWMol（Read-Write）を使います．
どちらも GraphMol.h で宣言されています。

ソースファイル

下のソースコードは，冒頭で紹介した GitHub にあるサンプルプログラムの "example1.cpp" の冒頭を改変したものです．
ページ冒頭の「考慮すべき点」で指摘されているように，生ポインタではなく共有ポインタにしました．
ここでは，これを作業ディレクトリ ~/wk に exp1.cc という名前で保存します．


#include <iostream>
#include <GraphMol/GraphMol.h>
#include <GraphMol/SmilesParse/SmilesParse.h>

int main(int argc, char **argv) {
 std::shared_ptr<RDKit::ROMol> mol1(RDKit::SmilesToMol("Cc1ccccc1"));  //分子の構築
 std::cout << "Number of atoms " << mol1->getNumAtoms() << std::endl;  //分子から情報所得

 return 0;
}

冒頭に示した C++ API のページの RDKit::ROMol Class Reference を見ると理解が深まります．
例えば，

SmilesToMol()：RDKit 名前空間の公開関数．Smiles 形式の文字列（ここでは "Cc1ccccc1"）を引数にして RWMol へのポインタを返す
getNumAtoms()： RDKit::ROMol クラスの公開関数．bool 型引数を取ることができて，true なら水素以外の原子数を，false なら水素を含めた原子数を unsigned int 型で返す
GraphMol.h は，RDGeneral/export.h，ROMol.h，および RWMol.h をインクルードしているだけのファイル

ビルド

コンパイル

~/wk$ g++ exp1.cc -c -I/usr/include/rdkit -std=c++17
でコンパイルできました．「考慮すべき点」で記したように，

-I/usr/include/rdkit でヘッダファイルがインストールされたディレクトリを指定します
C++17 を使いました

リンク

利用しているクラス名から ~/wk$ g++ exp1.o -o exp1 -lRDKitSmilesParse
でリンクできました（-lRDKitGraphMol は無くてもリンクできました）．
ここでは exp1 という実行ファイルを生成しています．

-l オプションで指定しているライブラリ名は，「考慮すべき点」で触たようにヘッダファイルから推定できそうです．

実行

実行ファイルのあるディレクトリで
~/wk$ ./exp1
とすると，下のように出力されました．
Smiles 形式から分子が構築できたようです．

Number of atoms 7

例２ RDKit::MolFileToMol 関数を使った分子の読み込み

このセクションでは，MDL MOL 形式の分子ファイルに記述された情報で RDKit::ROMol インスタンスを構築する練習をします．
ここでは，RDKit::MolFileToMol 関数の戻り値を RDKit::ROMol のコンストラクタ引数とします．

分子ファイル

読み込みに使ったファイルを示します．エタンを Mol 形式で記述したファイルです．
MOL 形式には，通常 1 個の分子が記述されます．
Builcule で作成しました．
ファイルは，作業ディレクトリ ~/wk に ethane.mol という名前で保存してあるとします．



 OpenBabel12152307473D

  8  7  0  0  0  0  0  0  0  0999 V2000
   -0.6290   -0.4450    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6290   -1.1030   -0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6290   -1.1030    0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6290    0.4450    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.6290    1.1030   -0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6290    1.1030    0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5600    0.2140    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.5600   -0.2140    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
  1  4  1  0  0  0  0
  1  7  1  0  0  0  0
  4  5  1  0  0  0  0
  4  6  1  0  0  0  0
  4  8  1  0  0  0  0
M  END

ソースファイル

例１と同じくサンプルプログラムの "example1.cpp" を改変したものです．
ここでは，これを作業ディレクトリ ~/wk に exp2.cc という名前で保存します．


#include <iostream>
#include <GraphMol/GraphMol.h>
#include <GraphMol/FileParsers/FileParsers.h>

int main(int argc, char **argv) {
 std::string mol_file = "./ethane.mol";  //ファイルパス
 std::shared_ptr<RDKit::ROMol> mol2(RDKit::MolFileToMol(mol_file));  //分子の構築

 //分子情報の出力
 std::cout << "Number of atoms(atomic number > 1) " << mol2->getNumAtoms(true) << std::endl;  //水素以外の原子数
 std::cout << "Number of atoms " << mol2->getNumAtoms(false) << std::endl;  //水素を含めた原子数
 std::cout << "Number of heavy atoms " << mol2->getNumHeavyAtoms() << std::endl; //水素以外の原子数
 std::cout << "Number of bonds(atomic number > 1) " << mol2->getNumBonds(true) << std::endl;  //水素が関与しない結合数
 std::cout << "Number of bonds " << mol2->getNumBonds(false) << std::endl;  //水素の関与も含めた総結合数

 return 0;
}

"example1.cpp" では，環境変数 "RDBASE" が設定してあって，RDBASE を起点にファイルパスを設定しています．
ここでは簡便法として，環境変数は設定せず，ファイルパスを相対パスで記述しています．

C++ API の RDKit::ROMol Class Reference で目にした関数を使ってみました．

getNumHeavyAtoms()
getNumBonds(bool onlyHeavy=1)

C++ API によると，MolFileToMol() 関数は，RDKit 名前空間にある公開関数です．
関数の説明箇所を引用します．どの程度 sanitize されるか，どの程度 strict なのか，執筆時点では解りません．

RDKIT_FILEPARSERS_EXPORT RWMol * RDKit::MolFileToMol  (
  const std::string &   fName,
  bool   sanitize = true,
  bool   removeHs = true,
  bool   strictParsing = true 
 )

Parameters
    fName - string containing the file name
    sanitize - toggles sanitization and stereochemistry perception of the molecule
    removeHs - toggles removal of Hs from the molecule. H removal is only done if the molecule is sanitized
    strictParsing - if set to false, the parser is more lax about correctness of the contents.

ビルド

今回は -lRDKitGraphMol が必要でした．

~/wk$ g++ exp2.cc -c -I/usr/include/rdkit -std=c++17
~/wk$ g++ exp2.o -o exp2 -lRDKitGraphMol -lRDKitFileParsers

実行

~/wk$ ./exp2
とすると，下のように出力されました．
ファイルが正しく読み込まれたようです．

Number of atoms(atomic number > 1) 2
Number of atoms 8
Number of heavy atoms 2
Number of bonds(atomic number > 1) 1
Number of bonds 7

例3 RDKit::SDMolSupplier クラスを使った複数の分子の読み込み

このセクションでは，SDF 形式の分子ファイルを開く練習をします．
SDF 形式は，複数の MOL 形式の分子を一つのファイルにまとめたものです．
ここでは，ファイルパスを引数にして RDKit::SDMolSupplier インスタンスを作成し，個々の分子情報を while ループで取得します．

分子ファイル

下が練習に使った SDF ファイルです．
アンモニアと水分子を別々に MOL 形式で作成し，テキストエディタで 1 ファイルにまとめたものです．
複数の分子を区別するために，"$$$$" を区切り行とします．
ファイルは，作業ディレクトリ ~/wk に mols.sdf という名前で保存してあるとします．


ammonia
 OpenBabel12152314473D

  4  3  0  0  0  0  0  0  0  0999 V2000
    0.0000    0.0000    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000   -0.6580   -0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000   -0.6580    0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.9310    0.6580    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
  1  4  1  0  0  0  0
M  END
$$$$
H2O
 OpenBabel12152314523D

  3  2  0  0  0  0  0  0  0  0999 V2000
    0.3100    0.0000    0.3100 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.3100   -0.6580   -0.6210 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6210    0.6580    0.3100 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
M  END

ソースファイル

サンプルプログラムの "example2.cpp" の一部を改変したものです．
RDKit::SDMolSupplier というクラスを使っています．
ここでは，これを作業ディレクトリ ~/wk に exp3.cc という名前で保存します．


#include <iostream>
#include <GraphMol/GraphMol.h>
#include <GraphMol/FileParsers/MolSupplier.h>

int main(int argc, char **argv) {
 RDKit::SDMolSupplier mol_supplier("./mols.sdf");  //引数は SDF ファイルのパス

 //mol_supplier から 1 個ずつ分子を取得して原子数を表示
 std::unique_ptr<RDKit::ROMol> mol;
 while (!mol_supplier.atEnd()) {
  mol.reset(mol_supplier.next());
  std::cout << mol->getProp<std::string>("_Name") << " has "
    << mol->getNumAtoms(false) << " atoms." << std::endl;
 }

 return 0;
}

RDKit::SDMolSupplier には 3 種類のコンストラクタが記されていました．
引数なし，const std::string & を引数とするもの，および std::istream * を引数とするものです．
ここでは const std::string & 引数を使っています．

また，RDKit::SDMolSupplier には RDKit::ROMol * を返す関数が 2 個書いてありました．
ここでは ROMol * next () を利用しています．
もう一つは，ROMol * operator[] (unsigned int idx) です．配列でおなじみの [idx] の形で使います．

ビルド

~/wk$ g++ exp3.cc -c -I/usr/include/rdkit -std=c++17
~/wk$ g++ exp3.o -o exp3 -lRDKitGraphMol -lRDKitFileParsers -lRDKitRDGeneral

実行

~/wk$ ./exp3
とすると，下のように出力されました．
ファイルが正しく読み込まれたようです．

ammonia has 4 atoms.
H2O has 3 atoms.

例4 表記法のテスト

サンプルプログラムの "example3.cpp" では，いくつかの表記法と相互変換が例示されています．
このセクションでは，それらを抜粋して学びます．

キラルな分子の例

"example3.cpp" では，キラルな分子の読み取りと出力もしています．
ここでは，Builcule で MOL 形式の L-Ala を作成し，使ってみます．

L-Ala
 OpenBabel12162317323D

 13 12  0  0  1  0  0  0  0  0999 V2000
   -1.4260   -0.8400   -0.1880 N   0  0  0  0  0  0  0  0  0  0  0  0
   -0.1930   -0.0410   -0.3070 C   0  0  2  0  0  0  0  0  0  0  0  0
    1.0500   -0.9380   -0.1120 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.2240   -0.2590   -0.1170 O   0  0  0  0  0  0  0  0  0  0  0  0
    1.0240   -2.1590   -0.0600 O   0  0  0  0  0  0  0  0  0  0  0  0
   -0.2750    1.0870    0.7140 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.3030   -1.7150   -0.7030 H   0  0  0  0  0  0  0  0  0  0  0  0
   -2.2120   -0.3490   -0.6070 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.1370    0.3630   -1.3250 H   0  0  0  0  0  0  0  0  0  0  0  0
    2.1420    0.7040   -0.2380 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6040    1.7370    0.6770 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.3460    0.6950    1.7360 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1520    1.7170    0.5300 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  7  1  0  0  0  0
  1  8  1  0  0  0  0
  2  3  1  0  0  0  0
  2  6  1  0  0  0  0
  2  9  1  6  0  0  0
  3  4  1  0  0  0  0
  3  5  2  0  0  0  0
  4 10  1  0  0  0  0
  6 11  1  0  0  0  0
  6 12  1  0  0  0  0
  6 13  1  0  0  0  0
M  END

ソースコード

下にコードを示します．
前半ではキラル分子を開いて出力するだけなので，RDKit::ROMol を使っています．
後半では，ケクレ化と名前の設定という編集をおこなうので，RDKit::RWMol を使っています．

#include <iostream>
#include <GraphMol/FileParsers/FileParsers.h>
#include <GraphMol/SmilesParse/SmilesWrite.h>
#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/MolOps.h>

int main(int argc, char **argv) {
 //ファイルを開く．mol0⇔アラニン
 std::string mol_file = "./ala.mol";
 std::shared_ptr<RDKit::ROMol> mol0(RDKit::MolFileToMol(mol_file));
 std::cout << "キラリティの表示/非表示" << std::endl;
 //第二引数を true にすると
 std::cout << RDKit::MolToSmiles(*mol0, true) << std::endl;
 //第二引数を false にすると
 std::cout << RDKit::MolToSmiles(*mol0, false) << std::endl;
 std::cout << std::endl;

 std::cout << "芳香環の表記" << std::endl;
 //大文字小文字で記述．mol1⇔芳香環のピリジン
 std::shared_ptr<RDKit::RWMol> mol1(RDKit::SmilesToMol("c1cccnc1"));
 std::cout << RDKit::MolToSmiles(*mol1) << std::endl;
 //二重結合を記述．mol1⇔ケクレ化したピリジン
 std::shared_ptr<RDKit::RWMol> mol2(RDKit::SmilesToMol("C1=CC=CN=C1"));
 std::cout << RDKit::MolToSmiles(*mol2) << std::endl;
 RDKit::MolOps::Kekulize(*mol2);
 std::cout << RDKit::MolToSmiles(*mol2) << std::endl;
 std::cout << std::endl;

 std::cout << "ブロック化" << std::endl;
 mol2->setProp("_Name", "ケクレ表記のピリジン");
 std::cout << RDKit::MolToMolBlock(*mol2) << std::endl;

 return 0;
}

ビルド

~/wk$ g++ exp4.cc -c -I/usr/include/rdkit -std=c++17
~/wk$ g++ exp4.o -o exp4 -lRDKitSmilesParse -lRDKitFileParsers -lRDKitGraphMol -lRDKitRDGeneral

実行結果

~/wk$ ./exp4
とすると，下のように出力されました．

キラリティの表示/非表示
C[C@H](N)C(=O)O
CC(N)C(=O)O

芳香環の表記
c1ccncc1
c1ccncc1
C1=CC=NC=C1

ブロック化
ケクレ表記のピリジン
     RDKit          2D

  6  6  0  0  0  0  0  0  0  0999 V2000
    1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500    1.2990    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0
  2  3  1  0
  3  4  2  0
  4  5  1  0
  5  6  2  0
  6  1  1  0
M  END

参考書の検索

Amazon の「コンピュータ・IT」本カテゴリーでの，「python」での検索結果です
Amazon の「科学・テクノロジー」本カテゴリーでの，「ケモインフォマティクス」の検索結果です
Amazon の「科学・テクノロジー」本カテゴリーでの，「有機化学」の検索結果です