inst/doc/NPMVignette.R

In NetPathMiner: NetPathMiner for Biological Network Construction, Path Mining and Visualization

## ---- echo=FALSE, results="hide", warning=FALSE-------------------------------
suppressPackageStartupMessages({library('NetPathMiner')})

## ----Load_package, echo=TRUE, eval=TRUE, results="hide"-----------------------
library(NetPathMiner)

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  graph <- KGML2igraph(filename = file)
#  graph <- SBML2igraph(filename = file)

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  library(rBiopaxParser)
#  biopax = readBiopax(file)
#  graph <- BioPAX2igraph(biopax = biopax)

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  graph <- KGML2igraph(filename = c(file1, file2))

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  graph <- KGML2igraph(filename = ".")

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  # Extract all MIRIAM identifiers from an SBML file.
#  graph <- SBML2igraph(filename = file, miriam = "all")
#  
#  # Extract only miram.go identifiers from a BioPAX file.
#  graph <- BioPAX2igraph(biopax = biopax, miriam = "go")

## ---- echo=FALSE, eval=TRUE, results="hide"-----------------------------------
file <- file.path(find.package("NetPathMiner"), "extdata", "hsa00860.xml")

## ---- echo=TRUE, eval=FALSE, results="hide"-----------------------------------
#  graph <- KGML2igraph(filename = file, parse.as = "signaling")
#  
#  graph <- KGML2igraph(filename = file, parse.as = "signaling",
#  	expand.complexes = TRUE)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
data("ex_sbml")
graph <- ex_sbml
graph

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
head( V(graph) )

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
head( E(graph) )

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
head( V(graph)[ reactions ] )

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
V(graph)[ "reaction_71850" ]$attr

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
getAttrNames(graph)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
getAttrStatus(graph, pattern = "^miriam.")

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  require("RCurl")
#  # Fetch uniprot annotation
#  graph <- fetchAttribute(graph, organism = "Homo sapiens", target.attr = "miriam.ncbigene" , source.attr = "miriam.uniprot")
#  
#  # Fetch ChEBI annotation.
#  graph <- fetchAttribute(graph, target.attr = "miriam.chebi", source.attr = "miriam.kegg.compound")

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
rgraph <- makeReactionNetwork(graph, simplify=FALSE)
rgraph

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  rgraph <- simplifyReactionNetwork(rgraph)
#  rgraph <- makeReactionNetwork(graph, simplify=TRUE)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
# Expand complexes of gene network.
ggraph <- expandComplexes(rgraph, v.attr = "miriam.uniprot",
		keep.parent.attr= c("^pathway", "^compartment"))

# Convert reaction network to gene network.
ggraph <- makeGeneNetwork(rgraph)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
data(ex_microarray)


## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  # Assign weights to edges.
#  if(require("RCurl") && url.exists( NPMdefaults("bridge.web") ))
#  	rgraph <- fetchAttribute(rgraph, organism = "Homo sapiens",
#  						target.attr = "miriam.affy.probeset",
#  						source.attr = "miriam.uniprot")

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  # This requires an internet connection, and RCurl and ALL packages to be present.
#  # Instead, we will actually use a processed ALL data, where features are converted
#  # to miriam.uniprot annotation. (Next chunk)
#  
#  library(ALL)
#  data(ALL)
#  rgraph <- assignEdgeWeights(microarray = exprs(ALL), graph = rgraph,
#  weight.method = "cor", use.attr="miriam.affy.probeset", y=ALL$mol.bio, bootstrap = FALSE)

## ---- echo=FALSE, eval=TRUE---------------------------------------------------
# This is what is evaluated.
data(ex_microarray)
rgraph <- assignEdgeWeights(microarray = ex_microarray, graph = rgraph,
weight.method = "cor", use.attr="miriam.uniprot", y=colnames(ex_microarray), bootstrap = FALSE)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
rgraph$y.labels
head( E(rgraph)$edge.weights )

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
ranked.p <- pathRanker(rgraph, method = "prob.shortest.path",
	K = 25, minPathSize = 6)

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  pathsample <- samplePaths(rgraph, max.path.length = vcount(rgraph),
#  num.samples = 1000, num.warmup = 10)
#  
#  ranked.p <- pathRanker(rgraph, method = "pvalue",
#  sampledpaths = pathsample ,alpha=0.1)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
# Get paths as edge IDs.
eids <- getPathsAsEIDs(paths = ranked.p, graph = rgraph)

## ---- echo=TRUE, eval=TRUE, results="hide"------------------------------------
# Convert paths to other networks.
eids <- getPathsAsEIDs(paths = ranked.p, graph = ggraph)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
# Clustering.
ybinpaths <- pathsToBinary(ranked.p)
p.cluster <- pathCluster(ybinpaths, M = 2)

## ---- fig=TRUE, pdf=TRUE, echo=TRUE, eval=TRUE--------------------------------
plotClusters(ybinpaths, p.cluster)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
p.class <- pathClassifier(ybinpaths, target.class = "BCR/ABL", M = 2)

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  plotClassifierROC(p.class)

## ---- fig=TRUE, pdf=TRUE, echo=TRUE, eval=TRUE--------------------------------
plotClusters(ybinpaths, p.class)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
plotNetwork(rgraph, vertex.color="compartment.name")

## ---- fig=TRUE, pdf=TRUE, echo=TRUE, eval=FALSE-------------------------------
#  plotPaths(ranked.p, rgraph)
#  
#  # With clusters
#  plotPaths(ranked.p, graph, path.clusters=p.class)

## ---- fig=TRUE, pdf=TRUE, echo=TRUE, eval=TRUE--------------------------------
plotAllNetworks(ranked.p, metabolic.net = graph, reaction.net = rgraph,
		path.clusters=p.class, vertex.label = "", vertex.size = 4)

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  layout.c <- clusterVertexByAttr(rgraph, "pathway", cluster.strength = 3)
#  v.color <- colorVertexByAttr(rgraph, "pathway")
#  plotPaths(ranked.p , rgraph, clusters=p.class,
#  	layout = layout.c, vertex.color = v.color)

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  plotCytoscapeGML(graph, file="example.gml", layout = layout.c,
#  				vertex.size = 5, vertex.color = v.color)

## ---- echo=TRUE, eval=TRUE, results="hide"------------------------------------
getGeneSets(graph, use.attr="compartment", gene.attr="miriam.uniprot")

## ---- echo=TRUE, eval=TRUE, results="hide"------------------------------------
getGeneSetNetworks(graph, use.attr="compartment")

## ---- echo=TRUE, eval=FALSE---------------------------------------------------
#  graphNEL <- toGraphNEL(graph, export.attr="^miriam.")