This is an R wrapper for Pavlidis Lab’s ermineJ. A tool for gene set enrichment analysis with multifunctionality correction.

Table of Contents

Installation

ermineR requries 64 bit version of java to function. If you are a Mac user make sure you have the java SDK.

After java is installed you can install ermineR by doing

devtools::install_packages('PavlidisLab/ermineR')

If ermineR cannot find your java home by itself. Use either install rJava or use Sys.setenv(JAVA_HOME=javaHome) to point ermineR to the right path.

Usage

See documentation for ora, roc, gsr, precRecall and corr to see how to use them.

An explanation of what each method does is given. We recommend users start with the precRecall (for gene ranking-based enrichment analysis) or ora (for hit-list over-representation analysis).

Replicable go terms

GO terms are updated frequently so results can differ between versions. The default option of all ermineR functions is to get the latest GO version however this means you may get different results when you repeat the experiment later. If you want to use a specific version of GO, ermineR provides functions to deal with that.

  • goToday: Downloads the latest version of go to a path you provide
  • getGoDates: Lists all dates where a go version is available, from the most recent to oldest
  • goAtDate: Given a valid date, downloads the Go version from a specific date to a file path you provide

To use a specific version of GO, make sure to set geneSetDescription argument of all ermineR functions to the file path where you saved the go terms

Examples

Use GSR with gene scores

Here we will use a mock scores file located in our tests directory. The score file is specifically created to be enriched in genes with the term GO:0051082.

scores = read.table("tests/testthat/testFiles/pValues", header=T, row.names = 1)
head(scores)
##                pvalue
## 206190_at   0.3163401
## 208385_at   0.5186824
## 65086_at    0.6620389
## 202281_at   0.4068895
## 211622_s_at 0.9128846
## 219257_s_at 0.2936740

This scores file only includes scores for 118 genes. The file was generated using GPL96’s probesets so that is the annotation we’ll be using. Any gene that is not reperesented by the score file will be ignored.

gsrOut = gsr(annotation = 'GPL96',
                 scores = scores,
                 scoreColumn = 1,
                 iterations = 10000,
                 bigIsBetter = FALSE,
                 logTrans = TRUE)

head(gsrOut$results) %>% knitr::kable()
Name ID NumProbes NumGenes RawScore Pval CorrectedPvalue MFPvalue CorrectedMFPvalue Multifunctionality Same as GeneMembers
protein folding GO:0006457 44 26 3.091833 0.00e+00 0.000000 0.0000000 0.0000000 0.219 NA AIP
cellular component biogenesis GO:0044085 59 27 2.548777 0.00e+00 0.000000 0.0001000 0.0033000 0.769 NA AAMP
unfolded protein binding GO:0051082 61 33 3.306600 0.00e+00 0.000000 0.0000000 0.0000000 0.175 NA AAMP
protein binding GO:0005515 114 63 1.922271 5.94e-05 0.001959 0.0001511 0.0039900 0.564 NA AAMP
cellular component assembly GO:0022607 58 26 2.513885 1.00e-04 0.002640 0.0001000 0.0044000 0.808 NA CALR
intracellular organelle part GO:0044446 99 52 1.980901 1.00e-04 0.002200 0.0019580 0.0198779 0.302 NA AIP

Use Precision Recall with gene scores

We will use the same scores file from the example above

precRecallOut = precRecall(annotation = 'GPL96',
                           scores = scores,
                           scoreColumn = 1,
                           iterations = 10000,
                           bigIsBetter = FALSE,
                           logTrans = TRUE)

