PaxtoolsR Tutorial

The paxtoolsr package exposes a number of the algorithms and functions provided by the Paxtools Java library and Pathway Commons webservice allowing them to be used in R.

Overview

BioPAX, Paxtools, Pathway Commons, and the Simple Interaction Format

The Biological Pathway Exchange (BioPAX) format is a community-driven standard language to represent biological pathways at the molecular and cellular level and to facilitate the exchange of pathway data. BioPAX can represent metabolic and signaling pathways, molecular and genetic interactions and gene regulation networks. Using BioPAX, millions of interactions, organized into thousands of pathways, from many organisms are available from a growing number of databases. This large amount of pathway data in a computable form will support visualization, analysis and biological discovery. The BioPAX format using syntax for data exchange based on the OWL (Web Ontology Language) that aids pathway data integration; classes in the BioPAX ontology are described here. Ontologies are formal systems for knowledge representation allowing machine-readability of pathway data; one well-known example of a biological ontology is the Gene Ontology for biological terms.

Paxtools is a Java libary that allows users to interact with biological pathways represented in the BioPAX language. Pathway Commons is a resource that integrates biological pathway information for a number of public pathway databases, including Reactome, PantherDB, HumanCyc, etc. that are represented using the BioPAX language.

NOTE: BioPAX can encode very detailed information about biological processes. Analysis of this data, however, can be complicated as one needs to consider a wide array of n-ary relationships, different states of entities and generics. An alternative approach is to derive higher order relations based on a set of templates to define a simple binary network between biological entities and use conventional graph algorithms to analyze it. For many users of this package, the binary representation termed the Simple Interaction Format (SIF) will be the main entry point to the usage of BioPAX data. Conversion of BioPAX datasets to the SIF format is done through a series of simplification rules.

Limitations

The Paxtools Java library produces that full model of a given BioPAX data set that can be searched via code. The paxtoolsr provides a limited set of functionality mainly to produce SIF representations of networks that can be analyzed in R.

Basics

Installation

biocLite("paxtoolsr")

Getting Started

Load paxtoolsr package:

library(paxtoolsr)

A list of all accessible vignettes and methods is available with the following command:

help.search("paxtoolsr")

For help on any paxtoolsr package functions, use one of the following command formats:

help(graphPc)
?graphPc

Common Function Return Types

paxtoolsr return two main types of values data.frame and XMLInternalDocument. Data.frames are table like data structures. XMLInternalDocument is a representation provided by the XML package and this data structure form is returned for functions that search or return raw BioPAX results. An XMLInternalDocument can be used as the input for any function requiring a BioPAX file.

Handling BioPAX OWL Files

paxtoolsr provides several functions for handling BioPAX OWL files. paxtoolsr provides several functions for handling BioPAX OWL files: merging, validation, conversion to other formats. Many databases with protein-protein interactions and pathway information export the BioPAX format and BioPAX files; databases that support the BioPAX format can be found on PathGuide, a resource for pathway information.

Merging BioPAX Files

We illustrate how to merge two BioPAX files. Only entities that share IDs will be merged; no additional merging occurs on cross-references. The merging occurs as described further in the Java library documentation. Throughout this tutorial we use the system.file() command to access sample BioPAX files included with the paxtoolsr package. Merging may result in warning messages caused as a result of redundant actions being checked against by the Java library; these messages may be ignored.

file1 <- system.file("extdata", "raf_map_kinase_cascade_reactome.owl", package = "paxtoolsr")
file2 <- system.file("extdata", "circadian_clock.owl", package = "paxtoolsr")

mergedFile <- mergeBiopax(file1, file2)

Here we summarize information about one of the BioPAX files provide in the paxtoolsr package. The summarize() function produces a counts for various BioPAX classes and can be used to filter through BioPAX files matching particular characteristics. In the example below, we show that the merged file contains the sum of the Catalysis elements from the original two BioPAX files. This can be used iterate over and to identify files with particular properties quickly or to summarize across the files from a set.

s1 <- summarize(file1)
s2 <- summarize(file2)
s3 <- summarize(mergedFile)

s1$Catalysis
s2$Catalysis
s3$Catalysis

## [1] "5"
## [1] "18"
## [1] "23"

Validating BioPAX Files

To validate BioPAX paxtoolsr the types of validation performed are described in the BioPAX Validator publication by Rodchenkov I, et al.

errorLog <- validate(system.file("extdata", "raf_map_kinase_cascade_reactome.owl", 
    package = "paxtoolsr"), onlyErrors = TRUE)

Converting BioPAX Files to Other Formats

It is often useful to convert BioPAX into other formats. Currently, paxtoolsr supports conversion to Gene Set Enrichment Analysis (GSEA, .gmt), Systems Biology Graphical Notation (SBGN, .sbgn), Simple Interaction Format.

Simple Interaction Format (SIF) Network

The basic SIF format includes a three columns: PARTICIPANT_A, INTERACTION_TYPE, and PARTICIPANT_B; possible INTERACTION_TYPEs are described here.

sif <- toSif(system.file("extdata", "biopax3-short-metabolic-pathway.owl", package = "paxtoolsr"))

SIF representations of networks are returned as data.frame objects. SIF representations can be readily be visualized in network analysis tools, such as Cytoscape, which can be interfaced with through the R package, RCytoscape.

head(sif)

##               PARTICIPANT_A INTERACTION_TYPE              PARTICIPANT_B
## 1  Adenosine 5'-diphosphate  used-to-produce  Adenosine 5'-triphosphate
## 2  Adenosine 5'-diphosphate      reacts-with beta-D-glucose 6-phosphate
## 3  Adenosine 5'-diphosphate  used-to-produce             beta-D-glucose
## 4 Adenosine 5'-triphosphate  used-to-produce   Adenosine 5'-diphosphate
## 5 Adenosine 5'-triphosphate  used-to-produce beta-D-glucose 6-phosphate
## 6 Adenosine 5'-triphosphate      reacts-with             beta-D-glucose

Extended Simple Interaction Format (SIF) Network

Often analysis requires additional items of information, this could be the literature references related to a resource or the name of the data source where an interaction was derived. This information can be retrieved as part of an extended SIF network. A BioPAX dataset can be converted to extended SIF network.

# Select additional node and edge properties
nodeProps <- c("EntityReference/name", "EntityReference/xref")
edgeProps <- "Interaction/dataSource/displayName"

results <- toSifnx(inputFile = system.file("extdata", "dna_replication.owl", 
    package = "paxtoolsr"), nodeProps = nodeProps, edgeProps = edgeProps)

