複数のパラメータが相互作用しあっている状態で最適なパラメーターセットを探索する際に、あまり精度の高いフィードバックが得られないと最適解の近傍の密度が高くなり、さらに準最適解があちこちに散らばるようなことがおこる。

それはPCAとかICAで空間を適当にとってあげることが出来て、特に第一主成分と第二主成分を平面にとってプロットしてやると、散布図として表現できる。

さてこれを、グラフにしてやればCytoscapeで扱えて便利だろうと思ったことがあったのでやり方を考えてみた。

アルゴリズム

1. 任意の点から非類似度Tの距離にある点の数をすべての点について数える
2. 最も点を含むものが多い点をrootとし、含まれる点との間で親子関係をつくる
3. 2の点を除いて1を行い新たな親子を決める。その際に親から非類似度T内にある既存の子のうち一番近いものを結ぶ（階層ができる）
4. これをすべての点がなくなるまで繰り返す。

例

以下の様な散布図を考える。

総当りで近傍の数をかぞえると、次の5つの点が含まれる空間が一番密度が高いことが分かる。

中心と含まれる点に対して親子関係をつくる(P1, C1-C4)。

これらを除いて、近傍の数をかぞえると次に密度が高いのが右側で3つの点を含んでいる。

これらに関しても先ほどと同様に親子関係をつくっておくが、新しくできた親に関しては既存の子とつなげられないかをチェックする。

するとC2 -> P2をつなげるのが良いことがわかる。

同様に行う（図の左側）。

これもC4 -> P3とつなげることができる。

最終的に全部つなげることができた。

実際に階層型の木を描くと以下のようになる。

単純だけど、全体の構造を保ったままCytoscapeで眺められていいかなーと思うんだけど。

もっといい方法あったら教えてもらえると嬉しいです。

まぜてもらうことになりました。

多様化する分子設計用ソフトウエア

僕の思うLBDDの行き着くところは、当たり前のように量子化学計算が行われていて軌道の相互作用として結合論が語られる未来だと思っているので、まぁそんな感じの話に持って行きたいなーと思っています。

FMOをやっても、表面的には相関エネルギーっていう数値しか出てこないけど、軌道の係数とか理解出来ないとドラッグデザインに生かせないし、ケミストもMOに対する直感的な理解とか必須だと思うしね。

みたいな話を、informaticsとかSBDDと絡めてお話できればいいと思うんだけど、突然心変わりしてHaskellの文脈でLBDDを捉えてみるとかPDCAサイクルはドラムンベースでいうところのベースラインそのものであるとかそういう謎のストーリーを作るかもしれない。

昨日は製薬業界の集まりがあって、他社のヒトと少し話す機会があって、LBDDどういう感じですかねと言われて、量子化学計算に真面目に取り組んでますよーっていう話から、ケミストは電子吸引基とか供与基とか言う割に計算して確かめようとしないんですよねー（そうですよねー）っていう流れになったので、そちらのケミストは計算するんですかねー？って聞いたら、するヒトはするし、しないヒトはしないっていう答えが返ってきた。想定通り。

もう少し知りたいのは、ちゃんと計算するケミストは、結局ただの計算マニアで結局普通のヒトなのか、それとも論理的で優秀な傾向が強いのか？ということかな。

さて、pygamessをpybelに対応させたので、より簡潔に書けるようになった。

>>> from pygamess import Gamess
>>> import pybel
>>> g = Gamess()
>>> g.run_type('optimize')
>>> mol = pybel.readstring('smi','C')
>>> mol.make3D()
>>> optimized_mol = g.run(mol)
>>> optimized_mol.energy
-37.0895866208

やっぱコンストラクタに引数をわたせるようにするべきだよなぁ。

初めてのPython 第3版
Mark Lutz
オライリージャパン / 4830円 ( 2009-02-26 )

量子化学計算したいならこれかな。

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

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

Gamessのpythonラッパーであるpygamessを0.2.0にアップデートした。

主な変更点

• pybelに対応した
• テストをpyvowsに移行した

pybelも使えるようになったので、よりpythonicにコードを書けるようになっていい感じ。 構造さえ用意すれば、こんなに簡単に量子化学計算ができて、量子化学計算にありがちなごちゃごちゃしたインプットファイルの作成から解放される。

>>> import pybel, pygamess
>>> g = pygamess.Gamess()
>>> mol = pybel.readfile("mol", "examples/ethane.mol").next()
>>> nmol = g.run(mol)
>>> nmol.energy
-78.30530748

それから、テストをpyvowsに変えてから快適です。自分にはBDDはあっているなぁ。

TODO

ドキュメントをきちんと書く

リリースしている安定版はインストールできなかったので、svnのtrunk(r3983)を入れた。それからbrewを使っているのでsetup.pyをちょっと変える必要があった。

/opt/localはmacports用の設定だと思うので/usr/X11に変更しないとGL/gl.hがないとかそんなエラーを吐くはず。

