#' @title Group Comparison of Network Properties
#'
#' @description Calculate and compare network properties for microbial networks
#'   using Jaccard's index, the Rand index, the Graphlet Correlation Distance, 
#'   and permutation tests.
#'
#' @param x object of class \code{microNetProps} (returned by 
#'   \code{\link[NetCoMi]{netAnalyze}}).
#' @param permTest logical. If \code{TRUE}, a permutation test is conducted 
#'   to test centrality measures and global network properties for group
#'   differences. Defaults to \code{FALSE}. May lead to a considerably increased
#'   execution time!
#' @param jaccQuant numeric value between 0 and 1 specifying the quantile 
#'   used as threshold to identify the most central nodes for each centrality
#'   measure. The resulting sets of nodes are used to calculate Jaccard's index
#'   (see details). Default is 0.75.
#' @param lnormFit logical indicating whether a log-normal distribution should
#'   be fitted to the calculated centrality values for determining Jaccard's
#'   index (see details). If \code{NULL} (default), the value is adopted
#'   from the input, i.e., equals the method used for determining hub nodes.
#' @param testRand logical. If \code{TRUE}, a permutation test is conducted 
#'   for the adjusted Rand index (with H0: ARI = 0). Execution time may be 
#'   increased for large networks.
#' @param nPermRand integer giving the number of permutations used for testing 
#'   the adjusted Rand index for being significantly different from zero. 
#'   Ignored if \code{testRand = FALSE}. Defaults to 1000L.
#' @param gcd logical. If \code{TRUE} (default), the Graphlet Correlation 
#'   Distance (GCD) is computed.
#' @param gcdOrb numeric vector with integers from 0 to 14 defining the orbits 
#'   used for calculating the GCD. Minimum length is 2. 
#'   Defaults to c(0, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11), 
#'   thus excluding redundant orbits such as the orbit o3.
#' @param verbose logical. If \code{TRUE} (default), status messages are shown.
#' @param nPerm integer giving the number of permutations if 
#'   \code{permTest = TRUE}. Default is 1000L.
#' @param adjust character indicating the method used for multiple testing
#'   adjustment of the permutation p-values. Possible values are \code{"lfdr"}
#'   (default) for local false discovery rate correction (via
#'   \code{\link[fdrtool]{fdrtool}}), \code{"adaptBH"} for the adaptive
#'   Benjamini-Hochberg method \cite{(Benjamini and Hochberg, 2000)},  or one of
#'   the methods provided by \code{\link[stats]{p.adjust}} (see
#'   \code{p.adjust.methods()}).
#' @param trueNullMethod character indicating the method used for estimating the
#'   proportion of true null hypotheses from a vector of p-values. Used for the
#'   adaptive Benjamini-Hochberg method for multiple testing adjustment (chosen
#'   by \code{adjust = "adaptBH"}). Accepts the provided options of the
#'   \code{method} argument of \code{\link[limma]{propTrueNull}}:
#'   \code{"convest"}(default), \code{"lfdr"}, \code{"mean"}, and \code{"hist"}.
#'   Can alternatively be \code{"farco"} for the "iterative plug-in method"
#'   proposed by \cite{Farcomeni (2007)}.
#' @param cores integer indicating the number of CPU cores used for
#'   permutation tests. If cores > 1, the tests are performed in parallel.
#'   Is limited to the number of available CPU cores determined by
#'   \code{\link[parallel]{detectCores}}. Defaults to 1L (no parallelization).
#' @param logFile character string naming the log file to which the current 
#'   iteration number is written (if permutation tests are performed). Defaults
#'   to \code{NULL} so that no log file is generated.
#' @param seed integer giving a seed for reproducibility of the results.
#' @param fileLoadAssoPerm character giving the name or path (without file 
#'   extension)  of the file containing the "permuted" association/dissimilarity 
#'   matrices that was generated by setting \code{storeAssoPerm} to 
#'   \code{TRUE}. Only used for permutation tests. If \code{NULL}, no 
#'   existing associations are used. 
#' @param fileLoadCountsPerm character giving the name or path (without file 
#'   extension)  of the file containing the "permuted" count matrices that was 
#'   generated by setting \code{storeCountsPerm} to \code{TRUE}. 
#'   Only used for permutation tests, and if \code{fileLoadAssoPerm = NULL}. 
#'   If \code{NULL}, no existing count matrices are used.
#' @param storeAssoPerm logical indicating whether the association/dissimilarity 
#'   matrices for the permuted data should be saved to a file.
#'   The file name is given via \code{fileStoreAssoPerm}. If \code{TRUE}, 
#'   the computed "permutation" association/dissimilarity matrices can be reused
#'   via \code{fileLoadAssoPerm} to save runtime. Defaults to \code{FALSE}.
#'   Ignored if \code{fileLoadAssoPerm} is not \code{NULL}.
#' @param fileStoreAssoPerm character giving the name of a file to which the 
#'   matrix with associations/dissimilarities of the permuted data is saved. 
#'   Can also be a path.
#' @param storeCountsPerm logical indicating whether the permuted count matrices
#'   should be saved to an external file. Defaults to \code{FALSE}.
#'   Ignored if \code{fileLoadCountsPerm} is not \code{NULL}.
#' @param fileStoreCountsPerm character vector with two elements giving the 
#'   names of two files storing the permuted count matrices belonging to the 
#'   two groups.
#' @param returnPermProps logical. If \code{TRUE}, the global properties and 
#'   their absolute differences for the permuted data are returned.
#' @param returnPermCentr logical. If \code{TRUE}, the centralities and 
#'   their absolute differences for the permuted data are returned.
#' @param assoPerm only needed for output generated with NetCoMi v1.0.1! A list 
#'   with two elements used for the permutation procedure.
#'   Each entry must contain association matrices for \code{"nPerm"}
#'   permutations. This can be the \code{"assoPerm"} value as part of the
#'   output either returned by \code{diffnet} or \code{\link{netCompare}}. 
#' @param dissPerm only needed for output generated with NetCoMi v1.0.1! 
#'   Usage analog to \code{assoPerm} if a dissimilarity measure has been used 
#'   for network construction.
#' @details
#'   \strong{Permutation procedure:}\cr
#'   Used for testing centrality measures and global network properties for
#'   group differences.\cr
#'   The null hypothesis of the tests is defined as
#'   \deqn{H_0: c1_i - c2_i = 0,} where \eqn{c1_i} and
#'   \eqn{c2_i} denote the centrality values of taxon i in group 1 and 2,
#'   respectively.\cr
#'   To generate a sampling distribution of the differences under \eqn{H_0},
#'   the group labels are randomly reassigned to the samples while the group
#'   sizes are kept. The associations are then re-estimated for each permuted
#'   data set. The p-values are calculated as the proportion of
#'   "permutation-differences" being larger than or equal to the observed 
#'   difference. In non-exact tests, a pseudo-count is added to the numerator 
#'   and denominator to avoid p-values of zero. Several methods for adjusting 
#'   the p-values for multiplicity are available.
#'
#'   \strong{Jaccard's index:}\cr
#'   Jaccard's index expresses for each centrality measure how equal the sets of
#'   most central nodes are among the two networks.\cr
#'   These sets are defined as nodes with a centrality value above a defined
#'   quantile (via \code{jaccQuant}) either of the empirical distribution of the
#'   centrality values (\code{lnormFit = FALSE}) or of a fitted log-normal
#'   distribution (\code{lnormFit = TRUE}).\cr
#'   The index ranges from 0 to 1, where 1 means the sets of most central nodes
#'   are exactly equal in both networks and 0 indicates that the
#'   most central nodes are completely different.\cr
#'   The index is calculated as suggested by \cite{Real and Vargas (1996)}.
#'   \cr\cr
#'   \strong{Rand index:}\cr
#'   The Rand index is used to express whether the determined clusterings are
#'   equal in both groups. The adjusted Rand index (ARI) ranges from -1 to 1,
#'   where 1 indicates that the two clusterings are exactly equal. The expected
#'   index value for two random clusterings is 0. The implemented test procedure
#'   is in accordance with the explanations in \cite{Qannari et al. (2014)},
#'   where a p-value below the alpha levels means that ARI is significantly
#'   higher than expected for two random clusterings.
#'   \cr\cr
#'   \strong{Graphlet Correlation Distance:}\cr
#'   A graphlet-based distance measure, which is defined as the Euclidean 
#'   distance of the upper triangle values of the Graphlet Correlation 
#'   Matrices (GCM) of two networks \cite{(Yaveroglu et al., 2014)}.
#'   The GCM of a network is a matrix with Spearman's correlations between the 
#'   network's node orbits \cite{(Hocevar and Demsar, 2016)}.
#'   See \code{\link{calcGCD}} for details.
#'   
#' @return Object of class \code{microNetComp} with the following
#'   elements:\cr
#'   \tabular{ll}{
#'   \code{jaccDeg,jaccBetw,jaccClose,jaccEigen}\tab Values of Jaccard's index
#'   for the centrality measures\cr
#'   \code{jaccHub}\tab Jaccard index for the sets of hub nodes\cr
#'   \code{randInd}\tab Adjusted Rand index\cr
#'   \code{randIndLCC}\tab Adjusted Rand index for the largest connected 
#'   component (LCC)\cr
#'   \code{gcd}\tab Graphlet Correlation Distance (object of class \code{gcd} 
#'   returned by \code{\link{calcGCD}})\cr
#'   \code{gcdLCC}\tab Graphlet Correlation Distance for the LCC\cr
#'   \code{properties}\tab List with calculated network properties\cr
#'   \code{propertiesLCC}\tab List with calculated network properties of the 
#'   LCC\cr
#'   \code{diffGlobal}\tab Vectors with differences of global properties\cr
#'   \code{diffGlobalLCC}\tab Vectors with differences of global properties for
#'   the LCC\cr
#'   \code{diffCent}\tab Vectors with differences of the centrality values\cr
#'   \code{countMatrices}\tab The two count matrices returned
#'   by \code{netConstruct}\cr
#'   \code{assoMatrices}\tab The two association matrices returned
#'   by \code{netConstruct}\cr
#'   \code{dissMatrices}\tab The two dissimilarity matrices returned
#'   by \code{netConstruct}\cr
#'   \code{adjaMatrices}\tab The two adjacency matrices returned
#'   by \code{netConstruct}\cr
#'   \code{groups}\tab Group names returned by \code{netConstruct}\cr
#'   \code{paramsProperties}\tab Parameters used for network analysis
#'   }
#'   \strong{Additional output if permutation tests are conducted:}
#'  \tabular{ll}{
#'  \code{pvalDiffGlobal}\tab P-values of the tests for differential global
#'  properties\cr
#'  \code{pvalDiffGlobalLCC}\tab P-values of the tests for differential
#'  global properties in the LCC\cr
#'  \code{pvalDiffCentr}\tab P-values of the tests for differential centrality
#'  values\cr
#'  \code{pvalDiffCentrAdjust}\tab Adjusted p-values of the tests for
#'  differential centrality values\cr
#'  \code{permDiffGlobal}\tab \code{nPerm} x 10 matrix containing the absolute 
#'  differences of the ten global network properties (computed for the whole 
#'  network) for all \code{nPerm} permutations\cr
#'  \code{permDiffGlobalLCC}\tab \code{nPerm} x 11 matrix containing the 
#'  absolute differences of the eleven global network properties (computed for 
#'  the LCC) for all \code{nPerm} permutations\cr
#'  \code{permDiffCentr}\tab List with absolute differences of the four 
#'  centrality measures for all \code{nPerm} permutations. Each list contains 
#'  a \code{nPerm} x \code{nNodes} matrix.
#'  }
#' @examples
#' 
#' # Load data sets from American Gut Project (from SpiecEasi package)
#' data("amgut2.filt.phy")
#'
#' # Split data into two groups: with and without seasonal allergies
#' amgut_season_yes <- phyloseq::subset_samples(amgut2.filt.phy, 
#'                                       SEASONAL_ALLERGIES == "yes")
#' amgut_season_no <- phyloseq::subset_samples(amgut2.filt.phy, 
#'                                      SEASONAL_ALLERGIES == "no")
#'
#' amgut_season_yes
#' amgut_season_no
#' 
#' # Filter the 121 samples (sample size of the smaller group) with highest 
#' # frequency to make the sample sizes equal and thus ensure comparability.
#' n_yes <- phyloseq::nsamples(amgut_season_yes)
#'
#' # Network construction
#' amgut_net <- netConstruct(data = amgut_season_yes,
#'                           data2 = amgut_season_no,
#'                           measure = "pearson",
#'                           filtSamp = "highestFreq",
#'                           filtSampPar = list(highestFreq = n_yes),
#'                           filtTax = "highestVar",
#'                           filtTaxPar = list(highestVar = 30),
#'                           zeroMethod = "pseudoZO", normMethod = "clr")
#'
#' # Network analysis
#' # Note: Please zoom into the GCM plot or open a new window using:
#' # x11(width = 10, height = 10)
#' amgut_props <- netAnalyze(amgut_net, clustMethod = "cluster_fast_greedy")
#'
#' # Network plot
#' plot(amgut_props,
#'      sameLayout = TRUE,
#'      title1 = "Seasonal allergies",
#'      title2 = "No seasonal allergies")
#'
#' #--------------------------
#' # Network comparison
#'
#' # Without permutation tests
#' amgut_comp1 <- netCompare(amgut_props, permTest = FALSE)
#' summary(amgut_comp1)
#'
#' \donttest{
#'   # With permutation tests (with only 100 permutations to decrease runtime)
#'   amgut_comp2 <- netCompare(amgut_props,
#'                             permTest = TRUE,
#'                             nPerm = 100L,
#'                             cores = 1L,
#'                             storeCountsPerm = TRUE,
#'                             fileStoreCountsPerm = c("countsPerm1",
#'                                                     "countsPerm2"),
#'                             storeAssoPerm = TRUE,
#'                             fileStoreAssoPerm = "assoPerm",
#'                             seed = 123456)
#'
#' # Rerun with a different adjustment method ...
#' # ... using the stored permutation count matrices
#' amgut_comp3 <- netCompare(amgut_props, adjust = "BH",
#'                           permTest = TRUE, nPerm = 100L,
#'                           fileLoadCountsPerm = c("countsPerm1",
#'                                                  "countsPerm2"),
#'                           seed = 123456)
#' 
#' # ... using the stored permutation association matrices
#' amgut_comp4 <- netCompare(amgut_props, adjust = "BH",
#'                           permTest = TRUE, nPerm = 100L, 
#'                           fileLoadAssoPerm = "assoPerm",
#'                           seed = 123456)
#'   
#' # amgut_comp3 and amgut_comp4 should be equal
#' all.equal(amgut_comp3$adjaMatrices, amgut_comp4$adjaMatrices)
#' all.equal(amgut_comp3$properties, amgut_comp4$properties)
#' 
#' summary(amgut_comp2)
#' summary(amgut_comp3)
#' summary(amgut_comp4)
#' 
#' #--------------------------
#' # Use 'createAssoPerm' to create "permuted" count and association matrices
#' createAssoPerm(amgut_props, nPerm = 100, 
#'                computeAsso = TRUE,
#'                fileStoreAssoPerm = "assoPerm",
#'                storeCountsPerm = TRUE, 
#'                fileStoreCountsPerm = c("countsPerm1", "countsPerm2"),
#'                append = FALSE, seed = 123456)
#' 
#' amgut_comp5 <- netCompare(amgut_props, permTest = TRUE, nPerm = 100L, 
#'                           fileLoadAssoPerm = "assoPerm")
#' 
#' all.equal(amgut_comp3$properties, amgut_comp5$properties)
#' 
#' summary(amgut_comp5)
#' }
#'
#' @seealso \code{\link{summary.microNetComp}}, \code{\link{netConstruct}},
#'   \code{\link{netAnalyze}}
#' @references 
#'   \insertRef{benjamini2000adaptive}{NetCoMi} \cr\cr
#'   \insertRef{farcomeni2007some}{NetCoMi} \cr\cr
#'   \insertRef{gill2010statistical}{NetCoMi} \cr\cr
#'   \insertRef{hocevar2016computation}{NetCoMi}\cr\cr
#'   \insertRef{qannari2014significance}{NetCoMi}\cr\cr
#'   \insertRef{real1996probabilistic}{NetCoMi}\cr\cr
#'   \insertRef{yaveroglu2014revealing}{NetCoMi}
#' @import foreach doSNOW filematrix
#' @importFrom stats sd
#' @importFrom WGCNA randIndex
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export