The results object is a list with two entries: nodes and edges. nodes will be a data.frame where each row corresponds to a biological entity, an EntityReference, and will contain any user-selected node properties as additional columns. Similarly, edges will be a data.frame with a SIF extended with any user-selected properties for an Interaction as additional columns. Information on possible properties for an EntityReference or Interaction is available through the BioPAX ontology. It is also possible to download a pre-computed extended SIF representation for the entire Pathway Commons database that includes information about the data sources for interactions and identifiers for nodes; refer to documentation of the method for more details about the returned entries.

NOTE: downloadPc may take several minutes to complete. It is suggested that the results of this command be saved locally rather than using this command frequently.

results <- downloadPc(format = "SIFNX")

Searching Pathway Commons

Networks can also be loaded using Pathway Commons rather than from local BioPAX files. First, we show how Pathway Commons can be searched.

searchResults <- searchPc("Q06609")

All functions that query Pathway Commons include a flag verbose that allows users to see the query URL sent to Pathway Commons for debugging purposes.

searchResults <- searchPc("Q06609", verbose = TRUE)

## URL:  http://purl.org/pc2/5/search.xml?q=Q06609&page=0

Pathway Commons search results are returned as an XML object.

str(searchResults)

## Classes 'XMLInternalDocument', 'XMLAbstractDocument' <externalptr>

These results can be filtered using the XML package using XPath expressions; examples of XPath expressions and syntax. The examples here shows how to pull the name for the pathway and the URI that contains information about the pathway in the BioPAX format.

xpathSApply(searchResults, "/searchResponse/searchHit/name", xmlValue)[1]

## [1] "DNA repair protein RAD51 homolog 1"

xpathSApply(searchResults, "/searchResponse/searchHit/pathway", xmlValue)[1]

## [1] "http://identifiers.org/reactome/REACT_111183.2"

Alternatively, these XML results can be converted to data.frames using the XML and plyr libraries.

library(plyr)
searchResultsDf <- ldply(xmlToList(searchResults), data.frame)

This type of searching can be used to locallly save BioPAX files retrieved from Pathway Commons.

## Use an XPath expression to extract the results of interest. In this case,
## the URIs (IDs) for the pathways from the results
searchResults <- xpathSApply(searchResults, "/searchResponse/searchHit/uri", 
    xmlValue)

## Generate temporary file to save content into
biopaxFile <- tempfile()

## Extract the URI for the first pathway in the search results and save into
## a file
uri <- searchResults[1]
saveXML(getPc(uri, "BIOPAX"), biopaxFile)

Extracting Information from BioPAX Datasets Using traverse()

The traverse function allows the extraction of specific entries from BioPAX records. traverse() functionality should be available for any uniprot.org or purl.org URI.

# Get the protein Uniprot ID for a given HGNC gene symbol
id <- idMapping("AKT1")

# Covert the Uniprot ID to a URI that would be found in Pathway Commons
uri <- paste0("http://identifiers.org/uniprot/", id)

# Get URIs for only the ModificationFeatures of the protein
xml <- traverse(uri = uri, path = "ProteinReference/entityFeature:ModificationFeature")

# Extract all the URIs
uris <- xpathSApply(xml, "//value/text()", xmlValue)

# For the first URI get the modification position and type
tmpXml <- traverse(uri = uris[1], path = "ModificationFeature/featureLocation:SequenceSite/sequencePosition")
cat("Modification Position: ", xpathSApply(tmpXml, "//value/text()", xmlValue))

## Modification Position:  129

tmpXml <- traverse(uri = uris[1], path = "ModificationFeature/modificationType/term")
cat("Modification Type: ", xpathSApply(tmpXml, "//value/text()", xmlValue))

## Modification Type:  O-phospho-L-serine

Common Data Visualization Pathways and Network Analysis

Visualizing SIF Interactions from Pathway Commons using R Graph Libraries

A common use case for paxtoolsr to retrieve a network or sub-network from a pathway derived from a BioPAX file or a Pathway Commons query. Next, we show how to visualize subnetworks loaded from BioPAX files and retrieved using the Pathway Commons webservice. To do this, we use the igraph R graph library because it has simple methods for loading edgelists, analyzing the networks, and visualizing these networks.

Next, we show how subnetworks queried from Pathway Commons can be visualized directly in R using the igraph library. Alternatively, these graphical plots can be made using Cytoscape either by exporting the SIF format and then importing the SIF format into Cytoscape or by using the RCytoscape package to work with Cytoscape directly from R.

library(igraph)

We load the network from a BioPAX file using the SIF format:

sif <- toSif(system.file("extdata", "biopax3-short-metabolic-pathway.owl", package = "paxtoolsr"))

# graph.edgelist requires a matrix
g <- graph.edgelist(as.matrix(sif[, c(1, 3)]), directed = FALSE)
plot(g, layout = layout.fruchterman.reingold)

Pathway Commons Graph Query

Next, we show how to perform graph search using Pathway Commons useful for finding connections and neighborhoods of elements. The example below shows the extraction of a sub-network connecting a set of proteins that is then filtered for a specific interaction type: “controls-state-change-of”. State change here indicates a modification of another molecule (e.g. post-translational modifications). This interaction type is directed, and it is read as “A controls a state change of B”.

genes <- c("AKT1", "IRS1", "MTOR", "IGF1R")
t1 <- graphPc(source = genes, kind = "PATHSBETWEEN", format = "BIOPAX")
t2 <- toSif(t1)
t3 <- t2[which(t2[, 2] == "controls-state-change-of"), ]

g <- graph.edgelist(as.matrix(t3[, c(1, 3)]), directed = FALSE)
plot(g, layout = layout.fruchterman.reingold)

Overlaying Experimental Data on Pathway Commons Networks

Often, it is useful not only to visualize a biological pathway, but also to overlay a given network with some form of biological data, such as gene expression values for genes in the network.

library(RColorBrewer)

# Generate a color palette that goes from white to red that contains 10
# colors
numColors <- 10
colors <- colorRampPalette(brewer.pal(9, "Reds"))(numColors)

# Generate values that could represent some experimental values
values <- runif(length(V(g)$name))

# Scale values to generate indicies from the color palette
xrange <- range(values)
newrange <- c(1, numColors)

factor <- (newrange[2] - newrange[1])/(xrange[2] - xrange[1])
scaledValues <- newrange[1] + (values - xrange[1]) * factor
indicies <- as.integer(scaledValues)

# Color the nodes based using the indicies and the color palette created
# above
g <- set.vertex.attribute(g, "color", value = colors[indicies])

