pythonの組み込みウェブサーバーは便利で

import SimpleHTTPServer
SimpleHTTPServer.test()

という打鍵をpythonのインタラクティブシェルを起動してやっているのだけどPython - Swiss Army Knife Web Serverのコメントにあるような、

python -m SimpleHTTPServer

なんてできることを初めて知った。ということで、次の一行を.bashrcに書いておけば

alias pyhttpd='python -m SimpleHTTPServer'

いつでもpython組み込みのサーバーが起動して快適だ。例えば、画像突っ込んだディレクトリの内容なんかはさくっと見れていい感じ。

-mで呼び出すとモジュールの__main__が呼ばれるのね。モジュールの中身をみたら

if __name__ == '__main__':
    test()

って書いてあった。なるほど。(デーモンじゃないけど分かりやすいのでpyhttpdってコマンドで使うことにした。)

NMFを追っかけてたらMetagenes and molecular pattern discovery using matrix factorizationという論文を見つけたので、週末はこの論文を読みながら色々やってみた。NMFの便利なところは元の特徴(この論文の場合は遺伝子発現量)からより少ない任意の特徴量(論文中ではmetagene)に変換できるところであり、さらにそのままクラスターの分割に利用できる。

たとえば2つのmetageneで表現した場合、より発現量の大きいmetageneで分割すれば2つのクラスに分けられる。(QSARだったらdescriptorからmeta discriptorが導かれてそれに基づいてクラス分類ができるでしょう)

続いて、重要なのがクラスの安定性である。要するに最適なクラスタの数はいくつなのかということである。これに対して、この論文ではConsensus Clusteringというリサンプリングと隣接行列(connectivity matrix)を利用する方法をモディファイした方法を使っている。

ここで隣接行列はi番目のサンプルとj番目のサンプルが同じクラスなら1、それ以外なら0である。この行列のn回の平均値をコンセンサスマトリックスとする。コンセンサスマトリックスの値は0-1の間をとり、サンプルi,jが常に同じクラスになる場合は1、常に異なるクラスの場合は0である。フラフラするばあいはその間の値をとる。元の論文のconsensus clusteringアルゴリズムはデータの80%をランダムサンプリングして評価するのに対し、NMFの場合、初期値の行列はランダムな数字にしているので適当にループ(n)をまわすだけでよい。

NMFの実装は集合知プログラミングのものを用い、コンセンサスマトリックスを評価するコードをnumpyを使ってかいた。

import nmf
from numpy import *

def consensus(a,kstart,kend,nloop):
    """ calculate consensus matrix
    """
    (n,m)     = a.shape
    consensus = zeros((kend+1,m,m))
    conn      = zeros((m,m))

    #i = 0
    for j in range(kstart,kend+1):
        connac = zeros((m,m))

        for l in range(nloop):
            #i += 1
            (w,h) = nmf.factorize(a,pc=j)

            conn   = nmfconnectivity(h)
            connac = connac + conn

        consensus[j] = connac / float(nloop)
    return consensus

def nmfconnectivity(h):
    """ calculate connective matrix
    """
    (k,m) = h.shape

    ar = []
    for i in range(m):
        max_i = 0
        max_v = 0
        for index,v in enumerate(h[:,i]):
            if v > max_v :
                max_v = v
                max_i = index
        ar.append(max_i)

    mat1 = tile(matrix(ar),(m,1))
    mat2 = tile(matrix(ar).T,(1,m))

    return array(mat1 == mat2, dtype=int)

これをcc.pyという名前で保存しテスト用のセットを適当に用意して実行。

from numpy import *
from pylab import *
import cc

kstart = 2
kend = 4

testarray = array([[0,0,1,1,1,0,0,0,0],
                   [0,0,1,1,0,0,0,0,0],
                   [1,1,0,0,1,1,0,0,1],
                   [1,1,0,0,0,0,0,1,0],
                   [1,1,0,1,0,0,0,0,0],
                   [0,0,0,0,0,1,1,1,1],
                   [1,0,0,0,0,1,1,1,0],
                   [0,0,0,0,0,0,1,1,1],
                   [0,0,0,1,0,1,1,1,1]
                   ])

cons = cc.consensus(testarray.T,kstart,kend+1,20)


for i in range(kstart,kend+1):
    pcolormesh(cons[i])
    if i == kstart:
        colorbar()
    savefig('ccr' + str(i) + '.png')

見ればすぐ分かるが３つにクラスタリングできそうなマトリックス。

結果

2クラスタで分けた場合。

3クラスタで分けた場合

4クラスタで分けた場合

もうちょっと実際のデータでやってみないとあれだなぁ。それとConsensus ClusteringのCDFプロットみたいなのが欲しいところ。

参考書籍

集合知プログラミング
Toby Segaran
オライリージャパン / ￥ 3,570 ()
通常24時間以内に発送

参考論文

• Consensus Clustering
• Metagenes and molecular pattern discovery using matrix factorization
• matlabのコードなどを参考にした

