perlmol / Drkcore

PerlMolでSMARTS

Mon, 23 Oct 2006 20:00:19 +0919

Chemistry::File::SMARTS Version 0.22

基本的な事柄

implicit,explicit
implicit:明示的に指定されていないbondとか(水素)原子も含んだ完全な形で。
smiles記法では（明らかに省略できる）水素原子は省略するのでexplicit,implicitという区分けがある。MOLもそうだけど。

http://cheminformatics.seesaa.net/article/23253094.html

などを参照するとイメージはわかるはず。

patern match check

パターンマッチの確認は以下のコードでおこなった。

#!/usr/bin/perl

use warnings;
use strict;
use Chemistry::File::SMARTS;
use Chemistry::File::SMILES;
use Chemistry::Ring 'aromatize_mol';

my ($SMARTS,$SMILES) = @ARGV;
print "SMARTS: $SMARTS\n";
print "SMILES: $SMILES\n";

my $pattern = Chemistry::Pattern->parse($SMARTS, format => 'smarts');
my $react = Chemistry::Mol->parse($SMILES, format => 'smiles');
aromatize_mol($react);

my $count = 0;

while ($pattern->match($react)){
    $count++;
}

print $count,"\n";

Atomic Primitives

*
全ての原子にマッチ
a,A
Aromatic(a) or Aliphatic(A)。 PerlMolではChemistry::Ringでaromatize_molしとかないと認識しない。
D<n>,X<n>,x<n>
結合原子数(n)を有するAtomにマッチ。明示的(D)か、内在的(X)か。つまり、Dは明示的にHとの結合が指定されていない限りこれをカウントしない。xはリングのコネクティビティだが、PerlMolではサポートされてないっぽい。
H<n>,h<n>
明示的なHをカウントする(H)か内在的なものも含めるか(h)
R<n>,r<n>
n原子からなるSmallest Set of Smallest Rings (SSSR)のメンバーにマッチ(R)。PerlMolだとRかどうかだけしか判断していないので<n>は関係ないかも。リングサイズが最小のSSSRのメンバーにマッチ(r)。芳香環かそうでないかは判断せず、n員環のメンバーであればマッチ

こんな感じ

$ ./smartstest [r6] 'C1CCCCC1CCC' # 6
$ ./smartstest [r6] 'c1ccccc1CCC' # 6
$ ./smartstest [r6] 'c1cc1CCC'    # 0
$ ./smartstest [r3] 'c1cc1CCC'    # 3

v<n> バレンス。いわゆる手を伸ばせる数(Nは3で、Cは4)
<n>,+<n>,#n チャージの有無と原子番号
@,@@,@<c><n>,@<c><n>? キラルはPerlMolではサポートされてない。

Bond Primitives

-,=,#,: 一重結合(-)、二重結合(=)、三重結合(#)、芳香属性の結合(:) aromatizeしてるのに二重結合がカウントされるのがわからん。局在化した芳香環の二重結合だけ読んでるっぽいので、SMARTS書くとき注意

下の例

$ ./smartstest "c:c" 'c1ccccc1'
SMARTS: c:c
SMILES: c1ccccc1
6
$ ./smartstest "c=c" 'c1ccccc1'
SMARTS: c=c
SMILES: c1ccccc1
3
$ ./smartstest "c-c" 'c1ccccc1'
SMARTS: c-c
SMILES: c1ccccc1
0

~,@ ワイルドカード(~)、リング結合全て(@)
/,,/?,? PerlMolでサポートされてるかどうか調べてない(とりあえず使わないので後回し)

Logical Operator

!,&,,(commma),;

否定(!),論理積(&,; ただし、&は最初に評価される),論理和(,)

Component-level grouping

PerlMolではサポートされていない

こういうのをblogに突っ込むのは無理がある。素直にWiki+Markdown入れようと思ったヨ。

Perlクイズを読み終えた

Thu, 03 Aug 2006 19:57:55 +0919

一ヶ月くらい前から暇をみてはPerlクイズを解いていたのだが、ここ2,3日は集中してすすめることができた。というわけで、本日終了。

最後のほうはパズルばかりでかなり手ごわいクイズが並んだが、一通りやり遂げた。

結城浩のPerlクイズ
結城浩
ソフトバンククリエイティブ / ?円 ( 2002-08 )

特に面白かったのはここらへん

初盤ではハッシュの操作とか正規表現なんかの問題が多かったが、特にソート関係は勉強になった。

あとは、「あらかじめ行数がわかっていなくてもランダムに行を表示する方法は特にイテレータ使ってる場合にランダムに一個選びたい場合にイイ。

chemichoっていうプログラムを書いていたときに、perlmolの部分構造マッチのメソッドがイテレータで、ヒット数を返すメソッドがなかったため、マッチした部分の中からランダムで一箇所選ぶためにわざわざ一回スキャンしてヒット数を調べてから、int(ヒット数*rand)番目を再度イテレータ使って取得するという二回同じイテレータをまわすプログラムだったという。