156         EXT = "/usr/X11"
157         inc_dirs=["ov/src",
158                   "layer0","layer1","layer2",
159                   "layer3","layer4","layer5",
160                   EXT+"/include",
161                   EXT+"/include/GL",
162                   EXT+"/include/freetype2",
163                   "modules/cealign/src",
164                   "modules/cealign/src/tnt",
165                   "generated/include",
166                   "generated/src",
167                   ]

それからPmwはソースをダウンロードしてきて入れた。なくても動くと思うが、import errorが出るので気持ち悪かったから後から入れておいた。

８章のマルコフモデルの実装までやった。ちなみに二単語のプレフィックスによるモデルを作ってた。ちなみにNodeで作る人工無脳は一単語ですね。こっちはリアルタイム化させているので面白いです。

9章はググらせるだけなのでこれでおしまいでいいやと。

• kzfm / koi-suru-program

恋するプログラム―Rubyでつくる人工無脳
秋山 智俊
毎日コミュニケーションズ / ?円 ( 2005-04 )

いい本だと思うんだけどなかなか復刊しないですね。

coffeescriptでディープコピーをやる方法を探していたらCoffeeScript Cookbookを見つけたので、Smooth CoffeeScriptとあわせて読んでみる予定。

化合物の合成ストーリーはケミスト言語による文章構築

ところで、人口無能はつくってて楽しいですね。とか書いていて、ふとケミストの合成の流れもマルコフ性があるんじゃないかなぁと思った。思ったというより確信に近い。プレフィックスに構造とイシュー(活性改善とか代謝よくするとか)を入れてやれば、次に修飾したがるところって容易に想像できるよね。実際ベースラインとしてのファーマコフォアをたもちながらっていうよりは、直近の結果にものすごい左右されやすいという性格があってだなぁ（略）。

例えば、ある文章中の単語をランダムに並べ換えて、それをもとの文章に戻すっていう問題を考えた場合マルコフ性をうまく利用すればどの位うまく解けるんだろうかね？

そういう風に考えると、特許化合物のリストっていうのはランダムな単語リストと同じになるわけで、合成のマルコフ性を利用して合成の流れを再構築しなおすってのは面白いかもしれんなぁと思った。

そのうちもう少し良く考えてみようっと。

pybel(openbabel)で立体構造の立ち上げを行うのは簡単だ。

import pybel
mol = pybel.readstring("smi", "c1ccccc1C")
mol.make3D()
mol.write("mol")

コードを眺めていたら構造をたちあげるメソッドは

• make3d
• localopt
• globalopt # ただし2.3.1ではコメントアウトされてた

make3dもglobaloptもlocaloptのstep数を変えて呼び出している。

localoptメソッドはデフォルトではmmff94パラメータを使って最急降下法で構造最適化をしている。

Perlでコードを書いていてCytoscape用にファイルを出力したかったんだがGraph::GMLは読み込み専用だし、Graph::XGMMLはCytoscapeじゃ読み込めないし、、、

Graph::XGMMLのソースコードを読んでみたら、XML::Writerを使えばいいみたいなので書いた。

sub write_xgmml {
  open my $output, '>', 'output.xgmml';
  my $xml  = XML::Writer->new(OUTPUT=>$output);
  $xml->xmlDecl('UTF-8');
  $xml->startTag('graph',
                directed=>"1",
                label=>"Sample1",
                'xmlns:dc'=>"http://purl.org/dc/elements/1.1/",
                'xmlns:xlink'=>"http://www.w3.org/1999/xlink",
                'xmlns:rdf'=>"http://www.w3.org/1999/02/22-rdf-syntax-ns#",
                'xmlns:cy'=>"http://www.cytoscape.org",
                'xmlns'=>"http://www.cs.rpi.edu/XGMML"
               );

  my $i = 1;
  my %dict;

  #add node
  foreach my $key ( keys %$ids ){
    $dict{$key} = $i;
    $xml->startTag('node',
                  id => $i,
                  label => $key,
                 );
    $xml->emptyTag('att',
                  type => 'string',
                  name => 'canonicalName',
                  value => $key
                 );

    $xml->emptyTag('att',
                  type => 'real',
                  name => 'property',
                  value => $ids->{$key}->{'property'} ? $ids->{$key}->{'property'} : 0.0
                 );

    $xml->endTag('node');
    $i++;
  }

  #add edge
  for my $e (@$spt) {
    $xml->startTag('edge',
                  label  => "$e->[0] to $e->[1]",
                  source => $dict{$e->[0]},
                  target => $dict{$e->[1]}
                 );
    $xml->emptyTag('att',
                  type=>'string',
                  name=>'canonicalName',
                  value=> "$e->[0] to $e->[1]"
                 );
    $xml->endTag('edge');
  }

  $xml->endTag('graph');
  $xml->end;
  close($output);
}

基本的にstartTagとendTagで包んでemptyTagで属性を入れてくだけですね。あとCytoscapeのIDって整数じゃないとダメだったような気がするので、そうなるように書いたんだけど実際どうだったかは忘れた。

アジャイル開発を創薬系にどう取り込んでいくかっていう観点で。

