R/BioconductorのGOstats packageをもちいたGene Ontology(GO)やKEGGのenrichment解析

GO enrichment analysis

R/Bioconductorの GOstats packageをもちいてGene Ontology (GO) enrichment analysisをおこなう。
-1*log10(p-value)の棒グラフや、各アノテーションに該当した遺伝子のリストも一緒に取得する。

準備

以下の2種類の遺伝子リストをEntrez GeneのID番号（=NCBI Gene, LocusLink）で用意する。

• EntrezGeneID_interest.txt : 自分の注目する遺伝子セット
• EntrezGeneID_background.txt : その遺伝子セットをみつけてきたときのおおもとの遺伝子全体のリスト

Gene SymbolなどからBioMartなどでEntrez Gene IDに変換すると楽かも。

Entrez Geneのリスト読み込み

BGFILE <- "EntrezGeneID_background.txt"
INFILE <- "EntrezGeneID_interest.txt"

geneIds <- as.character(read.table(INFILE,header=FALSE)[,1])
univIds <- as.character(read.table(BGFILE,header=FALSE)[,1])

閾値などの設定

ANNOTATION <- "org.Mm.eg.db"
ONTOLOGY <- "BP"
PVALUE<-0.01
TOP <- 30
OUTNAME <- "myinterest"

ANNOTATIONはマウスなら"org.Mm.eg.db"だしヒトなら"org.Hs.eg.db"
ONTOLOGYは"BP","MF","CC"
（それぞれbiological process, molecular function, cellular component）

ライブラリの読み込み

library("GOstats")
library("GO.db")
library("annotate")
library("org.Mm.eg.db")

関数を作成

GOHYPERG <- function(geneIds,univIds,ANNOTATION,ONTOLOGY,PVALUE,TOP,OUTNAME) {
    param <- new(
　     "GOHyperGParams",
　     geneIds=geneIds,
　     universeGeneIds=univIds,
　     annotation= ANNOTATION,
　     ontology=ONTOLOGY
　 )
　 GOHypG <- hyperGTest(param)
　 DETECTGO <- sigCategories(
　     GOHypG,
　     pvalue <- PVALUE
　 )
　 if (length(DETECTGO)==0) {
　     RESULT_GOHYP <- data.frame(
　         GOid="",
　         pvalue="",
　         oddratio="",
　         expected="",
　         obserbed="",
　         background="",
　         numberMappedGenes=paste(geneMappedCount(GOHypG),length(geneIds),sep="/"),
　         numberMappedBkgGenes=paste(universeMappedCount(GOHypG),length(univIds),sep="/"),
　         GOterm="-",
　         observedGenes="-"
　     )
　 } else {
　     RESULT_GOHYP <- data.frame(
　         GOid=DETECTGO,
　         pvalue=pvalues(GOHypG)[DETECTGO],
　         oddratio=oddsRatios(GOHypG)[DETECTGO],
　         expected=expectedCounts(GOHypG)[DETECTGO],
　         obserbed=geneCounts(GOHypG)[DETECTGO],
　         background=universeCounts(GOHypG)[DETECTGO],
　         numberMappedGenes=paste(geneMappedCount(GOHypG),length(geneIds),sep="/"),
　         numberMappedBkgGenes=paste(universeMappedCount(GOHypG),length(univIds),sep="/"),
　         GOterm=sapply(
　             DETECTGO,
　             function(goIDs){
　                 sapply(
　                     goIDs,
　                     function(goID){Term(get(goID,GOTERM))}
　                 )
　             }
　         ),
　         observedGenes=sapply(
　             DETECTGO,
　             function(X){
　                 paste(
　                     sapply(
　                         geneIdsByCategory(GOHypG,catids = NULL)[[X]],
　                         function(XXX){paste(get(XXX,org.Mm.egSYMBOL),"(",XXX,")",sep="")}
　                     ),
　                     collapse=","
　                 )
　             }
　         )
　     )
　     
    　 BARHEIGHT <- -1 * log(RESULT_GOHYP[1:TOP,2],10)[TOP:1]
　     BARNAME <- RESULT_GOHYP$GOterm[TOP:1]
　 　 SAVEPDFNAME <- paste("RESULT_GOHYP_",ONTOLOGY,"_",OUTNAME,".pdf",sep="")
　 　 pdf(SAVEPDFNAME)
　 　 par(oma=c(2,2,2,2),mar=c(6,16,1,1))
　 　 barplot(
　         BARHEIGHT,
　         width=1,
　         space=NULL,
　         horiz=TRUE,
　         las=1,
　         main=paste(strsplit(SAVEFILENAME,".txt")[[1]][1]),
　         cex.main=.7,
　         sub=paste("(TOP",TOP," GOs)",sep=""),
　         cex.sub=.7,
　         names.arg=BARNAME,
　         cex.names=.6,
　         xlab="-1 * log10(p-value)",
　         cex.axis = 1
　     )
　     dev.off()
　 }
　         
　 SAVEFILENAME <- paste("RESULT_GOHYP_",ONTOLOGY,"_",OUTNAME,".txt",sep="")
　 write.table(RESULT_GOHYP,SAVEFILENAME,sep="\t",eol="\n",quote=FALSE,row.names=FALSE)
　 
　 return(OUTNAME)
　 
} ## end of func. GOHYPERG