head(precRecallOut$results) %>% knitr::kable()
Name ID NumProbes NumGenes RawScore Pval CorrectedPvalue MFPvalue CorrectedMFPvalue Multifunctionality Same as GeneMembers
protein folding GO:0006457 44 26 0.7209345 0e+00 0.00000 0.0000 0.0000000 0.219 NA AIP
unfolded protein binding GO:0051082 61 33 1.0000000 0e+00 0.00000 0.0000 0.0000000 0.175 NA AAMP
protein binding GO:0005515 114 63 0.9259590 1e-04 0.00440 0.0002 0.0088000 0.564 NA AAMP
extracellular region GO:0005576 55 30 0.6050500 1e-04 0.00330 0.0023 0.0337333 0.264 NA AAMP
intracellular GO:0005622 129 70 0.9568881 1e-04 0.00264 0.0011 0.0363000 0.652 NA AAMP
intracellular part GO:0044424 129 70 0.9568881 1e-04 0.00220 0.0011 0.0290400 0.653 NA AAMP

Use ORA with a hitlist

library(dplyr)


# genes for GO:0051082
hitlist = c("AAMP", "AFG3L2", "AHSP", "AIP", "AIPL1", "APCS", "BBS12", 
            "CALR", "CALR3", "CANX", "CCDC115", "CCT2", "CCT3", "CCT4", "CCT5", 
            "CCT6A", "CCT6B", "CCT7", "CCT8", "CCT8L1P", "CCT8L2", "CDC37", 
            "CDC37L1", "CHAF1A", "CHAF1B", "CLGN", "CLN3", "CLPX", "CRYAA", 
            "CRYAB", "DNAJA1", "DNAJA2", "DNAJA3", "DNAJA4", "DNAJB1", "DNAJB11", 
            "DNAJB13", "DNAJB2", "DNAJB4", "DNAJB5", "DNAJB6", "DNAJB8", 
            "DNAJC4", "DZIP3", "ERLEC1", "ERO1B", "FYCO1", "GRPEL1", "GRPEL2", 
            "GRXCR2", "HEATR3", "HSP90AA1", "HSP90AA2P", "HSP90AA4P", "HSP90AA5P", 
            "HSP90AB1", "HSP90AB2P", "HSP90AB3P", "HSP90AB4P", "HSP90B1", 
            "HSP90B2P", "HSPA1A", "HSPA1B", "HSPA1L", "HSPA2", "HSPA5", "HSPA6", 
            "HSPA8", "HSPA9", "HSPB6", "HSPD1", "HSPE1", "HTRA2", "LMAN1", 
            "MDN1", "MKKS", "NAP1L4", "NDUFAF1", "NPM1", "NUDC", "NUDCD2", 
            "NUDCD3", "PDRG1", "PET100", "PFDN1", "PFDN2", "PFDN4", "PFDN5", 
            "PFDN6", "PIKFYVE", "PPIA", "PPIB", "PTGES3", "RP2", "RUVBL2", 
            "SCAP", "SCG5", "SERPINH1", "SHQ1", "SIL1", "SPG7", "SRSF10", 
            "SRSF12", "ST13", "SYVN1", "TAPBP", "TCP1", "TMEM67", "TOMM20", 
            "TOR1A", "TRAP1", "TTC1", "TUBB4B", "UGGT1", "ZFYVE21")
oraOut = ora(annotation = 'Generic_human',
             hitlist = hitlist)

head(oraOut$results) %>% knitr::kable()
Name ID NumProbes NumGenes RawScore Pval CorrectedPvalue MFPvalue CorrectedMFPvalue Multifunctionality Same as GeneMembers
unfolded protein binding GO:0051082 123 123 114 0 0 0 0 0.641 NA AAMP
chaperone-mediated protein folding GO:0061077 81 81 37 0 0 0 0 0.679 NA BAG1
protein binding involved in protein folding GO:0044183 43 43 31 0 0 0 0 0.641 NA BBS12
‘de novo’ protein folding GO:0006458 54 54 31 0 0 0 0 0.707 NA BAG1
chaperone binding GO:0051087 96 96 23 0 0 0 0 0.758 NA AHSA1
response to topologically incorrect protein GO:0035966 182 182 26 0 0 0 0 0.903 NA ACADVL

Using your own GO annotations

If you want to use your own GO annotations instead of getting files provided by Pavlidis Lab, you can use makeAnnotation after turning your annotations into a list. See the example below

