バイオ・ノート: 2008

2008年6月20日金曜日

go-perlサンプルコードを読んでみた

まずは、昨日走らせたサンプルルコードを読んでみる。


tomo$ less go-perl.pl
    #!/usr/bin/perl

    use GO::Parser;

    my $parser = new GO::Parser({handler=>'obj'}); # create parser object
    $parser->parse("gene_ontology.obo"); # parse file -> objects
    my $graph = $parser->handler->graph; # get L object
    my $term = $graph->get_term("GO:0001303"); # fetch a term by ID
    printf "Got term: %s %sn", $term->acc, $term->name;
    my $ancestor_terms = $graph->get_recursive_parent_terms($term->acc);
    foreach my $anc_term (@$ancestor_terms) {
      printf " Ancestor term: %s %sn", $anc_term->acc, $anc_term->name;
    }

GO::Parserというのがこのモジュールの名前らしい。


    my $parser = new GO::Parser({handler=>'obj'}); # create parser object
    $parser->parse("gene_ontology.obo"); # parse file -> objects

ここで、GO::parserオブジェクトをGO object modelをhandlerとして作成。そして、"gene_ontology.obo"ファイルをパースして、gene ontologyのデータをGO::Parserオブジェクトに取り込む。
handlerとして選択可能なモデルは、ここにあるGO::Handler::*すべてのようだ。
パース可能なファイルフォーマットは、こちらに記載されている。このサンプルでは、各Gene ontologyの階層関係を記述したgene_ontology.oboファイルを読み込んでいる。


my $graph = $parser->handler->graph; # get L object

ここでは、GO::Parserオブジェクトのhandlerメソッドを呼び出し、得られたGO::Handler::objのgraphメソッドを呼び出して、GO::Model::Graphオブジェクトを得ている。
このGO::Model::Graphオブジェクトに、ファイルから読み出した情報が格納されていることになる。
GO::Model::Graphオブジェクトのメソッドは、こちらに列記されている。


my $term = $graph->get_term("GO:0001303"); # fetch a term by ID

ここでは、GO::Model::Graphオブジェクトのget_termメソッドを使用して、GO IDからGO::Model::Termオブジェクトを獲得し、


printf "Got term: %s %sn", $term->acc, $term->name;

ここで、GO::Model::Termオブジェクトの要素をそれぞれ呼び出して表示する。

要するに、GO::Model::GraphオブジェクトとGO::Model::Termオブジェクトの内容を調べれば、一通り扱えそうだ。

go-perlのインストール

環境：
MacBook(Mac OS X 10.4.11)

go-perlの入手先
http://search.cpan.org/~cmungall/go-perl-0.09/

インストールはCPAN shellから行いました。

sudo perl -MCPAN -e shell
install GO::Parser

途中でData::Stagをインストールするかどうかを聞いてくるので、yesと答えます。
特にエラー無くインストールは終了。

動作確認には、go-perlのdocumentに記載されているサンプルコードを実行してみました。


    tomo$ less go-perl.pl
    #!/usr/bin/perl

    use GO::Parser;

    my $parser = new GO::Parser({handler=>'obj'}); # create parser object
    $parser->parse("gene_ontology.obo"); # parse file -> objects
    my $graph = $parser->handler->graph; # get L object
    my $term = $graph->get_term("GO:0001303"); # fetch a term by ID
    printf "Got term: %s %sn", $term->acc, $term->name;
    my $ancestor_terms = $graph->get_recursive_parent_terms($term->acc);
    foreach my $anc_term (@$ancestor_terms) {
      printf " Ancestor term: %s %sn", $anc_term->acc, $anc_term->name;
    }

こちらから、"gene_ontology.obo"ファイルをダウンロードして来て、上記go-perl.plと同じディレクトリに置いておきます。
http://www.geneontology.org/ontology/gene_ontology.obo