# get.vertex.attribute(h, 'color')

plot(g, layout = layout.fruchterman.reingold)

Network Statistics

Often it is useful to produce statistics on a network. Here we show how to determine SIF network statistics and statistics on BioPAX files.

SIF Network Statistics

Once Pathway Commons and BioPAX networks are loaded as graphs using the igraph R package, it is possible to analyze these networks. Here we show how to get common statistics for the current network retrieved from Pathway Commons:

# Degree for each node in the igraph network
degree(g)

##    AKT1  AKT1S1  DEPTOR    IRS1   MLST8    MTOR   RPTOR    AKT2    AKT3 
##      11       2       2      11       6       7       3       9       5 
##    IGF1     CRK    IRS2    IRS4   IGF1R    IGF2   IGF2R     INS    INSR 
##       4       2       7       7       4       3       3       3       3 
##   INSRR MAPKAP1    PRR5  RICTOR 
##       3       3       1       3

# Number of nodes
length(V(g)$name)

## [1] 22

# Clustering coefficient
transitivity(g)

## [1] 0.04265403

# Network density
graph.density(g)

## [1] 0.2207792

# Network diameter
diameter(g)

## [1] 4

Another common task determine paths between nodes in a network.

# Get the first shortest path between two nodes
tmp <- get.shortest.paths(g, from = "IRS1", to = "MTOR")

# igraph seems to return different objects on Linux versus OS X for
# get.shortest.paths()
if (class(tmp[[1]]) == "list") {
    path <- tmp[[1]][[1]]  # Linux
} else {
    path <- tmp[[1]]  # OS X
}

# Convert from indicies to vertex names
V(g)$name[path]

## [1] "IRS1" "MTOR"

Gene Set Enrichment Analysis with Pathway Commons

The processing of the microarray data is taken from the following webpage: Bioconductor Tutorial on Microarray Processing and Gene Set Analysis with for grabbing gene sets from a Pathway Commons pathway and using same data as in the example, but stored in the estrogen R package.

To access microarray data sets, users should consider retrieving data from the NCBI Gene Expression Omnibus (GEO) using the GEOQuery package.

The first thing we'll do is load up the necessary packages.

library(paxtoolsr)  # To retrieve data from Pathway Commons
library(limma)  # Contains geneSetTest
library(affy)  # To load microarray data
library(hgu95av2)  # Annotation packages for the hgu95av2 platform
library(hgu95av2cdf)
library(XML)  # To parse XML files

We then retrieve a pathway of interest using the the Pathway Commons search functionality and convert this pathway to a gene set.

# Generate a Gene Set Search Pathway Commons for 'glycolysis'-related
# pathways
searchResults <- searchPc(q = "glycolysis", type = "pathway")

## Use an XPath expression to extract the results of interest. In this case,
## the URIs (IDs) for the pathways from the results
searchResults <- xpathSApply(searchResults, "/searchResponse/searchHit/uri", 
    xmlValue)

## Generate temporary files to save content into
biopaxFile <- tempfile()
gseaFile <- tempfile()

## Extract the URI for the first pathway in the search results and save into
## a file
uri <- searchResults[1]
saveXML(getPc(uri, "BIOPAX"), biopaxFile)

## Generate a gene set for the BioPAX pathway with gene symbols
tmp <- toGSEA(biopaxFile, gseaFile, "HGNC Symbol", FALSE)
geneSet <- tmp$geneSet

Finally, we process the microarray data and apply the gene set entrichment analysis.

# Process Microarray Data Load/process the estrogen microarray data
estrogenDataDir <- system.file("extdata", package = "estrogen")
targets <- readTargets(file.path(estrogenDataDir, "estrogen.txt"), sep = "")

abatch <- ReadAffy(filenames = file.path(estrogenDataDir, targets$filename))
eset <- rma(abatch)

## Background correcting
## Normalizing
## Calculating Expression

f <- paste(targets$estrogen, targets$time.h, sep = "")
f <- factor(f)
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)

fit <- lmFit(eset, design)

cont.matrix <- makeContrasts(E10 = "present10-absent10", E48 = "present48-absent48", 
    Time = "absent48-absent10", levels = design)

fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

## Map the gene symbols to the probe IDs
symbol <- unlist(as.list(hgu95av2SYMBOL))

### Check that the gene symbols are on the microarray platform
genesOnChip <- match(geneSet, symbol)
genesOnChip  # CHECK FOR ERROR HERE 

##  [1]   331  7748  3445  4057   614   834    NA  7635    NA  6589  3730
## [12] 11306  8934    NA  3583    NA    NA  7585

genesOnChip <- genesOnChip[!is.na(genesOnChip)]

### Grab the probe IDs for the genes present
genesOnChip <- names(symbol[genesOnChip])
genesOnChip <- match(genesOnChip, rownames(fit2$t))
genesOnChip <- genesOnChip[!is.na(genesOnChip)]

## Run the Gene Set Test from the limma Package
geneSetTest(genesOnChip, fit2$t[, 1], "two.sided")

## [1] 0.2391206

ID Mapping

Functions and results from paxtoolsr functions can be used in conjunction with the ID mapping functions of the biomaRt Bioconductor package; users should check the the biomaRt package documentation for a list of possible ID types.

sif <- toSif(system.file("extdata", "raf_map_kinase_cascade_reactome.owl", package = "paxtoolsr"))

# Generate a mapping between the HGNC symbols in the SIF to the Uniprot IDs
library(biomaRt)
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

hgnc_symbol <- c(sif$PARTICIPANT_A, sif$PARTICIPANT_B)
output <- getBM(attributes = c("hgnc_symbol", "uniprot_sptrembl"), filters = "hgnc_symbol", 
    values = hgnc_symbol, mart = ensembl)

# Remove blank entries
output <- output[output[, 2] != "", ]

head(output)

##   hgnc_symbol uniprot_sptrembl
## 1        RAF1           L7RRS6
## 3        RAF1           B4E1N6
## 4        RAF1           H7C155
## 5      MAP2K2           B3KS97
## 6      MAP2K2           G5E9C7
## 8        CDK1           E5RIU6

Troubleshooting

File Paths

Use properly delimited and full paths (do not use relative paths, such as ../directory/file or ~/directory/file) to files should be used with the paxtoolsr package.

toSif("/directory/file")
# or
toSif("X:\\directory\\file")

Memory Limits: Specify JVM Maximum Heap Size

