R/netCompare.R
In stefpeschel/NetCoMi: Network Construction and Comparison for Microbiome Data
Defines functions
netCompare
Documented in
netCompare
#' @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)
}
stefpeschel/NetCoMi documentation built on Nov. 12, 2024, 7:12 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions netCompare
Documented in netCompare
#' @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)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.