これって無駄だなと思っていたので、上のやり方はエレガントじゃ！と感動した。

openbabelのperlラッパー

Wed, 07 Jun 2006 23:32:51 +0919

openbabelにはobfitっていうツールがあることをケムインフォマティクスに虚空投げで知った。なかなかイカス。今まで、「babelって調子いいフォーマットコンバーターだよな」ぐらいの認識だったヨ。

でもobfitはmol,molで重ね合わせるコマンドライン向きのプログラムなんで、普通にライブラリをMCSでぶん回すには、シェルスクリプト書いたりしなきゃならんし、シェルスクリプトでまわすとin,outでテンポラリのファイルがたくさんできるからやだナァ。

と思ったのでドキュメント漁りをしてみたら、Chemistry::OpenBabelっていうperlラッパーがあった。うーんCPANでみたことないなぁ、、、と思ったらCPANにはあがってなくて、展開したアーカイブのscripts/perl/で配布してるみたい。

ソースを眺めたら、SWIG使っているのね。普通にインストールはOK。C++の関数をそのままPerlでラップしてるみたいなので、obfitのソース眺めながら、そのままperlで書きなおしてやれば、perlスクリプトとして使えるかな?って感じなのであとでやる。

それよりもSWIG覚えたいかも。実用perlプログラミングの初版には詳しく書いてあったはずだけど、第二版ではばっさり削られてInline::Cとかになってた。こんなん覚えてくるとperlでjavaライブラリのCDKとか呼びたくなるにちがいないんだろうけど、Inline::Javaって使えるのかな？これもあとで調べる。

実用Perlプログラミング
スリラムスリニバサン
オライリー・ジャパン / ?円 ( 1998-11 )

ただ、API見る感じだと、Chemistry::OpenBabeで分子のモディファイとかリアクション関連いじるのやりにくくなかろうか？やっぱPerlMolははずせないわなと思ったが、Readmeに、そのうちPerlMolとopenbabelのグルーモジュール作るよん！って書いてあった。ナイス、楽しみ。

Future development will include a "glue module" to integrate with PerlMol and other Perl chemistry projects.

で、さらにケムインフォマティクスに虚空投げでOELibという記述が何度か出てきたので、もしやOpenEyeかと思ったら、やっぱそうだった。OpenEyeのライブラリなんだったら別にPerlMol使わなくても、化合物いじり出来るんではないの？とか思ったが、FAQに歴史的な経緯が書いてあった。だから、やっぱ素直にPerlMol使ったほうがいいのでしょう。OEChemはOELibを書き直してるってことはデータの内部構造結構違うんだろうな。

というわけで、ツールキットみたなかではいまんとこOEChemが一番良く出来てる気がしてる。

In silico Library Enumeration of Synthetically Feasible Libraries

Thu, 18 May 2006 21:54:42 +0919

みんな、似たようなこと考えるのね。

Cheminfostream: In silico Library Enumeration of Synthetically Feasible Libraries

“The core functionality to enumerate focused de novo type libraries in Oracle for a number of organic transformations is already available”

一応、Feasibleだったとしてもだ、、、

Sparseすぎ
ActiveCompoundあったとしても、活性相関探るために近傍のcompsつくるために手動かさなきゃならんので意外にrapidじゃない。H2Lの部分は速くならないってこと。
造り過ぎ
楽しいので、ついつい作りすぎてしまうネ。100,000　x 100,000で100億でしょ。さらに多段階合成なんかでもう一回かますともう大変。作るほうはドンドコドンドコ量産されるのでとっても嬉しいんだけど。

でも、大変なのはこっからだ。

どうやって絞込みすんの?

精度のいい予測つかって精度99.0%のモデルでも1億化合物残るね。Rule of X(5,3)だとかいっても、そもそも考慮されていたりするから、意外に減らせなかったり。そうすっと結構泣く。

調子にのって作りすぎなきゃよかった

みたいな。みんな若気の至り汁が出て困ったことがあるはず。

Diels-Alder Reaction (PerlMol)

Sat, 11 Feb 2006 17:15:53 +0919

PerlMolで化学反応のルーチンってどう書けばええの？と聞かれたことはないが、こんな感じで僕は書いてマスヨってことで、ひとつDiels-Alder反応をさせてみる。

適当に用意したこんな反応だ。

#!/usr/bin/perl

use strict;
use warnings;
use Chemistry::File::SMILES;
use Chemistry::File::SMARTS;
use Chemistry::Ring 'aromatize_mol';

#反応させる基質を用意
my $react1 = Chemistry::Mol->parse('ClC=C dienophile', format => 'smiles');
my $react2 = Chemistry::Mol->parse('OC=C-C=C diene', format => 'smiles');

