CRANのNewsletterが更新されていたので、一昨日の帰りの新幹線の中で読んでたら、fingerprintっていうパッケージが追加されてた。

CRAN - Package fingerprint

This package contains functions to manipulate binary fingerprints of arbitrary length. A fingerprint is represented by an object of S4 class 'fingerprint' which is internally represented a vector of integers, such that each element represents the position in the fingerprint that is set to 1. The bitwise logical functions in R are overridden so that they can be used directly with 'fingerprint' objects. Distances metrics are also available. Fingerprints can be converted to Euclidean vectors (i.e., points on the unit hypersphere)i and can also be folded using XOR. Arbitrary fingerprint formats can be handled via line handlers. Currently handlers are provided for CDK, MOE and BCI fingerprint data.

CDKとかMOE,BCI用のハンドラがついてるようです。で、Similarity metricは何使えんのかな？とマニュアル読んでみると、

• tanimoto
• euclidean
• mt
• dice

の4つ。このなかで、MT係数って聞いたことないナァ。ちょっと調べてみたけどわからなかったので直接ソースを読んでみることにした。はじめてのパッケージイントロスペクション。

で読んでみたら、あー、Rのパッケージ読むの単なる食わず嫌いだったかも。思ったよりややこしいことはしてませんな、これからは積極的に読んだほうがいいかもと思った。

というわけで、ちょっとまとめておいた。前提としてフィンガープリントの共通のビットが立つかどうかを下の表であらわしておく。

tanimoto

まずはわかりやすいタニモト係数から。

 dist <- c / (a+b+c)

dは考慮されない、つまりビットの立たないもの同士は類似とみなされない。

euclidean

dist <- sqrt((d+c) / (a+b+c+d))

ユークリッド距離はお互いビットが立ってないという情報も似ていると判断されるわけですな。

dice

dist <- c / (.5*a + .5*b + c)

diceの距離も基本的にはタニモトと一緒。僕はあまり使わない。

mt

size <- length(fp1)

t1 <- c/(size-d)
t0 <- d/(size-c)
phat <- ((size-d) + c)/(2*size)
dist <- (2-phat)*t1/3 + (1+phat)*t0/3

t1はc/(a+b+c)、t0はd/(a+b+d)とかける。phatは全フィンガープリントでのビットの立っている割合（ビットの密度）。ビットの密度を考慮してるのはわかったんだが、t1,t0の分母の意味がちょっと理解できてないかも(aとbの情報を考慮してんだろうな)。おそらくビットの密度に偏りがある場合にできるだけ補正できるような距離なのかなとか思った。タニモトにビットが立たない場合の情報の扱いを加えた感じか。

と、ここまで考えてやっと気付いた。

MTってもしかしてModified Tanimoto?

ビンゴだった。

昔作ったコードが気に入らなかったので書き直したら、そこそこ動くようになったので、Chemicho-0.01.tar.gzとしてダウンロードできるようにしておいた。これはPerlMolを使って実装された化合物にランダムに変異を入れていくツールで、例えば、H2Lでヒット化合物に修飾を軽く加えてみたりとかしてくれる、ケミストチックな小人さんに育つことを期待してるわけだ。

• AddCarbon(メチル基付加)
• Aromatize(ジエンを芳香環化)
• DeAromatize(芳香環をジエン化)
• DecreaseBondOrder(結合次数を減らす)
• GenerateRing(環化)
• Halogenize(メチルをハロゲンに)
• IncreaseBondOrder(結合次数を増やす)

設定はYAMLで。

Generationは世代数で、cycleは世代で生成される化合物の数。ちなみに第0世代は親。例えばcycle3で第三世代まで発生させると3^3で27化合物で4世代目には3^4で81化合物できる。ランダムな置換の頻度はweightで変えられるようになってる。

reactionはSMARTS使ってこんな感じで反応させてる。

いまんとこ、変異用のモジュールが少ないのでバリエーションが少ないのだけど、これから増やす予定。あと、cycle間で重複を許しているのでこの部分は改善する。

TODO

• コンバーター作成
• c1ccccc1 をc1cscc1に変える
• 極性基の導入、削除
• エステル化
• 芳香環に極性基の導入
• その他
• 重複チェック
• オブジェクトをsmilesにしてsqliteで管理する
• QSAR的なツールと絡めて、変異の頻度を自由に変えられるように。（例えばclogpがでかくなったら極性導入の頻度を高くしたり、MWTが大きくなりすぎたら、大胆に削るとか）
• wikiかなんかにちゃんとしたドキュメントとか書く。

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の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

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

知らんかった。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が出来上がるな。

plaggerで planet chemoinformaticsをつくってみた。

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法で解いていると。感覚的にはそうなんだろなぁと思うが、ここら辺の理解度はちょっとあやしい。

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+CDKとRpyを組み合わせて色々できると面白いナァと思う。

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

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

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

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

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

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