OpenBabelで分子の配座探索

OpenBabelでエネルギー計算ができるようになって、さらに便利になりました。そこで、超シンプルな配座探索プログラム(obconfopt)を作ってみました。

基本的には、ランダムに配座を発生させ、最安定な配座を探すというランダム探索法を採用しています。

　Step 1. 初期配座の入力
　↓
　Step 2. torsion角をランダムに設定し、新しい配座の生成。
　↓
　Step 3. Energy Minimization
　↓
　Step 4. 前に発生させた配座と異なっていれば保持する。
　↓
　Step 5. 終了判定

終了するまでは、Step2からStep5を繰り返します。
ここで、Step2はobrotateを、Step3はobenergyとobminimizeを、そしてStep4はobfitを参考にすることができます。簡単ですので、興味のある方は、作ってみてください。

例としてペンタンを使ってテストしました。


実用レベルで使うには、環構造を考慮したり、ヒューリスティックなアルゴリズムを採用したりとまだまだ改良する部分は多いですが、全くのゼロからobconfoptに相当するプログラムを書くのは、相当大変ですので、OpenBablの便利さを再認識しました。

OBAtomIteratorの使い方

OpenBabelのOBAtomIeratorを使って、分子の中心座標を計算してみます。

vector3 calcMass(OBMol& mol){

　vector3 center(0,0,0);
　OBAtomIterator i;
　OBAtom *atom;

　for(atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)){
　　center += atom->GetVector();
　}

　center /= (float) mol.NumAtoms();

　return center;
}

マクロ定義がされていますので、for文を次のように書くこともできます。

　FOR_ATOMS_OF_MOL(atom,mol)
　{
　　center += atom->GetVector();
　}

ライブラリによって、原子が配列の0から格納されていたり、1から格納されていたりするので、Iteratorを使うことによって、コーディングミスを減らすことができると思います。

logPとTPSAの計算

OpenBabel2.1.0になって、エネルギー計算を筆頭に有益な機能がいくつか追加されているようです。

今回は、logPとTPSAを計算してみます。

ob_exp1.cpp


私の環境では、コンパイルを次のようにして行いました。
　
　$ g++ -I/usr/local/include/openbabel-2.0 -L/usr/local/lib ob_exp1.cpp -lopenbabel

実行

　$ ./a.out
LogP 1.3922
TPSA 20.23

あと、PDBをfasta形式で出力できるようになってました。

　$ babel -ipdb input.pdb -ofsa output.fsa

単純な変換ですが、便利だと思います。

Open Babel GUI

OpenBabelのWindows GUI版を使ってみました。なかなか使い勝手もよく便利だと思います。そういえば、BabelWinというGUI版のbabelもありましたね。私はコマンドライン版に慣れきってしまっているため使うことはほとんどないですが、GUI版は、バイナリで配布されているため、インストール作業も必要なく便利なのは確かだと思います。





人気ブログランキング(クリックして応援してね)

PythonでOpenBabel その2

PythonからOpenBabelのクラスを利用して、簡易版obgrepを作成したいと思います。

simple_obgrep.pyはsdf形式のファイルを入力し、部分構造としてベンゼン(c1ccccc1)を含んでいる分子を検出するシンプルなプログラムです。

実行：
　$ python simple_obgrep.py > output.sdf

simple_obgrep.py:




人気ブログランキング(クリックして応援してね)

OpenBabelでsimilarity search

OpenBabelでは、3種類のfingerprintを利用してsimilarity searchを行うことができます。fingerprintの比較はTanimoto coefficientで行われています。

それでは、実際にsimilarity searchをしてみます。

　$ babel 'query.smi' 'mydata.sdf' -ofpt -xfFP2

　結果：
　> Tanimoto from first mol = 0.166667
　> Tanimoto from first mol = 0.285714
　> Tanimoto from first mol = 0.27027
　・
　・
　> Tanimoto from first mol = 0.0802469
　> Tanimoto from first mol = 0.149425
　> Tanimoto from first mol = 0.122449
　39 molecules converted

ここでは、query.smiとmydata.sdfに含まれる全ての分子との間のTanimoto coefficientが計算されています。-xfFP2はfingerprintにFP2を利用することを示しています。fingerprintの種類は以下のコマンドで知ることができます。

　> babel -F
　FP2 -- Indexes linear fragments up to 7 atoms.
　FP3 -- SMARTS patterns specified in the file patterns.txt
　FP4 -- SMARTS patterns specified in the file SMARTS_InteLigand.txt