By default paxtoolsr uses a maximum heap size limit of 512MB. For large BioPAX files, this limit may be insufficient. The code below shows how to change this limit and observe that the change was made.

NOTE: This limit cannot be changed once the virtual machine has been initialized by loading the library, so the memory heap size limit must be changed beforehand.

options(java.parameters = "-Xmx1024m")

library(paxtoolsr)

# Megabyte size
mbSize <- 1048576

runtime <- .jcall("java/lang/Runtime", "Ljava/lang/Runtime;", "getRuntime")
maxMemory <- .jcall(runtime, "J", "maxMemory")
maxMemoryMb <- maxMemory/mbSize
cat("Max Memory: ", maxMemoryMb, "\n")

Session Information

sessionInfo()

## R Under development (unstable) (2014-10-07 r66723)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] biomaRt_2.23.0      hgu95av2cdf_2.15.0  hgu95av2_2.2.0     
##  [4] affy_1.45.0         Biobase_2.27.0      BiocGenerics_0.13.0
##  [7] limma_3.23.1        RColorBrewer_1.0-5  igraph_0.7.1       
## [10] paxtoolsr_1.1.4     plyr_1.8.1          rjson_0.2.15       
## [13] RCurl_1.95-4.3      bitops_1.0-6        XML_3.98-1.1       
## [16] rJava_0.9-6         BiocStyle_1.5.2     knitr_1.7          
## 
## loaded via a namespace (and not attached):
##  [1] AnnotationDbi_1.29.1  BiocInstaller_1.17.1  DBI_0.3.1            
##  [4] GenomeInfoDb_1.3.6    IRanges_2.1.4         RSQLite_1.0.0        
##  [7] Rcpp_0.11.3           S4Vectors_0.5.2       affyio_1.35.0        
## [10] evaluate_0.5.5        formatR_1.0           preprocessCore_1.29.0
## [13] stats4_3.2.0          stringr_0.6.2         tools_3.2.0          
## [16] zlibbioc_1.13.0

References

• Cerami EG, Gross BE, Demir E, Rodchenkov I, Babur O, Anwar N, Schultz N, Bader GD, Sander C. Pathway Commons, a web resource for biological pathway data. Nucleic Acids Res. 2011 Jan;39(Database issue):D685-90. doi: 10.1093/nar/gkq1039. Epub 2010 Nov 10.
• Rodchenkov I, Demir E, Sander C, Bader GD. The BioPAX Validator. Bioinformatics. 2013 Oct 15;29(20):2659-60. doi: 10.1093/bioinformatics/btt452. Epub 2013 Aug 5.

Pathway Commons Graph Query

Next, we show how to perform graph search using Pathway Commons useful for finding connections and neighborhoods of elements. The example below shows the extraction of a sub-network connecting a set of proteins that is then filtered for a specific interaction type: “controls-state-change-of”. State change here indicates a modification of another molecule (e.g. post-translational modifications). This interaction type is directed, and it is read as “A controls a state change of B”.

genes <- c("AKT1", "IRS1", "MTOR", "IGF1R")
t1 <- graphPc(source = genes, kind = "PATHSBETWEEN", format = "BIOPAX")
t2 <- toSif(t1)
t3 <- t2[which(t2[, 2] == "controls-state-change-of"), ]

g <- graph.edgelist(as.matrix(t3[, c(1, 3)]), directed = FALSE)
plot(g, layout = layout.fruchterman.reingold)

Overlaying Experimental Data on Pathway Commons Networks

Often, it is useful not only to visualize a biological pathway, but also to overlay a given network with some form of biological data, such as gene expression values for genes in the network.

library(RColorBrewer)

# Generate a color palette that goes from white to red that contains 10
# colors
numColors <- 10
colors <- colorRampPalette(brewer.pal(9, "Reds"))(numColors)

# Generate values that could represent some experimental values
values <- runif(length(V(g)$name))

# Scale values to generate indicies from the color palette
xrange <- range(values)
newrange <- c(1, numColors)

factor <- (newrange[2] - newrange[1])/(xrange[2] - xrange[1])
scaledValues <- newrange[1] + (values - xrange[1]) * factor
indicies <- as.integer(scaledValues)

# Color the nodes based using the indicies and the color palette created
# above
g <- set.vertex.attribute(g, "color", value = colors[indicies])

# get.vertex.attribute(h, 'color')

plot(g, layout = layout.fruchterman.reingold)

Network Statistics

Often it is useful to produce statistics on a network. Here we show how to determine SIF network statistics and statistics on BioPAX files.

SIF Network Statistics

Once Pathway Commons and BioPAX networks are loaded as graphs using the igraph R package, it is possible to analyze these networks. Here we show how to get common statistics for the current network retrieved from Pathway Commons:

# Degree for each node in the igraph network
degree(g)

##    AKT1  AKT1S1  DEPTOR    IRS1   MLST8    MTOR   RPTOR    AKT2    AKT3 
##      11       2       2      11       6       7       3       9       5 
##    IGF1     CRK    IRS2    IRS4   IGF1R    IGF2   IGF2R     INS    INSR 
##       4       2       7       7       4       3       3       3       3 
##   INSRR MAPKAP1    PRR5  RICTOR 
##       3       3       1       3

# Number of nodes
length(V(g)$name)

## [1] 22

# Clustering coefficient
transitivity(g)

## [1] 0.04265403

# Network density
graph.density(g)

## [1] 0.2207792

# Network diameter
diameter(g)

## [1] 4

Another common task determine paths between nodes in a network.

# Get the first shortest path between two nodes
tmp <- get.shortest.paths(g, from = "IRS1", to = "MTOR")

# igraph seems to return different objects on Linux versus OS X for
# get.shortest.paths()
if (class(tmp[[1]]) == "list") {
    path <- tmp[[1]][[1]]  # Linux
} else {
    path <- tmp[[1]]  # OS X
}

# Convert from indicies to vertex names
V(g)$name[path]

## [1] "IRS1" "MTOR"

Gene Set Enrichment Analysis with Pathway Commons

The processing of the microarray data is taken from the following webpage: Bioconductor Tutorial on Microarray Processing and Gene Set Analysis with for grabbing gene sets from a Pathway Commons pathway and using same data as in the example, but stored in the estrogen R package.

To access microarray data sets, users should consider retrieving data from the NCBI Gene Expression Omnibus (GEO) using the GEOQuery package.

The first thing we'll do is load up the necessary packages.

library(paxtoolsr)  # To retrieve data from Pathway Commons
library(limma)  # Contains geneSetTest
library(affy)  # To load microarray data
library(hgu95av2)  # Annotation packages for the hgu95av2 platform
library(hgu95av2cdf)
library(XML)  # To parse XML files