アジャイルサムライ−達人開発者への道−
Jonathan Rasmusson
オーム社 / 2730円 ( 2011-07-16 )

第V部のアジャイルなプログラミングに関しては昔考えたことがあるので、II−IV部あたりを。III部の計画づくりはうまく取り込むのが難しそうな気がするけどIV部の運営は色々参考になる。

まず、ソフトウェア開発と初期創薬開発の似ている部分異なる部分を書いておくが、特に異なる部分が重要。

異なる点

• 工学的にコントロールできるわけではなく、発見に頼る部分がある。つまり作ってみて評価してみないとわからない部分が大きい
• 成功したかしないかの二択。9割完成したとかはない（こういうのは失敗とみなされる）
• ベロシティがうまく見積もれない。ブレークスルードリブンなので非連続の進捗になりやすい
• あっちを動かすとこっちがおかしくなるというなかでバランスをとリながらベストを探るという多次元最適化戦略になることが多い

似ている点

• 何人もの異なるエキスパートの協調作業を強いられる。役割はきちんよ分かれているが、優秀な人材は幾つかのフィールドをまたぐことができる（インフラもフロントエンドもできるプログラマみたいな感じ）
• 自分のジョブしかできないヒトはスコープが狭いし、メタな視点に立てない。学習意欲のない人の生産性が低いもの一緒

みんなをバスに乗せる

スタートを切る前からだめになってしまうプロジェクトの主な理由は次の二点

• 答えるべき問いに答えられない
• 手強い質問をする勇気を持てない

まぁ、創薬プロジェクトでもよくありますね、認識しておかなければいけない問題を敢えて目をつぶって進めてた結果、激ハマリっていうパターンが。

手強い質問を先にすませてしまうためにインセプションデッキというツールが使える。

インセプションデッキは10の手ごわい質問と問題から構成されており、いずれの課題もプロジェクトを開始する前に聞いておかないとまずい質問ばかりだ

エレベーターピッチ

効能

• 明快になる
• チームの意識を顧客に向けさせる
• 核心を捉える

テンプレート

• [潜在的なニーズを満たしたり抱えている課題を解決したい]したい
• [対象顧客]向けの、
• [プロダクト名]というプロダクトは、
• [プロダクトのカテゴリ]です。
• これは[重要な利点、対価に見合う説得力のある理由]ができ、
• [代替手段の最右翼]とは違って、
• [差別化の決定的な特徴]が備わっている。

やらないことリスト

「やる」、「やらない」、「あとで決める」をきちんと見える化しておくツール。

創薬プロジェクトで「やらない」にタスクを入れるのはなかなか難しいかも。なので、「やる」と「後で決める」を明確に分けられればいいかな。

アジャイルなプロジェクト運営

ジャストインタイム分析の利点

• 最新かつ最も充実した情報に基づいて分析できる
• プロジェクトが進むにつれて分析がうまくなっていく
• 手戻りが大量に発生することを避けられる

最後の手戻りを避けられるのは大きな利点だな。それからリリースボードとストーリーボードは便利そうだ。

アジャイルな計画づくり

これが一番納得出来なかった。ベロシティが見積もれないようなプロジェクトだと与えられた期間を頑張るしかないんじゃないかなぁと思う。

それから名前が付いてるのかどうかわからんけど、探索的なプロジェクトだと締め切り間際になると集中力を発揮して何とかやり切るパワーみたいなのってありますよね。そういうのも重要だと思うんだよね。

というわけで、この部に関してはモヤモヤ感が残っているのでまた後で書くかもしれない。

2.3.1がリリースされたようです。

個人的に興味があるのは

• PNG files from Open Babel contain molecular information and can be read to give the MDL Molfile.
• Pybel now uses the built-in 2D depiction, and no longer needs OASA.

とABINITのフォーマットに対応したあたりかな。

あと、openbabel-python.iをいじってたので、ここをいじった場合のコンパイルのオプションをメモっておく。swigが有効になるようにしないといけないのに気付かなくてハマった。

cmake ../openbabel-2.3.1 -DPYTHON_BINDINGS=ON -DEIGEN2_INCLUDE_DIR=/usr/local/tmp/eigen-eigen-2.0.12 -DRUN_SWIG=ON

OBGenericからOBOrbitalDataへのキャストをできるようにして、vectorの設定もしたので、手元のpythonバインディングでは

orb = toOrbitalData(mol.GetData(openbabel.ElectronicData))
orb.GetAlphaOrbitals()[orb.GetAlphaHOMO()-1].GetEnergy()

とやるとHOMOのエネルギー(eV)を得られるようになっている。

追記12.01.28

homebrewでいれたpythonで使いたい場合optionで指示する

$ cmake ../openbabel-2.3.1 -DPYTHON_BINDINGS=ON \
  -DPYTHON_LIBRARY=/usr/local/lib/libpython2.7.dylib -DPYTHON_EXECUTABLE=/usr/local/bin/python \
  -DEIGEN2_INCLUDE_DIR=/Users/kzfm/openbabel/eigen-eigen-2.0.17 -DRUN_SWIG=ON