実行結果は以下の通り、


    tomo$ perl go-perl.pl
    Got term: GO:0001303 nucleolar fragmentation during replicative aging
     Ancestor term: GO:0001302 replicative cell aging
     Ancestor term: GO:0007569 cell aging
     Ancestor term: GO:0007568 aging
     Ancestor term: GO:0032502 developmental process
     Ancestor term: GO:0008150 biological_process
     Ancestor term: GO:0007576 nucleolar fragmentation
     Ancestor term: GO:0007000 nucleolus organization and biogenesis
     Ancestor term: GO:0006997 nuclear organization and biogenesis
     Ancestor term: GO:0006996 organelle organization and biogenesis
     Ancestor term: GO:0016043 cellular component organization and biogenesis
     Ancestor term: GO:0009987 cellular process
     Ancestor term: GO:0008150 biological_process
     Ancestor term: GO:0007569 cell aging
     Ancestor term: GO:0007568 aging
     Ancestor term: GO:0032502 developmental process
     Ancestor term: GO:0008150 biological_process

2008年6月1日日曜日

GOに属する遺伝子群の抽出(6)

結局、evidence codeがIEAのGOはすべて除外して登録してみた。
トランスクリプトとの対応で、10297レコード。だいぶ近づいて来た。多分トランスクリプトとタンパク質との重なりが少しあるのであろう。タンパク質コードとの対応表を作ってみる事にする。

GOに属する遺伝子群の抽出(5)

やはり、GOはGene Productに対してつけられているので、具体的にはタンパク質に対してつけられているはず。よって、トランスクリプトとの対応の方が正しいようだ。
いまは、トランスクリプトから引いたGOに対して、親を検索し、親まで含めたすべてのGOをトランスクリプトに関連づけるようにしているが、その際に元となったGOのエビデンスコードをそこから引いた親すべてに付与するようにしてテーブルに格納してみる。

GOに属する遺伝子群の抽出(4)

トランスクリプトでは無く遺伝子とGOを対応させてみたが、やはりレコード数が多い。
考えられることといえば、Humanのみのデータを抽出するか、evidence_codeのIEAを除外してみるかだな。もともとのトランスクリプトはHumanのものだから、それらからひいたGOの親も、Humanのものだと思われるし、そもそもGOのヒエラルキーの上部は、種間で共通にしてあるはず。
とりあえず、今度はevidence codeをテーブルに突っ込んでみて、SQLで検索し様子を見てみることとする。

2008年5月31日土曜日

GOに属する遺伝子群の抽出(3)

という事で、ENSEMBL GENE IDをキーとして、GOをテーブルに格納する事にした。そして、別にENSEMBL Transcript IDとENSEMBL Gene IDの対応をつけるテーブルを作成し、テーブルを結合してENSEMBL Transcript IDからGO番号を取り出すようにしてみる。

GOに属する遺伝子群の抽出(2)

データベース上でGO treeをたどるのは、非常に面倒だ。幸いGOデータベースには、graph_pathというテーブルがあり、GO番号を入れてやることにより簡単に親のGO番号を取得することができる。
そこで、EBSEMBL transcript IDとGO番号を対応付けたデータベースを以下のように作り直すことにした。

GOに関してはBiological Processのみを抽出
オリジナルのテーブルから引いてきたGO番号を使用して、graph_pathテーブルに検索をかけて親に当たるGO番号をすべて抽出
それらすべてをENSEMBLE transcript IDと対応付ける新しいデータベースを作成

これで、新しいデータベースには、親のGO番号を含めたすべてのGO番号がENSEMBL Transcript IDと対応付けられているはずなので、あとはENSEMBLE Transcript IDのリストを渡して、各GO番号の出現頻度を計算するスクリプトを書けばいいことになる。
オリジナルのテーブルは224679レコードだったが、果たして新しいテーブルはどの程度の大きさになるだろうか？

1時間くらいスクリプトが走りっぱなし（非力なマシンなんです）で、トータル586984レコード。心配していたほど巨大にはならなかった。ただ、確認のため、GO番号がGO:0008150でレコード数を検索してみると、なんと27380レコードもある。
完全にGOを網羅しているテーブルであれば、結果は9069になるはず。対応もとをトランスクリプトームにしたために、おそらく遺伝子で見るとredundantになっているのだろう。また、Biological Processだけを抽出したはずなので、distinct ENSEMBL Transcript IDで検索すると、27380にならなければいけないはずだが、何故か27646と返ってきた。んー。どうしたものか。

