今日は子守で会社を休んでいるのでPRMLのサンプリングの章でもじっくり読んでやるかという有意義そうな予定を組んでいたのだけど、なぜかmacbookにblast入れたりしていた。

しかもportで入れたblastが古くてあれやなぁと。調べたら今時はBLAST+(ブラプラって読むの？)らしくて、Biopythonのほうも対応しているらしいのでこっち入れなあかんやろと、何が目的だっけ？そもそもblast入れる目的なんだっけ的なよくあるパターンに。

NCBIのダウンロードサイトからMac OSX用のをダウンロードすればよい。

/usr/loca/ncbi/binにインストールされるのでパスを切っておく。

コマンドは引数も含めていろいろ変更されていて、blastallのpオプションで指定していたのがそのままコマンドになっている。formatdbがmakeblastdbになっていたりとか。

引数も短縮形じゃなくて、意味がわかるようなものに変更されている。

blastp -query test.fasta -db pdbaa

あと、ホームディレクトリでblast実行すると以下のエラーが出るんだけど、最初どこに問題があるんだかわからなかった。

  what():  NCBI C++ Exception:
    "/am/ncbiapdata/release/blast/src/2.2.23/IntelMAC-universal/c++/GCC401-
ReleaseMT--IntelMAC-universal/../src/objtools/blast/seqdb_reader
/seqdbimpl.cpp", line 412: Error: OID not found

Abort trap

で、dbのファイルを絶対パスで指定してやると解決した。

blastp -query test.fasta -db /Users/kzfm/blast/db/pdbaa

ただしデータベースが見つからないときのエラーは

BLAST Database error: No alias or index file found for protein database

こんなんだからなぁ。

なんなんだろうね。バグ？

Bioinformatics Programming Using Python
Mitchell L. Model
Oreilly & Associates Inc / 5075円 ( 2009-12-23 )

python2.5に入れる場合はNumericがsourceforgeになかったり、MxTextToolsの古いバージョンを入れたりとか、いろいろとあれなので、ドキュメント読みましょうってことです。

で、RCSBからgz圧縮されたpdbファイルをとりにいってBioPythonで扱うサンプル

import urllib2, StringIO, gzip
from Bio.PDB.PDBParser import PDBParser

def fetch_pdb(id):
    url = 'http://www.rcsb.org/pdb/files/%s.pdb.gz' % id
    content = urllib2.urlopen(url).read()
    sf = StringIO.StringIO(content)
    return gzip.GzipFile(fileobj=sf)

if __name__ == "__main__":
    p=PDBParser(PERMISSIVE=1)
    s=p.get_structure("1bgw", fetch_pdb("1bgw"))

    for model in s.get_list():
        for chain in model.get_list():
            for residue in chain.get_list():
                if residue.has_id("CA"):
                    ca_atom=residue["CA"]
                    print ca_atom.get_coord()

ところで、製薬系のドラッグデザインっていうかchemoinformaticsなSoftwareはpythonで拡張できるものが多いんだけれども、その割にはpythonでガリガリ書くよっていうヒトはあんまし見かけたことないんですよね。

Bioinformatics Programming in Python: A Practical Course for Beginners という書籍が出るらしい。

Bioinformatics Programming in Python: A Practical Course for Beginners
Ruediger-Marcus Flaig
Wiley-VCH / ￥ 7,525 ()
通常1～3週間以内に発送

TOC見る限りはあんま欲しいとは思わないけど、他にもpythonでbioinformaticsという本が2冊くらいは出るらしいので、そっちのほうに期待。

Fetching PDB files remotely in pure Python codeを見つけて、 数行でpdbのファイルをフェッチできるなんてやるなRCSBとpythonのコンボとか思ったんだけど、コメントに「生でフェッチすんのは環境に悪いでよ」とか書いてあった。

確かにPDB界ではgzip圧縮したpdbファイルをとってきてローカルで展開すんのが昔からのナラワシだよなと、試しにgzとかZをくっつけてwgetしてみるとその通りのファイルがダウンロードできた。をを。

で、おーこれは、数行スクリプトにgzipモジュールかませばいいだけなんじゃなかろうかと思て書き書きしてみたが、これってファイルしかとれないのかな?

zlibならどうだといじってみたけど、

>>> import urllib
>>> import zlib
>>> url = 'http://www.rcsb.org/pdb/files/1ab6.pdb.gz'
>>> text = urllib.urlopen(url).read()
>>> zlib.decompress(text)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
zlib.error: Error -3 while decompressing data: incorrect header check

なんかヘッダがーっておこられた。やっぱgzipでなんとかすべきなのか? ともうちっと追いかけたらGzipFileのほうを使えばよいらしいことに辿りついた。

>>> import gzip, StringIO, urllib
>>> url = 'http://www.rcsb.org/pdb/files/1ab6.pdb.gz'
>>> content = urllib.urlopen(url).read()
>>> sf = StringIO.StringIO(content)
>>> dec = gzip.GzipFile(fileobj=sf)
>>> data = dec.read()

ちと長いがこれで読めた。というわけで、

>>> import urllib, StringIO, gzip
>>> def fetch_pdb(id):
...     url = 'http://www.rcsb.org/pdb/files/%s.pdb.gz' % id
...     content = urllib.urlopen(url).read()
...     sf = StringIO.StringIO(content)
...     return gzip.GzipFile(fileobj=sf).read()

ちょっと長くなったけど帯域にやさしい。

最初サンプルみたときにはちょっと理解できなさそうだったので保留にしてたけど、bioperlとかpythonの話が面白そうだったのでポチッた。

ビューティフルコード
Brian Kernighan,Jon Bentley,まつもとゆきひろ
オライリージャパン / ￥ 3,990 ()
在庫あり。

そういえばpython magazineの最新号にbiopythonの特集が載ってたので買って読んだんだけど、なんていうか入門的な内容でいまいち楽しめなかった。

むしろpytableのほうが面白かった。pylabなんかと連携させて使うのがよいらしい。自分でオリジナルな計算かけたときのフォーマットとして使ってみようかな。

perl,python,Rでオームな解析をするための本。ツールやデータベースの説明が主なので、ハウツーな感じではなくて、インフォマティクス側からみたバイオロジーとかケミストリーのサービスとかツールのレビューに近い感じかも。

オープンソースで始めるゲノム・プロテオーム・メタボローム解析
樋口 千洋
オーム社 / ￥ 6,090 ()
在庫あり。

個人的には5章のケモインフォのとこと、9章のRの章が面白かった。

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

AMMOS: Automated Molecular Mechanics Optimization tool for in silico Screening

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

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 ()
在庫あり。

参考論文

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

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

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

集合知プログラミング
Toby Segaran
オライリージャパン / ￥ 3,570 ()
在庫あり。

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

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