CDKで型があわねーとかいうエラー

Depth-FirstでみっけたStructure-CDKが使いたいということでcdk-20050826.jarからcdk-20060714.jarにあげたらこの前つくったモジュールが動かなくなった。

sdg.setMolecule(mol);
sdg.generateCoordinates();
mol = sdg.getMolecule();

のところで、

Smi2Mol_899b.java:49: 互換性のない型
検出値  : org.openscience.cdk.interfaces.IMolecule
期待値  : org.openscience.cdk.Molecule
        mol = sdg.getMolecule();
                             ^
エラー 1 

もげー、型が違うってなんでじゃ?

planet chemoinformatics をちょっと修正

planet chemoinformaticsのRSS Auto Discoveryが微妙におかしい(/が抜けてた)ので修正。これで、livedoor readerとかbloglinesに登録できる。ま、直接rssを放り込んでも良いんだけど。僕はwebでほとんどみないので、デフォルトのテンプレートをいじって綺麗にしたいという欲求が沸いてこない(今のとこ)。だからRSSリーダーで読んだほうが良いと思うヨ。ちなみに、一日3回更新するように設定してます。

あと、docking.orgも追加した。rssが日付と本文を含まないので、Filter::EntryFullTextを書いてみた。

author: kzfm
handle: http://docking.org/.*?article\.pl\?
extract: (<TABLE WIDTH="100%" BORDER="0" CELLPADDING="5" CELLSPACING="0">
<TR><TD BGCOLOR="#FFFFFF">.*?                  </TD></TR></TABLE>)
.*?NAME="returnto" VALUE="//docking.org/article.pl\?sid=(\d\d/\d\d/\d\d)/\d+">
extract_capture: body date
extract_date_format: %y/%m/%d

なかなか快適になってきた。

openbabelって簡易DBっぽくつかえるのね

知らんかった。Zincの drug-likeのsmiを使ってみた。200万件くらいあるので、結構でかめ。

babel druglike.smi -ofs

でインデックス作成する(結構かかる、athron3000+mem512Mで2時間くらい)。 インデックスが作成できたら、SMARTSで検索してみる。