2008年5月30日金曜日

GOに属する遺伝子群の抽出

とりあえず、Gene Ontologyの解析から始めることにした。
単純に各GO numberに属する遺伝子群をnon redundantに抽出するのはそんなに難しくない。
が、どうやら登録されているGOというのは、ヒエラルキーの下の方のGOみたいだ。ということは、tree構造をたどっていって、それぞれ親に当たるGOに対しても処理をしていく必要がありそうだ。

Gene Ontologyデータベース

gene ontologyデータベースをダウンロードして、ローカルのmysqlに再構築したのは良いけれど、その構造がよくつかめていなかった。
http://www.geneontology.org/images/go-database-ER-diagram.png
ここに、データベースの構造図が載っているが、複雑なデータベース設計なんて見た事無いので、よく理解できなかった。しかし、
http://www.geneontology.org/GO.database.schema.shtml
ここの記述とにらめっこしながら、結合をたどっていくと、何だか少しずつ分かって来たようだ。
はじめは、termテーブルは、GOとその説明だけが登録されていると誤解していたのだが、relationship_type_idの記述などなどもtermテーブルに登録されていたようだ。

各データベースID間の対応表

今回は、もとのデータがENSEMBLのtranscript IDで与えられているので、それらからKEGG, NCBI, Uniprotへの対応表を作成した。
ENSEMBLからNCBIのRNA accession id, Protein accession id, Gene id, refseq idに関しては、ENSEMBLサイトのBiomartでそれぞれデータをダウンロードし、それをperlでローカルのデータベースに登録した。そのまま取り込んでもいいんだけど、重複したデータをチェックしながら取り込みたかったので、その辺をスクリプトにして処理をした。一応、これらの対応表は完成。ヒトのデータだけだけど、どの程度カバーしているのかは不明。一応レコードとしては3万ちょっとはあるけれど、トランスクリプトームと考えると、心もとないかな？そうでもないかな？

2008年5月27日火曜日

タンパク質相互作用データベース

Cytoscapeのネットワークを作成するには、タンパク質相互作用データベースからデータを引っ張ってくる必要がある。Reactome, intActなどなどあるが、intActのデータがCytoscapeに取り込みやすい形式のようだ。ただ、ここで出てくるタンパク質のIDは、Swiss-port emblのも。また対応表を作らないと行けないのか？
何かとデータベースが分散していると煩雑である。

intActのデータはこちら。
ftp://ftp.ebi.ac.uk/pub/databases/intact/current/psimitab/intact.zip

2008年5月24日土曜日

統計解析法

具体的に思うようにデータセットを得ることができたとして、どういった統計学的な指標をもって判断すればいいのかまだ良く分からない。ある一群の遺伝子が特定のgene ontologyに偏った分布をしていることを示すには、Fisher's exact testを用いている論文が多いようだ。

よく理解していないが、あるgene ontologyに含まれる遺伝子数と含まれない遺伝子数、対象の遺伝子群とそれ以外の遺伝子群で2x2表を作り、Fisher's exact testを行えばいいのであろうか？この場合の母集団は、すべての遺伝子数になるのかな？それとも、すべての遺伝子にontologyが付いているわけではないので、gene ontologyが付いている遺伝子のnon-redundantなサブセットになるのかな？ただ、これをgene ontologyごとに繰り返したら、いわゆる多重比較の問題になってしまう。しかし、gene ontologyすべてを分割表の一方に持ってきて、2 x n表を作ってしまうと、ひとつの遺伝子が複数のgene ontologyグループに属しているので、各列が独立とはいえない気もする。 http://gostat.wehi.edu.au/example.html ここのアルゴリズムの解析を読むと、対照とする遺伝子群は、今回の解析の場合、アノテーションがついているすべての遺伝子群になりそうです。で、多重比較の問題は、Benjamini and Hochberg correction controls the false discovery rateを使用しているものが多いようですね。

