## ----include=FALSE-------------------------------------------------------
library(knitr)
options(width=64,digits=2)
opts_chunk$set(size="small")
opts_chunk$set(tidy=TRUE,tidy.opts=list(width.cutoff=50,keep.blank.line=TRUE))
opts_knit$set(eval.after='fig.cap')
# for a package vignette, you do want to echo.
# opts_chunk$set(echo=FALSE,warning=FALSE,message=FALSE)
opts_chunk$set(warning=FALSE,message=FALSE)
opts_chunk$set(cache=TRUE,cache.path="cache/mmnet")

## ----install-pkg, eval=FALSE---------------------------------------------
## ## install release version of mmnet
## source("http://bioconductor.org/biocLite.R")
## biocLite("mmnet")
## 
## ##install the latest development version
## useDevel()
## biocLite("mmnet")

## ----load-pkg,eval=TRUE, include=FALSE-----------------------------------
library(mmnet)

## ----load-pkg2,eval=FALSE------------------------------------------------
## library(mmnet)

## ----sample-download, eval=FALSE-----------------------------------------
## download.file("ftp://ftp.metagenomics.anl.gov/projects/10/4440616.3/raw/507.fna.gz",
##               destfile = "obesesample.fna.gz")
## download.file("ftp://ftp.metagenomics.anl.gov/projects/10/4440823.3/raw/687.fna.gz",
##               destfile = "leansample.fna.gz")

## ----login-mgrast, eval = FALSE, echo=TRUE-------------------------------
## ## login on MG-RAST
##  login.info <- loginMgrast(user="mmnet",userpwd="mmnet")

## ----upload-mgrast, eval=FALSE,echo=TRUE---------------------------------
## ## select the sample data to upload
## seq.file <- c("obesesample.fna.gz", "leansample.fna.gz")
## ## upload sequence data to MG-RAST
## metagenome.id <- lapply(seq.file, uploadMgrast, login.info = login.info)

## ----submit-mgrast, eval=FALSE, echo=TRUE--------------------------------
## ## submit MGRAST job
##   metagenome.id <- lapply(seq.file, submitMgrastJob, login.info = login.info)
##   show(metagenome.id)

## ----check-metagenome, eval=FALSE,echo=TRUE------------------------------
## ## check MGRAST project status
##   metagenome.status <- lapply(metagenome.id, checkMgrastMetagenome, login.info = login.info)
## ## apparently, status of completed annotation metagenome is TRUE
##   show(metagenome.status)

## ----getannotation,eval=FALSE, echo=TRUE---------------------------------
## ## private data
##   private.annotation <- lapply(metagenome.id, getMgrastAnnotation, login.info=login.info)
## ## public annotation data, does not require login.info
##   public.annotation <- lapply(metagenome.id, getMgrastAnnotation)

## ----data-load-----------------------------------------------------------
data(anno)
summary(anno)

## ----MG-RAST-mmnet, eval=FALSE-------------------------------------------
## ## first login on MG-RAST
## login.info <- loginMgrast("mmnet", "mmnet")
## ## prepare the metagenomic sequence for annotation
## seq <- "obesesample.fna.gz"
## ## mgrast annotation
## uploadMgrast(login.info, seq)
## metagenome.id2 <- submitMgrastJob(login.info, seqfile = basename(seq))
## while (TRUE) {
##   status <- checkMgrastMetagenome(metagenome.id = metagenome.id2)
##   if (status)
##     break
##   Sys.sleep(5)
##   cat("In annotation, please waiting...")
##   flush.console()
## }
## ## if annotation profile is public,take login.info as NULL
## anno2 <- getMgrastAnnotation(metagenome.id2, login.info=login.info)

## ----estimate, echo=TRUE-------------------------------------------------
 mmnet.abund <- estimateAbundance(anno)
 show(mmnet.abund)

## ----MG-RAST-BIOM, fig.align='center',fig.cap='enzymatic genes abundance comparison between MG-RAST and mmnet package',dev='pdf',fig.show='hold',out.width='.7\\linewidth', out.height='.7\\linewidth'----
## download BIOM functional abundance profile of the two sample metagenome from MG-RAST
if(require(RCurl)){
  function.api <- "http://api.metagenomics.anl.gov/1/matrix/function"
  mgrast.abund <- read_biom(getForm(function.api,
                             .params=list(id="mgm4440616.3",id="mgm4440823.3",
                                          source="KO",result_type="abundance")))
}
## obtain the intersect ko abundance of MG-RAST and esimatiAbundance
intersect.ko <- intersect(rownames(mgrast.abund),rownames(mmnet.abund))
## compare the two by taking one metagenome
mgrast.abund1 <- biom_data(mgrast.abund)[,1][intersect.ko]
mmnet.abund1 <- biom_data(mmnet.abund)[,1][intersect.ko]
if(require(ggplot2)){
  p <- qplot(mgrast.abund1,mmnet.abund1) + 
            geom_abline(slope=1, intercept=0,color="blue") +
            ylim(0, 400) + xlim(0, 400)
  print(p)
}

