_in house_でケモインフォ関連のプログラミングをするのに、二次元構造の立ち上げツール（またはライブラリ）と、構造描画ツールは必須だが、いまのとこCDKを使うという選択肢しかない。

実用PerlプログラミングにInline::Javaを使うとjavaのコードをperlに組み込めると書いてあったので、CDKを組み込んでみたヨ。

実用 Perlプログラミング 第2版
Simon Cozens
オライリージャパン / 3360円 ( 2006-03-01 )

まぁ、要するに、
Q.そこまでしてperlで処理したいのか？
A. もちろんですヨ、奥さん
ということなんだが。

さて、インストールはj2sdkの場所を指定しないといけないのでcpan -iだとこける。そのため地道に、buildディレクトリで再インストール。

$ cd .cpan/build/Inline-Java-0.51
$ perl Makefile.PL J2SDK=/usr/java/jdk1.5.0_08/
$ make java
$ make
$ make test
$ make install

ただ、make javaで

注:入力ファイルの操作のうち、未チェックまたは安全ではないものがあります。
注:詳細については、-Xlint:unchecked オプションを指定して再コンパイルしてください。

とかいう警告が出るが、testは通るのでまぁ気にしないということで。

で、コードを書いてみる。

参考にしたjavaのコード

• http://cheminfo.informatics.indiana.edu/~rguha/code/java/
• http://cheminformatics.seesaa.net/article/20455721.html

:::java use Inline Java => <<'END_OF_JAVA_CODE' ; import java.io.StringWriter; import java.io.IOException; import org.openscience.cdk.Molecule; import org.openscience.cdk.exception.CDKException; import org.openscience.cdk.io.MDLWriter; import org.openscience.cdk.smiles.SmilesParser; import org.openscience.cdk.layout.StructureDiagramGenerator; import org.openscience.cdk.tools.HydrogenAdder;

class smi2mol{

    public smi2mol(){
    }

    public String convert(String smi){
    Molecule mol = null;

    try {
        SmilesParser sp = new SmilesParser();
        mol = sp.parseSmiles(smi);
    }catch (Exception ise){
        ise.printStackTrace();
    }

    StringWriter w = null;
    MDLWriter  mw = null;

    try {
        w = new StringWriter();
        mw = new MDLWriter(w);
    } catch (Exception e) {
        e.printStackTrace();
    }

    HydrogenAdder ha = new HydrogenAdder();
    try{
        ha.addExplicitHydrogensToSatisfyValency(mol);
    } catch (Exception e){
        e.printStackTrace();
    }

    StructureDiagramGenerator sdg = new StructureDiagramGenerator();
        try{
        sdg.setMolecule(mol);
        sdg.generateCoordinates();
        mol = sdg.getMolecule();
        mw.write(mol);  
    } catch (Exception e){
        e.printStackTrace();
    }

    return w.toString();
    }
}

END_OF_JAVA_CODE

my $alu = new smi2mol();    
my $smiles_string = shift;
print $alu->convert($smiles_string);

実行してみる

::: sh $ perl test.pl "C(=O)O"

  CDK

  5  4  0  0  0  0  0  0  0  0999 V2000
    1.2990   -0.7500    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000    0.0000    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    2.5981   -0.0000    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    1.2990   -2.2500    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
    3.8971   -0.7500    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
  2  1  2  0  0  0  0 
  3  1  1  0  0  0  0 
  1  4  1  0  0  0  0 
  3  5  1  0  0  0  0 
M  END

水素が付加されて、座標も計算されたMOLファイルが出力される。

javaって例外を返す可能性がある場合try,catch構文入れないと強制的にエラーになるのね。try,catch構文は読みやすいけど、強制されんのはちょっとやだナァと。単にフォーマットをコンバートするだけのコードなのに結構長くなっちゃったし。

::: perl use CDK; use Chemistry::Mol;

my $mol = CDK::smi2sdf($molecule);
print $mol;

みたいにモジュールにしてuseして使えばCDKを意識しなくてよくなる。

コメントでhaskellでフィボナッチの別解をもらったが、パッと見てもなぜこれでフィボナッチ数列ができるのかいまいちわからんかったので、よく読んでみたヨ。

まずは型チェック。 見ればわかるとおりfibはリストを返す関数ですな。つまり、fibを実行するとフィボナッチ数がどんどこどんどことそれこそとめどなく溢れてくる(はず)

fib :: Num a => [a]
fib = 0:1:zipWith (+) fib (tail fib)

これを　fib.hsという名前で保存してghciで実行してみる。

$ghci
*Main> :l fib.hs
*Main> fib

をを～。凄い勢いでフィボナッチ数列がどんどこど(以下略)。ちなみにlengthでリストの長さをはかってみると、プロンプトがかえってこないゾ。さすが無限リスト。

zipWithの型とtailの型はこんな感じ。

Prelude> :t zipWith
zipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
Prelude> :t tail
tail :: [a] -> [a]

zipWithは高階関数で、+関数をfibとtail fibに適用して新しいリストをつくるわけですな。tailは先頭の要素を取り除いたリストを返す関数だから、、、うーんやってることはわかるんだけど、イマイチイメージできん。

