RパッケージのGOseq。日本語でも、GOseqの使い方の説明はネット上に散見されます。ただ、多くはヒトなどゲノム情報が手に入る生物種を対象にしていて、いわゆる非モデル生物の場合の説明は見かけません。マニュアルを見ても、コードの例までは書いてません。
そこで、非モデル生物の場合のenrichment analysisのコード例を以下に書いておきます。
まずBioconductorからGOseqパッケージをインストールします。そして以下を実行。始めにデータセットを作ります。非モデル生物の場合は、GOをあらかじめBlast2GOなどのソフトで取得しておく必要があります(ここでは適当に設定)。
library (goseq)
map_test <- list( gene1=c("GO:0050790","GO:0004197"), gene2=NA, gene3=NA, gene4=c("GO:0005975","GO:0004553"),gene5=NA,gene6= c("GO:0005975") , gene7=c("GO:0005975"), gene8=NA, gene9=c("GO:0004553"), gene10= c("GO:0004553") )
genelist_test <- data.frame( DEgenes=c(1,0,1,0,1,0,1,0,1,1), bias.data = c(500,6110,2000,410,600, 1000, 4500,1200,400,8010) )
## DEgenesは発現変動したかどうか(1なら変動)を示す
## bias.dataは遺伝子長
rownames(genelist_test)<-c( "gene1", "gene2", "gene3", "gene4", "gene5", "gene6","gene7", "gene8", "gene9", "gene10" )
次に長さの補正。
pwf_test <- nullp(genelist_test[,1], bias.data=genelist_test[,2])
rownames( pwf_test) <- rownames(genelist_test)
最後にgoseq関数でenrichment analysisをおこなう。 ここでuse_genes_without_catの引数をTRUEに設定することが、ゲノム情報のない生物での解析のポイント。
result <- goseq (pwf = pwf_test, gene2cat = map_test , use_genes_without_cat = TRUE, method = "Wallenius")
なおmethodにWallenius分布以外を用いれば、特に遺伝子長は関係ないので、nullp関数の下りはなくてもenrichment解析ができます(pwf_test <- nullp (...)のコードを飛ばして、最後にgoseq (pwf = genelist, ...)と書く)。
結果。
参考にしたページ:Using goseq on non-model organism - how to define genome?