Bioconductorから、
goTools
rat2302
topGO
の各パッケージをダウンロード、インストールします。Rを起動して、コンソールに以下をコピーペースト（最初の一回のみやればOK）。
#------ ここから ------
source("http://bioconductor.org/biocLite.R")#おまじない
biocLite("goTools")#goToolsパッケージのインストール
biocLite("hgu133a")#hgu133aパッケージのインストール
biocLite("rat2302")#rat2302パッケージのインストール
biocLite("topGO", dependencies=TRUE)#'topGO' およびその他必要となる‘XML’, ‘RBGL’, ‘graph’, ‘SparseM’, ‘ALL’,
#‘hgu95av2.db’, ‘multtest’, ‘Rgraphviz’, ‘hgu133a.db’の9つのパッケージがインストールされます
#------ ここまで ------
パイチャート（パッケージ：goTools）
抽出した遺伝子セットの中にある遺伝子の機能をざっと分類して、パイチャートで表したい。
-goToolsパッケージに含まれているサンプルデータを使って練習-
1. 準備
Rのコンソールに以下のコマンドを打ち込む（コピーペースト）と、パッケージ"goTools"とサンプルデータが読み込まれます。
#------ ここから ------
library(goTools) #パッケージの読み込み
data(probeID) #サンプルデータ（probeID）の読み込み
#------ ここまで ------
なお、ここで使うontoPlot()関数は、環境によってはうまく動作しなかったり、計算機パワーの問題で時間がかかったりすることがあります。そんな時のために、計算結果を格納したオブジェクトを用意しました。
res
res2
res3
res4
res5
以上をこの項で練習を行う前に、上記各リンクをクリックし、「保存」で作業ディレクトリ（デスクトップなど）にダウンロードしておいてください。
パッケージに含まれているサンプルデータの場合は作業ディレクトリは指定しなくてもいいですが、後半の自前データを想定した項目ではあらかじめ作業ディレクトリを指定する必要があります。この段階で例えば「デスクトップ」などを作業ディレクトリにしておき、上記オブジェクトをダウンロードしておくとよいと思います。
次に、コンソールに
#------ ここから ------
ls() #読み込まれているオブジェクトを表示
#------ ここまで ------
と入力（コピーペースト）します。
[1] "affylist" "operonlist"
と表示されていればパッケージのインストール、読み込み、サンプルデータの読み込み全てがうまくいっています。
実際に自分のデータに応用するためには、サンプルデータがどのような形式で記述されているか
（どのようなオブジェクトとしてRに読み込まれているか）を知る必要があります。
"affylist" "operonlist"は「リスト（ベクトルや行列など異なる構造のデータを集めて 1 個のオブジェクトにしたもの）」です。
データを眺めてみるには、
#------ ここから ------
affylist #"affylist"の内容を表示
#------ ここまで ------
と入力（コピーペースト）します。こんな感じ(affylist.txt)に出力されるはずです。L1、L2、L3はそれぞれ100個、90個、70個のAffy_IDからなる、list形式のデータ（オブジェクト）であることがわかります。この後に出てくる、"ontoCompare()"という関数は、入力にlist形式のオブジェクトが必要なようです。listは、affylistの中身であるL1、L2、L3のように異なる長さのベクトルを同時に扱えるので便利です。
2. サンプルデータを、ontoCompare()関数を使ってデフォルトでセットされているGO Term毎に棒グラフ表示
#------ ここから ------
par(mar = c(10, 4, 5, 3)) #グラフを描画する領域を指定
res <- ontoCompare(affylist, probeType = "hgu133a", plot = FALSE) #"affylist"の内容を読み込んで、デフォルトでセットされているEnd Nodeに属するAffy_IDの割合（L1の場合は100個、L2は90個、L3は70個をそれぞれ1とした割合）を棒グラフで表示(plot = FALSEとしているので表示されない）
ontoPlot(res) #棒グラフを表示
#------ ここまで ------
いざという時のためにオブジェクトを用意しておきました。上記にあまりに時間がかかる、とにかく結果を見たい、というときには以下をコピペしてください。
#------ ここから ------
load("090814_affylist_res") #あらかじめ計算しておいたオブジェクトres2の読み込み
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(mar = c(10, 4, 5, 3)) #グラフを描画する領域を指定
ontoPlot(res) #棒グラフを表示
#------ ここまで ------
3. サンプルデータを、デフォルトでセットされているGO Term毎にontoCompare()関数を使ってパイチャート表示
#------ ここから ------
res2 <- ontoCompare(affylist["L1"], probeType = "hgu133a", method = "TIDS", plot = FALSE) #"Affylist"の中の"L1"のみの内容を読み込んで、デフォルトでセットされているEnd Nodeに属するAffy_IDの割合をres2に格納(plot = FALSEなのでグラフは表示されない）
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(mar = c(5, 8, 5, 8)) #グラフを描画する領域を指定
ontoPlot(res2, cex = 0.7) #パイチャートで表示（読み込んだlistの中身がベクトル1つだけの場合はパイチャート、2つ以上のベクトルの場合は棒グラフで表示される。cexは文字サイズの指定で、この場合は0.7）
#------ ここまで ------
いざという時のためにオブジェクトを用意しておきました。上記にあまりに時間がかかる、とにかく結果を見たい、というときには以下をコピペしてください。
#------ ここから ------
load("090814_affylist_L1_res2") #あらかじめ計算しておいたオブジェクトres2の読み込み
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(mar = c(1, 2, 1, 2)) #グラフを描画する領域を指定
ontoPlot(res2, cex=0.7) #パイチャートで表示（文字サイズ0.7）
#------ ここまで ------
goToolsのPDFマニュアル
-テキストファイル形式のサンプルデータを使って練習-
1. 準備
サンプルデータ1 (BATtop300.txt)
サンプルデータ2 (WATtop300.txt)
サンプルデータ3 (Livertop300.txt)
2008年の私の論文より。それぞれ24時間絶食ラットの褐色脂肪組織、白色脂肪組織、肝臓において、コントロールに比べて発現上昇した上位300遺伝子のリストです。
バックグラウンドデータ1 (all.txt)
Rat Genome 230 2.0に搭載されている全probe set IDのリストです。
以上をこの項で練習を行う前に、作業ディレクトリに(ファイルを右クリックして「対象をファイルに保存」で)ダウンロードしておいてください。
2. サンプルデータ1を、デフォルトでセットされているGO Term毎にontoCompare()関数を使ってパイチャート表示
実際に自分のデータに応用するためには、データをリスト形式で読み込む必要があります。私はここで結構苦労しました。まず失敗した例を示します。最初は、とりあえずread.table()関数でファイルを読み込んで、list形式に後付けで変換してやればよいと思い、以下のように
#------ ここから ------
BAT <- as.list(NULL) #BAT(list形式）の初期化
data <- read.table("BATtop300.txt", header = T) #データの読み込み
BAT <- as.list(data) #データの型をlistに変換、BATに格納
BAT #BATの中身を表示
page(BAT) #別ページでBATを表示
#------ ここまで ------
やろうとして見事にこけました。読み込めていそうなのですが、中身を上記のaffylistやoperonlistと比較すると微妙に違う（"Affy_ID"が""でくくられていない、要素が文字として認識されていない）ことがわかります。これでontoCompare()を使おうとすると、エラーが出て走りません。
いろいろ調べていくうちに、データの読み込みを別の方法でやるとうまくいくことに気がつきました（Webサイト「R-Tips」の「46.落穂ひろい」を参照しました）。
#------ ここから ------
BAT <- scan("BATtop300.txt", list(Probe_ID=""),skip=1) #データの読み込み、listとして読み込んで、"Probe_ID"を列ラベル、1行目（要素としての"Probe_ID"）は読み込まない
BAT #BATの中身を表示
#------ ここまで ------
これだと、affylistなどの形式と同じように読み込めているように見えます。サンプルデータ1の中身はベクトルが一つだけですので、ontoCompare()関数によってパイチャートが出力されるはずです。パイチャートを描かせてみましょう。ただし、以下はパソコンのスペックによっては非常に時間がかかる可能性があります。私の現在使っているiMac (3.06 Ghz, Intel Core 2 Duo, 4GBメモリ, OSX 10.5.8）だと十分強で計算が終わりますが、以前使用していたWindowsマシン だと一晩かかりました。
#------ ここから ------
res3 <- ontoCompare(BAT["Probe_ID"], probeType = "rat2302", goType="BP") #"BAT"の中の"Probe_ID"のみの内容を読み込んで、GOのBiological Process 以下のEnd Nodeに属するAffy_IDの割合をresに格納(明示的に指定しないと plot = FALSEなのでグラフは表示されない）
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(mar = c(1, 2, 1, 2)) #グラフを描画する領域を指定
ontoPlot(res3, cex=0.9) #パイチャートで表示（文字サイズ0.9）
#------ ここまで ------
いざという時のためにオブジェクトを用意しておきました。上記にあまりに時間がかかる、とにかく結果を見たい、というときには以下をコピペしてください。
#------ ここから ------
load("090814_BATtop300_res3") #あらかじめ計算しておいたオブジェクトres3の読み込み
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(mar = c(1, 2, 1, 2)) #グラフを描画する領域を指定
ontoPlot(res3, cex=0.9) #パイチャートで表示（文字サイズ0.9）
#------ ここまで ------
3. 複数サンプルデータを、カスタマイズしたGO Term毎にontoCompare()関数を使って棒グラフ表示（1）
もちろん、自分の得たあるデータセットについてどのような機能を持った遺伝子がどのくらい含まれているかについて、パイチャートで表すのもよいのですが、他のデータセットと比較したい、あるいはバックグラウンドデータと比較したい、なんていうこともあると思います。また、デフォルトの分類だけでなく、自分である程度機能分類を指定したい、ということもあると思います。そんなケースにも対応可能です。やってみましょう。
#------ ここから ------
hoge <- scan("all.txt", list(all=""),skip=1) #バックグラウンドデータの読み込み
hoge2 <- scan("Livertop300.txt", list(top_300=""),skip=1) #肝臓データの読み込み
hoge4 <- scan("BATtop300.txt", list(top_300=""),skip=1) #褐色脂肪組織データの読み込み
hoge5 <- scan("WATtop300.txt", list(top_300=""),skip=1) #白色脂肪組織データの読み込み
hoge6 <- c(list(hoge[[1]]),list(hoge2[[1]]),list(hoge4[[1]]),list(hoge5[[1]])) #list()の()の中身はリストではなくベクターにする必要があるみたい。
myname <- c("all", "Liver", "BAT","WAT") # ラベルをつける。（R-tips 26:names属性と要素のラベル参照）
names(hoge6) <- myname #付けたラベルをオブジェクト"hoge6"のラベルとして使用
#試しにendnodeを以下のセットでやってみる
#GO:0006810 : transport ( 15525 )
#GO:0006629 : lipid metabolism ( 3599 )
#GO:0008152 : metabolism ( 51001 )
#GO:0009058 : biosynthesis ( 14115 )
#GO:0005975 : carbohydrate metabolism ( 3823 )
#GO:0006259 : DNA metabolism ( 5122 )
#GO:0019538 : protein metabolism ( 16428 )
#GO:0016070 : RNA metabolism ( 3603 )
BPendnode <- c("GO:0006810","GO:0006629","GO:0008152","GO:0009058","GO:0005975","GO:0006259","GO:0019538","GO:0016070") #上記のGO termをカスタマイズしたEndNodeとしてセット
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(xpd=T) #描画パラメータの設定（作図領域におさまるように）
par(cex=0.8, mai=c(1.5,0.7,0.7,0.5)) #描画パラメータの設定（ラベルの字の大きさ、マージン（下左上右））
res4 <- ontoCompare(hoge6, probeType = "rat2302", endnode = BPendnode, goType="BP",plot =F) #カスタマイズしたEndoNoddeを用いる。listの数が2個以上なので棒グラフになるはず
ontoPlot(res4) #結果の（棒）グラフを表示
#------ ここまで ------
いざという時のためにオブジェクトを用意しておきました。上記にあまりに時間がかかる、とにかく結果を見たい、というときには以下をコピペしてください。
#------ ここから ------
load("090814_Barplot_res4") #あらかじめ計算しておいたオブジェクトres4の読み込み
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(xpd=T) #描画パラメータの設定（作図領域におさまるように）
par(cex=0.8, mai=c(1.5,0.7,0.7,0.5)) #描画パラメータの設定（ラベルの字の大きさ、マージン（下左上右））
ontoPlot(res4) #グラフを表示
#------ ここまで ------
4. 複数サンプルデータを、カスタマイズしたGO Term毎にontoCompare()関数を使って棒グラフ表示（2）
実際に論文を書いた時は、先行するR-Y. Li et al,BBRC 344(2006)562-570の論文の図を参考にしました。その論文で使われていた分類をもとに、EndNodeをカスタマイズしてみます。
#------ ここから ------
#endnodeのセットを、R-Y. Li et al,BBRC 344(2006)562-570の論文に近づけてみる
#GO:0008150 : biological_process ( 144962 )
#GO:0051301 : cell division ( 1060 )
#GO:0007154 : cell communication ( 12527 )
#GO:0006928 : cell motility ( 1752 )
#GO:0007582 : physiological process ( 94562 )
#GO:0006955 : immune response ( 1522 )
#GO:0006350 : transcription ( 8833 )
#GO:0006412 : protein biosynthesis ( 6116 )
#GO:0008152 : metabolism ( 51001 )
BPendnode2 <- c("GO:0008150","GO:0051301","GO:0007154","GO:0006928","GO:0007582","GO:0006955","GO:0006350","GO:0006412","GO:0008152") #上記論文のセットに近いカテゴリ
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(xpd=T) #描画パラメータの設定（作図領域におさまるように）
par(cex=0.8, mai=c(1.5,0.7,0.7,0.5)) #描画パラメータの設定（ラベルの字の大きさ、マージン（下左上右））
res5 <- ontoCompare(hoge6, probeType = "rat2302", endnode = BPendnode2, goType="BP",plot =F) #カスタマイズしたEndoNoddeを用いる。listの数が2個以上なので棒グラフになるはず
ontoPlot(res5) #結果の（棒）グラフを表示
#------ ここまで ------
いざという時のためにオブジェクトを用意しておきました。上記にあまりに時間がかかる、とにかく結果を見たい、というときには以下をコピペしてください。
#------ ここから ------
load("090814_Barplot_res5") #あらかじめ計算しておいたオブジェクトres4の読み込み
windows(height=9, width=13) #作図領域の指定（高さ9インチ、幅13インチ）
par(xpd=T) #描画パラメータの設定（作図領域におさまるように）
par(cex=0.8, mai=c(1.5,0.7,0.7,0.5)) #描画パラメータの設定（ラベルの字の大きさ、マージン（下左上右））
ontoPlot(res5) #グラフを表示
#------ ここまで ------
topGO（工事中）
topGOは、あるGO termについて、「自分がアップロードした遺伝子セット中に占めるそのカテゴリに含まれる遺伝子の割合」と、
「バックグラウンド（今回の場合はGeneChip上の全プローブセット）中に占めるそのカテゴリに含まれる遺伝子の割合」とを比較し、
どれほど有意に割合が高まっているか（濃縮されているか）を計算し、グラフ表示してくれます。いわゆるGene Enrichment Analysisを行うためのパッケージです。
分子ネットワークを可視化する強力なツールである"Cytoscape"の、プラグインである"BiNGO"とよく似た出力が返ってきます。
（準備）Rの他に、Graphvizというソフト（フリー）が必要。この辺（R単独でできないという点）がちょっとめんどくさいですね。Rを起動する前に、
Graphvizのサイトの、Downloadから
graphviz-2.24.msi（現時点のWindows版stableの最新バージョン）
をデスクトップにダウンロードする。デスクトップの上記アイコンをダブルクリックし、インストールする。設定等は全てデフォルトでOK。
スタートメニュー＞コントロールパネル＞システム＞詳細設定＞環境変数 をチェック、"システム環境変数"の中をスクロールダウンして"Path"を表示、クリックして編集、
"C:\Program Files\Graphviz2.24\bin"; が最後にくっついていたら、セミコロンを残してクォーテーションマークをはずす。最終的に、
・・・;C:\Program Files\Graphviz 2.24\bin;
という記述になればOK。
のはずだったのですが、Windows VistaでR-2.9.1という環境では上記最新のgraphvizをインストールしてもエラーが出て最後のグラフが出力されません。graphvizの方の問題のようです。
-ALLパッケージに含まれているサンプルデータを使って練習-
#------ ここから ------
library(topGO)
library(ALL)
data(ALL)
ls()
BPterms <- ls(GOBPTerm)
str(BPterms)
affyLib <- paste(annotation(ALL), "db", sep = ".")
library(package = affyLib, character.only = TRUE)
library(genefilter)
f1 <- pOverA(0.25, log2(100))
f2 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1, f2)
eset <- ALL[genefilter(ALL, ff), ]
ls()
geneNames <- featureNames(eset)
length (geneNames)
myInterestedGenes <- sample(geneNames, 100)
geneList <- factor(as.integer(geneNames %in% myInterestedGenes))
names(geneList) <- geneNames
str(geneList)
GOdata <- new("topGOdata", ontology = "MF", allGenes = geneList, annot = annFUN.db, affyLib = affyLib)
test.stat <- new("classicCount", testStatistic = GOFisherTest,name = "Fisher test")
resultFis <- getSigGroups(GOdata, test.stat)
showSigOfNodes(GOdata, score(resultFis), firstTerms = 5,useInfo = "all")
#------ ここまで ------
リンク
Rでマイクロアレイデータ解析（門田さん）
R本家のサイト
Bioconductor
CRAN（筑波大学のミラーサイト）
R & Bioconductor Manual
です。