bioinformatics / Drkcore

PDBデータに化合物情報を追記して一つのファイルにしたい

Thu, 01 Dec 2011 21:29:27 +0900

mol形式(sdf形式)のデータだと化合物の区切りが$$$$なので、化合物を追加したい場合は何も考えずにファイルに追記するだけでいいのでよいですね。

PDB形式のデータにsdf形式の化合物情報をマージしたいんだけど、いい方法ないかなぁと調べてみたところmol2でOKだった。

両方mol2形式にして

cat compounds.mol2 >> protein_data.mol2

ってやればマージできる。

どういう用途を想定しているかっていうと

ある適当な部分構造(substructure)を持っている化合物の複合体結晶構造に、同じsubstructureをもつ別の化合物のコンフォマーを発生しつつ複合体のsubstructureの座標でalignする

つまり

obconformerでコンフォマーを発生
複合体結晶構造のリガンドの部分構造を使ってobfitでコンフォマーをアライン
一つのファイルにまとめてドッキングモデル完成

みたいなことをやりたかったわけです。こういうのはファイルが2つに分かれてるとユーザーのヒトとか使いにくいしどういう計算したのかわからなくなっちゃうからね。

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

Half Sphere Exposureという指標

Thu, 24 Nov 2011 19:32:42 +0900

biopythonのMLに「蛋白内部に埋没している残基をどうやってけいさんすんの？」っていう質問が流れてて、HSEっていう指標が実装されているのを知った。

HSEってのはCalphaとCbetaのベクトルと直交する平面で球を切ってUpとDownの半球のことで、その中に他の残基のCalphaとCbetaが幾つあるか数えるという単純なCNっていう指標で溶媒接触表面積の代わりに使えるらしい。

この指標って例えば（潜在的な）リガンド結合部位の予測に使えたりするんだろうか？

PPI阻害剤なんかのターゲット部位予測に使えたら面白いかもねと思った。

Python for Bioinformatics (Chapman & Hall/CRC Mathematical & Computational Biology)
Sebastian Bassi
Chapman and Hall/CRC / 5857円 ( 2009-10-07 )

pubmedのidをdoiに変換する

Fri, 13 May 2011 20:36:20 +0900

pubmedのidからdoiを調べたい。

BeautifulSoupでXMLをパースするのが良いのだが、ソース見たらげんなりした(preってなんやねん)。というわけで、MEDLINE形式のデータから正規表現でdoiを抜き出してます。

#!/usr/bin/env python

import urllib2
import re,sys

def pmid2doi(pmid):
   url = "http://www.ncbi.nlm.nih.gov/pubmed/%s?dopt=MEDLINE" % pmid
   r = re.compile('AID - (10.\d+/.+?) \[doi\]')
   response = urllib2.urlopen(url)
   if response.code == 200:
       s = response.read()
       m = r.search(s)
       return m.group(1)
   else:
       return "error: %d" % response.code

if __name__ == '__main__':
   if len(sys.argv) == 2:
       print pmid2doi(sys.argv[1])
   else:
       print "usage: %s [pmid]" % sys.argv[0]

コマンドラインから使う場合には

$ pmid2doi.py 20053000
10.1021/ci900416a

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

Rによるバイオインフォマティクスデータ解析　第2版 -Bioconductorを用いたゲノムスケールのデータマイニング

Wed, 19 Jan 2011 05:28:53 +0900

献本ありがとうございます

Rによるバイオインフォマティクスデータ解析　第2版－Bioconductorを用いたゲノムスケールのデータマイニング－
樋口千洋
共立出版 / 4830円 ( 2011-01-08 )

内容はバイオインフォマティクスに限らずに割と広い内容をカバーした感じで、クックブックと逆引きの中間的なスタイルと言えば良いのだろうか？

Rのインストールから基本的な操作は（大体どの本にもあるように）載っていて

データマイニングとしては

PCA
ICA
PLS
MDS
SPE
k-means,Fuzzy cmeans
spectral clustering
NMF
SOM
decision tree
kNN
SVM
RF
LASSO
MARS

がサンプルコードとともに簡潔に説明されている。

