化学系だったら一度はコードを書いたことのあるヒュッケル法による軌道エネルギー計算はPython/Numpy/Pybelを使うと200行ぐらいで書けるヨというエントリをみつけた。

FORTRANだとこんな感じ。

smi23dっていう3D Coordinate Generatorが出てた。

おー素晴らしいよーマジで

後で触ってみよう。これが使えると色々できるようになりますな。

職場でいつものようにGUIをチマチマいじってたら、隣で塩を取り除くコマンドがないっちゅよと困っていた。openbabelにあるんじゃねーの？どれどれと探してみたけどなかった。さらに、いやあるだろふつーはとかいってググってもcgiみたいなもんしかみつからなかった。

smilesだったら.でsplitしてゴニョればOKなんだけど、座標情報落っこちるから、perlmolとかopenbabelで実装して任意のフォーマットにできるようなコマンドが欲しい。

それを仮にwashmolとすれば

washmol -i input.sdf [-o out.smi]

みたいに。デフォルトではインプットと同じフォーマットで返すような、調子いいコマンドが欲しくなった。

僕は塩を取り除くのはwashと言ってるんだけど、desaltとか色々と呼び方があんのね。

さて、bioinformaticsとかchemoinformaticsとか言われてるようなアレは、化学とか生物学に対して情報学的観点からアプローチしたりするわけです。こういう点から見ると、薬ってもんは蛋白質にうまいことはまる鍵みたいなもんで、蛋白質の穴にいい感じではまるような化合物を(コンピュータを駆使して)設計していくのが(コンピュテーショナルな)ドラッグデザイン(CADD)という分野だヨ。

で、蛋白質というものは複数のアミノ酸から構成され(100から数百)そしてアミノ酸は20種類程度存在し、それぞれ数十の原子から構成されているわけだ。要するに一対多の階層構造をとる。

protein -> amino-acid -> atom

みたいな。

でそれぞれprotein aminoacid(aa) atomみたいなクラスを用意すれば

$atom1 = Atom->new({name => 'C1', type => 'C'});
...
$aa1 = AminoAcid->new({name => 'GLY', atoms => [$atom1, $atom2, ]})
...
$protein = Protein->new({name=> 'ProteinA', aminoacids => [$aa1,$aa2,$aa3...]});

みたいにそれぞれ配列に突っ込めば蛋白質を表現できて、あるatomオブジェクトを与えられた場合にそれがどのアミノ酸に属しているのか調べるのに

for my $aminoacid ($protein->aminoacids){
    for my $atom ($aminoacid->atoms){
        return $aminoacid if $qatom == $atom;
    }
}

みたいにdepthfirstで探索していけばいいんだろうけど、ちょっと探索効率が悪いので、atomオブジェクトに $atom->{parent} = $parent_aminoacidみたいに親のアミノ酸オブジェクト返すような属性追加したんだけど、これだと構造が複雑になってなんか気持ち悪い。

他にうまいやり方ってあんのかなと思ったお盆の夏2007。

OSRAがええよとか言ってたわけだが、そもそも画像のメタデータに構造情報埋め込めばいいやん的発想なのがこのエントリ。

use Image::ExifTool qw(:Public);
use XXX;

$info = ImageInfo('rosiglitazone.png');

XXX $info;

で実行

$ wget http://depth-first.com/demo/20070801/rosiglitazone.png
$ perl etest.pl rosiglitazone.png 
---
BitDepth: 8
ColorType: RGB with Alpha
Compression: Deflate/Inflate
Directory: .
ExifToolVersion: 6.90
FileModifyDate: 2007:08:01 21:18:16
FileName: rosiglitazone.png
FileSize: 8 kB
FileType: PNG
Filter: Adaptive
ImageHeight: 109
ImageSize: 327x109
ImageWidth: 327
Interlace: Noninterlaced
MIMEType: image/png
PixelUnits: Unknown
PixelsPerUnitX: 1
PixelsPerUnitY: 1
molfile: |-
  name
  params
  comments
   25 27  0  0  0  0  0  0  0  0  0 V2000
      1.6910   -6.1636    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
      2.5571   -6.6636    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
      3.4231   -6.1636    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
      3.4231   -5.1636    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
      2.5571   -4.6636    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
      1.6910   -5.1636    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
      4.2891   -4.6636    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    ...
   23 19  1  0  0  0  0
   22 24  2  0  0  0  0
   20 25  2  0  0  0  0
  M  END