CytoscapeのBiNGOというアノテーション解析プラグインのデータを、Rのfishre.testの結果と比較したら、一致する結果を得たので、おそらくfisher.test自体のやり方は間違っていないようです。で、同様にRでenjamini and Hochberg correctionができるかどうか調べているところ。
http://sekhon.berkeley.edu/stats/html/p.adjust.html
これを使えばできそうな感じです。

2008年5月23日金曜日

生物学的ネットワークの解析

まえのエントリーで紹介した論文が引いていたレビュー論文
もうちょっと詳しいことが載っていたので、まとめて見ます。

ネットワークトポロジーを表現する4つの指標

Average degree (K)
Clustering coefficient (C)
Average path length (L)
diameter (D)

トポロジカル・ディストリビューションを表現する4つの指標

degree distribution p(k)
degree distribution of cluster coefficient C(k)
shortest path distribution SP(i)
topological coefficient distribution TC(k)

いずれも、RのiGraphで計算できそうです。

Shihua Zhang et al., “Discovering functions and revealing mechanisms at molecular level from biological networks,” Proteomics 7, no. 16 (August 2007): 2856-69, doi:10.1002/pmic.200700095.

Protein Protein interactionのトポロジカルな特徴を示す指標

degree
clustering coefficient :
betweenness centrality
closeness centrality

でそれぞれの説明。

degree

ノードから出ているエッジの数。degreeが高いノードは、ハブと呼ばれ、多くのほかのたんぱく質と機能的な結合があることを示す。

clustering coefficient

隣接したノード間で、考えられるエッジの数に対する、実際に存在するエッジの割合

betweenness centrality

ノード間を結ぶ最短経路のうち、特定のノードを通過する数
betweenness centralityが高いノードは、そのネットワークの「ボトルネック」と呼ばれる

closeness centrality

その他のすべてのノードへの最短経路の数の逆数

その他の指標

in-degree, out-degree

あるネットワーク内でサブネットワークがモジュールを形成しているかどうかの指標。in-degreeは、特定のノードがもつコネクションのうちサブネットワーク内のコネクションの数。out=degreeは、特定のノードがもつコネクションのうち、サブネットワークに含まれないコネクションの数。
この場合のコネクションとは、エッジのことかな？

Chun-Wei Hsu, Hsueh-Fen Juan, and Hsuan-Cheng Huang, “Characterization of microRNA-regulated protein-protein interaction network,” Proteomics 8, no. 10 (May 19, 2008): 1975-1979, doi:10.1002/pmic.200701004.

2008年5月22日木曜日

Cytoscapeチュートリアル（勝手日本語訳）-基本編-始めに

Cytoscape チュートリアル　はじめに

Cytoscapeは、生物学ネットワークデータを描写・解析するためのオープンソース・ソフトウェアです。Cytoscapeの基本的な機能には、グラフの自動レイアウト、発現や遺伝子オントロジーなど、そのほかのデータとネットワークデータの統合、ノードやエッジの属性に応じて視覚的な特徴づけを行うことなどがあります。くわえて、cytoscapeはオープンなプラグイン環境を提供しており、だれでも自由にプラグインを書いて機能を拡張したり、すでにあるサードパーティ製のプラグイン・ライブラリをメンテナンスしたりすることが出来ます。現在、文献検索やモジュール検索などのプラグインが利用できます。

このページでは、Cytoscapeとプラグインについてのチュートリアルを提供します。最初のパート（基本編）では、Cytoscapeの基本的機能について説明します。Java Web Startを用いると、お手持ちのコンピュータにCytoscapeをインストールすることなしに、このチュートリアルを実行することが出来ます。しかし、別のページで述べた様にJava Web Startにはセキュリティ上の問題があることには注意してください。Java Web Start はWindowsもしくはMacOSXであれば、インストールすれば動作するはずです。しかし、Linuxでは追加のセットアップが必要ですので、こちらをご覧ください。

後半のチュートリアル（発展編）では、特定のプラグインを用いたCytoscapeのより高度な使い方について説明します。すべてのプラグインは自由に利用可能ですが、そのほとんどが使用許諾が必要です。中には再配布が制限されているものもあります。よって、発展編のチュートリアルではJava Web Startは利用できません。そこで、あなたのお手持ちのコンピュータにCytoscapeおよびそのプラグインをインストールする方法を記載しています。