We then retrieve a pathway of interest using the the Pathway Commons search functionality and convert this pathway to a gene set.

# Generate a Gene Set Search Pathway Commons for 'glycolysis'-related
# pathways
searchResults <- searchPc(q = "glycolysis", type = "pathway")

## Use an XPath expression to extract the results of interest. In this case,
## the URIs (IDs) for the pathways from the results
searchResults <- xpathSApply(searchResults, "/searchResponse/searchHit/uri", 
    xmlValue)

## Generate temporary files to save content into
biopaxFile <- tempfile()
gseaFile <- tempfile()

## Extract the URI for the first pathway in the search results and save into
## a file
uri <- searchResults[1]
saveXML(getPc(uri, "BIOPAX"), biopaxFile)

## Generate a gene set for the BioPAX pathway with gene symbols
tmp <- toGSEA(biopaxFile, gseaFile, "HGNC Symbol", FALSE)
geneSet <- tmp$geneSet

Finally, we process the microarray data and apply the gene set entrichment analysis.

# Process Microarray Data Load/process the estrogen microarray data
estrogenDataDir <- system.file("extdata", package = "estrogen")
targets <- readTargets(file.path(estrogenDataDir, "estrogen.txt"), sep = "")

abatch <- ReadAffy(filenames = file.path(estrogenDataDir, targets$filename))
eset <- rma(abatch)

## Background correcting
## Normalizing
## Calculating Expression

f <- paste(targets$estrogen, targets$time.h, sep = "")
f <- factor(f)
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)

fit <- lmFit(eset, design)

cont.matrix <- makeContrasts(E10 = "present10-absent10", E48 = "present48-absent48", 
    Time = "absent48-absent10", levels = design)

fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

## Map the gene symbols to the probe IDs
symbol <- unlist(as.list(hgu95av2SYMBOL))

### Check that the gene symbols are on the microarray platform
genesOnChip <- match(geneSet, symbol)
genesOnChip  # CHECK FOR ERROR HERE 

##  [1]   331 11456  4897  7748  3445  4057  4643   614   144  2901   834
## [12]    NA  7635  5840  6589  3730   935 11306  8934  3583    NA    NA
## [23]    NA  7585

genesOnChip <- genesOnChip[!is.na(genesOnChip)]

### Grab the probe IDs for the genes present
genesOnChip <- names(symbol[genesOnChip])
genesOnChip <- match(genesOnChip, rownames(fit2$t))
genesOnChip <- genesOnChip[!is.na(genesOnChip)]

## Run the Gene Set Test from the limma Package
geneSetTest(genesOnChip, fit2$t[, 1], "two.sided")

## [1] 0.08986801

ID Mapping

Functions and results from paxtoolsr functions can be used in conjunction with the ID mapping functions of the biomaRt Bioconductor package; users should check the the biomaRt package documentation for a list of possible ID types.

sif <- toSif(system.file("extdata", "raf_map_kinase_cascade_reactome.owl", package = "paxtoolsr"))

# Generate a mapping between the HGNC symbols in the SIF to the Uniprot IDs
library(biomaRt)
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

hgnc_symbol <- c(sif$PARTICIPANT_A, sif$PARTICIPANT_B)
output <- getBM(attributes = c("hgnc_symbol", "uniprot_sptrembl"), filters = "hgnc_symbol", 
    values = hgnc_symbol, mart = ensembl)

# Remove blank entries
output <- output[output[, 2] != "", ]

head(output)

##   hgnc_symbol uniprot_sptrembl
## 1        RAF1           L7RRS6
## 3        RAF1           B4E1N6
## 4        RAF1           H7C155
## 5      MAP2K2           B3KS97
## 6      MAP2K2           G5E9C7
## 8        CDK1           E5RIU6

Troubleshooting

File Paths

Use properly delimited and full paths (do not use relative paths, such as ../directory/file or ~/directory/file) to files should be used with the paxtoolsr package.

toSif("/directory/file")
# or
toSif("X:\\directory\\file")

Memory Limits: Specify JVM Maximum Heap Size

By default paxtoolsr uses a maximum heap size limit of 512MB. For large BioPAX files, this limit may be insufficient. The code below shows how to change this limit and observe that the change was made.

NOTE: This limit cannot be changed once the virtual machine has been initialized by loading the library, so the memory heap size limit must be changed beforehand.

options(java.parameters = "-Xmx1024m")

library(paxtoolsr)

# Megabyte size
mbSize <- 1048576

runtime <- .jcall("java/lang/Runtime", "Ljava/lang/Runtime;", "getRuntime")
maxMemory <- .jcall(runtime, "J", "maxMemory")
maxMemoryMb <- maxMemory/mbSize
cat("Max Memory: ", maxMemoryMb, "\n")

Session Information

sessionInfo()

## R Under development (unstable) (2014-11-23 r67046)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
##  [3] LC_TIME=en_US.UTF-8           LC_COLLATE=C                 
##  [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
##  [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
##  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
## [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] biomaRt_2.23.5      hgu95av2cdf_2.15.0  hgu95av2_2.2.0     
##  [4] affy_1.45.0         Biobase_2.27.0      BiocGenerics_0.13.3
##  [7] limma_3.23.3        RColorBrewer_1.1-2  igraph_0.7.1       
## [10] paxtoolsr_1.1.5     plyr_1.8.1          rjson_0.2.15       
## [13] RCurl_1.95-4.5      bitops_1.0-6        XML_3.98-1.1       
## [16] rJava_0.9-6         BiocStyle_1.5.3     knitr_1.8          
## 
## loaded via a namespace (and not attached):
##  [1] AnnotationDbi_1.29.11 BiocInstaller_1.17.1  DBI_0.3.1            
##  [4] GenomeInfoDb_1.3.10   IRanges_2.1.33        RSQLite_1.0.0        
##  [7] Rcpp_0.11.3           S4Vectors_0.5.14      affyio_1.35.0        
## [10] evaluate_0.5.5        formatR_1.0           preprocessCore_1.29.0
## [13] stats4_3.2.0          stringr_0.6.2         tools_3.2.0          
## [16] zlibbioc_1.13.0

References

