PerlMolでSMARTS

Chemistry::File::SMARTS Version 0.22

基本的な事柄

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

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

patern match check

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#!/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クイズを読み終えた

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

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

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


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

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

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

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

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

openbabelのperlラッパー

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って使えるのかな?これもあとで調べる。

ProductName 実用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

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

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)

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

dareact

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

#!/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で確認してみると

daproduct

ちゃんと環まいとるヨ。

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

  • 基質の用意と反応する部分構造のパターン定義
  • 生成物のオブジェクトを作成し基質を放り込む(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)

当初、CMLからSVGに変換するのはxsltでも用意すればと考えていた。だが、xsltイマイチ分からん、というより覚えるのがだるい。というわけで、「だらだら流れる(逐次処理する)んだったらSAXでも使おうかナァ」なんてノリでとりかかってみた。 しかし、bondのatomRefがatomidで、座標をatomエレメントから呼んでこないといけないので、やっぱSAXもメンドイかなってことで最終的にDOMりましたヨ。

perlのXML::DOMSVGモジュール使用。

とりあえず出来たSVGをconvertでpngにしてみた。

SVGはココだが一回セーブしてからでないとビューアーがうまく立ち上がらないような。

あとは、アトムラベルのとことかダブルトリプルボンドの処理とかだけど、丁寧に書いていけばよいな。それから、スケーリングもちゃんと書かないといけないなぁ。もうちょっと、ましなスクリプトになったら落とせるようにしておくヨ。

しかし、重要な問題が残っている(まぁコレに気付いてやる気レス化したため、プログラム書くのが延び延びになったんだが、、、)

それは、

openbabelPerlMolも座標を自動で発生してくれない!

ということだ。ちゅうわけで、SMILESからCMLにすると泣く。まぁ座標の問題はそのうち対応するみたいなんで期待してる。

Perlmolのsmilesには名前が必須みたいだ

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');

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

ケミッチョ

さて、ここで、ケミッチョという言葉を定義しよう(なんつって)。

ケミッチョとは、「それなりに生物学的等価性とか、(できれば)化学反応論に沿った形での合成化合物アイデアをたたき出す何か(でもいわゆる知的ではない)」と定義します。

そして、ケミッチョのperlでの実装をchemichoとしてみた。つまり、図のように、まろやかに変異させていきつつ、中間状態をいくつか保持しておいて、サイクルをまわすようなコードを書くと。ちなみに、図はperlmolで反応を記述しやすいように部分構造検索が反応のトリガーになるようにしてます。

chemicho

変異の重みは適当に制御できるようにYAMLかXMLかなんかのコンフィグレーションファイル形式にすればいいでしょう。ってことで、暇を見てはコード書き書きしてます。

一応、創薬的な偉さのイメージは↓な感じになると思う。

ケミッチョ(ランダム) << (メディシナル)ケミスト << 偉い仔(いわゆる神がかりな)

とりあえず、ランダムにアイデア生成する人工無能っぽいモノをつくってみたい(ネーミングも含めて) そして、グーグル界でオンリーワンを狙え!みたいなノリで、思いついたのが半年以上前のこと。で、エントリ書くのサボっているうちにグーグルでのオンリーワンはもろくも崩れさった模様。

あと独りブレストにはやっぱwemaがいい!図はwemaで書いてます。

debian(sarge)のperl5.8.4はperlmolが動かん

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をインストール

colinuxはtelnetとftpの設定してしまえば、後はpoderosaみたいなターミナルソフトでアクセスだ(colinuxのコンソールは使いづらい) あとキーボードを日本語対応にしたり、タイムゾーンを変更したりはココ参照だ。

とりあえず、せっかくdebian入れたことだし、sargeにアップグレードしておく。/etc/apt/source.list変えて、コマンドポチポチするだけダ。apt-getかなり便利と思った。

さて、Linux側の環境を整えたら、BioperlPerlMolをインストールするヨ。

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

PerlMolはCPAN経由で入れる。

perl -MCPAN -e 'install PerlMol'

BioperlはCPAN経由だとやたらとはまって口が半開きになって、なんでだよチクショーなんてむかつきまくってイカリンコ汁が漏れるので、素直にapt-getで(それが大人スタイル)。

Package: bioperl

素直にしこしこapt-get installしましょう。blastとかclustalwなんかも一緒に入れておくと良いですぞ。

© kzfm 2003-2016