これらのチュートリアルと平行してCytoscapeマニュアルを参照することも出来ます。

以下のウェブブラウザーがこれらのチュートリアルではサポートされています。

Netscape 5.0 or later
Internet Explorer 5.0 or later
Mozilla 6.0 or later

ご意見・ご感想は、Cytoscape-discuss メーリングリスト（日本語ディスカッショングループ）までお寄せください。

Cytoscape チュートリアル　基本編

チュートリアル 1: 始めてみよう

チュートリアル 2: フィルタ & エディタ

チュートリアル 3: 外部データの取得

チュートリアル 4: 発現解析

ん～。ブラウザ名がなんとも懐かしいですね。

CytoscapeをJavaWebStart経由で実行した際のセキュリティに関する重要な情報

Cytoscape JavaWebStartセキュリティ情報

CytoscapeをJavaWebStart経由で実行した際のセキュリティに関する重要な情報

JavaWebStartのリンクをクリックすると、Cytoscapeがダウンロードされ、あなたのコンピュータ上で実行されます。この仕組みは、あなたのお手元のコンピュータで走るJava appletに近いものですが、重要な点で異なっているところがあります。

最も重要な点は、Cytosacpeは、その他の多くのウェブ・アプリケーションと異なり、FULL PERMISSIONSで動作するように設定されています。もともと、Cytoscapeはコマンドラインから実行するように書かれており、それゆえネットワークファイルを読み書き出来るなど、ローカル・コンピュータにアクセスする様に作られているからです。このことは、CytoscapeをJavaWebStart経由で実行する際にも、このプログラムをあなたのお手元のコンピュータにダウンロード、インストール、実行した場合と同様のセキュリティ上のリスクを伴うことを意味しています。もしセキュリティを心配するのであれば、このことが本当にあなたが行おうとしていることかどうか、もう一度確認してください。Cytoscapeを実行することで、あなたは制限と責任に関する免責事項を含む使用許諾に同意したものとみなされます。こちらから使用許諾のコピーを見ることが出来ます。

WebStartリンクをクリックすると、まずJavaWebStartは、アプリケーションを実行する前にセキュリティ上のリスクに関する警告メッセージを表示します。

続けてアプリケーションを起動する場合には"Start"を、中止する場合には"Exit"をクリックしてください。もし"Start"をクリックしてアプリケーションを起動した場合、今後同じアプリケーションを実行する場合には、あなたがお手元のコンピュータのアーカイブからこのアプリケーションを消去しなければ、このメッセージが再び出現することはありません。（下記参照）

Cytoscapeでは、ユーザーによる外部で開発されたプラグインの読み込みが可能です。これらのプラグインは、ローカルコンピュータへの読み書きを含めた、すべてのCytoscapeコア・プラットフォームの権限を継承します。それぞれのプラグインは、それぞれの使用許諾があります。Cytoscapeチームでは、コア・プラットフォームで開発配布されたプラグイン以外のいかなるプラグインに対してもその安全性と動作に関して保障するものでもありませんし、責任を追うものでもありません。

JavaWebStartは、アプリケーションをあなたのお手元のコンピュータ上に保存します（WebStartアーカイブ)。ですので、次回からはダウンロードなしで実行することが出来ます。また、Java Web Start executableそのものを実行することにより、ダウンロードしたアプリケーションを管理することが出来ます。これらは通常Javawsと呼ばれており、あなたのJavaもしくはJava Web Startアプリケーションフォルダにあるはずです。Java Web Startはまた、自動的にウェブ版Cytoscapeのアップデートをチェックしダウンロードします。Cytoscape Web Startアプリケーションは、インターネットへの接続が得られなければ動作しないでしょう。

通常のCytoscapeセッションで、設定情報を保存するためにあなたのローカルコンピュータに複数のファイルを保存します。vizmap.propsは、あなたのホームディレクトリに保存され、visual mappingに関する設定が記録されています。これらのファイルは、ユーザーへの確認を得ずに穂存されます。そして、これらのファイルを消去しても、Cytoscapeの機能には影響しません。