実行

GOHYPERG(geneIds,univIds,ANNOTATION,ONTOLOGY,PVALUE,TOP,OUTNAME)

"BP","MF","CC"全部やりたいなら

for (ONTOLOGY in c("BP","MF","CC")) {
    GOHYPERG(geneIds,univIds,ANNOTATION,ONTOLOGY,PVALUE,TOP,OUTNAME)
}

KEGGでenrichment analysis

GOのかわりにKEGGでもやれる。"BP","MF","CC"のかわりに"over","under"になる

ライブラリ読み込み

library("KEGG.db")

関数を作成

KEGGHYPERG <- function(geneIds,univIds,ANNOTATION,OVERUNDER,PVALUE,TOP,OUTNAME) {
    keggparam <- new(
        "KEGGHyperGParams",
        geneIds=geneIds,
        universeGeneIds=univIds,
        annotation= ANNOTATION,
        pvalueCutoff=PVALUE,
        testDirection=OVERUNDER
    )
    KEGGHypG <- hyperGTest(keggparam)
    DETECTKEGG <- sigCategories(KEGGHypG,pvalue<-PVALUE)
    if (length(DETECTKEGG)==0) {
        RESULT_KEGGHYP <- data.frame(
            KEGGid="",
            pvalue="",
            oddratio="",
            expected="",
            obserbed="",
            background="",
            numberMappedGenes=paste(geneMappedCount(KEGGHypG),length(geneIds),sep="/"),
            numberMappedBkgGenes=paste(universeMappedCount(KEGGHypG),length(univIds),sep="/"),
            KEGGterm="-",
            observedGenes="-"
        )
    } else {
        RESULT_KEGGHYP <- data.frame(
            KEGGid=DETECTKEGG,
            pvalue=pvalues(KEGGHypG)[DETECTKEGG],
            oddratio=oddsRatios(KEGGHypG)[DETECTKEGG],
            expected=expectedCounts(KEGGHypG)[DETECTKEGG],
            obserbed=geneCounts(KEGGHypG)[DETECTKEGG],
            background=universeCounts(KEGGHypG)[DETECTKEGG],
            numberMappedGenes=paste(geneMappedCount(KEGGHypG),length(geneIds),sep="/"),
            numberMappedBkgGenes=paste(universeMappedCount(KEGGHypG),length(univIds),sep="/"),
            KEGGterm=sapply(DETECTKEGG,function(keggIDs){sapply(keggIDs,function(keggID){get(keggID,KEGGPATHID2NAME)})}),
            observedGenes=sapply(
                DETECTKEGG,
                function(X){paste(sapply(
                    geneIdsByCategory(KEGGHypG, catids = NULL)[[X]],
                    function(XXX){paste(get(XXX,org.Mm.egSYMBOL),"(",XXX,")",sep="")}),collapse=",")}
            )
        )
        SAVEFILENAME <- paste("RESULT_KEGGHYP_",OVERUNDER,OUTNAME,".txt",sep="")
        write.table(RESULT_KEGGHYP,SAVEFILENAME,sep="\t",eol="\n",quote=FALSE,row.names=FALSE)
        BARHEIGHT <- -1 * log(RESULT_KEGGHYP[1:TOP,2],10)[TOP:1]
        BARNAME <- RESULT_KEGGHYP$KEGGterm[TOP:1]
        SAVEPDFNAME <- paste("RESULT_KEGGHYP_",OVERUNDER,OUTNAME,".pdf",sep="")
        pdf(SAVEPDFNAME)
        par(oma=c(2,2,2,2),mar=c(6,16,1,1))
        barplot(
            BARHEIGHT,
            width=1,
            space=NULL,
            horiz=TRUE,
            las=1,
            main=paste(strsplit(SAVEFILENAME,".txt")[[1]][1]),
            cex.main=.7,
            sub=paste("(TOP",TOP," KEGGs)",sep=""),
            cex.sub=.7,
            names.arg=BARNAME,
            cex.names=.6,
            xlab="-1 * log10(p-value)",
            cex.axis = 1
        )
        dev.off()
    }
    return(OUTNAME)
} ## end of func. KEGGHYPERG

実行

for (OVERUNDER in c("over","under")) {
    KEGGHYPERG(geneIds,univIds,ANNOTATION,OVERUNDER,PVALUE,TOP,OUTNAME)
}