## ----IMGSample, eval=FALSE-----------------------------------------------
## ## Load the IMG/M sample data
## IMGFile <- system.file("extdata/IMGSample.tab",package="mmnet")
## IMGSample <- read.delim2(IMGFile,quote = "")
## ## Create BIOM file for network construction
## abundance <- IMGSample$Gene.Count
## abundance <- abundance/sum(abundance)
## abundance <- as.data.frame(abundance)
## KO <- IMGSample$KO.ID
## KO <- as.data.frame(gsub("KO:","",KO))
## biom.data <- make_biom(abundance, observation_metadata = KO)
## biom.data$type <- "enzymatic genes abundance"

## ----IMGSample2, eval=FALSE----------------------------------------------
## ## Construct and analyze SSN
## ssn <- constructSSN(biom.data)
## topologicalAnalyzeNet(ssn)

## ----initial-data,echo=TRUE----------------------------------------------
loadMetabolicData()
summary(RefDbcache)

## ----construct-network,echo=TRUE-----------------------------------------
  ssn <- constructSSN(mmnet.abund)
  g <- ssn[[1]]
  summary(g)
  abund <- get.vertex.attribute(g,"abundance",index=V(g))
  summary(abund)

## ----topologicalAnalyzeNet,echo=TRUE,fig.align='center',fig.cap='Topological metabolic network analysis, linking topological features and enzymatic gene abundances',dev='pdf',fig.show='hold',out.width='.7\\linewidth', out.height='.7\\linewidth'----
topo.net <- topologicalAnalyzeNet(g)
## network with topological features as attributes
topo.net

## ----differential,echo=TRUE,fig.align='center',fig.cap='Differential metabolic network analysis, enzymatic genes that are associated with specific state appear as colored nodes (red=enriched, green=depleted)',dev='pdf',fig.show='hold',out.width='.7\\linewidth', out.height='.7\\linewidth'----
state <- c("obese", "lean")
differential.net <- differentialAnalyzeNet(ssn, sample.state= state, method="OR", cutoff = c(0.5, 2))
summary(differential.net)

## ----showMetagenomicNet,echo=TRUE,fig.align='center',fig.cap='Visualization of State Specific Network  using \textit{showMetagenomicNet} with node size  proportional to log(abundance)',dev='pdf',fig.show='hold', out.width='.7\\linewidth', out.height='.7\\linewidth'----
## the reference network
# showMetagenomicNet(RefDbcache$network,mode="ref")
#the state specific metabolic network
showMetagenomicNet(g, mode="ssn", vertex.label = NA, edge.width = 0.3, edge.arrow.size = 0.1, edge.arrow.width = 0.1, layout = layout.fruchterman.reingold)

## ----RCytoscape, eval = FALSE, echo = TRUE-------------------------------
## if (require(RCytoscape)){
##   refnet <- ssn[[1]]
##   net <- igraph.to.graphNEL(refnet)
##   ## initialize the edge attibute
##   #edge.attr=list.edge.attributes(refnet)
##   #edge.attr.class = sapply(edge.attr, class)
##   #edge.attr.class[edge.attr.class=="character"]="char"
##   ## init node attributes
##   node.attr=list.vertex.attributes(refnet)
##   if (length(node.attr)){
##     node.attr.class = sapply(node.attr, class)
##     node.attr.class[node.attr.class=="character"]="char"
##     for (i in 1:length(node.attr))
##       net<- initNodeAttribute(net, attribute.name = node.attr[i],
##                             attribute.type = node.attr.class[i], default.value = "0")
##   }
##   ## our metagenomic network does not have edge attributes, set them all to 1
##   net <- initEdgeAttribute(net, attribute.name = "weight",
##                               attribute.type = "numeric", default.value = "1")
##   ## create a network window in Cytoscape
##   cw <- new.CytoscapeWindow ('net', graph = net, overwriteWindow = TRUE)
##   ## transmits the CytoscapeWindowClass's graph data, from R to Cytoscape, nodes,
##   ## edges, node and edge attributes
##   displayGraph (cw)
## }

## ----echo=FALSE----------------------------------------------------------
sessionInfo()

## ----closeConnetions-----------------------------------------------------
allCon <- showConnections()
socketCon <- as.integer(rownames(allCon)[allCon[, "class"] == "sockconn"])
sapply(socketCon, function(ii) close.connection(getConnection(ii)) )