2008年5月19日月曜日

Cytoscape

最終的なデータの描写には、Cytoscapeを利用しようと思っています。
こちらは、開発者のお一人である大野さんの立ち上げられている日本語のサイト。
こちらは、本家サイト。

実は、ENSEMBL IDからGOやPathwayなどのあのテーションの付加はこのソフト内から簡単にできるのですが、100遺伝子からなるリストを500組一度に処理したいので、cytoscapeに取り込めるようなファイル形式で結果を出力するスクリプトを自分で書く事にしました。
cytoscapeでもバッチ処理みたいな事はできるようですし、Javaで独自に拡張もできそうですが、cytoscapeの内部構造の勉強をする必要が出てきますので、今回は見合わせようかと思っています。

2008年5月18日日曜日

ENSEMBL IDからrefseq IDへの変換

私が対象としている、遺伝子はENSEMBL IDでリストになっています。Gene Ontologyの情報は、ENSEMBLのperl APIを使えば取得できる事は確認できたのですが、KEGGからパスウェイデータを取得するためには、ENSEMBL IDでは検索できません。

KEGG側から見てみると、KEGGパスウェイを検索するには、KEGG gene IDが必要なようです。
で、KEGG gene IDとNCBI-IDの対応表は、KEGGのFTPサイトに転がっています。
NCBI-IDとrefsqe-IDとの対応表は、
ftp://ftp.ncbi.nih.gov:/gene/DATA/gene2refseq.gz
ここの転がっていて、
最後に、refseq-IDとENSEMBL-IDの対応は、
http://jp.youtube.com/watch?v=hu1Sr1kKWUI
ここに従って、Biomartで得る事ができます。

これで、理論的には、
ENSEMBL-ID -> refseq-ID -> NCBI-ID -> KEGG gene id -> KEGG pathway
と行けそうです。
ただ、この辺をいつもネットワーク越しでアクセスしていたら大変そうなので、ローカルのデータベースに突っ込むのが良さそうです。

なんか、もっと簡単な方法がありそうなものですが、、、。

追記：
ENSEMBL perl APIで遊んでいたら、実はtranscriptのDBentryにrefseqDNAの情報が入っている事が分かりました。ここでいう、refseqDNAはいわゆるMN_*の形式ですから、そのままKEGG gene idに変換してパスウェイの検索に使えそうです。

さらに追記：
KEGGに落ちていたrefseq_idとKEGG gene_idとの対応表にあるrefseq_idとは、どうやらprotein_giのようだ。

このブログの目的です

今更ですが、遺伝子オントロジー解析、パスウェイ解析を行おうと、勉強をはじめました。

100程度の遺伝子の組が500組程度あり、それぞれの遺伝子群がどのような特徴を共有しているのかを、調べてしまおうという、私にしてみれば壮大で無謀な計画です。

このブログでは、その過程で得られたTIPSについて、おもに私の健忘録代わりに、書いていこうと思っています。

登録: 投稿 (Atom)

2008年6月20日金曜日

2008年6月1日日曜日

2008年5月31日土曜日

2008年5月30日金曜日

2008年5月27日火曜日

2008年5月24日土曜日

統計解析法

2008年5月23日金曜日

生物学的ネットワークの解析

ネットワークトポロジーを表現する4つの指標

トポロジカル・ディストリビューションを表現する4つの指標

Protein Protein interactionのトポロジカルな特徴を示す指標

degree

clustering coefficient

betweenness centrality

closeness centrality

その他の指標

in-degree, out-degree

2008年5月22日木曜日

Cytoscape チュートリアル はじめに

Cytoscape チュートリアル 基本編

コメント

CytoscapeをJavaWebStart経由で実行した際のセキュリティに関する重要な情報

2008年5月19日月曜日

2008年5月18日日曜日

ブログ アーカイブ

自己紹介

Cytoscape チュートリアル　はじめに

Cytoscape チュートリアル　基本編

ブログアーカイブ