8章はバイオ系データの解析、チップとか。odesolveを利用したシミュレーションのサンプルもあって、SBMLRは面白そうだなぁと思った。メカニズムがどうなっているのかはモデルと実験系の不一致をよく突き詰めて考えることでしかきちんとした理解は得られないと思っている。

最後のほうの章は統合環境、データベースの連携、サーバー構築あたりの話。

あと、twitteRを使って生体組織名をつぶやくとその組織の図が返ってくるという例が載ってた。ちなみにRでもPitがあります。

芳泉閣でゆったりしませんか？

Mon, 25 Oct 2010 05:50:12 +0900

熱海でだらだらと温泉につかりながら、コード書いたり書かなかったりという場を設けました。キーワードはchemoinformatics,bioinformatics,R,Python,Rubyで製薬企業でコード書いているヒト多めです。僕はFlaskいじるかopenbabelのMCS実装かなんかをする予定です(でもいい日本酒がゲットできたら飲んでるかも)

2010.10.29-30
芳泉閣@熱海
一泊二日で10500くらい

まだ何人か分の空きがありますので興味があれば私にメールかtwitterでメッセージを下さい。金曜から土曜にかけてという日程です。

今のとこの参加者

kzfm
zgmfx20a
bohohu
garuby
shirahakase
hiro_h

MacBookにBLAST+をいれた

Thu, 22 Apr 2010 15:10:16 +0900

今日は子守で会社を休んでいるので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実行すると以下のエラーが出るんだけど、最初どこに問題があるんだかわからなかった。

terminate called after throwing an instance of 'ncbi::CSeqDBException'
  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 / ￥ 5,915 ()
在庫あり。

ゲノム言語ATGC

Tue, 12 Jan 2010 11:24:59 +0900

プログラムとして実行できるfasta形式のプログラミング言語を作ってみた。いちおうチューリング完全(なはず)。

>HELLO_WORLD
ccggaccgcg gggcaccgcc ggcggaccgc cgccggaccg cccggcgacc gccgccccac
cgcccgccca ccgcggggga ccgccgcccc accgccgccg gaccgccgcc ggaccgccgg
cgcaccgcgg cgggaaccga cacatccata ccacagaacc caaaa

これはatgcというコマンドで解釈して実行します。

$ bin/atgc hello_world.fasta 
Hello world!

ゲノム的にGCリッチなほうがいいだろうということで0と1にg,cをそれぞれ割り当てて数字を表現するようにしてる。exitコマンドには終止コドンを割り当てようかとも思ったが、なんとなくaaaにしてみた(polyA)。

しかも（というか当たり前だけど）blastでホモロジーサーチがかけられるし、multifastaにしておけばソース管理もできるうえに、データベース化してインデックスはっておけば、NCBIのツール群でコマンド一発で取り出せる。

ただ今回作ったHELLO_WORLDの配列はblastnだといい感じにヒットしなくて悲しかったので、blastxかけたらブラックコットンウッドからなんかひっかかった。

>gb|ABK94795.1|  unknown [Populus trichocarpa]
Length=229

 Score = 33.5 bits (75),  Expect = 5.9
 Identities = 14/21 (66%), Positives = 15/21 (71%), Gaps = 0/21 (0%)
 Frame = +1

Query  100  GPPPDRRRTAAGTDTSIPQNP  162
            GPPPDRRRT  GT  S P +P
Sbjct  209  GPPPDRRRTRQGTTKSEPASP  229

コンパイラ側のソースはこんな感じでほとんどWhitespaceだ。

VMとかは特にいじってないのでEsotericの本を参照のこと。

Rubyで作る奇妙なプログラミング言語 ~Esoteric Language~
原悠
毎日コミュニケーションズ / ￥ 2,814 ()
在庫あり。

VMを使った中間言語方式の強力さを理解した。

python2.5でBioPythonを使う

Tue, 12 Jan 2010 11:23:48 +0900

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

Tue, 12 Jan 2010 11:23:21 +0900

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冊くらいは出るらしいので、そっちのほうに期待。

pythonでHTTPのgzipデータを読み込む

Tue, 12 Jan 2010 11:06:36 +0900

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()

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