Contents
Overview
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
source("http://bioconductor.org/biocLite.R")
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 BioPAX and Pathway Commons 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", "biopax3-short-metabolic-pathway.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] "2"
## [1] "7"
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)
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.
## Search Pathway Commons for 'glycolysis'-related pathways
searchResults <- searchPc(q = "glycolysis", type = "pathway")
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.
## Search Pathway Commons for 'glycolysis'-related pathways
searchResults <- searchPc(q = "glycolysis", type = "pathway", verbose = TRUE)
## URL: http://www.pathwaycommons.org/pc2/search.xml?q=glycolysis&page=0&type=pathway
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] "Glycolysis / Gluconeogenesis"
xpathSApply(searchResults, "/searchResponse/searchHit/pathway", xmlValue)[1]
## [1] "http://identifiers.org/kegg.pathway/hsa00020"
Alternatively, these XML results can be converted to data.frames using the XML and plyr libraries.
library(plyr)
searchResultsDf <- ldply(xmlToList(searchResults), data.frame)
# Simplified results
simplifiedSearchResultsDf <- searchResultsDf[, c("name", "uri", "biopaxClass")]
This type of searching can be used to locally 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 a URI for a pathway in the search results and save into a file
idx <- which(grepl("panther", simplifiedSearchResultsDf$uri) & grepl("glycolysis",
simplifiedSearchResultsDf$name, ignore.case = TRUE))
uri <- simplifiedSearchResultsDf$uri[idx]
saveXML(getPc(uri, "BIOPAX"), biopaxFile)
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. This can be used to extract the neighborhood of a single gene 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”.
gene <- "BDNF"
t1 <- graphPc(source = gene, kind = "neighborhood", format = "BINARY_SIF", verbose = TRUE)
## URL: http://www.pathwaycommons.org/pc2/graph?kind=neighborhood&source=BDNF&format=BINARY_SIF
t2 <- t1[which(t1[, 2] == "controls-state-change-of"), ]
g <- graph.edgelist(as.matrix(t2[, c(1, 3)]), directed = FALSE)
plot(g, layout = layout.fruchterman.reingold)
The example below shows the extraction of a sub-network connecting a set of proteins:
genes <- c("AKT1", "IRS1", "MTOR", "IGF1R")
t1 <- graphPc(source = genes, kind = "PATHSBETWEEN", format = "BINARY_SIF",
verbose = TRUE)
## URL: http://www.pathwaycommons.org/pc2/graph?kind=PATHSBETWEEN&source=AKT1&source=IRS1&source=MTOR&source=IGF1R&format=BINARY_SIF
t2 <- t1[which(t1[, 2] == "controls-state-change-of"), ]
g <- graph.edgelist(as.matrix(t2[, 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)
## ACVR1 AKT1 AKT2 AKT3 ACVR1B ACVR2A ACVR2B ACVRL1
## 3 187 176 172 3 3 3 3
## ADRBK1 ADRBK2 AKT1S1 DEPTOR IRS1 LAMTOR1 LAMTOR2 LAMTOR3
## 3 3 2 2 13 1 1 1
## LAMTOR4 LAMTOR5 MLST8 MTOR RHEB RPTOR RRAGA RRAGB
## 1 1 6 7 1 3 1 1
## RRAGC RRAGD ALK AMHR2 ARAF BMPR1A BMPR1B BMPR2
## 1 1 3 3 3 3 3 3
## BRAF BTK CAMK2A CAMK2B CAMK2D CAMK2G CCNH CDC7
## 3 3 3 3 3 3 3 3
## CDK1 CDK2 CDK3 CDK4 CDK5 CDK6 CDK7 CDK8
## 3 3 3 3 3 3 3 3
## CDK9 CHUK CSK CSNK1A1 CSNK1D CSNK1E CSNK1G1 CSNK1G2
## 3 3 3 3 3 3 3 3
## CSNK1G3 EGFR FGFR1 FGFR2 FGFR3 FGFR4 FLT1 FLT4
## 3 3 3 3 3 3 3 3
## FYN GFRA1 GRK4 GRK5 GRK6 GRK7 GSK3A GSK3B
## 3 3 3 3 3 3 3 3
## IGF1 CRK IRS2 IRS4 IGF1R IGF2 IGF2R IKBKB
## 4 3 8 7 7 3 3 3
## IKBKE IKBKG INS INSR INSRR IRAK1 IRAK2 IRAK4
## 3 3 3 6 3 3 3 3
## ITK JAK1 JAK2 JAK3 KDR KIT LCK LTK
## 3 3 3 3 3 3 3 3
## MAP2K1 MAP2K2 MAP2K3 MAP2K4 MAP2K5 MAP2K6 MAP2K7 MAP3K10
## 3 3 3 3 3 3 3 3
## MAP3K11 MAP3K12 MAP3K13 MAP3K14 MAP3K1 MAP3K2 MAP3K3 MAP3K5
## 3 3 3 3 3 3 3 3
## MAP3K7 MAP3K8 MAP3K9 MAP4K2 MAP4K3 MAPK10 MAPK11 MAPK12
## 3 3 3 3 3 3 3 3
## MAPK13 MAPK14 MAPK1 MAPK3 MAPK7 MAPK8 MAPK9 MAPKAP1
## 3 3 3 3 3 3 3 3
## MAPKAPK2 MAPKAPK3 MAPKAPK5 MET MKNK1 MNAT1 MOS NGFR
## 3 3 3 3 3 3 3 3
## NLK NTRK1 PAK1 PAK2 PAK4 PAK6 PAK7 PDGFRA
## 3 3 3 3 3 3 3 3
## PDGFRB PDPK1 PLK1 PLK3 PRKAA1 PRKAA2 PRKAB1 PRKAB2
## 3 3 3 3 3 3 3 3
## PRKACA PRKACB PRKACG PRKAG1 PRKAG2 PRKAG3 PRKAR1A PRKAR1B
## 3 3 3 3 3 3 3 3
## PRKAR2A PRKAR2B PRKCA PRKCB PRKCD PRKCE PRKCG PRKCH
## 3 3 3 3 3 3 3 3
## PRKCI PRKCQ PRKCZ PRKG1 PRKG2 PRR5 PTK2 PTK2B
## 3 3 3 3 3 3 3 3
## RAF1 RICTOR RIPK1 RIPK2 RIPK3 RPS6KA1 RPS6KA2 RPS6KA3
## 3 3 3 3 3 3 3 3
## RPS6KA6 RPS6KB1 RPS6KB2 RYK SRC SYK TEC TGFBR1
## 3 3 3 3 3 3 3 3
## TGFBR2 TYK2 WEE1 ZAP70
## 3 3 3 3
# Number of nodes
length(V(g)$name)
## [1] 196
# Clustering coefficient
transitivity(g)
## [1] 0.0008877991
# Network density
graph.density(g)
## [1] 0.02946102
# Network diameter
diameter(g)
## [1] 3
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)
## Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail'
## when loading 'hgu95av2cdf'
## Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head'
## when loading 'hgu95av2cdf'
library(XML) # To parse XML files
We then retrieve a pathway of interest using the the Pathway Commons search functionality.
# 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
## integer(0)
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] 2
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] != "", ]
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")
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.