patterns.txtとSMARTS_InteLigand.txtは、通常、/usr/local/share/openbabel/2.0.2などにインストールされていると思います。このファイルに独自にpatternを追加すればオリジナルのfingerprintで探索ができます。



人気ブログランキング(クリックして応援してね)

PythonでOpenBabel

OpenBabelではSWIGを用いて、そのほとんどクラスをPythonから利用できるようにしています(Perlからも利用できるみたいです)。

OpenBabelを普通にインストールしただけではPythonから利用できないので、ちょっとしたインストール作業が必要です。

1. openbabel-2.0.1.tar.gzを展開したディレクトリに移動します。

　＞cd $OpenBabelSrc/openbabel-2.0.2/scripts/python

2. ビルド、インストール作業を行います。

　＞python setup.py build
　＞python setup.py install

これで、PythonからOpenBabelのクラスが利用できるようになります。

3. 簡単なサンプルプログラムを動かしてみます。

　sample.py：
　import openbabel

　mol = openbabel.OBMol()
　conv = openbabel.OBConversion()
　conv.SetInAndOutFormats("smi","mdl")
　conv.ReadString(mol,"c1ccccc1O")
　print mol.NumAtoms()
　mol.AddHydrogens()
　print mol.NumAtoms()
　conv.WriteFile(mol,'phenol.mol')

実行:
　$python sample.py

出力：
　7　<-水素なしの時の原子数
　13　<-水素ありの時の原子数

出力ファイル(phenol.mol):

OpenBabel

13 13 0 0 0 0 0 0 0 0999 V2000
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 C 0 0 0 0 0
0.0000 0.0000 0.0000 O 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
0.0000 0.0000 0.0000 H 0 0 0 0 0
1 6 2 0 0 0
1 2 1 0 0 0
1 8 1 0 0 0
2 3 2 0 0 0
2 9 1 0 0 0
3 4 1 0 0 0
3 10 1 0 0 0
4 5 2 0 0 0
4 11 1 0 0 0
5 6 1 0 0 0
5 12 1 0 0 0
6 7 1 0 0 0
7 13 1 0 0 0
M END

OpenBabelでは、まだ2次元座標の生成機能がないこともあり、xyzが全て0となります。
ScriptからOpenBabelを利用できるとコンパイル作業もなく便利ですね。



人気ブログランキング(クリックして応援してね)

iBabel

iBabelというMac OS/Xで動作するOpenBabelのGUIがあるそうです。私の周りには、残念ながらMacがないので動作確認ができないのですが、ファイル変換だけでなく、部分構造検索、物性計算、さらには、Jmol, JChemPaint, Marvinを使った分子構造の入力/閲覧もできるみたいです。また、iBabelは、AppleScriptというスクリプト言語で実装されているようで、これもおもしろそう。

既存のソフトウェアのシンプルな組合せですが、簡単にsdfの中身を覗くことができ、便利なのではないでしょうか。



人気ブログランキング(クリックして応援してね)

OpenBabelの小技

mol/sdf形式中のbond typeは、次の様に定義されています。

1=Single, 2=Double,
3=Triple, 4=Aromatic,
5=Single or Double,
6=Single or Aromatic,
7=Double or Aromatic,
8=Any

例えば、ナフタレンの場合、全ての結合は4としてもOKとなります。したがって、4と指定されたナフタレンのmol形式ファイル(naft.mol)を分子表示ソフトで見ると、下図左のような表記になります。しかしながら,下図右のような表記を希望する方も多いと思います。このような場合、4のbond typeを適切に2と1のbond typeへと変換する必要がありますが、当然、手作業でやりたくありません。
Openbabelでは、単純にmolからmolへ変換することにより、上記作業を実現することができます。

　$ babel -imol naft.mol -omol naft2.mol



ネット上で公開されている化合物のsdf形式は、bond typeが4と指定されているものが多くありますので、この小技を使えばbond typeを1,2に戻せると思います(ただし、全てのケースにおいてうまく戻せるのかは分かりません)。それなりに使える機会があるのではないでしょうか。



人気ブログランキング(クリックして応援してね)

OELibでねじれ角の計算

6月4日の記事で、シクロヘキサンのねじれ角(C1-C2-C3-C4)の計算にOELibで作成したtorsionを利用しました。
そこで、今回は、OELibを用いた分子の入力からねじれ角の計算までをtorsion.cppを通して紹介したいと思います。

