備忘録 a record of inner life

やったことや考えたこと・本・論文・音楽の備忘録。 特に環境科学・生態毒性に関して。

非モデル生物のGO enrichment analysisをGOseqでおこなう

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, ...)と書く)。

 

結果。 

> result
category over_represented_pvalue under_represented_pvalue numDEInCat numInCat term ontology
1 GO:0004197 0.6636645      1.0000000       1    1    cysteine-type endopeptidase activity MF
4 GO:0050790 0.6636645      1.0000000       1    1    regulation of catalytic activity BP
2 GO:0004553 0.7150671      0.7978634       2    3    hydrolase activity, hydrolyzing O-glycosyl compounds MF
3 GO:0005975 0.9661791      0.3357600        1    3    carbohydrate metabolic process BP

 

 

参考にしたページ:Using goseq on non-model organism - how to define genome?