library('org.Hs.eg.db') # get go terms from bioconductor 
goAnnots = as.list(org.Hs.egGO)
goAnnots = goAnnots %>% lapply(names)
goAnnots %>% head
## $`1`
##  [1] "GO:0002576" "GO:0008150" "GO:0043312" "GO:0005576" "GO:0005576"
##  [6] "GO:0005615" "GO:0031012" "GO:0031093" "GO:0034774" "GO:0070062"
## [11] "GO:0072562" "GO:1904813" "GO:0003674"
## 
## $`2`
##  [1] "GO:0001869" "GO:0002576" "GO:0007597" "GO:0010951" "GO:0022617"
##  [6] "GO:0043547" "GO:0048863" "GO:0051056" "GO:0005576" "GO:0005576"
## [11] "GO:0005829" "GO:0031093" "GO:0070062" "GO:0072562" "GO:0002020"
## [16] "GO:0004867" "GO:0005096" "GO:0005102" "GO:0005515" "GO:0019838"
## [21] "GO:0019899" "GO:0019959" "GO:0019966" "GO:0043120" "GO:0048306"
## 
## $`3`
## NULL
## 
## $`9`
## [1] "GO:0006805" "GO:0005829" "GO:0004060"
## 
## $`10`
## [1] "GO:0006805" "GO:0005829" "GO:0004060" "GO:0005515"
## 
## $`11`
## NULL

The goAnnots object we created has go terms per entrez ID. Similar lists can be obtained from other species db packages in bioconductor and some array annotation packages. We will now use the makeAnnotation function to create our annotation file. This file will have the names of this list (entrez IDs) as gene identifiers so any score or hitlist file you provide should have the entrez IDs as well.

makeAnnotation only needs the list with gene identifiers and go terms to work. But if you want to have a complete annotation file you can also provide gene symbols and gene names. Gene names have no effect on the analysis. Gene symbols matter if you are providing custom gene sets and using “Option 2” or if same genes are represented by multiple gene identifiers (eg. probes). Gene symbols will also be returned in the GeneMembers column of the output. If they are not provided, gene IDs will also be used as gene symbols

Here we’ll set them both for good measure.

geneSymbols = as.list(org.Hs.egSYMBOL) %>% unlist
geneName = as.list(org.Hs.egGENENAME) %>% unlist

annotation = makeAnnotation(goAnnots,
                            symbol = geneSymbols,
                            name = geneName,
                            output = NULL, # you can choose to save the annotation to a file
                            return = TRUE) # if you only want to save it to a file, you don't need to return

Now that we have the annotation object, we can use it to run an analysis. We’ll try to generate a hitlist only composed of genes annotated with GO:0051082.

mockHitlist = goAnnots %>% sapply(function(x){'GO:0051082' %in% x}) %>% 
    {goAnnots[.]} %>% 
    names

mockHitlist %>% head
## [1] "14"  "325" "811" "821" "871" "908"
oraOut = ora(annotation = annotation,
             hitlist = mockHitlist)

head(oraOut$results) %>% knitr::kable()
Name ID NumProbes NumGenes RawScore Pval CorrectedPvalue MFPvalue CorrectedMFPvalue Multifunctionality Same as GeneMembers
unfolded protein binding GO:0051082 113 113 113 0 0 0 0 0.644 NA AAMP
protein binding involved in protein folding GO:0044183 28 28 23 0 0 0 0 0.543 NA BBS12
chaperone-mediated protein folding GO:0061077 60 60 26 0 0 0 0 0.674 NA BBS12
‘de novo’ protein folding GO:0006458 36 36 21 0 0 0 0 0.661 NA BBS12
chaperone binding GO:0051087 89 89 22 0 0 0 0 0.765 NA AHSA1
response to topologically incorrect protein GO:0035966 179 179 25 0 0 0 0 0.907 NA ACADVL

We can see GO:0051082 is the top scoring hit as expected.