torsion.cppのソースコードを下のインラインフレーム内に記載しております。まずは、OELibのヘッダファイルであるopenbabel/mol.hとopenbabel/obconversion.hをincludeし、namespaceを定義します。ちなみに、torsionのUsageは以下のとおりです。
Usage: torsion filename atom1 atom2 atom3 atom4
ここで、filenameは分子ファイル名を示し、atomXは、ねじれ角を構成する原子のIDを示します。
次に引数チェックのあと、分子ファイルの拡張子からファイル形式を判定しています。分子の読み込みは、ファイルから1分子づつfor文の無限ループ内で行っており、mol.Empty()がtrueとなればbreakします。プログラムのメインとなるねじれ角の計算はOBMolクラスのGetTorsionを用いて計算しています。
今回定義したシクロヘキサンのねじれ角は下図のようになります。


torsion.cpp


コンパイル：
$ g++ -I/usr/local/include/openbabel-2.0 -L/usr/local/lib -o torsion　torsion.cpp -llibopenbabel

実行：
$ torsion cyclo.sdf 1 2 3 4

出力：
conf-1 torsion_angle 304.814
conf-2 torsion_angle 291.955
..
..
conf-100 torsion_angle 329.156

無事にねじれ角が計算できると思います。




人気ブログランキング(クリックして応援してね)

バベる!5世：obrotate

obrotateは任意の二面角を任意の角度に回転させるプログラムです。
obfitと同様にまずSMARTパターンに一致する部分構造を検索し、
回転の軸とする原子を決定します。

例としてブタンの二面角を0度から360度まで30度刻みで
回転させたいと思います。

$ obrotate "CCCC" butane.mol 1 2 3 4 0 > butane0.mol
$ obrotate "CCCC" butane.mol 1 2 3 4 30 > butane30.mol
$ obrotate "CCCC" butane.mol 1 2 3 4 60 > butane60.mol
....
....
$ obrotate "CCCC" butane.mol 1 2 3 4 360 > butane360.mol

ここで、butane.molは水素原子の付加されたbutaneのmolファイルであり、"CCCC"というSMARTパターンでbutane.mol内を検索することを指示しています。1 2 3 4はbutane.mol内の原子IDを示しているのではなく、SMARTパターン上の原子を示しています。最後の引数が指定した二面角の角度となります。

生成されたmolファイルからアニメーションgifを作ってみました。
指定どおり回転していることが分かります。


任意の二面角の回転は、化合物のコンフォメーション解析やドッキングなど利用できる場面が多くあると思います。obrotateはOELibを使ってシンプルに実装されていますので、興味のある方はソースコードをチェックしてみてください。



人気ブログランキング(クリックして応援してね)

バベる!4世:obfit

obfitは、分子の重ね合わせを行うことができます。特徴は、任意に指定したSMARTパターンに一致する部分構造をベースとして、重ね合わせを行うところにあります。もっとも単純な実装の1つとして、各分子ごとに3つの原子を任意に指定し、これをベースに重ね合わせを行うというものが挙げられますが、自動化には不向きです。obfitのように部分構造検索と組み合わせれば、複数分子の重ね合わせを瞬時に得ることができ便利ですね。

obfitの使用例として、下の化合物Aに化合物Bを重ね合わせたいと思います。ここで、重ね合わせに用いるSMARTパターンは"n1c(ccn1C)"とします。



$ obfit "n1c(ccn1C)" compA.mol compB.mol > compBm.mol

compBm.molは、compA.molに重ね合わせを行ったcompB.molの座標をもっています。
次に、重ね合わせの結果を見てみます。



おお！見事に重なっていますね。

今回は、SMARTパターンとして、化合物Aと化合物Bの最大共通部分構造(Maximum Common Substructure:MCS)を用いています(目視で選択)。OpenBabelでは、まだMCSは実装されていないと思いますが、そのうち実装されると期待しています。分子の重ね合わせは、ファーマコフォアの推定や3D-QSARの前処理に必須ですので、利用価値が高いと思います。



人気ブログランキング(クリックして応援してね)

バベる!3世：obprop

あのダビスタの園部さんの作品である「カルチョビット」期待もこめて買いました。おもしろい！久しぶりに楽しめるゲームに出会えてラッキーです。

さて、今日はobpropについてメモしたいと思います。obpropは分子量と分子に含まれる環構造の数を計算するプログラムです。

$ obprop compA.mol
name a1
mol_weight 172.226
num_rings 2

ここでnum_ringsは、Smallest Set of Smallest Rings (SSSR)検出アルゴリズムにより算出されています。SSSRは芳香族環の検出の前処理など多くの場面で使う機会があると思います。シンプルなプログラムですが、それなりに使う場面があるかも？しかし、せめてRule of 5の判定ができる程度のpropertyが算出できるとうれしいですね。自分で作れということでしょうか。



人気ブログランキング(クリックして応援してね)

バベる!2世: obgrep