• Cerami EG, Gross BE, Demir E, Rodchenkov I, Babur O, Anwar N, Schultz N, Bader GD, Sander C. Pathway Commons, a web resource for biological pathway data. Nucleic Acids Res. 2011 Jan;39(Database issue):D685-90. doi: 10.1093/nar/gkq1039. Epub 2010 Nov 10.
• Rodchenkov I, Demir E, Sander C, Bader GD. The BioPAX Validator. Bioinformatics. 2013 Oct 15;29(20):2659-60. doi: 10.1093/bioinformatics/btt452. Epub 2013 Aug 5.

library(RColorBrewer)

# Generate a color palette that goes from white to red that contains 10
# colors
numColors <- 10
colors <- colorRampPalette(brewer.pal(9, "Reds"))(numColors)

# Generate values that could represent some experimental values
values <- runif(length(V(g)$name))

# Scale values to generate indicies from the color palette
xrange <- range(values)
newrange <- c(1, numColors)

factor <- (newrange[2] - newrange[1])/(xrange[2] - xrange[1])
scaledValues <- newrange[1] + (values - xrange[1]) * factor
indicies <- as.integer(scaledValues)

# Color the nodes based using the indicies and the color palette created
# above
g <- set.vertex.attribute(g, "color", value = colors[indicies])

# get.vertex.attribute(h, 'color')

plot(g, layout = layout.fruchterman.reingold)

# Degree for each node in the igraph network
degree(g)

##    AKT1  AKT1S1  DEPTOR    IRS1   MLST8    MTOR   RPTOR    AKT2    AKT3 
##      11       2       2      11       6       7       3       9       5 
##    IGF1     CRK    IRS2    IRS4   IGF1R    IGF2   IGF2R     INS    INSR 
##       4       2       7       7       4       3       3       3       3 
##   INSRR MAPKAP1    PRR5  RICTOR 
##       3       3       1       3

# Number of nodes
length(V(g)$name)

## [1] 22

# Clustering coefficient
transitivity(g)

## [1] 0.04265403

# Network density
graph.density(g)

## [1] 0.2207792

# Network diameter
diameter(g)

## [1] 4

# Get the first shortest path between two nodes
tmp <- get.shortest.paths(g, from = "IRS1", to = "MTOR")

# igraph seems to return different objects on Linux versus OS X for
# get.shortest.paths()
if (class(tmp[[1]]) == "list") {
    path <- tmp[[1]][[1]]  # Linux
} else {
    path <- tmp[[1]]  # OS X
}

# Convert from indicies to vertex names
V(g)$name[path]

## [1] "IRS1" "MTOR"

library(paxtoolsr)  # To retrieve data from Pathway Commons
library(limma)  # Contains geneSetTest
library(affy)  # To load microarray data
library(hgu95av2)  # Annotation packages for the hgu95av2 platform
library(hgu95av2cdf)
library(XML)  # To parse XML files

# Generate a Gene Set Search Pathway Commons for 'glycolysis'-related
# pathways
searchResults <- searchPc(q = "glycolysis", type = "pathway")

## Use an XPath expression to extract the results of interest. In this case,
## the URIs (IDs) for the pathways from the results
searchResults <- xpathSApply(searchResults, "/searchResponse/searchHit/uri", 
    xmlValue)

## Generate temporary files to save content into
biopaxFile <- tempfile()
gseaFile <- tempfile()

## Extract the URI for the first pathway in the search results and save into
## a file
uri <- searchResults[1]
saveXML(getPc(uri, "BIOPAX"), biopaxFile)

## Generate a gene set for the BioPAX pathway with gene symbols
tmp <- toGSEA(biopaxFile, gseaFile, "HGNC Symbol", FALSE)
geneSet <- tmp$geneSet

# Process Microarray Data Load/process the estrogen microarray data
estrogenDataDir <- system.file("extdata", package = "estrogen")
targets <- readTargets(file.path(estrogenDataDir, "estrogen.txt"), sep = "")

abatch <- ReadAffy(filenames = file.path(estrogenDataDir, targets$filename))
eset <- rma(abatch)

## Background correcting
## Normalizing
## Calculating Expression

f <- paste(targets$estrogen, targets$time.h, sep = "")
f <- factor(f)
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)

fit <- lmFit(eset, design)

cont.matrix <- makeContrasts(E10 = "present10-absent10", E48 = "present48-absent48", 
    Time = "absent48-absent10", levels = design)

fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

## Map the gene symbols to the probe IDs
symbol <- unlist(as.list(hgu95av2SYMBOL))

### Check that the gene symbols are on the microarray platform
genesOnChip <- match(geneSet, symbol)
genesOnChip  # CHECK FOR ERROR HERE 

##  [1]   331 11456  4897  7748  3445  4057  4643   614   144  2901   834
## [12]    NA  7635  5840  6589  3730   935 11306  8934  3583    NA    NA
## [23]    NA  7585

genesOnChip <- genesOnChip[!is.na(genesOnChip)]

### Grab the probe IDs for the genes present
genesOnChip <- names(symbol[genesOnChip])
genesOnChip <- match(genesOnChip, rownames(fit2$t))
genesOnChip <- genesOnChip[!is.na(genesOnChip)]

## Run the Gene Set Test from the limma Package
geneSetTest(genesOnChip, fit2$t[, 1], "two.sided")

## [1] 0.08986801

sif <- toSif(system.file("extdata", "raf_map_kinase_cascade_reactome.owl", package = "paxtoolsr"))

# Generate a mapping between the HGNC symbols in the SIF to the Uniprot IDs
library(biomaRt)
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

hgnc_symbol <- c(sif$PARTICIPANT_A, sif$PARTICIPANT_B)
output <- getBM(attributes = c("hgnc_symbol", "uniprot_sptrembl"), filters = "hgnc_symbol", 
    values = hgnc_symbol, mart = ensembl)

# Remove blank entries
output <- output[output[, 2] != "", ]

head(output)

##   hgnc_symbol uniprot_sptrembl
## 1        RAF1           L7RRS6
## 3        RAF1           B4E1N6
## 4        RAF1           H7C155
## 5      MAP2K2           G5E9C7
## 6      MAP2K2           B3KS97
## 8        CDK1           E5RIU6

toSif("/directory/file")
# or
toSif("X:\\directory\\file")

options(java.parameters = "-Xmx1024m")

library(paxtoolsr)

# Megabyte size
mbSize <- 1048576

runtime <- .jcall("java/lang/Runtime", "Ljava/lang/Runtime;", "getRuntime")
maxMemory <- .jcall(runtime, "J", "maxMemory")
maxMemoryMb <- maxMemory/mbSize
cat("Max Memory: ", maxMemoryMb, "\n")

sessionInfo()