#dieneとdienophileのパターンをSMARTSで定義
my $dienophile_pat = Chemistry::Pattern->parse("C=C", format => 'smarts');
my $diene_pat = Chemistry::Pattern->parse("C=C-C=C", format => 'smarts');

#まずおまじない
aromatize_mol($react1);
aromatize_mol($react2);

#生成物に名前をつけ,基質をくっつける。
my $name = $react1->name . "+" . $react2->name;
my $prod = Chemistry::Mol->new(name => $name);
$prod->combine($react1, $react2);

#dieneとdienophilieのパターンを探し、みつかったら各オブジェクトにマップ
$dienophile_pat->match($prod);
my @atom_map1 = $dienophile_pat->atom_map;
my @bond_map1 = $dienophile_pat->bond_map;

$diene_pat->match($prod);
my @atom_map2 = $diene_pat->atom_map;
my @bond_map2 = $diene_pat->bond_map;

#反応用の原子が揃っているかチェック
if ($atom_map1[0] && $atom_map1[1] && $atom_map2[0] && $atom_map2[3]) {

#結合の作成
    $prod->new_bond(atoms => [$atom_map1[0], $atom_map2[0]], order => '1');
    $prod->new_bond(atoms => [$atom_map1[1], $atom_map2[3]], order => '1');

#dieneとdienophileの結合次数を変える
    $bond_map1[0]->order(1);
    $bond_map2[0]->order(1);
    $bond_map2[1]->order(2);
    $bond_map2[2]->order(1);

my $smi = $prod->print(format => 'smiles', unique => 1, name => 1);

#SMILESで出力
    print "$smi\n";
}



早速実行してみると、

$dareac.pl 
OC1C=CCCC1Cl    dienophile+diene



jchempaintで確認してみると

ちゃんと環まいとるヨ。
さて、スクリプトをくどくどと書いたけど、処理の流れは単純でわかりやすいのですヨ。

基質の用意と反応する部分構造のパターン定義
生成物のオブジェクトを作成し基質を放り込む(combine)
パターンマッチさせたら、原子と結合に位置のマップ
マップされた結合と原子に対して処理
結合新しく作ったり、切り離したり
原子オブジェクトの削除とか

エッジとノードで表現するせいなんだと思うがマップされた結合と原子に対して処理するルーチンが結構煩雑になり気味で、デバッグは大抵ココに集中したりする。dieneの端っこ結合作ったら、pな軌道が隣り合うから間は二重結合になるだろって思うが    $bond_map2[1]->order(2)って明示しないと、ラジカルとして表現されちゃったりとか。sp2,sp3の情報入れとけば自動で判別するようにできるような気がするが難しいのだろうか?
あとは、サンプルなのではしょってますが、この反応ではpara体も生成させることは出来ます。

#結合の作成
    $prod->new_bond(atoms => [$atom_map1[0], $atom_map2[3]], order => '1');
    $prod->new_bond(atoms => [$atom_map1[1], $atom_map2[0]], order => '1');



って書けばpara体出来上がり。
というわけで、普通は二通りの生成物を提示するんだろうけど、
Diels-Alder反応
Diels-Alder付加体の二つの置換基はオルト位/パラ位を占めるように位置選択的に付加することが有機電子論から経験的に予測される(オルト-パラ則)。詳細にはFrontier軌道のローブ係数を求めて比較することで説明がなされる。すなわち、下図のようにHOMO/LUMOの係数が大きい点同士が重なるように付加する。

とあるように、軌道の計算してやればどっちが出来やすいか予測することは可能なので、反応処理ルーチンにPyQuanteみたいな量子化学計算モジュールを組み込んで判断させてやれれば素敵かなぁと思っている。



CML2SVG(Perl)
Thu, 02 Feb 2006 17:32:47 +0919
当初、CMLからSVGに変換するのはxsltでも用意すればと考えていた。だが、xsltイマイチ分からん、というより覚えるのがだるい。というわけで、「だらだら流れる(逐次処理する)んだったらSAXでも使おうかナァ」なんてノリでとりかかってみた。
しかし、bondのatomRefがatomidで、座標をatomエレメントから呼んでこないといけないので、やっぱSAXもメンドイかなってことで最終的にDOMりましたヨ。
perlのXML::DOMとSVGモジュール使用。
とりあえず出来たSVGをconvertでpngにしてみた。


SVGはココだが一回セーブしてからでないとビューアーがうまく立ち上がらないような。
あとは、アトムラベルのとことかダブルトリプルボンドの処理とかだけど、丁寧に書いていけばよいな。それから、スケーリングもちゃんと書かないといけないなぁ。もうちょっと、ましなスクリプトになったら落とせるようにしておくヨ。
しかし、重要な問題が残っている（まぁコレに気付いてやる気レス化したため、プログラム書くのが延び延びになったんだが、、、）
それは、
openbabelもPerlMolも座標を自動で発生してくれない!
ということだ。ちゅうわけで、SMILESからCMLにすると泣く。まぁ座標の問題はそのうち対応するみたいなんで期待してる。