...
  at etest.pl line 6

論文のpdfとかもこういう感じで画像を埋め込んでくれればちょっとはデータが取り出しやすくていいかも。

LLで話されるような内容とは対極に近い、Rの統計解析環境としての本。かなりニッチなところをついてきた感があるが、Rのユーザー層を考えるとこっちのほうが普通の使い方なのかも。

R Commanderハンドブック―A Basic-Statistics GUI for R
舟尾 暢男
九天社 / ￥ 3,360 (2007-08)
通常4～5日以内に発送

RにおけるR Commanderってelispに対するemacsみたいな感じなのかな。

あとRってpythonみたいに対話環境がデフォルトだし。

なんつうか不思議な言語ですな。

PerlMolにはseparateメソッドがあるので、SMARTSで適当にマッチさせて結合を壊してからseparateすると化合物の構造を分割できる。

芳香環の炭素原子とヘテロ原子の間を切って、分割した構造をsmilesで出力するスクリプトはこう書ける。

use warnings;
use strict;
use Chemistry::File::SMARTS;
use Chemistry::File::SMILES;
use Chemistry::Ring 'aromatize_mol';

my $SMILES = shift;
my $react = Chemistry::Mol->parse($SMILES, format => 'smiles');
my $SMARTS = 'c-A';
my $pattern = Chemistry::Pattern->parse($SMARTS, format => 'smarts');

aromatize_mol($react);
molsplit($react);

sub molsplit {
 my $mol = shift;

 if($pattern->match($mol)){
   my @nbd = $pattern->bond_map;
   my @nfg = $pattern->atom_map;

   $nbd[0]->delete;
   for my $fg (@nfg){
     my $hcnt = $fg->implicit_hydrogens();
     $fg->implicit_hydrogens(++$hcnt);
   }

   for my $frag ($mol->separate){
     molsplit($frag);
   }
 }else{
   print $mol->print(format => 'smiles', unique => 1, name => 1),"\n";
 }
}

試しにpioglitazoneを分割してみると

$ perl molsplit.pl "CCC1=CN=C(C=C1)CCOC2=CC=C(C=C2)CC3C(=O)NC(=O)S3"
CC      
c1ccncc1        
CCO     
c1ccccc1        
CC1SC(=O)NC1=O

Mining Drug Spaceで知ったOSRAっていう化学構造認識ソフトがかなりよさげ。

そもそも、なんでコンピューターで認識できる形式でデータを流通させないのか?という本質的な問題はあるんだけど、この業界は画像(しかも二次元)ってのが標準なんで、こういう認識系のツールは重要だったり。

早速pubchemなんかの画像を認識させてみるときちんと認識される。OSRAのサイトではパテントをgifに変換したのを用意しているけど、1ページに複数の構造があっても別のものとしてきちんと認識する。ただ、フリーハンドで書いたのはベンゼンみたいな単純なものでも認識しなかった。

自分用に集めてるパテントとかジャーナルのpdfなんかを自動的に構造抽出してデータベースに突っ込むようにするとよいかも。

POV-Ray でレイトレできるとやっぱ調子いい。

Avogadro

構造活性相関(2D,3D-QSAR)とかで、アッセイ系の精度が気になって仕方がない。なんか生データはもちろんだし、さらに実験者の手技まで気になって数字をまともに受け取れない。

と同時に、構造も厳密に求めたくなり、量子化学計算でいちいち構造最適化しないと安心できない。

そんな感じで深みにはまっていった結果、一つの確信を得た。

活性の不確かさΔActを厳密に求めようとすると、構造の不確かさΔConfが増えていく、まだその逆も然り

これを、SARの不確定性原理と名付けてみた。

僕らは、予測モデルの精度を上げるために、厳密に活性コンフォメーションを求めたいし、ノイズのないアッセイデータを使いたい。でも常に、どちらかのノイズが邪魔をして満足のいく予測モデルはつくれない。

というわけで、QSARなんてばらつきを許容するような手法でやってるぐらいがちょうどいいのかもねーという結論に達し、同時にDrugDesignへのモチベーションがやや失せてきた。

という夢を見た。