netCompare <- function(x,
                       permTest = FALSE,
                       jaccQuant = 0.75,
                       lnormFit = NULL,
                       testRand = TRUE,
                       nPermRand = 1000L,
                       gcd = TRUE,
                       gcdOrb = c(0, 2, 5, 7, 8, 10, 11, 6, 9, 4, 1),
                       verbose = TRUE,
                       nPerm = 1000L,
                       adjust = "adaptBH",
                       trueNullMethod = "convest",
                       cores = 1L,
                       logFile = NULL,
                       seed = NULL, 
                       fileLoadAssoPerm  = NULL,
                       fileLoadCountsPerm = NULL,
                       storeAssoPerm = FALSE,
                       fileStoreAssoPerm = "assoPerm",
                       storeCountsPerm = FALSE,
                       fileStoreCountsPerm = c("countsPerm1", "countsPerm2"),
                       returnPermProps = FALSE,
                       returnPermCentr = FALSE,
                       assoPerm = NULL, dissPerm = NULL) {
  
  # Check input arguments
  argsIn <- as.list(environment())
  
  if (verbose) message("Checking input arguments ... ", appendLF = FALSE)
  
  argsOut <- .checkArgsNetComp(argsIn)
  
  for (i in 1:length(argsOut)) {
    assign(names(argsOut)[i], argsOut[[i]])
  }
  
  if (verbose) message("Done.")
  
  #-----------------------------------------------------------------------------
  
  # Parameters used for network construction
  parNC <- x$paramsNetConstruct
  
  # Parameters used for network analysis
  parNA <- x$paramsProperties
  
  assoType <- x$input$assoType
  distNet <- ifelse(assoType == "dissimilarity", TRUE, FALSE)
  
  xgroups <- x$input$groups
  matchDesign <- x$input$matchDesign
  sampleSize <- x$input$sampleSize
  twoNets <- x$input$twoNets
  callNetConstr <- x$input$call
  
  if (is.null(lnormFit)) {
    lnormFit <- parNA$lnormFit
  }
  
  count1 <- x$input$countMat1
  count2 <- x$input$countMat2
  countsJoint <- x$input$countsJoint
  
  count_norm1 <- x$input$normCounts1
  count_norm2 <- x$input$normCounts2
  
  assoMat1 <- x$input$assoMat1
  assoMat2 <- x$input$assoMat2
  
  dissMat1 <- x$input$dissMat1
  dissMat2 <- x$input$dissMat2
  
  adja1 <- x$input$adjaMat1
  adja2 <- x$input$adjaMat2
  
  #---------------------------------------------------------------------------
  
  if (!twoNets) {
    stop("Network comparison not possible because a single network
                    has been constructed.")
  }
  
  if (all(adja1[upper.tri(adja1)] == 0)) {
    stop(paste0("Network properties cannot be compared because network of group'",
                xgroups[1], "' is empty."))
  }
  
  if (all(adja2[upper.tri(adja2)] == 0)) {
    stop(paste0("Network properties cannot be compared because network of group'",
                xgroups[2], "' is empty."))
  }
  
  #---------------------------------------------------------------------------
  # Calculate and compare network properties
  
  if (!is.null(seed)) set.seed(seed)
  
  if (permTest & verbose) message("Calculate network properties ... ",
                                  appendLF = FALSE)
  
  props <- .calcDiffProps(adja1 = adja1, adja2 = adja2,
                          dissMat1 = dissMat1, dissMat2 = dissMat2,
                          assoMat1 = assoMat1, assoMat2 = assoMat2,
                          avDissIgnoreInf = parNA$avDissIgnoreInf,
                          sPathNorm = parNA$sPathNorm,
                          sPathAlgo = parNA$sPathAlgo,
                          connectivity = parNA$connectivity,
                          normNatConnect = parNA$normNatConnect,
                          weighted = parNC$weighted,
                          clustMethod = parNA$clustMethod,
                          clustPar = parNA$clustPar,
                          clustPar2 = parNA$clustPar2,
                          weightClustCoef = parNA$weightClustCoef,
                          hubPar = parNA$hubPar,
                          hubQuant = parNA$hubQuant,
                          jaccQuant = jaccQuant,
                          lnormFit = lnormFit,
                          weightDeg = parNA$weightDeg,
                          normDeg = parNA$normDeg,
                          normBetw = parNA$normBetw,
                          normClose = parNA$normClose,
                          normEigen = parNA$normEigen,
                          centrLCC = parNA$centrLCC,
                          testRand = testRand,
                          nPermRand = nPermRand,
                          gcd = gcd, gcdOrb = gcdOrb)
  
  if (permTest & verbose) message("Done.")
  
  callArgs <- argsIn
  callArgs$x <- NULL
  
  output <- list(jaccDeg = props$jaccDeg,
                 jaccBetw = props$jaccBetw,
                 jaccClose = props$jaccClose,
                 jaccEigen = props$jaccEigen,
                 jaccHub = props$jaccHub,
                 randInd = props$randInd,
                 randIndLCC = props$randIndLCC,
                 gcd = props$gcd,
                 gcdLCC = props$gcdLCC,
                 diffGlobal = props$diffsGlobal,
                 diffGlobalLCC = props$diffsGlobalLCC,
                 diffCentr = props$diffsCentr,
                 properties = props$props,
                 propertiesLCC = props$propsLCC,
                 countMatrices = list(count1 = count1,
                                      count2 = count2),
                 assoMatrices = list(assoMat1 = assoMat1,
                                     assoMat2 = assoMat2),
                 dissMatrices = list(dissMat1 = dissMat1,
                                     dissMat2 = dissMat2),
                 adjaMatrices = list(adja1 = adja1,
                                     adja2 = adja2),
                 groups = list(group1 = xgroups[1],
                               group2 = xgroups[2]),
                 paramsProperties = x$paramsProperties,
                 call = match.call(), 
                 callArgs = callArgs)
  
  #-----------------------------------------------------------------------------
  # Generate test statistics for the permuted data
  
  if (permTest) {
    
    if (!is.null(seed)) set.seed(seed)
    
    if (assoType == "dissimilarity") {
      n1 <- ncol(count1)
      n2 <- ncol(count2)
      n <- n1 + n2
      xbind <- cbind(count1, count2)
      
    } else {
      if (parNC$jointPrepro) {
        n1 <- nrow(count_norm1)
        n2 <- nrow(count_norm2)
        
        xbind <- rbind(count_norm1, count_norm2)
        
      } else {
        n1 <- nrow(count1)
        n2 <- nrow(count2)
        
        xbind <- rbind(count1, count2)
      }
      
      n <- n1 + n2
    }
    
    
    nVar = ncol(adja1)
    
    #---------------------------------------------------------------------------
    # Load permutation association matrices
    
    if (!is.null(fileLoadAssoPerm)) {
      
      if (storeAssoPerm && identical(fileLoadAssoPerm, fileStoreAssoPerm)) {
        storeAssoPerm <- FALSE
      }
      
    } else if (!is.null(fileLoadCountsPerm)) {
      
      if (storeCountsPerm && identical(fileLoadCountsPerm, 
                                       fileStoreCountsPerm)) {
        storeCountsPerm <- FALSE
      }
      
    } else {
      perm_group_mat <- .getPermGroupMat(n1 = n1, n2 = n2, n = n,
                                         nPerm = nPerm,
                                         matchDesign = matchDesign)
    }
    
    #-----------------
    # Store permutation count matrices
    
    if (storeCountsPerm) {
      fmat_counts1 <- fm.create(filenamebase = fileStoreCountsPerm[1], 
                                nrow = (n1 * nPerm), ncol = nVar)
      fmat_counts2 <- fm.create(filenamebase = fileStoreCountsPerm[2], 
                                nrow = (n2 * nPerm), ncol = nVar)
      
      if (verbose) {
        message("Files '", 
                paste0(fileStoreCountsPerm[1], ".bmat, "),
                paste0(fileStoreCountsPerm[1], ".desc.txt, \n  "),
                paste0(fileStoreCountsPerm[2], ".bmat, and "),
                paste0(fileStoreCountsPerm[2], ".desc.txt created. "))
      }
      
      close(fmat_counts1)
      close(fmat_counts2)
    }
    
    #-----------------
    # Store permutation association matrices
    
    if (storeAssoPerm) {
      fmat = fm.create(filenamebase = fileStoreAssoPerm, 
                       nrow = (nVar * nPerm), ncol = (2 * nVar))
      
      if (verbose) {
        message("Files '", 
                paste0(fileStoreAssoPerm, ".bmat and "),
                paste0(fileStoreAssoPerm, ".desc.txt created. "))
      }
      
      close(fmat)
    }
    
    #---------------------------------------------------------------------------
    # Run permutation tests
    
    if (verbose) message("Execute permutation tests ... ")
    
    if (!is.null(seed)) {
      seeds <- sample.int(1e8, size = nPerm)
    } else {
      seeds <- NULL
    }
    
    if (cores > 1) {
      if (parallel::detectCores() < cores) cores <- parallel::detectCores()
      
      cl <- parallel::makeCluster(cores, outfile = "")
      
      doSNOW::registerDoSNOW(cl)
      
      '%do_or_dopar%' <- get('%dopar%')
      
    } else {
      '%do_or_dopar%' <- get('%do%')
    }
    
    if (verbose) {
      pb <- utils::txtProgressBar(0, nPerm, style=3)
      
      progress <- function(n) {
        utils::setTxtProgressBar(pb,n)
      }
      
      opts <- list(progress=progress)
    } else {
      opts <- list()
    }
    
    if (!is.null(logFile)) cat("", file=logFile, append=FALSE)
    
    p <- NULL
    
    propsPerm <- foreach(
      p = 1:nPerm,
      .packages = c("filematrix"),
      .export = c(
        ".calcAssociation",
        ".calcDiffProps",
        ".calcJaccard",
        ".calcProps",
        ".getOrbcounts",
        ".normCounts",
        ".scaleDiss",
        ".sparsify",
        ".transToDiss",
        ".transToSim",
        ".transToAdja",
        ".zeroTreat",
        "calcGCD",
        "calcGCM",
        "cclasso",
        "gcoda",
        "multAdjust"
      ),
      .options.snow = opts
      
    ) %do_or_dopar% {
      
      if (!is.null(logFile)) {
        cat(paste("iteration", p, "\n"), file = logFile, append = TRUE)
      }
      
      if (verbose) progress(p)
      
      if (!is.null(seed)) set.seed(seeds[p])
      
      if (!is.null(assoPerm)) {
        # Load permutation association matrices (old version)
        assoMat1.tmp <- assoPerm[[1]][[p]]
        assoMat2.tmp <- assoPerm[[2]][[p]]
        count1.tmp <- count2.tmp <- NULL
        
      } else if (!is.null(dissPerm)) {
        # Load permutation dissimilarity matrices (old version)
        assoMat1.tmp <- dissPerm[[1]][[p]]
        assoMat2.tmp <- dissPerm[[2]][[p]]
        count2.tmp <- count2.tmp <- NULL
        
      } else if (!is.null(fileLoadAssoPerm)) {
        fmat <- fm.open(filenamebase = fileLoadAssoPerm, readonly = TRUE)
        
        # Load permutation asso/diss matrices
        assoMat1.tmp <- fmat[(p - 1) * nVar + (1:nVar), 1:nVar]
        assoMat2.tmp <- fmat[(p - 1) * nVar + (1:nVar), nVar + (1:nVar)]
        
        dimnames(assoMat1.tmp) <- dimnames(assoMat1)
        dimnames(assoMat2.tmp) <- dimnames(assoMat2)
        
        count1.tmp <- count2.tmp <- NULL
        
        close(fmat)
        
      } else {
        if (!is.null(fileLoadCountsPerm)) {
          fmat_counts1 <- fm.open(filenamebase = fileLoadCountsPerm[1],
                                  readonly = TRUE)
          
          fmat_counts2 <- fm.open(filenamebase =  fileLoadCountsPerm[2],
                                  readonly = TRUE)
          
          # load permutation count matrices
          count1.tmp <- fmat_counts1[(p - 1) * n1 + (1:n1), 1:nVar]
          
          count2.tmp <- fmat_counts2[(p - 1) * n2 + (1:n2), 1:nVar]
          
          close(fmat_counts1)
          close(fmat_counts2)
          
        } else {
          # Generate permutation count matrices
          if (assoType == "dissimilarity") {
            count1.tmp <- xbind[, which(perm_group_mat[p,] == 1)]
            count2.tmp <- xbind[, which(perm_group_mat[p,] == 2)]
            
          } else {
            count1.tmp <- xbind[which(perm_group_mat[p,] == 1),]
            count2.tmp <- xbind[which(perm_group_mat[p,] == 2),]
          }
          
          if (!parNC$jointPrepro) {
            # Zero treatment and normalization necessary if two count matrices
            # were used for network construction or
            # dissimilarity network is created
            
            suppressMessages(
              count1.tmp <- .zeroTreat(countMat = count1.tmp,
                                       zeroMethod = parNC$zeroMethod,
                                       zeroParam = parNC$zeroPar,
                                       needfrac = parNC$needfrac,
                                       needint = parNC$needint,
                                       verbose = FALSE)
            )
            
            suppressMessages(
              count2.tmp <- .zeroTreat(countMat = count2.tmp,
                                       zeroMethod = parNC$zeroMethod,
                                       zeroParam = parNC$zeroPar,
                                       needfrac = parNC$needfrac,
                                       needint = parNC$needint,
                                       verbose = FALSE)
            )
            
            suppressMessages(
              count1.tmp <- .normCounts(countMat = count1.tmp,
                                        normMethod = parNC$normMethod,
                                        normParam = parNC$normPar,
                                        zeroMethod = parNC$zeroMethod,
                                        needfrac = parNC$needfrac,
                                        verbose = FALSE)
            )
            
            suppressMessages(
              count2.tmp <- .normCounts(countMat = count2.tmp,
                                        normMethod = parNC$normMethod,
                                        normParam = parNC$normPar,
                                        zeroMethod = parNC$zeroMethod,
                                        needfrac = parNC$needfrac,
                                        verbose = FALSE)
            )
          }
          
          if (storeCountsPerm) {
            # Store permutation count matrices
            fmat_counts1 <- fm.open(filenamebase = fileStoreCountsPerm[1])
            
            fmat_counts2 <- fm.open(filenamebase = fileStoreCountsPerm[2])
            
            fmat_counts1[(p - 1) * n1 + (1:n1),  1:nVar] <- count1.tmp
            fmat_counts2[(p - 1) * n2 + (1:n2),  1:nVar] <- count2.tmp
            
            close(fmat_counts1)
            close(fmat_counts2)
          }
        }
        
        assoMat1.tmp <- .calcAssociation(count1.tmp,
                                         measure = parNC$measure,
                                         measurePar = parNC$measurePar,
                                         verbose = FALSE)
        
        assoMat2.tmp <- .calcAssociation(count2.tmp,
                                         measure = parNC$measure,
                                         measurePar = parNC$measurePar,
                                         verbose = FALSE)
        
        if (storeAssoPerm) {
          fmat <- fm.open(filenamebase = fileStoreAssoPerm)
          
          fmat[(p - 1) * nVar + (1:nVar), 1:nVar] <- assoMat1.tmp
          fmat[(p - 1) * nVar + (1:nVar), nVar + (1:nVar)] <- assoMat2.tmp
          
          close(fmat)
        }
      }
      
      #----------------------------------------------------
      
      if (distNet && parNC$scaleDiss) {
        assoMat1.tmp <- .scaleDiss(assoMat1.tmp)
        assoMat2.tmp <- .scaleDiss(assoMat1.tmp)
      }
      
      #----------------------------------------------------
      # Sparsification and transformation
      
      if (parNC$sparsMethod == "softThreshold") {
        if (length(parNC$softThreshPower) < 2) {
          power1 <- power2 <- parNC$softThreshPower
        } else {
          power1 <- parNC$softThreshPower[1]
          power2 <- parNC$softThreshPower[2]
        }
        
        if (length(parNC$softThreshCut) < 2) {
          stcut1 <- stcut2 <- parNC$softThreshCut
        } else {
          stcut1 <- parNC$softThreshCut[1]
          stcut2 <- parNC$softThreshCut[2]
        }
      }
      
      sparsReslt <- .sparsify(assoMat = assoMat1.tmp,
                              countMat = count1.tmp,
                              sampleSize = n1,
                              measure = parNC$measure,
                              measurePar = parNC$measurePar,
                              assoType = assoType,
                              sparsMethod = parNC$sparsMethod,
                              thresh = parNC$thresh[1],
                              alpha = parNC$alpha[1],
                              adjust = parNC$adjust,
                              lfdrThresh = parNC$lfdrThresh,
                              trueNullMethod = parNC$trueNullMethod,
                              nboot = parNC$nboot,
                              softThreshType = parNC$softThreshType,
                              softThreshPower = power1,
                              softThreshCut = stcut1,
                              kNeighbor = parNC$kNeighbor,
                              knnMutual = parNC$knnMutual,
                              cores = 1L,
                              verbose = FALSE)
      
      if (distNet) {
        dissMat1.tmp <- sparsReslt$assoNew
        assoMat1.tmp <- NULL
        
      } else {
        assoMat1.tmp <- sparsReslt$assoNew
        dissMat1.tmp <- .transToDiss(x = assoMat1.tmp,
                                     dissFunc = parNC$dissFunc,
                                     dissFuncPar = parNC$dissFuncPar)
      }
      
      simMat1.tmp <- .transToSim(x = dissMat1.tmp,
                                 simFunc = parNC$simFunc,
                                 simFuncPar = parNC$simFuncPar)
      
      adja1.tmp <- .transToAdja(x = simMat1.tmp, weighted = parNC$weighted)
      
      # Network 2
      sparsReslt <- .sparsify(assoMat = assoMat2.tmp,
                              countMat = count2.tmp,
                              sampleSize = n2,
                              measure = parNC$measure,
                              measurePar = parNC$measurePar,
                              assoType = assoType,
                              sparsMethod = parNC$sparsMethod,
                              thresh = parNC$thresh[2],
                              alpha = parNC$alpha[2],
                              adjust = parNC$adjust,
                              lfdrThresh = parNC$lfdrThresh,
                              trueNullMethod = parNC$trueNullMethod,
                              nboot = parNC$nboot,
                              softThreshType = parNC$softThreshType,
                              softThreshPower = power2,
                              softThreshCut = stcut2,
                              kNeighbor = parNC$kNeighbor,
                              knnMutual = parNC$knnMutual,
                              cores = 1L,
                              verbose = FALSE)
      
      if (distNet) {
        dissMat2.tmp <- sparsReslt$assoNew
        assoMat2.tmp <- NULL
        
      } else {
        assoMat2.tmp <- sparsReslt$assoNew
        
        dissMat2.tmp <- .transToDiss(x = assoMat2.tmp,
                                     dissFunc = parNC$dissFunc,
                                     dissFuncPar = parNC$dissFuncPar)
      }
      
      simMat2.tmp <- .transToSim(x = dissMat2.tmp,
                                 simFunc = parNC$simFunc,
                                 simFuncPar = parNC$simFuncPar)
      
      adja2.tmp <- .transToAdja(x = simMat2.tmp, weighted = parNC$weighted)
      
      dimnames(adja1.tmp) <- dimnames(adja1)
      dimnames(adja2.tmp) <- dimnames(adja2)
      dimnames(assoMat1.tmp) <- dimnames(assoMat1)
      dimnames(assoMat2.tmp) <- dimnames(assoMat2)
      dimnames(dissMat1.tmp) <- dimnames(dissMat1)
      dimnames(dissMat2.tmp) <- dimnames(dissMat2)
      
      #----------------------------------------------------
      # Compute network properties
      prop.tmp <- .calcDiffProps(adja1 = adja1.tmp,
                                 adja2 = adja2.tmp,
                                 dissMat1 = dissMat1.tmp,
                                 dissMat2 = dissMat2.tmp,
                                 assoMat1 = assoMat1.tmp,
                                 assoMat2 = assoMat2.tmp,
                                 avDissIgnoreInf = parNA$avDissIgnoreInf,
                                 sPathNorm = parNA$sPathNorm,
                                 sPathAlgo = parNA$sPathAlgo,
                                 connectivity = parNA$connectivity,
                                 normNatConnect = parNA$normNatConnect,
                                 weighted = parNC$weighted,
                                 clustMethod = parNA$clustMethod,
                                 clustPar = parNA$clustPar,
                                 clustPar2 = parNA$clustPar2,
                                 weightClustCoef = parNA$weightClustCoef,
                                 hubPar = parNA$hubPar,
                                 hubQuant = parNA$hubQuant,
                                 jaccQuant = jaccQuant,
                                 lnormFit = lnormFit,
                                 weightDeg = parNA$weightDeg,
                                 normDeg = parNA$normDeg,
                                 normBetw = parNA$normBetw,
                                 normClose = parNA$normClose,
                                 normEigen = parNA$normEigen,
                                 centrLCC = parNA$centrLCC,
                                 nPermRand = nPermRand,
                                 testJacc = FALSE,
                                 testRand = FALSE,
                                 gcd = gcd,
                                 gcdOrb = gcdOrb)
      
      out <- list(diffsGlobal = NULL,
                  diffsGlobalLCC = NULL,
                  absDiffsCentr = NULL)
      
      for (i in c("diffavDiss",
                  "diffavPath",
                  "diffDensity",
                  "diffVertConnect",
                  "diffEdgeConnect",
                  "diffNatConnect",
                  "diffPEP",
                  "diffClustCoef",
                  "diffModul")) {
        
        out$diffsGlobal[[i]] <- prop.tmp$diffsGlobal[[i]]
        out$diffsGlobalLCC[[i]] <- prop.tmp$diffsGlobalLCC[[i]]
      }
      
      out$diffsGlobal[["diffnComp"]] <- 
        prop.tmp$diffsGlobal$diffnComp
      
      out$diffsGlobalLCC[["difflccSize"]] <- 
        prop.tmp$diffsGlobalLC$difflccSize
      
      out$diffsGlobalLCC[["difflccSizeRel"]] <- 
        prop.tmp$diffsGlobalLC$difflccSizeRel
      
      out$absDiffsCentr <- prop.tmp$absDiffsCentr
      
      out$props <- prop.tmp$props
      out$propsLCC <- prop.tmp$propsLCC
      
      out$gcd <- prop.tmp$gcd$gcd
      out$gcdLCC <- prop.tmp$gcdLCC$gcd
      
      out
    }
    
    if (verbose) {
      # close progress bar
      close(pb)  
      if (cores > 1) {
        # stop cluster
        message("Stopping socket cluster ... ", appendLF = FALSE)
      }
    }
    
    if (cores > 1) parallel::stopCluster(cl)
    
    if (verbose) {
      message("Done.")
    }
    
    
    if (verbose) {
      message("Calculating p-values ... ", appendLF = FALSE)
    }
    
    #---------------------------------------------------------------------------
    # Create matrices with differences
    
    # Matrices for differences in global properties
    
    globNames <- c("avDiss", "avPath", "Density", "VertConnect", "EdgeConnect",
                   "NatConnect", "PEP", "ClustCoef", "Modul")
    globNamesDiff <- paste0("diff", globNames)
    
    namesWhole <- c("diffnComp", globNamesDiff)
    namesLCC <- c("difflccSize", "difflccSizeRel",
                  globNamesDiff)
    
    permDiffsGlobal <- matrix(0, nrow = nPerm, ncol = length(namesWhole))
    permDiffsGlobal_lcc <- matrix(0, nrow = nPerm, ncol = length(namesLCC))
    
    colnames(permDiffsGlobal) <- namesWhole
    colnames(permDiffsGlobal_lcc) <- namesLCC
    
    #----------------------------------------------
    # Matrices for differences in centralities
    
    absDiffsPermDeg <- absDiffsPermBetw <-
      absDiffsPermClose <- absDiffsPermEigen <-
      permDeg1 <- permDeg2 <- permBetw1 <- permBetw2 <-
      permClose1 <- permClose2 <- permEigen1 <- permEigen2 <-
      matrix(0, nrow = nPerm, ncol = ncol(adja1))
    
    #----------------------------------------------
    for (b in 1:nPerm) {
      # Store differences in global properties
      for (i in 1:9) {
        permDiffsGlobal[b,globNamesDiff[i]] <-
          as.numeric(propsPerm[[b]]$diffsGlobal[globNamesDiff[i]])
        
        permDiffsGlobal_lcc[b,globNamesDiff[i]] <-
          as.numeric(propsPerm[[b]]$diffsGlobalLC[globNamesDiff[i]])
      }
      
      permDiffsGlobal[b, "diffnComp"] <-
        as.numeric(propsPerm[[b]]$diffsGlobal["diffnComp"])
      
      permDiffsGlobal_lcc[b, "difflccSize"] <- 
        as.numeric(propsPerm[[b]]$diffsGlobalLCC["difflccSize"])
      
      permDiffsGlobal_lcc[b, "difflccSizeRel"] <- 
        as.numeric(propsPerm[[b]]$diffsGlobalLCC["difflccSizeRel"])
      
      # Store differences in centralities
      absDiffsPermDeg[b,]   <- propsPerm[[b]]$absDiffsCentr$absDiffDeg
      absDiffsPermBetw[b,]  <- propsPerm[[b]]$absDiffsCentr$absDiffBetw
      absDiffsPermClose[b,] <- propsPerm[[b]]$absDiffsCentr$absDiffClose
      absDiffsPermEigen[b,] <- propsPerm[[b]]$absDiffsCentr$absDiffEigen
      
      permDeg1[b, ]  <- propsPerm[[b]]$props$deg1
      permDeg2[b, ]  <- propsPerm[[b]]$props$deg2
      permBetw1[b, ]  <- propsPerm[[b]]$props$betw1
      permBetw2[b, ]  <- propsPerm[[b]]$props$betw2
      permClose1[b, ] <- propsPerm[[b]]$props$close1
      permClose2[b, ] <- propsPerm[[b]]$props$close2
      permEigen1[b, ] <- propsPerm[[b]]$props$eigen1
      permEigen2[b, ] <- propsPerm[[b]]$props$eigen2
    }
    
    # Save as list
    permDiffsCentr <- list(degree = absDiffsPermDeg,
                           between = absDiffsPermBetw,
                           close = absDiffsPermClose,
                           eigen = absDiffsPermEigen)
    
    permPropsCentr1 <- list(degree = permDeg1, 
                            between = permBetw1,
                            close = permClose1,
                            eigen = permEigen1)
    
    permPropsCentr2 <- list(degree = permDeg2, 
                            between = permBetw2,
                            close = permClose2,
                            eigen = permEigen2)
    
    #---------------------------------------------------------------------------
    # Create matrices for global properties of the permuted data
    
    permPropsGlobal1 <- matrix(0, nrow = nPerm, ncol = 10)
    permPropsGlobal1_lcc <- matrix(0, nrow = nPerm, ncol = 11)
    
    permPropsGlobal2 <- matrix(0, nrow = nPerm, ncol = 10)
    permPropsGlobal2_lcc <- matrix(0, nrow = nPerm, ncol = 11)
    
    propnames <- c("avDiss", "avPath", "Density", "VertConnect", 
                   "EdgeConnect", "NatConnect", "PEP", "ClustCoef", 
                   "Modul")
    propnames2 <- c("avDiss", "avPath", "density", "vertConnect", 
                    "edgeConnect", "natConnect", "pep", "clustCoef", 
                    "modularity")
    
    colnames(permPropsGlobal1) <- c("nComp", propnames)
    colnames(permPropsGlobal2) <- c("nComp", propnames)
    
    colnames(permPropsGlobal1_lcc) <- c("lccSize", "lccSizeRel", propnames)
    colnames(permPropsGlobal2_lcc) <- c("lccSize", "lccSizeRel", propnames)
    
    for (i in 1:9) {
      for (b in 1:nPerm) {
        permPropsGlobal1[b,i] <- 
          as.numeric(propsPerm[[b]]$props[paste0(propnames2[i], 1)])
        
        permPropsGlobal2[b,i] <- 
          as.numeric(propsPerm[[b]]$props[paste0(propnames2[i], 2)])
        
        permPropsGlobal1_lcc[b,i] <- 
          as.numeric(propsPerm[[b]]$propsLCC[paste0(propnames2[i], 1)])
        
        permPropsGlobal2_lcc[b,i] <- 
          as.numeric(propsPerm[[b]]$propsLCC[paste0(propnames2[i], 2)])
      }
    }
    
    
    for (b in 1:nPerm) {
      permPropsGlobal1[b, "nComp"] <- 
        as.numeric(propsPerm[[b]]$props["nComp1"])
      
      permPropsGlobal2[b, "nComp"] <- 
        as.numeric(propsPerm[[b]]$props["nComp2"])
      
      permPropsGlobal1_lcc[b, "lccSize"] <- 
        as.numeric(propsPerm[[b]]$propsLCC["lccSize1"])
      
      permPropsGlobal2_lcc[b, "lccSize"] <- 
        as.numeric(propsPerm[[b]]$propsLCC["lccSize2"])
      
      permPropsGlobal1_lcc[b, "lccSizeRel"] <- 
        as.numeric(propsPerm[[b]]$propsLCC["lccSizeRel1"])
      
      permPropsGlobal2_lcc[b, "lccSizeRel"] <- 
        as.numeric(propsPerm[[b]]$propsLCC["lccSizeRel2"])
    }
    
    #---------------------------------------------------------------------------
    # Store GCD values
    
    if (gcd) {
      gcdPerm <- gcdPermLCC <- numeric(nPerm)
      
      for (b in 1:nPerm) {
        gcdPerm[b] <- propsPerm[[b]]$gcd
        gcdPermLCC[b] <- propsPerm[[b]]$gcdLCC
      }
    }
    
    #---------------------------------------------------------------------------
    # Compute p-values
    
    # Define whether permutation test is exact or approximative
    if (is.null(x$input$matchDesign)) {
      if (nPerm == choose(n, n1)) {
        exact <- TRUE
      } else {
        exact <- FALSE
      }
    } else {
      if (nPerm == .getMaxCombMatch(matchDesign = x$input$matchDesign, n = n)) {
        exact <- TRUE
      } else {
        exact <- FALSE
      }
    }
    
    #------------------------------
    # P-values for the global properties
    
    pvalnames <- paste0("pval", globNames)
    
    pvalDiffGlobal <-  pvalDiffGlobalLCC <- list()
    
    pvalDiffGlobal[["pvalnComp"]] <- 
      .calcPermPval(tstat = props$diffsGlobal[["diffnComp"]],
                    tstatPerm = permDiffsGlobal[, "diffnComp"],
                    nPerm = nPerm, exact = exact)
    
    
    pvalDiffGlobalLCC[["pvallccSize"]] <-
      .calcPermPval(tstat = props$diffsGlobalLCC[["difflccSize"]],
                    tstatPerm = permDiffsGlobal_lcc[, "difflccSize"],
                    nPerm = nPerm, exact = exact)
    
    pvalDiffGlobalLCC[["pvallccSizeRel"]] <-
      .calcPermPval(tstat = props$diffsGlobalLCC[["difflccSizeRel"]],
                    tstatPerm = permDiffsGlobal_lcc[, "difflccSizeRel"],
                    nPerm = nPerm, exact = exact)
    
    for (i in seq_along(pvalnames)) {
      pvalDiffGlobal[[pvalnames[i]]] <-
        .calcPermPval(tstat = props$diffsGlobal[[globNamesDiff[i]]],
                      tstatPerm = permDiffsGlobal[, globNamesDiff[i]],
                      nPerm = nPerm, exact = exact)
      
      pvalDiffGlobalLCC[[pvalnames[i]]] <-
        .calcPermPval(tstat = props$diffsGlobalLCC[[globNamesDiff[i]]],
                      tstatPerm = permDiffsGlobal_lcc[, globNamesDiff[i]],
                      nPerm = nPerm, exact = exact)
    }
    
    #------------------------------
    # P-values for the centralities
    
    pvalDiffDeg <- sapply(1:ncol(adja1), function(i) {
      .calcPermPval(tstat = props$absDiffs$absDiffDeg[i],
                    tstatPerm = absDiffsPermDeg[, i],
                    nPerm = nPerm, exact = exact)})
    
    pvalDiffBetw <- sapply(1:ncol(adja1), function(i) {
      .calcPermPval(tstat = props$absDiffs$absDiffBetw[i],
                    tstatPerm = absDiffsPermBetw[, i],
                    nPerm = nPerm, exact = exact)})
    
    pvalDiffClose <- sapply(1:ncol(adja1), function(i) {
      .calcPermPval(tstat = props$absDiffs$absDiffClose[i],
                    tstatPerm = absDiffsPermClose[, i],
                    nPerm = nPerm, exact = exact)})
    
    pvalDiffEigen <- sapply(1:ncol(adja1), function(i) {
      .calcPermPval(tstat = props$absDiffs$absDiffEigen[i],
                    tstatPerm = absDiffsPermEigen[, i],
                    nPerm = nPerm, exact = exact)})
    
    names(pvalDiffDeg) <- names(pvalDiffBetw) <- 
      names(pvalDiffClose) <- names(pvalDiffEigen) <- colnames(adja1)
    
    #------------------------------
    # P-value for the GCD
    if (gcd) {
      pvalGCD <- .calcPermPval(props$gcd$gcd, gcdPerm,
                               nPerm = nPerm, exact = exact)
      
      pvalGCDLCC <- .calcPermPval(props$gcdLCC$gcd, gcdPermLCC,
                                  nPerm = nPerm, exact = exact)
    }
    #------------------------------
    
    if (verbose) message("Done.")
    
    if (verbose & adjust != "none") {
      message("Adjust for multiple testing using '", adjust, "' ... ",
              appendLF = FALSE)
    }
    
    # Adjust for multiple testing
    
    if (adjust == "adaptBH" && !requireNamespace("limma", quietly = TRUE)) {
      
      message("Installing missing package 'limma' ...")
      
      if (!requireNamespace("BiocManager", quietly = TRUE)) {
        utils::install.packages("BiocManager")
      }
      
      BiocManager::install("limma", dependencies = TRUE)
      message("Done.")
      
      message("Check whether installed package can be loaded ...")
      requireNamespace("limma")
      message("Done.")
    }
    
    pAdjustDiffDeg <- multAdjust(pvals = pvalDiffDeg, adjust = adjust,
                                 trueNullMethod = trueNullMethod,
                                 verbose = FALSE)
    
    pAdjustDiffBetw <- multAdjust(pvals = pvalDiffBetw, adjust = adjust,
                                  trueNullMethod = trueNullMethod,
                                  verbose = FALSE)
    
    pAdjustDiffClose <- multAdjust(pvals = pvalDiffClose, adjust = adjust,
                                   trueNullMethod = trueNullMethod,
                                   verbose = FALSE)
    
    pAdjustDiffEigen <- multAdjust(pvals = pvalDiffEigen, adjust = adjust,
                                   trueNullMethod = trueNullMethod,
                                   verbose = FALSE)
    
    if (verbose & adjust != "none") message("Done.")
    
    output$pvalDiffGlobal <- pvalDiffGlobal 
    
    output$pvalDiffGlobalLCC <- pvalDiffGlobalLCC
    
    output$pvalDiffCentr <- list(pvalDiffDeg = pvalDiffDeg,
                                 pvalDiffBetw = pvalDiffBetw,
                                 pvalDiffClose = pvalDiffClose,
                                 pvalDiffEigen = pvalDiffEigen)
    
    output$pvalDiffCentrAdjust <- list(pAdjustDiffDeg = pAdjustDiffDeg,
                                       pAdjustDiffBetw = pAdjustDiffBetw,
                                       pAdjustDiffClose = pAdjustDiffClose,
                                       pAdjustDiffEigen = pAdjustDiffEigen)
    
    if (returnPermProps) {
      output$permPropsGlobal1 <- permPropsGlobal1
      output$permPropsGlobal2 <- permPropsGlobal2
      output$permPropsGlobalLCC1 <- permPropsGlobal1_lcc
      output$permPropsGlobalLCC2 <- permPropsGlobal2_lcc
      output$permDiffGlobal <- permDiffsGlobal
      output$permDiffGlobalLCC <- permDiffsGlobal_lcc
    }
    
    if (returnPermCentr) {
      output$permPropsCentr1 <- permPropsCentr1
      output$permPropsCentr2 <- permPropsCentr2
      output$permDiffCentr <- permDiffsCentr
    }
    
    if (gcd) {
      output$gcd$pval <- pvalGCD
      output$gcd$gcdPerm <- gcdPerm
      
      output$gcdLCC$pval <- pvalGCDLCC
      output$gcdLCC$gcdPerm <- gcdPermLCC
      
      # Reorder
      output$gcd <- output$gcd[c("gcd", "pval", "gcdPerm", "ocount1",
                                 "ocount2", "gcm1", "gcm2")]
      
      output$gcdLCC <- output$gcdLCC[c("gcd", "pval", "gcdPerm", "ocount1",
                                       "ocount2", "gcm1", "gcm2")]
    }
  }
  
  class(output) <- "microNetComp"
  return(output)
}