Perlmolのsmilesには名前が必須みたいだ
Thu, 19 Jan 2006 23:35:52 +0919
perlmolでparse_smiするとwarningが出るのでなんでかな？と思って調べたところ、名前を入れてやらないと駄目らしい。

Use of uninitialized value in concatenation (.) or string at /usr/lib/perl5/site_perl/5.8.3/Chemistry/File/SMILES.pm line 500.
Oc1ccccc1Cl 


ってエラーがでるので、

my $react = Chemistry::Mol->parse('c1ccccc1Cl Cl-benzene', format => 'smiles'); 


ってしてやれば解消される。


ケミッチョ
Tue, 17 Jan 2006 21:49:09 +0919
さて、ここで、ケミッチョという言葉を定義しよう（なんつって）。
ケミッチョとは、「それなりに生物学的等価性とか、(できれば)化学反応論に沿った形での合成化合物アイデアをたたき出す何か(でもいわゆる知的ではない)」と定義します。
そして、ケミッチョのperlでの実装をchemichoとしてみた。つまり、図のように、まろやかに変異させていきつつ、中間状態をいくつか保持しておいて、サイクルをまわすようなコードを書くと。ちなみに、図はperlmolで反応を記述しやすいように部分構造検索が反応のトリガーになるようにしてます。

変異の重みは適当に制御できるようにYAMLかXMLかなんかのコンフィグレーションファイル形式にすればいいでしょう。ってことで、暇を見てはコード書き書きしてます。
一応、創薬的な偉さのイメージは↓な感じになると思う。
ケミッチョ(ランダム) << (メディシナル)ケミスト << 偉い仔(いわゆる神がかりな)
とりあえず、ランダムにアイデア生成する人工無能っぽいモノをつくってみたい（ネーミングも含めて)
そして、グーグル界でオンリーワンを狙え！みたいなノリで、思いついたのが半年以上前のこと。で、エントリ書くのサボっているうちにグーグルでのオンリーワンはもろくも崩れさった模様。
あと独りブレストにはやっぱwemaがいい！図はwemaで書いてます。


debian(sarge)のperl5.8.4はperlmolが動かん
Sun, 15 Jan 2006 15:19:32 +0919
R4でコード書いていたら、以下のエラーが。

Weak references are not implemented in the version of perl at /usr/local/share/perl/5.8.4/Chemistry/Atom.pm line 51


どうも、Scalar::Utilのweakenが駄目みたいで、perlのバージョンのせいっぽい感じ。
なんで、Perlのバージョンあげるか落とすかしないといけないんだけど、apt-getいまいち把握していないから、いじるの怖い。
と、思ったら、
colinuxだからイメージファイルコピーして、バックアップとっておけば、やりたい放題なことに気づく。
これから、ソースから5.8.7入れます。


colinux(debian)にBioperlとPerlMolをインストール
Tue, 20 Dec 2005 12:29:15 +0919
colinuxはtelnetとftpの設定してしまえば、後はpoderosaみたいなターミナルソフトでアクセスだ(colinuxのコンソールは使いづらい)
あとキーボードを日本語対応にしたり、タイムゾーンを変更したりはココ参照だ。
とりあえず、せっかくdebian入れたことだし、sargeにアップグレードしておく。/etc/apt/source.list変えて、コマンドポチポチするだけダ。apt-getかなり便利と思った。
さて、Linux側の環境を整えたら、BioperlとPerlMolをインストールするヨ。

PerlMol
オブジェクト指向で化合物とか反応を扱う。酸クロとかの反応用のモジュール作っておくとさらに便利。ただし、蛋白質とかの巨大なものを扱うにはどうじゃろか？今ひとつな気もする。構造立ち上げとかしっかりしてくればかなりいいツールになるんじゃないかと思っている。
Bioperl
とりあえずお約束。フォーマット変換とかは便利。ただし、モジュールの質はまちまち。ドキュメントに不備がある場合、直接コードを読むことになるが、そういうモジュールに限って泣きたくなる度が高い。で、結局自分で書いちゃったみたいな。

PerlMolはCPAN経由で入れる。

perl -MCPAN -e 'install PerlMol'



BioperlはCPAN経由だとやたらとはまって口が半開きになって、なんでだよチクショーなんてむかつきまくってイカリンコ汁が漏れるので、素直にapt-getで（それが大人スタイル）。
Package: bioperl 
素直にしこしこapt-get installしましょう。blastとかclustalwなんかも一緒に入れておくと良いですぞ。