ちょっと描きたいものがあったので調べたらpcolormeshを使えばいいだけだった。

from numpy import *
from pylab import *
a = rand(100,100)
pcolormesh(a)
colorbar()
savefig('colour.png')

NARのWeb Server issueを眺めてたらBioNMFというサービスをみつけた(SOAPでアクセスもできるようだ)。

Nonnegative matrix factorization(NMF)に関しては集合知プログラミングに書いてあるし、コード量も多くないので自分で書いたほうが自由に使えて良いような気もするが。

集合知プログラミング
Toby Segaran
オライリージャパン / ￥ 3,570 ()
通常24時間以内に発送

ところで、参考にあったpdfによると、ユークリッド距離の他にKullback-Leibler divergenceでの求めかたも書いてあったんだけど、これってどういう場合に使うといいのかいまいちわからん。

なんか具体的な適用事例はないのかなぁ。

Inline::Pythonを使うだけ。お手軽

use Inline Python => <<'END';
import pybel
def draw_png(smiles,file):
    mymol = pybel.readstring('smi',smiles)
    mymol.draw(filename=file, show=False)
    return True
END

draw_png("CCCc1ccccc1OC","xxx.png");

TODO: あとでクックブックに追加しておく

というエントリをケモインフォクックブックにのせておきました。

AMMOSというinduced fitを考慮したドッキングツールがopen sourceででるそうです。

AMMOS: Automated Molecular Mechanics Optimization tool for in silico Screening

ダウンロードできるようになったら試してみよう。DUDみたいなドッキングシミュレーションのベンチマーク用のデータベースもそろってきてるし、そろそろ今時のインターフェースで使いやすいドッキングシミュレーションツールのプロジェクトなど立ち上がりそうな気もしますが(もうなにか動いてるのかな)。

openbabelだったらobgenをつかえばいい。

pythonバインディングからやりたかったので、obgen.cppをよんだらOBBuilderクラスのインスタンスを生成すればいいらしいのだが、pythonからは現状OBBuilderのインスタンスをつくれないらしい。

MLを検索すると似たようなやりとりがあって、ラッパーを使えばいいらしく、pybelのソース見ろやーって書いてあったのでみた。

OPSってのは新たに導入された簡易プラグインシステムらしい。

で、ケモインフォクックブックにコードのせといた。

IPアドレスに応じて色を変えられるのであれば、アクセスログをLEDで光らせればいいんではなかろうか？と思ったので。

参考

• Identicon
• Image::Identicon
• ピンポン玉をLEDのディフューザーに
• PyGainer

Gainerのほうは秋月で買ったカソードのフルカラーLEDを刺した。抵抗は全て330Ω（それしかなかった）でピンポン球はダイソーかどっかで買った6個入り100円のヤツを錐とドライバーでぐりぐりした。

という単純な配線。

アクセスログは自宅サーバーに直接アクセスしてIPアドレスを返すようなCGIを用意した。 今回初めてPyGainerを使ってみた。

import sha,struct,urllib2,time,sys
from time import sleep
from PyGainer import PyGainer
#from random import randint

p = '/dev/cu.usbserial-A2002mCa'
g = PyGainer()
config = { 'baudrate' : 38400, 'timeout' : 5 }

#def rand_rgb():
#    return(randint(0,255),randint(0,255),randint(0,255))

def get_rgb():
    ip = urllib2.urlopen('http://myserver/taillog.cgi').read()

    print ''.join(ip.split('.'))

    s = sha.new(ip).digest()
    code = struct.unpack('5L',s)[0]

    blue  = int((code >> 6) & 0x0ff)
    green = int((code >> 15) & 0x0ff)
    red   = int((code >> 24) & 0x0ff)

    return (blue,green,red)

if g.open(p, config) == 0:
    if g.reset() == True:
        print "RESET"
    else:
        print "ERROR"

    sleep(0.1)
    g.version()
    g.configuration(1)

    sleep(0.1)

    for i in range(100):
        (blue,green,red) = get_rgb()
        g.set_specified_analog_output_port(0, green)     # G
        g.set_specified_analog_output_port(1, blue )     # B
        g.set_specified_analog_output_port(3, red  )     # R
        sleep(3)

    g.set_specified_analog_output_port(0,0)
    g.set_specified_analog_output_port(1,0)
    g.set_specified_analog_output_port(3,0)

    g.close()

動かした。うしろでカチャカチャうるさいのは犬が晩ご飯を催促している音。

Identiconでいうところのパターンのようなものの表現は点滅で実装しようかと思ったが、それって分かりにくいよなと。もう少し面白そうなデバイスはなかろうか。

bioinformaticsの学習道路はそこそこ充実している。

でも、chemoinformaticsだってオープンソースで色々できるし、独学でも結構いけるよ！なんて思うんだけど情報が足りないなぁと思った。

というわけで、ちょっと小道でも。

ケモインフォクックブック

特にperl,pythonでちょっとした化学情報をいじるコードをクックブック形式で書いていけたらなぁなんて思ってる。

週に一回くらいは更新するつもりで。