In Bioconductor/GSEABase: Gene set enrichment data structures and methods

GeneSet

## FIXME: <adjMat> adjacency matrix -- color w. +/- 1
## FIXME: limma topTable --> GeneColorSet
## w. verbose=TRUE
suppressPackageStartupMessages({
    library(GSEABase)
    library(hgu95av2.db)
    library(GO.db)
})

A GeneSet stores a set of related gene identifiers. Important components of the gene set are a vector of identifiers, general descriptive information about the set, and information about how the gene set was constructed. To construct a gene set, use GeneSet . For example, to create a gene set from the identifiers in a subset of the sample ExpressionSet in the r Biocpkg('Biobase') package use

data(sample.ExpressionSet) # from Biobase
egs <- GeneSet(sample.ExpressionSet[201:250,], setName="Sample")
egs

Each gene set may have a name. The gene set contains r length(geneIds(egs)) identifiers ('genes') from the ExpressionSet. These are accessible using geneIds, e.g.,

head(geneIds(egs))

The gene set records that the identifiers are probeset names from the annotation package r Biocpkg('hgu95av2.db'), and that the source of the gene set was an ExpressionSet. Additional details are available:

details(egs)

The set identifier, set version, and creation date provide mechanisms for carefully curating gene sets. Additional information is automatically copied from the expression set used to create egs.

  ## FIXME: GeneSet(AnnotationIdentifier("hgu95av2")) --> non-empty
  ## FIXME: GeneSet(AnnotationIdentifier("hgu95av2"),
  ## collectionType=GOCollection()) filters on GOCollection (or KEGG)

View additional methods for creating gene sets with:

showMethods("GeneSet", inherited=FALSE)

These are described on the GeneSet help page.

The identifier type of gene sets created from expression sets is r as.character(class(geneIdType(egs))). Additional predefined identifiers are available:

names(slot(getClass("GeneIdentifierType"), "subclasses"))

It is possible to map between identifier types (provided the corresponding map in the annotation package exists):

mapIdentifiers(egs, EntrezIdentifier())

mapIdentifiers consults the gene set to determine that annotation (probeset) identifiers are to be converted to Entrez IDs based on the mapping in the r Biocpkg('hgu95av2.db') package. The function mapIdentifiers can automatically determine many of the common maps; it is also possible to provide as a third argument an environment containing an arbitrary map. Use the verbose argument of mapIdentifiers to be informed when the map between identifier types is not 1:1.

A gene set can be created with different types of identifier, e.g., to create a gene set with Entrez IDs, use

library(annotate)                       # getEG
eids <- unique(getEG(geneIds(egs), "hgu95av2"))
eids <- eids[!is.na(eids)]
GeneSet(EntrezIdentifier(), geneIds=as.character(eids))

The collectionType of a gene set provides additional information about a gene set. Available collection types are

names(slot(getClass("CollectionType"), "subclasses"))

Collection types are most important when creating gene sets. For instance, the GOCollection class allows creation of gene sets based on gene ontology (GO) terms. The following command creates a gene set from two terms, including all genes with a particular evidence code.

GeneSet(GOCollection(c("GO:0005488", "GO:0019825"),
                     evidenceCode="IDA"),
        geneIdType=EntrezIdentifier("org.Hs.eg.db"),
        setName="Sample GO Collection")

This creates a gene set by consulting an object in the r Biocpkg('GO.db') package. A gene set created from an expression set, and with collection type GOCollection consults the appropriate environment in the annotation package associated with the expression set.

Gene sets encoded in XML following the schema and conventions of the Broad Institute can be read into R as follows:

fl <- system.file("extdata", "Broad1.xml", package="GSEABase")
bgs <- GeneSet(BroadCollection(), urls=fl)
bgs

The example above uses a local file, but the source for the gene set might also be a web address accessible via http:// the protocol. The file name is added to the url of the gene set. The set name and category of the Broad collection indicate that the gene set contains gene symbols from band 24 of the q arm of chromosome 16. The probe sets in chip hgu95av2 corresponding to these symbols can be determined by mapping identifiers

bgs1 <- mapIdentifiers(bgs, AnnotationIdentifier("hgu95av2"))
bgs1

Subsetting creates sets with just the symbols identified. Subsetting can use indices or symbol names.

bgs[1:5]
bgs[c("GALNS", "LOC646365")]

Logical operations provide a convenient way to identify genes with particular properties. For instance, the intersection

egs & bgs1

is empty (note that the identifiers in the two sets were of the same type), indicating that none of the identifiers in egs are on 16q24. Additional operations on sets include union (performed with |) and difference (setdiff).