## R Under development (unstable) (2015-01-31 r67686)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.1 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
##  [3] LC_TIME=en_US.UTF-8           LC_COLLATE=C                 
##  [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
##  [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
##  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
## [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] biomaRt_2.23.5      hgu95av2cdf_2.15.0  hgu95av2_2.2.0     
##  [4] affy_1.45.2         Biobase_2.27.1      BiocGenerics_0.13.4
##  [7] limma_3.23.8        RColorBrewer_1.1-2  igraph_0.7.1       
## [10] paxtoolsr_1.1.7     plyr_1.8.1          rjson_0.2.15       
## [13] RCurl_1.95-4.5      bitops_1.0-6        XML_3.98-1.1       
## [16] rJava_0.9-6         BiocStyle_1.5.3     knitr_1.9          
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.11.4           BiocInstaller_1.17.5  formatR_1.0          
##  [4] GenomeInfoDb_1.3.12   tools_3.2.0           zlibbioc_1.13.1      
##  [7] digest_0.6.8          evaluate_0.5.5        RSQLite_1.0.0        
## [10] preprocessCore_1.29.0 DBI_0.3.1             yaml_2.1.13          
## [13] stringr_0.6.2         IRanges_2.1.38        S4Vectors_0.5.19     
## [16] stats4_3.2.0          AnnotationDbi_1.29.17 rmarkdown_0.5.1      
## [19] htmltools_0.2.6       affyio_1.35.0

# Generate a Gene Set Search Pathway Commons for 'glycolysis'-related
# pathways
searchResults <- searchPc(q = "glycolysis", type = "pathway")

## Use an XPath expression to extract the results of interest. In this case,
## the URIs (IDs) for the pathways from the results
searchResults <- xpathSApply(searchResults, "/searchResponse/searchHit/uri", 
    xmlValue)

## Generate temporary files to save content into
biopaxFile <- tempfile()

## Extract the URI for the first pathway in the search results and save into
## a file
uri <- searchResults[2]
saveXML(getPc(uri, "BIOPAX"), biopaxFile)

And then, we convert this pathway to a gene set.

## Generate temporary files to save content into
gseaFile <- tempfile()

## Generate a gene set for the BioPAX pathway with gene symbols NOTE: Not all
## search results are guaranteed to result in gene sets
tmp <- toGSEA(biopaxFile, gseaFile, "HGNC Symbol", FALSE)
geneSet <- tmp$geneSet

Finally, we process the microarray data and apply the gene set entrichment analysis.

# Process Microarray Data Load/process the estrogen microarray data
estrogenDataDir <- system.file("extdata", package = "estrogen")
targets <- readTargets(file.path(estrogenDataDir, "estrogen.txt"), sep = "")

abatch <- ReadAffy(filenames = file.path(estrogenDataDir, targets$filename))
eset <- rma(abatch)

## Background correcting
## Normalizing
## Calculating Expression

f <- paste(targets$estrogen, targets$time.h, sep = "")
f <- factor(f)
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)

fit <- lmFit(eset, design)

cont.matrix <- makeContrasts(E10 = "present10-absent10", E48 = "present48-absent48", 
    Time = "absent48-absent10", levels = design)

fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

## Map the gene symbols to the probe IDs
symbol <- unlist(as.list(hgu95av2SYMBOL))

### Check that the gene symbols are on the microarray platform
genesOnChip <- match(geneSet, symbol)
genesOnChip  # CHECK FOR ERROR HERE 

##  [1]   331  7748  3445  4057   614   834    NA  7635    NA  6589  3730
## [12] 11306  8934    NA  3583    NA    NA  7585

genesOnChip <- genesOnChip[!is.na(genesOnChip)]

### Grab the probe IDs for the genes present
genesOnChip <- names(symbol[genesOnChip])
genesOnChip <- match(genesOnChip, rownames(fit2$t))
genesOnChip <- genesOnChip[!is.na(genesOnChip)]

## Run the Gene Set Test from the limma Package
geneSetTest(genesOnChip, fit2$t[, 1], "two.sided")

## [1] 0.2391206

sif <- toSif(system.file("extdata", "raf_map_kinase_cascade_reactome.owl", package = "paxtoolsr"))

# Generate a mapping between the HGNC symbols in the SIF to the Uniprot IDs
library(biomaRt)
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

hgnc_symbol <- c(sif$PARTICIPANT_A, sif$PARTICIPANT_B)
output <- getBM(attributes = c("hgnc_symbol", "uniprot_sptrembl"), filters = "hgnc_symbol", 
    values = hgnc_symbol, mart = ensembl)

# Remove blank entries
output <- output[output[, 2] != "", ]

toSif("/directory/file")
# or
toSif("X:\\directory\\file")

options(java.parameters = "-Xmx1024m")

library(paxtoolsr)

# Megabyte size
mbSize <- 1048576

runtime <- .jcall("java/lang/Runtime", "Ljava/lang/Runtime;", "getRuntime")
maxMemory <- .jcall(runtime, "J", "maxMemory")
maxMemoryMb <- maxMemory/mbSize
cat("Max Memory: ", maxMemoryMb, "\n")

sessionInfo()

## R version 3.2.0 alpha (2015-03-20 r68043)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
##  [3] LC_TIME=en_US.UTF-8           LC_COLLATE=C                 
##  [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
##  [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
##  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
## [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] hgu95av2cdf_2.15.0   hgu95av2_2.2.0       affy_1.45.2         
##  [4] Biobase_2.27.3       BiocGenerics_0.13.10 limma_3.23.11       
##  [7] RColorBrewer_1.1-2   igraph_0.7.1         paxtoolsr_1.1.11    
## [10] plyr_1.8.1           rjson_0.2.15         RCurl_1.95-4.5      
## [13] bitops_1.0-6         XML_3.98-1.1         rJava_0.9-6         
## [16] BiocStyle_1.5.3      knitr_1.9           
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.11.5           BiocInstaller_1.17.6  formatR_1.1          
##  [4] GenomeInfoDb_1.3.16   tools_3.2.0           zlibbioc_1.13.3      
##  [7] digest_0.6.8          evaluate_0.5.5        RSQLite_1.0.0        
## [10] preprocessCore_1.29.0 DBI_0.3.1             yaml_2.1.13          
## [13] stringr_0.6.2         IRanges_2.1.43        S4Vectors_0.5.22     
## [16] stats4_3.2.0          AnnotationDbi_1.29.20 rmarkdown_0.5.1      
## [19] htmltools_0.2.6       affyio_1.35.0

# Degree for each node in the igraph network
degree(g)