と悩んでたら、Wikipediaに書き下されてる表をみつけて納得した。perlで同じことやろうとするとまず有限個の要素をもつ配列があってみたいなところから始まるから、こういう表現はできないな。というわけで、目からうろこっていうか、ちょっと感動した。

ここで考えたようなn-1番目とn番目のフィボナッチ数が与えられたときにn+1番目のフィボナッチ数を求めるということを考えるのではなくて、そもそもフィボナッチ数列とは何かを考えて、それを無限リストとして表現することを考えればよいんですよと。

となると、リスト内包表記も感動するような綺麗さダ。

fibs = 0 : 1 : [ a+b | a <- fibs | b <- tail fibs ]

Haskell - Wikipedia

GHCの並列リスト内包表記(GHCの拡張は特別なコマンドラインフラグによって有効にしなければならない。

とあるのが、なにが問題なんだかわからんので、とりあえず実行してやろう。

fib4.hs:1:15: Illegal parallel list comprehension: use -fglasgow-exts
Failed, modules loaded: none.

デフォルトだと、同じリストが複数でてくるような表記の仕方(この場合はfibs)は問題あるということなのかな

ghci -fglasgow-exts

とオプションを指定して実行することで、OK。をを～。凄い勢いでフィボナッチ数列がどんどこど(以下略)。

haskell素敵過ぎ

U隊長がボニピンベストがITMSで買えないんだよねーとか言うので、amazonで。

Every Single Day -Complete BONNIE PINK (1995-2006)-
BONNIE PINK
ワーナーミュージック・ジャパン / 2492円 ( 2006-07-26 )

僕はEvil and Flowersが好き。

プログラム工学VIを読み終えて、次はHaskeller への道など読んでます。

さて、Hello, World! プログラムでllってコマンド出てきた。

llってls -lのことか?と自分のFC5端末でllと打つとファイルのリストが表示される。whichすると確かにls -lにalias切ってあるし、、、、

/etc/profile.d/のファイルで設定されてる。

最近のパッケージはllって標準なのね、昔のはなかったよな～
とか思いながら、今もってる最古のOSでll打ってみた。

Kondara MNU/Linux release 1.2 (Ayaka)
Kernel 2.2.17-16ksmp-ata66 on a 2-processor i686
login:

$ which ll
alias ll='ls -l --color=tty'
        /bin/ls

うぉ、しっかりあるし、、、、。今までずっと無駄にls -lって打ってたことに気付いて凹んだ。ls -lがllなんて知らなかった。(ls -aはlaで使ってるけど)

profile.d見てみるとllでls -lのほかにl.ってのがあって,.で始まるファイルを表示するようにaliasきってあるようだ。

Statistics::R なんかいまいちとか言って、もっぱらRのヒストリファイルをいじってバッチ用に書き換えて使う日々が続いてたが、何を思ったのか、突然RSPerlを使いたくなったのでインストールしてみた。

まずは、Rのインストールから。rpmで楽々。

$ rpm -qa | grep R-
R-devel-2.3.1-1.fc5
R-2.3.1-1.fc5

RSPerlもコマンド一発で楽チン。

$ R CMD INSTALL  --configure-args='--with-in-perl' RSPerl_0.91-0.tar.gz

環境変数の設定スクリプトがあるので、/etc/profileに書いておく

$ vi /etc/profile
#for RSPerl
source /usr/lib/R/library/RSPerl/scripts/RSPerl.bsh

これで、OK。

でRSPerlのいいところは

• 配列をリファレンスで渡せる(Statics::Rとは違う)
• AUTOLOADでRの関数を呼び出せる。R::plotとかできる

ので、結構便利だ。

でもマトリックスとかデータフレームの扱いがちょっとめんどくさいのと、Rを起動するために、 &R::initR("--silent") とかやらないといけないのと、$r->plotとかやりたいのにR::plotと書かなきゃいけないので、長いコードだとごちゃごちゃしてきて見通しが悪くなる。

結局Rでバッチのほうがコードが読みやすかったりとか。やっぱRpyか。昨日ソフトのインストールに来てた業者の人も最近pythonの案件増えてますよ～とか言ってたし。

みんなのPython
柴田 淳
ソフトバンククリエイティブ / ?円 ( 2006-08-22 )

うー、欲しくなってきた。

Fedora core 5のjdkだとJchemPaintが立ち上がらないので、SUNのJDKをインストール

rpm -e jdk
rpm -ivh jdk-1_5_0_08-linux-i586.rpm

profileに環境変数書いておく

/etc/profile

#java
JAVA_HOME=/usr/java/jdk1.5.0_08
CLASSPATH=.:$JAVA_HOME/lib/tools.jar:/usr/local/javalib/cdk-20050826.jar:\
/usr/local/javalib/jchempaint-2.0.12.jar

クラスパスってcdkとかjchempaintのjarファイル指定する必要あるんだろうか？わからん。

/usr/bin/javac,java,jarのシンボリックリンクの張りなおしをして作業終了

懐かしの名曲のremixが出ていたので。

二天一流のほうは、今風っちゅうよりも、アブストっぽさが失われたっちゅうか、なんかよくわからん。It's A Jazz Thingはちょっと回転速めの低音厚めで好きな感じの仕上がり。オリジナルはあれはあれでユラユラ感がたまらんのだけど。
Ni Ten Ichi Ryu
It's A Jazz Thing

メタモ楽しみとか書いたけど、やっぱ、ムンベ。つまり、DJ AKIが一番の楽しみなのだ。

Structure-Based Drug Designに従事していると、トラディッショナルな目視での解釈に限界と疑問を感じる時が遅かれ早かれやってくるはず。

そうなると、解いた蛋白-リガンド結晶構造をGAMESSとかABINIT-MPなんかでFMO計算し始めるわけだ（そのうちただ単に計算したいという理由でPDBのデータなんかに手出したり）これにより、電子的な相互作用をエネルギーの大小で考えることができるようになるため、複合体に対する理解が深まるのでイイ。

でも、やがてこれだけでは満足できなくなり、より深い解釈を求め始めるようになる。つまり、生化学、構造生物学から、有機化学、量子化学へと足を踏み入れはじめるようになるわけだ。

ただし、このような場合に、勉強するのにちょうどいい本や資料がないなぁと困っていたところ、GAMESSについて書かれている本が出版されたというエントリをみっけたので、脊髄反射的にamazonで購入してしまった。

s2k's eye 2nd Ed.: 「すぐできる量子化学計算…」を買いました。

Q&A形式で、どこからでも読めるというのも好印象です。そして、何より、今までの成書ではあまり無かったGAMESSについての解説が多く入れられています(GaussianとGAMESSの解説書です)。応用計算については多くがGaussianについてですが、これはGAMESSというプログラムの性質上しょうがないところでしょう。

すぐできる 量子化学計算ビギナーズマニュアル (KS化学専門書)

講談社 / ?円 ( 2006-04-29 )

内容は、FAQ集という感じで、200pぐらいある割には半日かからずに読み終えた。もちろんGAMESSに関する記述も結構あるので、Gamess Input Examplesを補完する目的でも使える。まぁ、計算化学系の研究室にいれば、こういったノウハウとかこつみたいなものは、誰かに聞けばすぐに解決するようなものなんだが、、、、
製薬系って量子化学計算する人って思ったほどは多くないので、D章の計算実践編とF章の目的別対処法は参考になった。

このように、初学者が悩んだりつまずいたりしそうなところをQ&A形式で解説している本なので、実験グループで量子化学計算を担当された方を対象にしているという前書きは納得がいく。

ただし、Howto本ではないから、ザボとか一通り読んであって、Gaussian,GAMESS触ったことがあって、unix系のエディタでファイル操作が普通にできるヒト向けでしょう。あとは、より深く知るためのリファレンスが少なすぎる気がするなぁ。改訂版出すならそこら辺を充実してくれればありがたい。

Wikiあればいいのに、RjpWikiの量子化学計算版みたいなの。

メタモのチケットきた。

BLACK STROBE MARK FARINA が気になる。

昨日に引き続き3章を読み進めていたら、フィボナッチ数をキャッシュで高速化する例が。

Higher-Order Perl: Transforming Programs with Programs
Mark Jason Dominus
Morgan Kaufmann / 6402円 ( 2005-03-28 )

コードはこんな感じで、割と普通な。

my %cache;

sub fib{
    my $n = shift;
    unless (exists $cache{$n}){
    if($n < 2){
        $cache{$n} = 1;
    }else{   
        $cache{$n} = fib($n-1) + fib($n-2);
    }
    }
    return $cache{$n};
}

さて、昨日のタプルを使ったコードとどっちが速いのか比べてみた。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

#!/usr/bin/perl

use strict;
use warnings;
use Benchmark qw(cmpthese);

sub fibStep {
    my $n = shift;
    return [$n->[1], $n->[0] + $n->[1]];
}

sub fibPair {
    my ($n, $fs, $fp) = @_;
    return [0,1] if ($n == 0);
    return $fs->($fp->($n-1,$fs,$fp))
}

sub fastFib {
    my $n = shift;
    my $fs = \&fibStep;
    my $fp = \&fibPair;
    my $fib_num = $fp->($n,$fs,$fp);
    return $fib_num->[0];
}

my %cache;

sub fib{
    my $n = shift;
    unless (exists $cache{$n}){
    if($n < 2){
        $cache{$n} = 1;
    }else{   
        $cache{$n} = fib($n-1) + fib($n-2);
    }
    }
    return $cache{$n};
}

my $number = 50;
my $count = 10000000;

my $fibcache = \&fib;
my $fibfast = \&fib;

cmpthese($count, {
        'Cache' => $fibcache->($number),
        'Fast' => $fibfast->($number)
    });

とりあえず1000万回まわしてみたけどそれでもwarningがでるので、この際無視。

./fib.pl 
            (warning: too few iterations for a reliable count)
            (warning: too few iterations for a reliable count)
            Rate  Fast Cache
Fast  38461538/s    --  -27%
Cache 52631579/s   37%    --

キャッシュしたほうが速いというのがイマイチ理解できん。何でじゃろか？