Methods exist to directly subset expression sets using gene sets

sample.ExpressionSet[bgs,]

Remember that sample.ExpressionSet contains just r nrow(sample.ExpressionSet) probe sets, so the small size of the subset is not surprising. Note also that subsetting required mapping of symbol identifiers in bgs to AnnotationIdentifiers; this map used the annotation information in the expression set, and is not necessarily 1:1.

GeneColorSet

A GeneColorSet is a gene set with "coloring" to indicate how features of genes and phenotypes are associated. The following sample data describes how changes in expression levels of several genes (with Entrez and Symbol names) influence cisplatin resistance phenotype.

conn <- textConnection("
Entrez ID, Gene Symbol, Expression level, Phenotype response
##used to be MRP2
1244, ABCC2, Increase, Resistant
538, ATP7A, Increase, Resistant
540, ATP7B, Increase, Resistant
9961, MVP, Increase, Resistant
##the LRP below must be MVP
##LRP, Increase, Resistant - need to know which one
7507,XPA, Increase, Resistant
2067, ERCC1, Increase, Resistant
##TOP, Increase, Resistant  - need to know which one, notes say II
672, BRCA1, Increase, Resistant
3725, JUN, Increase, Resistant
#GCS, Increase, Resistant  - my notes say alpha-GCS - so which one?
##I only found gamma at PubMed as being related
2730, GCLM, Increase, Resistant")
tbl <- read.csv(conn, strip.white=TRUE, comment.char="#", stringsAsFactors = TRUE)
close(conn)
unlink(conn)

tbl

Note that three different aspects of data influence coloring: the phenotype under consideration (cisplatin resistance), whether expression responses refer to increasing or decreasing levels of gene expression, and whether the phenotypic response represents greater resistance or sensitivity to cisplatin. Here is the resulting gene color set:

gcs <- GeneColorSet(EntrezIdentifier(),
                    setName="A color set",
                    geneIds=as.character(tbl$Entrez.ID),
                    phenotype="Cisplatin resistance",
                    geneColor=tbl$Expression.level,
                    phenotypeColor=tbl$Phenotype.response)
gcs

Gene color sets can be used in the same way as gene sets, e.g., for subsetting expression sets (provided the map between identifiers is 1:1, so that coloring corresponding to identifiers can be determined). The coloring method allows access to the coloring information with a data frame interface; phenotype, geneColor and phenotypeColor provide additional accessors.

GeneSetCollection

A GeneSetCollection is a collection of gene sets. Sets in the collection must have distinct setNames, but can be a mix of GeneSet and GeneColorSet. Two convenient ways to create a gene set collection are by specifying a source of identifiers (e.g., an ExpressionSet or AnnotationIdentifier) and how the identifiers are to be induced into sets (e.g., by consulting the GO or KEGG ontologies):

gsc <- GeneSetCollection(sample.ExpressionSet[201:250,], setType=GOCollection())
gsc
gsc[["GO:0005737"]]

In this example, the annotation identifiers of the sample expression set are organized into gene sets based on their presence in GO pathways. Providing arguments such as evidenceCode to GOCollection act to select just those pathways satisfying the GO collection constraint:

GeneSetCollection(sample.ExpressionSet[201:300,],
                  setType=GOCollection(evidenceCode="IMP"))

Sets in the collection are named after the GO terms, and can be accessed by numeric index or name.

A file or url containing several gene sets defined by Broad XML can be used to create a gene set collection, e.g.,

   ## FIXME: BroadCollection default to paste("c", 1:4, sep="")
   ## FIXME: GeneSetCollection(BroadCollection(), urls=fl); filters on bcCategory
fl <- system.file("extdata", "Broad.xml", package="GSEABase")
gss <- getBroadSets(fl)
gss
names(gss)

Identifiers within a gene set collection can be mapped to a common type (provided maps are available) with, for example,

gsc <- mapIdentifiers(gsc, EntrezIdentifier())
gsc
gsc[["GO:0005737"]]

A convenient way to visualize a GeneSetCollection is with the r Biocpkg('ReportingTools') package.

## 'interesting' gene sets
idx <- sapply(gsc, function(x) length(geneIds(x))) > 2

library(ReportingTools)
gscReport <- HTMLReport(
    shortName="gsc_example",
    title="GSEABase Vignette GeneSetCollection",
    basePath=tempdir())
publish(gsc[idx], gscReport, annotation.db="org.Hs.eg")
url <- finish(gscReport)