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

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

1. GOに関してはBiological Processのみを抽出
2. オリジナルのテーブルから引いてきたGO番号を使用して、graph_pathテーブルに検索をかけて親に当たるGO番号をすべて抽出
3. それらすべてを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と返ってきた。んー。どうしたものか。

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万ちょっとはあるけれど、トランスクリプトームと考えると、心もとないかな？そうでもないかな？

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

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

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

統計解析法

統計解析法

具体的に思うようにデータセットを得ることができたとして、どういった統計学的な指標をもって判断すればいいのかまだ良く分からない。 ある一群の遺伝子が特定の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
これを使えばできそうな感じです。

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

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

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

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

1. Average degree (K)
   
2. Clustering coefficient (C)
3. Average path length (L)
4. diameter (D)
   

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

1. degree distribution p(k)
2. degree distribution of cluster coefficient C(k)
3. shortest path distribution SP(i)
4. topological coefficient distribution TC(k)

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

1. 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.