OpenBabelをインストールするとbabel以外にも以下の5つのプログラムがインストールされます。

obgrep : SMARTSを使った部分構造検索
obfit : 分子の重ね合わせ
obrotate : torsion angleの回転
obprop : 分子の物性計算
obchiral : キラル中心の検索

これらはどれもOpenBabelのAPIであるOELibを利用して作られたプログラムです。上記プログラムを全て自力で作ろうと思えば、それなりに大変な内容なのですが、OELibを使うと非常にシンプルに実装できることが分かると思います。また、OELibを使いこなすよい参考にもなります。

今回はobgrepについてメモしたいと思います。
上記ツール郡は、openbabel-2.0.1/tools以下にソースコードが置かれており、まずは自力でコンパイルしたいと思います。Makefileはこの場所に置かれていますので、これを使ってもいいのですが(見た目にちょっと複雑ですね)、よりシンプルな方法でコンパイルしたいと思います。

obgrep.cppを作業ディレクトリにコピーする。
$ cp (展開したディレクトリ)/openbabel-2.0.1/tools/obgrep.cpp .

obgrep.cppの以下の3行を変更する。
(変更前)
#include "mol.h"
#include "obconversion.h"
#include "parsmart.h"
(変更後)
#include "openbabel/mol.h"
#include "openbabel/obconversion.h"
#include "openbabel/parsmart.h"

コンパイルの実行
$ g++ -I/usr/local/include/openbabel-2.0 -L/usr/local/lib -o obgrep obgrep.cpp -llibopenbabel

無事にobgrep.exe(Cygwin上なので)ができました。
obgrepの使用法は以下のとおりです。

Usage: ./obgrep [options] "PATTERN"
Options:
-v Invert the matching, print non-matching molecules
-c Print the number of matched molecules
-f Full match, print matching-molecules when the number
of heavy atoms is equal to the number of PATTERN atoms
-n Only print the name of the molecules
-t NUM Print a molecule only if the PATTERN occurs NUM times inside the molecule.

例）1000化合物を含むdata.smiからフェノール基をもつ化合物の数を調べる。
$ ./obgrep -c "c1ccccc1[OH]" data.smi
67

67個の化合物が含まれていることが分かりました。検索にかかった時間はPentium4(2.6GHz,512MB)、Cygwin上で約20secです。研究室内の小規模な化学構造検索システムであればobgrepを利用して構築できるかもしれません。また、AtomTypeのアサインや、QSAR Descriptorとして用いる構造フラグメントの検索などにも利用できると思います。



人気ブログランキング(クリックして応援してね)

Windowsでバベる！

Babelといえば有名な化学構造のファイル変換ソフトウェアです。普段はもっぱらLinuxを使っているのですが、Windows上でも動かしたいとなと思いインストールしてみました。

今回はOpenBabelをソースからコンパイルしたいと思いますので、Cygwin上で作業を行います。Cygwinをインストールしていない場合は、まずはCygwinのインストールが必要です。

OpenBabelは以下のサイトより入手できます(現在のバージョンは2.0.1)。
http://openbabel.sourceforge.net/

openbabel-2.0.1.tar.gzを作業ディレクトリにコピーします。
$ cp openbabel-2.0.1.tar.gz (作業ディレクトリ)

作業ディレクトリに移動して、展開します。
$ tar zxvf openbabel-2.0.1.tar.gz

コンパイル・インストール作業です。
$ cd openbabel-2.0.1
$ ./configure
$ make
$ make install

無事にインストールできました。正直Cygwin上で一発でインストールできるとは思ってなかったので、Cygwin恐るべしっという感じです。

Babelの使い方はいたってシンプルです。
例)mol形式から水素原子を付加してPDB形式へ変換
$ babel -imol test.mol -opdb test.pdb -h

変換可能なファイル一覧は
$ babel -H
で取得できます。

mol2形式に変換するとGasteiger-Marsili法で計算された部分電荷が得られます。
$ babel -imol test.mol -omol2 test.mol2 -h
部分電荷はドッキングソフトウェアで利用したりできますので、役にたつと思います。

実は、以上のようなファイル変換はOpenBabelの一面にすぎません。
OpenBabelは、高機能なAPIを提供しており、これを利用することにより、SMILESによる部分構造検索や、分子構造の重ね合わせなど様々なことができるのです。しかも、OpenBabelはオープンソースですので、ケムインフォマティクスのエッセンスの詰まったソースコードを実際に見ることができ、学習するには最適な教材にもなります。具体的なプログラミングの仕方は、ぼちぼちメモしたいと思っています。



人気ブログランキング(クリックして応援してね)