2011/12/01 21:29:27
mol形式(sdf形式)のデータだと化合物の区切りが$$$$なので、化合物を追加したい場合は何も考えずにファイルに追記するだけでいいのでよいですね。
PDB形式のデータにsdf形式の化合物情報をマージしたいんだけど、いい方法ないかなぁと調べてみたところmol2でOKだった。
両方mol2形式にして
cat compounds.mol2 >> protein_data.mol2
ってやればマージできる。
どういう用途を想定しているかっていうと
ある適当な部分構造(substructure)を持っている化合物の複合体結晶構造に、同じsubstructureをもつ別の化合物のコンフォマーを発生しつつ複合体のsubstructureの座標でalignする
つまり
- obconformerでコンフォマーを発生
- 複合体結晶構造のリガンドの部分構造を使ってobfitでコンフォマーをアライン
- 一つのファイルにまとめてドッキングモデル完成
みたいなことをやりたかったわけです。こういうのはファイルが2つに分かれてるとユーザーのヒトとか使いにくいしどういう計算したのかわからなくなっちゃうからね。
2011/11/24 19:32:42
biopythonのMLに「蛋白内部に埋没している残基をどうやってけいさんすんの？」っていう質問が流れてて、HSEっていう指標が実装されているのを知った。
HSEってのはCalphaとCbetaのベクトルと直交する平面で球を切ってUpとDownの半球のことで、その中に他の残基のCalphaとCbetaが幾つあるか数えるという単純なCNっていう指標で溶媒接触表面積の代わりに使えるらしい。
この指標って例えば（潜在的な）リガンド結合部位の予測に使えたりするんだろうか？
PPI阻害剤なんかのターゲット部位予測に使えたら面白いかもねと思った。
2011/05/13 20:36:20
pubmedのidからdoiを調べたい。
BeautifulSoupでXMLをパースするのが良いのだが、ソース見たらげんなりした(preってなんやねん)。というわけで、MEDLINE形式のデータから正規表現でdoiを抜き出してます。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21 | #!/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
2011/01/19 05:28:53
献本ありがとうございます
内容はバイオインフォマティクスに限らずに割と広い内容をカバーした感じで、クックブックと逆引きの中間的なスタイルと言えば良いのだろうか？
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があります。
2010/10/25 05:50:12
熱海でだらだらと温泉につかりながら、コード書いたり書かなかったりという場を設けました。キーワードはchemoinformatics,bioinformatics,R,Python,Rubyで製薬企業でコード書いているヒト多めです。僕はFlaskいじるかopenbabelのMCS実装かなんかをする予定です(でもいい日本酒がゲットできたら飲んでるかも)
- 2010.10.29-30
- 芳泉閣@熱海
- 一泊二日で10500くらい
まだ何人か分の空きがありますので興味があれば私にメールかtwitterでメッセージを下さい。金曜から土曜にかけてという日程です。
今のとこの参加者
- kzfm
- zgmfx20a
- bohohu
- garuby
- shirahakase
- hiro_h
2010/04/22 15:10:16
今日は子守で会社を休んでいるので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
こんなんだからなぁ。
なんなんだろうね。バグ？
2010/01/12 11:24:59
プログラムとして実行できる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の本を参照のこと。
VMを使った中間言語方式の強力さを理解した。
2010/01/12 11:23:48
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でガリガリ書くよっていうヒトはあんまし見かけたことないんですよね。
2010/01/12 11:23:21
Bioinformatics Programming in Python: A Practical Course for Beginners という書籍が出るらしい。
TOC見る限りはあんま欲しいとは思わないけど、他にもpythonでbioinformaticsという本が2冊くらいは出るらしいので、そっちのほうに期待。
2010/01/12 11:06:36
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()
ちょっと長くなったけど帯域にやさしい。
Bioinformatics Programming Using Python
Python for Bioinformatics (Chapman & Hall/CRC Mathematical & Computational Biology)
Rによるバイオインフォマティクスデータ解析 第2版 －Bioconductorを用いたゲノムスケールのデータマイニング－
Bioinformatics Programming Using Python
Rubyで作る奇妙なプログラミング言語 ~Esoteric Language~
Bioinformatics Programming in Python: A Practical Course for Beginners