## ----setup, echo = FALSE, eval=TRUE, message=FALSE----------------------------
library(BiocStyle)
knitr::opts_knit$set(
  upload.fun = NULL,
  base.url = NULL) # use local files for images
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#"
)
# Override BiocStyle plot hook to avoid css_align issues
knitr::knit_hooks$set(plot = function(x, options) {
  paste0('![', basename(x), '](', x, ')')
})
runchunks = if (Sys.getenv("FORCE_VIGNETTE_REBUILD", "FALSE") == "TRUE") TRUE else FALSE

cache_file <- '../vignette_cache/SpiecEasi.RData'
if (!runchunks && file.exists(cache_file)) {
  load(cache_file)
  # If we loaded trimmed objects, reassign them to original names
  if (exists("se.mb.trimmed")) {
    se.mb.amgut <- se.mb.trimmed
    rm(se.mb.trimmed)
  }
  if (exists("se.gl.trimmed")) {
    se.gl.amgut <- se.gl.trimmed
    rm(se.gl.trimmed)
  }
  if (exists("se.trimmed")) {
    se <- se.trimmed
    rm(se.trimmed)
  }
  message("Loaded cached data from SpiecEasi.RData")
} else {
  if (!runchunks) {
    message("Cache file SpiecEasi.RData not found - building from scratch")
  }
  runchunks <- TRUE
}
saveout   = runchunks


## ----eval=FALSE---------------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("SpiecEasi")

## ----eval=FALSE---------------------------------------------------------------
# if (!require("remotes", quietly = TRUE))
#   install.packages("remotes")
# remotes::install_github("zdk123/SpiecEasi")

## ----eval=TRUE----------------------------------------------------------------
sessionInfo()

## ----eval=TRUE----------------------------------------------------------------
library(SpiecEasi)
data(amgut1.filt)

## ----eval=runchunks-----------------------------------------------------------
# depths <- rowSums(amgut1.filt)
# amgut1.filt.n  <- t(apply(amgut1.filt, 1, norm_to_total))
# amgut1.filt.cs <- round(amgut1.filt.n * min(depths))
# 
# d <- ncol(amgut1.filt.cs)
# n <- nrow(amgut1.filt.cs)
# e <- d

## ----eval=runchunks-----------------------------------------------------------
# set.seed(10010)
# graph <- SpiecEasi::make_graph('cluster', d, e)
# Prec  <- graph2prec(graph)
# Cor   <- cov2cor(prec2cov(Prec))
# 
# X <- synth_comm_from_counts(amgut1.filt.cs, mar=2, distr='zinegbin', Sigma=Cor, n=n)

## ----eval=runchunks-----------------------------------------------------------
# se <- spiec.easi(X, method='mb', lambda.min.ratio=1e-2, nlambda=15)

## ----eval=TRUE, fig.width=5, fig.height=5-------------------------------------
huge::huge.roc(se$est$path, graph, verbose=FALSE)
stars.pr(getOptMerge(se), graph, verbose=FALSE)
# stars selected final network under: getRefit(se)

## ----eval=runchunks, message=FALSE, warning=FALSE-----------------------------
# se.mb.amgut <- spiec.easi(amgut1.filt, method='mb', lambda.min.ratio=1e-2,
#                           nlambda=10, pulsar.params=list(rep.num=20))
# se.gl.amgut <- spiec.easi(amgut1.filt, method='glasso', lambda.min.ratio=1e-2,
#                           nlambda=10, pulsar.params=list(rep.num=20))
# sparcc.amgut <- sparcc(amgut1.filt)

## ----eval=TRUE----------------------------------------------------------------
## Define arbitrary threshold for SparCC correlation matrix for the graph
sparcc.graph <- abs(sparcc.amgut$Cor) >= 0.3
diag(sparcc.graph) <- 0
library(Matrix)
sparcc.graph <- Matrix(sparcc.graph, sparse=TRUE)
## Create igraph objects
ig.mb     <- adj2igraph(getRefit(se.mb.amgut))
ig.gl     <- adj2igraph(getRefit(se.gl.amgut))
ig.sparcc <- adj2igraph(sparcc.graph)

## ----eval=TRUE, fig.width=15, fig.height=7, message=FALSE---------------------
library(igraph)
library(Matrix)
## set size of vertex proportional to clr-mean
vsize    <- rowMeans(clr(amgut1.filt, 1))+6
am.coord <- layout.fruchterman.reingold(ig.mb)

par(mfrow=c(1,3))
plot(ig.mb, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="MB")
plot(ig.gl, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="glasso")
plot(ig.sparcc, layout=am.coord, vertex.size=vsize, vertex.label=NA, main="sparcc")

## ----eval=TRUE, fig.width=8, fig.height=5-------------------------------------
# Create edge lists for weight comparison
secor  <- cov2cor(getOptCov(se.gl.amgut))
sebeta <- symBeta(getOptBeta(se.mb.amgut), mode='maxabs')
elist.gl     <- summary(triu(secor*getRefit(se.gl.amgut), k=1))
elist.mb     <- summary(sebeta)
elist.sparcc <- summary(sparcc.graph*sparcc.amgut$Cor)

# Plot edge weight distributions
hist(elist.sparcc[,3], main='', xlab='edge weights')
hist(elist.mb[,3], add=TRUE, col='forestgreen')
hist(elist.gl[,3], add=TRUE, col='red')

## ----eval=TRUE, fig.width=9, fig.height=6-------------------------------------
dd.gl     <- degree.distribution(ig.gl)
dd.mb     <- degree.distribution(ig.mb)
dd.sparcc <- degree.distribution(ig.sparcc)

plot(0:(length(dd.sparcc)-1), dd.sparcc, ylim=c(0,.35), type='b',
      ylab="Frequency", xlab="Degree", main="Degree Distributions")
points(0:(length(dd.gl)-1), dd.gl, col="red" , type='b')
points(0:(length(dd.mb)-1), dd.mb, col="forestgreen", type='b')
legend("topright", c("MB", "glasso", "sparcc"),
        col=c("forestgreen", "red", "black"), pch=1, lty=1)

## ----echo = FALSE, eval=TRUE, message=FALSE-----------------------------------
# Save only the expensive computations if they exist and we're rebuilding
if (runchunks && exists("se.mb.amgut") && exists("se.gl.amgut")) {
  cache_file <- '../vignette_cache/SpiecEasi.RData'
  tryCatch({
    # Load trimming function and trim objects to reduce size
    source('../vignette_cache/trim_spiec_easi.R')
    se.mb.trimmed <- trim_spiec_easi(se.mb.amgut)
    se.gl.trimmed <- trim_spiec_easi(se.gl.amgut)
    se.trimmed <- trim_spiec_easi(se)
    
    # Save the trimmed expensive spiec.easi results plus other needed variables
    save(se.mb.trimmed, se.gl.trimmed, se.trimmed, sparcc.amgut, graph, X, Prec, Cor, 
         depths, amgut1.filt.cs, d, n, e, file=cache_file)
    message("Saved trimmed cache file: ", cache_file, " in directory: ", getwd())
  }, error = function(e) {
    message("Failed to save cache file: ", e$message)
  })
}