$ babel druglike.fs -osmi -s 'N#Cc1cccc(F)c1C#N'
2 candidates from fingerprint search phase
C(#N)c1c(c(c(c(c1C#N)F)F)F)F    ZINC00155214
C(#N)c1c(c(c(c(c1C#N)F)F)N)F    ZINC01507036
2 molecules converted

おお、7秒くらいで結果を返すので結構速い。
っていうか、babelのwikiみたら書いてあった。

babelはIDで検索はできないようだが、これはsmilesをsqliteにでも突っ込んでおけばいいので、Catalyst使えばお家いじり用のcompounds_dbが出来上がるな。

obfitのソースを眺めてみた

openbabelの分子重ね合わせツールのobfitをperlのモジュールにできんもんかな?とソース(C++)を眺めていたら、QTRFIT というアルゴリズム(実装はC++)を使っているようだ。

結局最大の重ねあわせを求めるのを固有値問題に落とし込んでるようだ。

QTRFIT

The eigenvector associated with the largest eigenvalue gives the desired quaternion. To generate a rotation matrix we consider the effects of the quaternion on the cartesian unit vectors

で、それをJacobi法で解いていると。感覚的にはそうなんだろなぁと思うが、ここら辺の理解度はちょっとあやしい。

RNews volume 6/3はオモシロげ

CRANでRNewsの最新版が読めるようになってた。今号はなかなか興味深い記事が多かった。全部読んでみたけど特に印象深かったのはこの3つ

Fitting dose-response curves from bioassays and toxicity testing

毒性試験との用量依存性曲線を描くためのコツというかそんな感じの内容。

酵素アッセイとかだと普通にシグモイド曲線のフィッティングをすればいいが、ED50なんかを求める系だと、hormesisがおきるので、特に低用量下でブランクよりも高い値がでる場合がある。そういう場合は linear-logistic modelを導入するとよいらしい。

The pls package

僕はよく使うので、multi-response modelsのとこと、varidation周りの話題が参考になった。

Generating, Using and Visualizing Molecular Information in R

RからCDKを使う話。SJavaをつかってRから記述子を抽出したりとか画像を表示したりとか。僕はRから他の言語を使おうとはあまり思わないほうなので参考になったという感じ。むしろ、jython+CDKRpyを組み合わせて色々できると面白いナァと思う。

いまさらながらSOAP::Lite

クライアントじゃなくてサーバーのほうをSOAP::Lite for Perlを参考に勉強中。久々の仕事絡みのプログラミングだ。

うちはケモインフォ関連の記述子発生ソフトや3次元構造立ち上げツールが複数のサーバーに分散されているという環境なので、rsh全開力技のスクリプトを書かないとフローが処理できない。でも、最近のツールにはSOAPじゃないと使えないのがあったり、あんまガシガシ書くとメンテが大変だったりするので、よく使うスクリプトはモジュール化しつつ、SOAPに対応することにした。

実はSOAPでサーバー側のスクリプト書くのは初めてで、色々悩むとこが多いなぁという感じ。特に記述子発生ツールはsmiles一化合物ごとに処理すると通信のオーバーヘッドが馬鹿になんないんじゃないの?だから束ねて計算させたほうがいいのかなとか、mopacとかgamessとか構造最適化みたいに意外に時間がかかる処理はどうしようかな?と。要するに非同期通信させないと駄目なんじゃないということなんですが。

で、いろいろ探し回ったところ、非同期にメソッドを実行するには?という遺伝研メソッドを発見。ほお、イメージとしては大分昔のNCBIがBLASTをリフレッシュ使って実装してたようなので捉えればいいのかな?

あとはSOAP::Liteでオブジェクトアクセスしたときにエラーハンドリングはどうやるんだろうか?やっぱ、関数型のやり取りのほうが、SOAP::Lite以外のクライアントとのやり取りがラクなんだろうなとか感覚的に思うことはあるんだけど、実際に色々作ってみて使ってみないとわからんな。

Inline::JavaでCDKのニ次元構造立ち上げをperlから

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

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

ProductName 実用 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のコード

:::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を意識しなくてよくなる。

Computer Applications in Pharmaceutical Research And Development

Structure-Property Relationshipsが最近面白げなもんで、今日の朝の演題はなかなか。QSPRとかやりだすと、その先にはパスウエイとかchip解析とかの連携が待っているんだろうなぁとか思ってただけに、新鮮だった。

Catalyst使った論文は常にチェックしているので3D-QSAR関連 の論文はおさえてたんだけど、Kohonen mapsの論文 は知らなかった。明日読んでみよう。

あんま、人でサーチかけたりとか更新チェックしてなかったけど、そういうこともしとくと新しいアイデアに出会うチャンスも増えるかな?とか思った。

あと、公演の最後で下の本を紹介してた。

amazon.comのほうだと中身が少し見れる たので、目次を眺めるとちょっと読んでみたい気が。

ソリューションガイドというものに駄文をのせてもらった

きわめてワタクシゴトなのですが、CBIという学会の年会で配布されるソリューションガイドといふものに、駄文を載せていただく機会を頂きました(ありがたや)

Yet Another Drug Discoveryとかいうタイトルで、(自分のパートは)かなり好き勝手に書かせてもらいましたが、目を通す機会などあったら、感想をこそっとでもズバッとでもきかせていただけるとありがたいかなと。

と、ふと思ったのでありますヨ。

ネタ的には、イントラにRSSリーダーを入れてあれとかこれとかRSS化してみようとか、pubmedに特化したdel.icio.usを作ってみたとか、まぁそんな感じです。

そして最近は、feed on feedsよりもユーザビリティ高そうなTiny Tiny RSSでも入れてみようかと思っている。