##    AKT1  AKT1S1  DEPTOR    IRS1   MLST8    MTOR   RPTOR    AKT2    AKT3 
##      11       2       2      13       6       7       3       9       5 
##    IGF1     CRK    IRS2    IRS4   IGF1R    IGF2   IGF2R     INS    INSR 
##       4       3       8       7       4       3       3       3       3 
##   INSRR MAPKAP1    PRR5  RICTOR 
##       3       3       1       3

# Number of nodes
length(V(g)$name)

## [1] 22

# Clustering coefficient
transitivity(g)

## [1] 0.1481481

# Network density
graph.density(g)

## [1] 0.2294372

# Network diameter
diameter(g)

## [1] 4

# Get the first shortest path between two nodes
tmp <- get.shortest.paths(g, from = "IRS1", to = "MTOR")

# igraph seems to return different objects on Linux versus OS X for
# get.shortest.paths()
if (class(tmp[[1]]) == "list") {
    path <- tmp[[1]][[1]]  # Linux
} else {
    path <- tmp[[1]]  # OS X
}

# Convert from indicies to vertex names
V(g)$name[path]

## [1] "IRS1" "MTOR"

library(paxtoolsr)  # To retrieve data from Pathway Commons
library(limma)  # Contains geneSetTest
library(affy)  # To load microarray data
library(hgu95av2)  # Annotation packages for the hgu95av2 platform
library(hgu95av2cdf)
library(XML)  # To parse XML files

# Generate a Gene Set Search Pathway Commons for 'glycolysis'-related
# pathways
searchResults <- searchPc(q = "glycolysis", type = "pathway")

## Use an XPath expression to extract the results of interest. In this case,
## the URIs (IDs) for the pathways from the results
searchResults <- xpathSApply(searchResults, "/searchResponse/searchHit/uri", 
    xmlValue)

## Generate temporary files to save content into
biopaxFile <- tempfile()

## Extract the URI for the first pathway in the search results and save into
## a file
uri <- searchResults[2]
saveXML(getPc(uri, "BIOPAX"), biopaxFile)

## Generate temporary files to save content into
gseaFile <- tempfile()

## Generate a gene set for the BioPAX pathway with gene symbols NOTE: Not all
## search results are guaranteed to result in gene sets
tmp <- toGSEA(biopaxFile, gseaFile, "HGNC Symbol", FALSE)
geneSet <- tmp$geneSet

# Process Microarray Data Load/process the estrogen microarray data
estrogenDataDir <- system.file("extdata", package = "estrogen")
targets <- readTargets(file.path(estrogenDataDir, "estrogen.txt"), sep = "")

abatch <- ReadAffy(filenames = file.path(estrogenDataDir, targets$filename))
eset <- rma(abatch)

## Background correcting
## Normalizing
## Calculating Expression

f <- paste(targets$estrogen, targets$time.h, sep = "")
f <- factor(f)
design <- model.matrix(~0 + f)
colnames(design) <- levels(f)

fit <- lmFit(eset, design)

cont.matrix <- makeContrasts(E10 = "present10-absent10", E48 = "present48-absent48", 
    Time = "absent48-absent10", levels = design)

fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

## Map the gene symbols to the probe IDs
symbol <- unlist(as.list(hgu95av2SYMBOL))

### Check that the gene symbols are on the microarray platform
genesOnChip <- match(geneSet, symbol)
genesOnChip  # CHECK FOR ERROR HERE 

##  [1]   331  7748  3445  4057   614   834    NA  7635    NA  6589  3730
## [12] 11306  8934    NA  3583    NA    NA  7585

genesOnChip <- genesOnChip[!is.na(genesOnChip)]

### Grab the probe IDs for the genes present
genesOnChip <- names(symbol[genesOnChip])
genesOnChip <- match(genesOnChip, rownames(fit2$t))
genesOnChip <- genesOnChip[!is.na(genesOnChip)]

## Run the Gene Set Test from the limma Package
geneSetTest(genesOnChip, fit2$t[, 1], "two.sided")

## [1] 0.2391206

sif <- toSif(system.file("extdata", "raf_map_kinase_cascade_reactome.owl", package = "paxtoolsr"))

# Generate a mapping between the HGNC symbols in the SIF to the Uniprot IDs
library(biomaRt)
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

hgnc_symbol <- c(sif$PARTICIPANT_A, sif$PARTICIPANT_B)
output <- getBM(attributes = c("hgnc_symbol", "uniprot_sptrembl"), filters = "hgnc_symbol", 
    values = hgnc_symbol, mart = ensembl)

# Remove blank entries
output <- output[output[, 2] != "", ]

toSif("/directory/file")
# or
toSif("X:\\directory\\file")

options(java.parameters = "-Xmx1024m")

library(paxtoolsr)

# Megabyte size
mbSize <- 1048576

runtime <- .jcall("java/lang/Runtime", "Ljava/lang/Runtime;", "getRuntime")
maxMemory <- .jcall(runtime, "J", "maxMemory")
maxMemoryMb <- maxMemory/mbSize
cat("Max Memory: ", maxMemoryMb, "\n")

sessionInfo()

## R version 3.2.0 (2015-04-16)
## Platform: x86_64-unknown-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
##  [3] LC_TIME=en_US.UTF-8           LC_COLLATE=C                 
##  [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
##  [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
##  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
## [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] hgu95av2cdf_2.16.0  hgu95av2_2.2.0      affy_1.46.0        
##  [4] Biobase_2.28.0      BiocGenerics_0.14.0 limma_3.24.5       
##  [7] RColorBrewer_1.1-2  igraph_0.7.1        paxtoolsr_1.2.1    
## [10] plyr_1.8.2          rjson_0.2.15        RCurl_1.95-4.6     
## [13] bitops_1.0-6        XML_3.98-1.2        rJava_0.9-6        
## [16] BiocStyle_1.6.0     knitr_1.10.5       
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.11.6           BiocInstaller_1.18.2  formatR_1.2          
##  [4] GenomeInfoDb_1.4.0    tools_3.2.0           zlibbioc_1.14.0      
##  [7] digest_0.6.8          evaluate_0.7          RSQLite_1.0.0        
## [10] preprocessCore_1.30.0 DBI_0.3.1             yaml_2.1.13          
## [13] stringr_1.0.0         IRanges_2.2.3         S4Vectors_0.6.0      
## [16] stats4_3.2.0          AnnotationDbi_1.30.1  rmarkdown_0.6.1      
## [19] magrittr_1.5          htmltools_0.2.6       stringi_0.4-1        
## [22] affyio_1.36.0