chemicho2png

chemichoの出力をpngにしたのでPioglitazone で試してみた。

設定はこんな感じで。

Generation: 10
Cycle: 1
Prefix: Piogli

Plugins:
  - module: AddCarbon
    weight: 1
  - module: Aromatize
    weight: 0.8
  - module: Bulky
    weight: 1
  - module: DeAromatize
    weight: 0.2
  - module: DecreaseBondOrder
    weight: 1
  - module: GenerateRing
    weight: 0.5
  - module: Halogenize
    weight: 1
  - module: Hydroxylation
    weight: 2
  - module: HeteroCyclization
    weight: 2
  - module: IncreaseBondOrder
    weight: 2
  - module: IntroduceAcid
    weight: 1

6番目の出力から10番目の出力はこんな感じ。

piogli5 piogli6

piogli7 piogli8

piogli9 piogli10

あんまバルキーな変異はよくないな。もうちょいヘテロを導入したりとか、ヘテロ原子増やさずに自由度を減らすような変異をいれたほうがいいかも。芳香環を壊しちゃうのも気になるナァ。

ポスター出します

またもや、個人的ネタで恐縮ですが、今週末のスポッポファイアーのユーザー会でポスター出します(多分)。

しかも名前だけで発表はしないでうろうろします。単なる質問要員としてノミネート。

発表自体は、ITインフラ系構築ネタになる予定。明日から2人で作り始めます。Ajaxとかバリバリ取り入れてるので、そこらへんの話になるかなぁ。か、いわゆる(インフォマティクス的に)テックレイヤーの人達に如何にデータトラッキングの重要性とか、一元管理を浸透させたかとかそういう苦労話メインになるかなぁ。明日、オーバービュー書きながらどっちに転ぶか決めまする。

ホントは、スポッポファイア meets Ajaxとかいうネタがやれればよかったんだけど、手動かす暇なくて。僕の頭の中で腐りかけてル。

gSpan

テキストいじったり、化合物検索したりとグラフ構造いじることが多いのでgSpan をちゃんと覚えようと。

gspan-del.icio.us

DFS-CODEとか枝刈りのあたりとか断片的には理解できるけど、全体としてイメージできてないな。

ちゃんと理解するために、あとでperlで書いてみよう。で、 Chemruby使ってこんなこともやれるみたいなので、同じことをPerlMol使ってできないかなとか思った。gSpanPPで。

MT距離ってのがいまいちつかめてない

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のパッケージ読むの単なる食わず嫌いだったかも。思ったよりややこしいことはしてませんな、これからは積極的に読んだほうがいいかもと思った。

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

fp1 bitありfp1 bitなし
fp2 bitありCB
fp2 bitなしAD

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

昔作ったコードが気に入らなかったので書き直したら、そこそこ動くようになったので、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かなんかにちゃんとしたドキュメントとか書く。

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