#' @title Plot Method for microNetProps Objects
#'
#' @description Plotting objects of class \code{microNetProps}.
#'
#' @aliases plot plot.microNetProps
#' @usage \method{plot}{microNetProps}(x,
#'   layout = "spring",
#'   sameLayout = FALSE,
#'   layoutGroup = "union",
#'   repulsion = 1,
#'   groupNames = NULL,
#'   groupsChanged = FALSE,
#'   labels = NULL,
#'   shortenLabels = "none",
#'   labelLength = 6L,
#'   labelPattern = c(5, "'", 3, "'", 3),
#'   charToRm = NULL,
#'   labelScale = TRUE,
#'   labelFont = 1,
#'   labelFile = NULL,
#'
#'   # Nodes:
#'   nodeFilter = "none",
#'   nodeFilterPar = NULL,
#'   rmSingles = "none",
#'   nodeSize = "fix",
#'   normPar = NULL,
#'   nodeSizeSpread = 4,
#'   nodeColor = "cluster",
#'   colorVec = NULL,
#'   featVecCol = NULL,
#'   sameFeatCol = TRUE,
#'   sameClustCol = TRUE,
#'   sameColThresh = 2L,
#'   nodeShape = NULL,
#'   featVecShape = NULL,
#'   nodeTransp = 60,
#'   borderWidth = 1,
#'   borderCol = "gray80",
#'
#'   # Hubs:
#'   highlightHubs = TRUE,
#'   hubTransp = NULL,
#'   hubLabelFont = NULL,
#'   hubBorderWidth = NULL,
#'   hubBorderCol = "black",
#'
#'   # Edges:
#'   edgeFilter = "none",
#'   edgeFilterPar = NULL,
#'   edgeInvisFilter = "none",
#'   edgeInvisPar = NULL,
#'   edgeWidth = 1,
#'   negDiffCol = TRUE,
#'   posCol = NULL,
#'   negCol = NULL,
#'   cut = NULL,
#'   edgeTranspLow = 0,
#'   edgeTranspHigh = 0,
#'
#'   # Additional arguments:
#'   cexNodes = 1,
#'   cexHubs = 1.2,
#'   cexLabels = 1,
#'   cexHubLabels = NULL,
#'   cexTitle = 1.2,
#'   showTitle = NULL,
#'   title1 = NULL,
#'   title2 = NULL,
#'   mar = c(1, 3, 3, 3),
#'   doPlot = TRUE,
#'   ...)
#'
#' @param x  object of class \code{microNetProps}
#' @param layout  indicates the layout used for defining node positions. Can be
#'   a character with one of the layouts provided by
#'   \code{\link[qgraph]{qgraph}}: \code{"spring"} (default), \code{"circle"},
#'   or \code{"groups"}. Alternatively, the layouts provided by igraph (see
#'   \code{\link[igraph:layout_]{layout\_}} are accepted (must be given as
#'   character, e.g. \code{"layout_with_fr"}). Can also be a matrix with row
#'   number equal to the number of nodes and two columns corresponding to the x
#'   and y coordinate.
#' @param sameLayout logical. Indicates whether the same layout should be used
#'   for both networks. Ignored if \code{x} contains only one network. See
#'   argument \code{layoutGroup}.
#' @param layoutGroup  numeric or character. Indicates the group, where the 
#'   layout is taken from if argument \code{sameLayout} is \code{TRUE}. 
#'   The layout is computed for group 1 (and adopted for group 2) if set to "1" 
#'   and it is computed for group 2 if set to "2". Can alternatively be set to
#'   "union" (default) to compute a union of both layouts, where the nodes are 
#'   placed as optimal as possible equally for both networks. 
#' @param repulsion positive numeric value indicating the strength of repulsive
#'   forces in the "spring" layout. Nodes are placed closer together for smaller
#'   values and further apart for higher values. See the \code{repulsion}
#'   argument of \code{\link[qgraph]{qgraph}}.
#' @param groupNames character vector with two entries naming the groups to
#'   which the networks belong. Defaults to the group names returned by 
#'   \code{\link{netConstruct}}: If the data set is split according to a group
#'   variable, its factor levels (in increasing order) are used. Ignored if
#'   arguments \code{title1} and \code{title2} are set or if a single network is
#'   plotted.
#' @param groupsChanged logical. Indicates the order in which the networks are
#'   plotted. If \code{TRUE}, the order is exchanged. See details. Defaults to
#'   \code{FALSE}.
#' @param labels defines the node labels. Can be a named character vector, which 
#'   is used for both groups (then, the adjacency matrices in \code{x} must 
#'   contain the same variables). 
#'   Can also be a list with two named vectors (names must match the row/column 
#'   names of the adjacency matrices). If \code{FALSE}, no labels are plotted. 
#'   Defaults to the row/column names of the adjacency matrices.
#' @param shortenLabels character indicating how to shorten node labels. 
#'   Ignored if node labels are defined via \code{labels}. NetCoMi's function 
#'   \code{\link{editLabels}()} is used for label editing.
#'   Available options are:
#'   \describe{
#'   \item{\code{"intelligent"}}{Elements of \code{charToRm} are removed,
#'   labels are shortened to length \code{labelLength}, and duplicates are
#'   removed using \code{labelPattern}.}
#'   \item{\code{"simple"}}{Elements of \code{charToRm} are  removed and labels
#'   are shortened to length \code{labelLength}.}
#'   \item{\code{"none"}}{Default. Original dimnames of the adjacency matrices 
#'   are used.} }
#' @param labelLength integer defining the length to which labels shall
#'   be shortened if \code{shortenLabels} is set to \code{"simple"} or 
#'   \code{"intelligent"}. Defaults to 6.
#' @param labelPattern vector of three or five elements, which is used if 
#'   argument \code{shortenLabels} is set to \code{"intelligent"}. 
#'   If cutting a label to length \code{labelLength} leads to duplicates, 
#'   the label is shortened according to \code{labelPattern}, 
#'   where the first entry gives the length of the first part, 
#'   the second entry is used a separator, and the third entry
#'   is the length of the third part. If \code{labelPattern} has five elements 
#'   and the shortened labels are still not unique, 
#'   the fourth element serves as further separator, and the fifth element gives
#'   the length of the last label part. Defaults to c(5, "'", 3, "'", 3).
#'   If the data contains, for example, three bacteria "Streptococcus1",
#'   "Streptococcus2" and "Streptomyces", they are by default shortened to 
#'   "Strep'coc'1", "Strep'coc'2", and "Strep'myc".
#' @param charToRm vector with characters to remove from node names. Ignored if
#'   labels are given via \code{labels}.
#' @param labelScale logical. If \code{TRUE}, node labels are scaled according 
#'   to node size
#' @param labelFont integer defining the font of node labels. Defaults to 1.
#' @param labelFile optional character of the form "<file name>.txt" naming a 
#'   file where the original and renamed node labels are stored. The file is 
#'   stored into the current working directory.
#' @param nodeFilter character indicating whether and how nodes should be
#'   filtered. Possible values are:
#'   \describe{
#'   \item{\code{"none"}}{Default. All nodes are plotted.}
#'   \item{\code{"highestConnect"}}{x nodes with highest connectivity (sum of
#'   edge weights) are plotted.}
#'   \item{\code{"highestDegree"}, \code{"highestBetween"}, 
#'   \code{"highestClose"}, \code{"highestEigen"}}{x nodes with highest
#'   degree/betweenness/closeness/eigenvector centrality are plotted.}
#'   \item{\code{"clustTaxon"}}{Only nodes belonging the same cluster as to
#'   variables that are given as character vector via  \code{nodeFilterPar}.}
#'   \item{\code{"clustMin"}}{Plotted are only nodes belonging to clusters with
#'   a minimum number of nodes of x.}
#'   \item{\code{"names"}}{Character vector with variable names to be 
#'   plotted}}\cr
#'   Necessary parameters (e.g. "x") are given via the argument
#'   \code{nodeFilterPar}.
#' @param nodeFilterPar parameters needed for the filtering method defined by
#'   \code{nodeFilter}.
#' @param rmSingles character value indicating how to handle unconnected nodes.
#'   Possible values are \code{"all"} (all single nodes are deleted),
#'   \code{"inboth"} (only nodes that are unconnected in both networks are
#'   removed) or \code{"none"} (default; no nodes are removed). Cannot be set to
#'   \code{"all"}, if the same layout is used for both networks.
#' @param nodeSize  character indicating how node sizes should be determined.
#'   Possible values are: \describe{
#'   \item{\code{"fix"}}{Default. All nodes have same size (hub size can be
#'   defined separately via \code{cexHubs}).}
#'   \item{\code{"degree"}, \code{"betweenness"}, \code{"closeness"},
#'   \code{"eigenvector"}}{Size scaled according to node's centrality}
#'   \item{\code{"counts"}}{Size scaled according to the sum of counts (of
#'   microbes or samples, depending on what nodes express).}
#'   \item{\code{"normCounts"}}{Size scaled according to the sum of normalized
#'   counts (of microbes or samples), which are exported by
#'   \code{netConstruct}.}
#'   \item{\code{"TSS", "fractions", "CSS", "COM", "rarefy", "VST", "clr", 
#'   "mclr"}}{Size scaled according to the sum of normalized
#'   counts. Available are the same options as for \code{normMethod} in
#'   \code{\link{netConstruct}}. Parameters are set via \code{normPar}.}
#'   }
#' @param normPar list with parameters passed to the function for normalization
#'   if \code{nodeSize} is set to a normalization method. Used analogously to
#'   \code{normPar} of \code{\link{netConstruct}()}.
#' @param nodeSizeSpread positive numeric value indicating the spread of node
#'   sizes. The smaller the value, the more similar are the node sizes. Node 
#'   sizes are calculated by: (x - min(x)) / (max(x) - min(x)) * nodeSizeSpread 
#'   + cexNodes. For \code{nodeSizeSpread = 4} (default) and 
#'   \code{cexNodes = 1}, node sizes range from 1 to 5.
#' @param nodeColor a character specifying the node colors. Possible values are
#'   \code{"cluster"} (colors according to determined clusters),
#'   \code{"feature"} (colors according to node's features defined by
#'   \code{featVecCol}), \code{"colorVec"} (the vector \code{colorVec}). For the
#'   former two cases, the colors can be specified via \code{colorVec}. If
#'   \code{colorVec} is not defined, the \code{rainbow} function from
#'   \code{grDevices} package  is used. Also accepted is a character value
#'   defining a color, which is used for all nodes. If \code{NULL}, "grey40" is
#'   used for all nodes.
#' @param colorVec a vector or list with two vectors used to specify node 
#'   colors. Different usage depending on the "nodeColor" argument:
#'   \describe{
#'   \item{\code{nodeColor = "cluster"}}{\code{colorVec} must be a vector. 
#'   Depending on the \code{sameClustCol} argument, the colors are used only in 
#'   one or both networks. If the vector is not long enough, a warning is 
#'   returned and colors from \code{rainbow()} are used for the remaining 
#'   clusters.}
#'   \item{\code{nodeColor = "feature"}}{Defines a color for each level of 
#'   \code{featVecCol}. Can be a list with two vectors used for the two networks 
#'   (for a single network, only the first element is used) or a vector, which 
#'   is used for both groups if two networks are plotted.}
#'   \item{\code{nodeColor = "colorVec"}}{\code{colorVec} defines a color for 
#'   each node implying that its names must match the node's names (which is 
#'   also ensured if names match the colnames of the original count matrix). 
#'   Can be a list with two vectors used for the two networks (for a single 
#'   network, only the first element is used) or a vector, which is used for 
#'   both groups if two networks are plotted.}}
#' @param featVecCol  a vector with a feature for each node. Used for coloring
#'   nodes if \code{nodeColor} is set to \code{"feature"}. Is coerced to a
#'   factor. If \code{colorVec} is given, its length must be larger than or
#'   equal to the number of feature levels.
#' @param sameFeatCol logical indicating whether the same color should be used 
#'   for same features in both networks (only used if two networks are plotted, 
#'   \code{nodeColor = "feature"}, and no color vector/list is given (via 
#'   \code{featVecCol})).
#' @param sameClustCol if TRUE (default) and two networks are plotted, clusters
#'   having at least \code{sameColThresh} nodes in common have the same color.
#'   Only used if \code{nodeColor} is set to \code{"cluster"}.
#' @param sameColThresh indicates how many nodes a cluster must have in common
#'   in the two groups to have the same color. See argument \code{sameClustCol}.
#'   Defaults to 2.
#' @param nodeShape character vector specifying node shapes. Possible values
#'   are \code{"circle"} (default), \code{"square"}, \code{"triangle"}, and
#'   \code{"diamond"}. If \code{featVecShape} is not \code{NULL}, the length of
#'   \code{nodeShape} must equal the number of factor levels given by
#'   \code{featVecShape}. Then, each shape is assigned to one factor level (in
#'   increasing order). If \code{featVecShape} is \code{NULL}, the first shape
#'   is used for all nodes. See the example.
#' @param featVecShape a vector with a feature for each node. If not \code{NULL},
#'   a different node shape is used for each feature. Is coerced to factor mode.
#'   The maximum number of factor levels is 4 corresponding to the four possible
#'   shapes defined via \code{nodeShape}.
#' @param nodeTransp an integer between 0 and 100 indicating the transparency of
#'   node colors. 0 means no transparency, 100 means full transparency. Defaults
#'   to 60.
#' @param borderWidth numeric specifying the width of node borders. Defaults to
#'   1.
#' @param borderCol character specifying the color of node borders. Defaults to
#'   "gray80"
#' @param highlightHubs logical indicating if hubs should be highlighted. If
#'   \code{TRUE}, the following features can be defined separately for hubs:
#'   transparency (by \code{hubTransp}), label font (by \code{hubLabelFont}),
#'   border width (by \code{hubBorderWidth}), and border color (by
#'   \code{hubBorderCol}).
#' @param hubTransp numeric between 0 and 100 specifying the color transparency
#'   of hub nodes. See argument \code{nodeTransp}. Defaults to
#'   \code{0.5*nodeTransp}. Ignored if \code{highlightHubs} is \code{FALSE}.
#' @param hubLabelFont integer specifying the label font of hub nodes. Defaults
#'   to \code{2*labelFont}. Ignored if \code{highlightHubs} is \code{FALSE}.
#' @param hubBorderWidth numeric specifying the border width of hub nodes.
#'   Defaults to \code{2*borderWidth}. Ignored if \code{highlightHubs} is
#'   \code{FALSE}.
#' @param hubBorderCol character specifying the border color of hub nodes.
#'   Defaults to \code{"black"}. Ignored if \code{highlightHubs} is 
#'   \code{FALSE}.
#' @param edgeFilter character specifying how edges are filtered. 
#'   Possible values are:
#'   \describe{
#'   \item{\code{"none"}}{Default. All edges are plotted.}
#'   \item{\code{"threshold"}}{For association networks, only edges 
#'   corresponding to an absolute association >= x are plotted. For 
#'   dissimilarity networks, only edges corresponding to a dissimilarity <= x 
#'   are plotted. The behavior is similar to that of sparsification via threshold 
#'   in netConstruct().}
#'   \item{\code{"highestWeight"}}{The first x edges with highest edge weight 
#'   are plotted.} 
#'   }
#'   x is defined by \code{edgeFilterPar}, respectively.
#' @param edgeFilterPar numeric specifying the "x" in \code{edgeFilter}.
#' @param edgeInvisFilter similar to \code{edgeFilter} but the edges are removed
#'   only after computing the layout so that edge removal does not influence the
#'   layout. Defaults to \code{"none"}.
#' @param edgeInvisPar numeric specifying the "x" in \code{edgeInvisFilter}.
#' @param edgeWidth numeric specifying the edge width. See argument
#'   \code{"edge.width"} of \code{\link[qgraph]{qgraph}}.
#' @param negDiffCol logical indicating if edges with a negative corresponding
#'   association should be colored different. If \code{TRUE} (default), argument
#'   \code{posCol} is used for edges with positive association and \code{negCol}
#'   for those with negative association. If \code{FALSE} and for dissimilarity
#'   networks, only \code{posCol} is used.
#' @param posCol vector (character or numeric) with one or two elements
#'   specifying the color of edges with positive weight and also for edges with
#'   negative weight if \code{negDiffCol} is set to \code{FALSE}. The first
#'   element is used for edges with weight below \code{cut} and the second for
#'   edges with weight above \code{cut}. If a single value is given, it is used
#'   for both cases. Defaults to \code{c("#009900", "darkgreen")}.
#' @param negCol vector (character or numeric) with one or two elements
#'   specifying the color of edges with negative weight. The first
#'   element is used for edges with absolute weight below \code{cut} and the
#'   second for edges with absolute weight above \code{cut}. If a single value
#'   is given, it is used for both cases. Ignored if \code{negDiffCol} is
#'   \code{FALSE}. Defaults to \code{c("red", "#BF0000")}.
#' @param cut defines the \code{"cut"} parameter of 
#'   \code{\link[qgraph]{qgraph}}. Can
#'   be either a numeric value (is used for both groups if two networks are
#'   plotted) or vector of length two. The default is set analogous to that in
#'   \code{\link[qgraph]{qgraph}}: "0 for graphs with less then 20 nodes. For
#'   larger graphs the cut value is automatically chosen to be equal to the
#'   maximum of the 75th quantile of absolute edge strengths or the edge
#'   strength corresponding to 2n-th edge strength (n being the number of
#'   nodes.)" If two networks are plotted, the mean of the two determined cut 
#'   parameters is used so that edge thicknesses are comparable.
#' @param edgeTranspLow numeric value between 0 and 100 specifying the
#'   transparency of edges with weight below \code{cut}. The higher this value,
#'   the higher the transparency.
#' @param edgeTranspHigh analogous to \code{edgeTranspLow}, but used for edges
#'   with weight ABOVE \code{cut}.
#' @param cexNodes numeric scaling node sizes. Defaults to 1.
#' @param cexHubs numeric scaling hub sizes. Only used if \code{nodeSize} is set
#'   to \code{"hubs"}.
#' @param cexLabels numeric scaling node labels. Defaults to 1. If set to 0, no
#'   node labels are plotted. 
#' @param cexHubLabels numeric scaling the node labels of hub nodes. Equals 
#'   \code{cexLabels} by default. Ignored, if \code{highlightHubs = FALSE}.
#' @param cexTitle numeric scaling title(s). Defaults to 1.2.
#' @param showTitle if \code{TRUE}, a title is shown for each network, which is
#'   either defined via \code{groupNames}, or \code{title1} and \code{title2}.
#'   Defaults to \code{TRUE} if two networks are plotted and \code{FALSE} for
#'   a single network.
#' @param title1 character giving a title for the first network.
#' @param title2 character giving a title for the second network (if existing).
#' @param mar a numeric vector of the form c(bottom, left, top, right) defining
#'   the plot margins. Works similar to the \code{mar} argument in
#'   \code{\link[graphics]{par}}. Defaults to c(1,3,3,3).
#' @param doPlot logical. If \code{FALSE}, the network plot is suppressed. 
#'   Useful for saving the output (e.g., the layout) without plotting.
#' @param ... further arguments being passed to \code{\link[qgraph]{qgraph}}, 
#'   which is used for network plotting.
#'
#' @return Returns (invisibly) a list with the following elements: \tabular{ll}{
#'   \code{q1,q2}\tab the qgraph object(s)\cr
#'   \code{layout}\tab layout(s) specifying node positions\cr
#'   \code{nodecolor}\tab one or two vectors with node colors (one for each 
#'   group) \cr
#'   \code{labels}\tab one or two vectors with node labels}
#'
#' @examples
#' # Load data sets from American Gut Project (from SpiecEasi package)
#' data("amgut1.filt")
#' 
#' # Network construction
#' amgut_net <- netConstruct(amgut1.filt, measure = "pearson",
#'                           filtTax = "highestVar",
#'                           filtTaxPar = list(highestVar = 50),
#'                           zeroMethod = "pseudoZO", normMethod = "clr",
#'                           sparsMethod = "threshold", thresh = 0.3)
#' 
#' # Network analysis
#' amgut_props <- netAnalyze(amgut_net)
#' 
#' ### Network plots ###
#' # Clusters are used for node coloring: 
#' plot(amgut_props, 
#'      nodeColor = "cluster")
#' 
#' # Remove singletons
#' plot(amgut_props, 
#'      nodeColor = "cluster", 
#'      rmSingles = TRUE)
#' 
#' # A higher repulsion places nodes with high edge weight closer together
#' plot(amgut_props, 
#'      nodeColor = "cluster", 
#'      rmSingles = TRUE, 
#'      repulsion = 1.2)
#' 
#' # A feature vector is used for node coloring
#' # (this could be a vector with phylum names of the ASVs)
#' set.seed(123456)
#' featVec <- sample(1:5, nrow(amgut1.filt), replace = TRUE)
#' 
#' # Names must be equal to ASV names
#' names(featVec) <- colnames(amgut1.filt)
#' 
#' plot(amgut_props, 
#'      rmSingles = TRUE, 
#'      nodeColor = "feature", 
#'      featVecCol = featVec,
#'      colorVec = heat.colors(5))
#' 
#' # Use a further feature vector for node shapes
#' shapeVec <- sample(1:3, ncol(amgut1.filt), replace = TRUE)
#' names(shapeVec) <- colnames(amgut1.filt)
#' 
#' plot(amgut_props,
#'      rmSingles = TRUE,
#'      nodeColor = "feature",
#'      featVecCol = featVec,
#'      colorVec = heat.colors(5), 
#'      nodeShape = c("circle", "square", "diamond"),
#'      featVecShape = shapeVec, 
#'      highlightHubs = FALSE)
#'
#' @importFrom  qgraph qgraph
#' @importFrom graphics par title
#' @importFrom grDevices rainbow
#' @importFrom stats quantile
#' @importFrom WGCNA labels2colors
#' @importFrom utils write.table
#' @method plot microNetProps
#' @export

plot.microNetProps <- function(x,
                               layout = "spring",
                               sameLayout = FALSE,
                               layoutGroup = "union",
                               repulsion = 1,
                               groupNames = NULL,
                               groupsChanged = FALSE,
                               labels = NULL,
                               shortenLabels = "none",
                               labelLength = 6L,
                               labelPattern = c(5, "'", 3, "'", 3),
                               charToRm = NULL,
                               labelScale = TRUE,
                               labelFont = 1,
                               labelFile = NULL,
                               #nodes
                               nodeFilter = "none",
                               nodeFilterPar = NULL,
                               rmSingles = "none",
                               nodeSize = "fix",
                               normPar = NULL,
                               nodeSizeSpread = 4,
                               nodeColor = "cluster",
                               colorVec = NULL,
                               featVecCol = NULL,
                               sameFeatCol = TRUE,
                               sameClustCol = TRUE,
                               sameColThresh = 2L,
                               nodeShape = NULL,
                               featVecShape = NULL,
                               nodeTransp = 60,
                               borderWidth = 1,
                               borderCol = "gray80",
                               #hubs
                               highlightHubs = TRUE,
                               hubTransp = NULL,
                               hubLabelFont = NULL,
                               hubBorderWidth = NULL,
                               hubBorderCol = "black",
                               #edges
                               edgeFilter = "none",
                               edgeFilterPar = NULL,
                               edgeInvisFilter = "none",
                               edgeInvisPar = NULL,
                               edgeWidth = 1,
                               negDiffCol = TRUE,
                               posCol = NULL,
                               negCol = NULL,
                               cut = NULL,
                               edgeTranspLow = 0,
                               edgeTranspHigh = 0,
                               cexNodes = 1,
                               cexHubs = 1.2,
                               cexLabels = 1,
                               cexHubLabels = NULL,
                               cexTitle = 1.2,
                               showTitle = NULL,
                               title1 = NULL,
                               title2 = NULL,
                               mar = c(1, 3, 3, 3),
                               doPlot = TRUE,
                               ...) {

  argsIn <- as.list(environment())
  
  argsOut <- .checkArgsPlotMNP(argsIn)

  for (i in 1:length(argsOut)) {
    assign(names(argsOut)[i], argsOut[[i]])
  }
  
  xgroups <- x$input$groups
  isempty1 <- x$isempty$isempty1
  isempty2 <- x$isempty$isempty2
  
  # if twoNets==TRUE, two networks are plotted
  twoNets <- x$input$twoNets
  
  if (twoNets & !is.null(xgroups)) {
    deletespace <- function(name) {
      if (grepl(" ", strsplit(name, "")[[1]])[1]) {
        gname <- strsplit(name, "")[[1]]
        name <- paste(gname[-1], collapse = "")
      }
      name
    }
    xgroups[1] <- deletespace(xgroups[1])
    xgroups[2] <- deletespace(xgroups[2])
  }
  
  opar <- par()
  
  adja1_orig <- x$input$adjaMat1
  
  #=============================================================================
  # Edge and node filtering
  
  dissMat1 <- if (is.null(x$input$dissEst1)) {
    NULL
  } else {
    x$input$dissMat1
  }
  
  adja1 <- .filterEdges(adja = adja1_orig, 
                        assoEst = x$input$assoMat1,
                        dissEst = dissMat1,
                        edgeFilter = edgeFilter,
                        edgeFilterPar = edgeFilterPar)
  
  keep1 <- .filterNodes(adja = adja1, nodeFilter = nodeFilter,
                        nodeFilterPar = nodeFilterPar, layout = layout,
                        degree = x$centralities$degree1,
                        between = x$centralities$between1,
                        close = x$centralities$close1,
                        eigen = x$centralities$eigenv1,
                        cluster = x$clustering$clust1)
  
  if (twoNets) {
    
    adja2_orig <- x$input$adjaMat2
    
    dissMat2 <- if (is.null(x$input$dissEst2)) {
      NULL
    } else {
      x$input$dissMat2
    }
    
    adja2 <- .filterEdges(adja = adja2_orig, 
                          assoEst = x$input$assoMat2,
                          dissEst = dissMat2,
                          edgeFilter = edgeFilter,
                          edgeFilterPar = edgeFilterPar)
    
    keep2 <- .filterNodes(adja = adja2_orig, nodeFilter = nodeFilter,
                          nodeFilterPar = nodeFilterPar, layout = layout,
                          degree = x$centralities$degree2,
                          between = x$centralities$between2,
                          close = x$centralities$close2,
                          eigen = x$centralities$eigenv2,
                          cluster = x$clustering$clust2)
    
    if (length(keep1) == 0 & length(keep2) == 0) {
      stop("No nodes remaining in both networks after node filtering.")
    }
    
    if (sameLayout) {
      keep.tmp <- union(keep1, keep2)
      keep <- colnames(adja1)[colnames(adja1) %in% keep.tmp]
      adja1 <- adja1[keep, keep]
      adja2 <- adja2[keep, keep]
    } else {
      adja1 <- adja1[keep1, keep1]
      adja2 <- adja2[keep2, keep2]
    }
    
  } else {
    if (length(keep1) == 0) {
      stop("No nodes remaining after node filtering.")
    }
    adja1 <- adja1[keep1, keep1]
    adja2 <- NULL
    adja2_orig <- NULL
  }
  
  
  if (rmSingles != "none") {
    adja1.tmp <- adja1
    
    diag(adja1.tmp) <- 0
    
    zeros1 <- sapply(1:nrow(adja1.tmp), function(i) { 
      all(adja1.tmp[i, ] == 0) })
    
    names(zeros1) <- colnames(adja1.tmp)
    
    if (any(zeros1 == TRUE)) {
      torm1 <- which(zeros1 == TRUE)
    } else {
      torm1 <- NULL
    }
    
    if (twoNets) {
      adja2.tmp <- adja2
      
      diag(adja2.tmp) <- 0
      
      zeros2 <- sapply(1:nrow(adja2.tmp), function(i) { 
        all(adja2.tmp[i, ] == 0) })
      
      names(zeros2) <- colnames(adja2.tmp)
      
      if (any(zeros1 == TRUE)) {
        torm2 <- which(zeros2 == TRUE)
      } else {
        torm2 <- NULL
      }
      
      torm_both <- intersect(torm1, torm2)
      
      if (rmSingles == "inboth" & length(torm_both) != 0) {
        adja1 <- adja1.tmp[-torm_both, -torm_both]
        adja2 <- adja2.tmp[-torm_both, -torm_both]
        
      } else if (rmSingles == "all") {
        if (sameLayout) {
          if (length(torm_both) != 0) {
            warning('Argument "rmSingles" set to "inboth"', 
                    ' due to the same layout.')
            adja1 <- adja1.tmp[-torm_both, -torm_both]
            adja2 <- adja2.tmp[-torm_both, -torm_both]
          }
        } else {
          if (length(torm1) != 0) adja1 <- adja1.tmp[-torm1, -torm1]
          if (length(torm2) != 0) adja2 <- adja2.tmp[-torm2, -torm2]
          
        }
      }
      
      isempty1 <- all(adja1[upper.tri(adja1)] == 0)
      isempty2 <- all(adja2[upper.tri(adja2)] == 0)
      
    } else {
      if (length(torm1) != 0) adja1 <- adja1.tmp[-torm1, -torm1]
      isempty1 <- all(adja1[upper.tri(adja1)] == 0)
      isempty2 <- NULL
      if (isempty1) stop("Network is empty and cannot be plotted.")
    }
  }
  
  kept1 <- which(colnames(adja1_orig) %in% colnames(adja1))
  kept2 <- which(colnames(adja2_orig) %in% colnames(adja2))
  
  #==========================================================================
  # Node colors
  
  if (nodeColor == "cluster") {
    
    if (!is.null(colorVec)) {
      
      if (is.list(colorVec)) {
        stop("'colorVec' must be a vector if clusters are used for node colors.
             Consider setting 'sameColThresh' to a high value.")
      }
      
      stopifnot(is.vector(colorVec))
    }
    
    clust1 <- x$clustering$clust1
    clust2 <- x$clustering$clust2
    
    if (!(is.null(clust1) & is.null(clust2))) {
      clustcolors <- .getClustCols(clust1 = clust1, clust2 = clust2,
                                    adja1 = adja1, adja2 = adja2,
                                    kept1 = kept1, kept2 = kept2,
                                    isempty1 = isempty1, isempty2 = isempty2,
                                    colorVec = colorVec,
                                    sameClustCol = sameClustCol,
                                    sameColThresh = sameColThresh,
                                    twoNets = twoNets)
      nodecol1 <- clustcolors$clustcol1
      nodecol2 <- clustcolors$clustcol2
      
    } else {
      message('No clusterings determined.')
      nodecol1 <- rep("grey40", ncol(adja1))
      if (twoNets) {
        nodecol2 <- rep("grey40", ncol(adja2))
      } else {
        nodecol2 <- NULL
      }
    }
    
  } else if (nodeColor == "feature") {
    featVecCol <- as.factor(featVecCol)
    stopifnot(all(colnames(adja1) %in% names(featVecCol)))
    
    if (is.null(colorVec)) {
      if (sameFeatCol) {
        colorVec1 <- colorVec2 <- grDevices::rainbow(length(levels(featVecCol)))
      } else {
        cols <- grDevices::rainbow(2 * length(levels(featVecCol)))
        colorVec1 <- cols[seq.int(1L, length(cols), 2L)]
        colorVec2 <- cols[seq.int(2L, length(cols), 2L)]
      }
      
    } else {
      if (is.list(colorVec)) {
        stopifnot(length(colorVec) == 2)
        colorVec1 <- colorVec[[1]]
        colorVec2 <- colorVec[[2]]
        
        if (length(colorVec1) < length(levels(featVecCol)) || 
            length(colorVec2) < length(levels(featVecCol))) {
          stop("Length of color vector(s) (argument 'colorVec') shorter than
               number of levels of 'featVecCol'.")
        }
        
      } else {
        if (length(colorVec) < length(levels(featVecCol))) {
          stop(paste("Argument 'featVecCol' has", length(levels(featVecCol)),
                     "levels but 'colorVec' has only length ", 
                     length(colorVec)))
        }
        
        colorVec1 <- colorVec2 <- colorVec
      }
    }
    
    feature1 <- featVecCol[colnames(adja1)]
    nodecol1 <- character(length(feature1))
    
    for (i in seq_along(levels(feature1))) {
      nodecol1[feature1 == levels(feature1)[i]] <- colorVec1[i]
    }
    
    nodecol2 <- NULL
    
    if (twoNets) {
      stopifnot(all(colnames(adja2) %in% names(featVecCol)))
      
      feature2 <- featVecCol[colnames(adja2)]
      nodecol2 <- character(length(feature2))
      
      for (i in seq_along(levels(feature2))) {
        nodecol2[feature2 == levels(feature2)[i]] <- colorVec2[i]
      }
    }
    
  } else if (nodeColor == "colorVec") {
    
    if (is.list(colorVec)) {
      stopifnot(length(colorVec) == 2)
      colorVec1 <- colorVec[[1]]
      colorVec2 <- colorVec[[2]]
      
    } else {
      colorVec1 <- colorVec2 <- colorVec
    }
    
    stopifnot(all(colnames(adja1) %in% names(colorVec1)))
    nodecol1 <- colorVec1[colnames(adja1)]
    
    nodecol2 <- NULL
    
    if (twoNets) {
      stopifnot(all(colnames(adja2) %in% names(colorVec2)))
      nodecol2 <- colorVec2[colnames(adja2)]
    }
    
  } else if (is.character(nodeColor)) {
    nodecol1 <- rep(nodeColor, ncol(adja1))
    nodecol2 <- NULL
    if (twoNets) {
      nodecol2 <- rep(nodeColor, ncol(adja2))
    }
    
  } else {
    nodecol1 <- rep("grey40", ncol(adja1))
    if (twoNets) {
      nodecol2 <- rep("grey40", ncol(adja2))
    }
    
  }
  
  if (nodeTransp > 0) {
    nodecol1 <- colToTransp(nodecol1, nodeTransp)
    if (twoNets) {
      nodecol2 <- colToTransp(nodecol2, nodeTransp)
    }
  }
  
  #==========================================================================
  # Node shapes
  
  if (is.null(featVecShape)) {
    if (is.null(nodeShape)) {
      nodeShape1 <- nodeShape2 <- "circle"
    } else {
      nodeShape2 <- nodeShape2 <- nodeShape[1]
    }
    
  } else {
    stopifnot(all(colnames(adja1) %in% names(featVecShape)))
    
    shape_fact <- as.factor(featVecShape)
    nlevel <- length(levels(shape_fact))
    
    if (nlevel > 4) {
      stop("Argument 'featVecShape' may at most have four factor levels.")
    }
    
    if (is.null(nodeShape)) {
      nodeShape <- c("circle", "square", "triangle", "diamond")
      
    } else {
      if (length(nodeShape) < nlevel) {
        stop(paste("Argument 'nodeShape' must have length", nlevel,
                   "because 'featVecShape' has", nlevel, "levels."))
      }
    }
    
    shapeVec <- as.character(featVecShape)
    names(shapeVec) <- names(featVecShape)
    
    for (i in 1:nlevel) {
      shapeVec[featVecShape == levels(shape_fact)[i]] <- nodeShape[i]
    }
    
    nodeShape1 <- shapeVec[colnames(adja1)]
    
    if (twoNets) {
      stopifnot(all(colnames(adja2) %in% names(featVecShape)))
      nodeShape2 <- shapeVec[colnames(adja2)]
    }
  }
  
  #==========================================================================
  # Node sizes
  
  if (!is.numeric(nodeSize)) {
    if (is.null(x$input$countMat1)) {
      counts.tmp <- x$input$countsJoint
    } else {
      counts.tmp <- x$input$countMat1
    }
    nodeSize1 <- .getNodeSize(nodeSize = nodeSize, normPar = normPar,
                               nodeSizeSpread = nodeSizeSpread,
                               adja = adja1, countMat = counts.tmp,
                               normCounts = x$input$normCounts1,
                               assoType = x$input$assoType, kept = kept1,
                               cexNodes = cexNodes, cexHubs = cexHubs,
                               hubs = x$hubs$hubs1, 
                               highlightHubs = highlightHubs,
                               degree = x$centralities$degree1,
                               between = x$centralities$between1,
                               close = x$centralities$close1,
                               eigen = x$centralities$eigenv1)
    if (twoNets) {
      if (is.null(x$input$countMat2)) {
        counts.tmp <- x$input$countsJoint
      } else {
        counts.tmp <- x$input$countMat2
      }
      nodeSize2 <- .getNodeSize(nodeSize = nodeSize, normPar = normPar,
                                 nodeSizeSpread = nodeSizeSpread,
                                 adja = adja2, countMat = counts.tmp,
                                 normCounts = x$input$normCounts2, 
                                 assoType = x$input$assoType, kept = kept2,
                                 cexNodes = cexNodes, cexHubs = cexHubs,
                                 hubs = x$hubs$hubs2, 
                                 highlightHubs = highlightHubs,
                                 degree = x$centralities$degree2,
                                 between = x$centralities$between2,
                                 close = x$centralities$close2,
                                 eigen = x$centralities$eigenv2)
    }
    rm(counts.tmp)
  }
  
  #===============================================
  # Border colors and fonts
  
  border1 <- rep(borderCol, nrow(adja1))
  labelFont1 <- rep(labelFont, nrow(adja1))
  borderWidth1 <- rep(borderWidth, nrow(adja1))
  
  if (highlightHubs) {
    if (is.null(hubLabelFont)) {
      hubLabelFont <- labelFont * 2
    }
    if (is.null(hubBorderWidth)) {
      hubBorderWidth <- borderWidth * 2
    }
    if (is.null(hubTransp)) {
      hubTransp <- nodeTransp / 2
    }
    
    hubs1 <- x$hubs$hubs1
    hubs1 <- hubs1[hubs1 %in% colnames(adja1)]
    
    border1[match(hubs1, rownames(adja1))] <- hubBorderCol
    labelFont1[match(hubs1, rownames(adja1))] <- hubLabelFont
    borderWidth1[match(hubs1, rownames(adja1))] <- hubBorderWidth
    
    if (nodeTransp != hubTransp) {
      hubcol1 <- colToTransp(nodecol1, hubTransp)
      nodecol1[rownames(adja1) %in% hubs1] <- 
        hubcol1[rownames(adja1) %in% hubs1]
    }
    
  }
  
  if (twoNets) {
    border2 <- rep(borderCol, nrow(adja2))
    labelFont2 <- rep(labelFont, nrow(adja2))
    borderWidth2 <- rep(borderWidth, nrow(adja2))
    
    if (highlightHubs) {
      hubs2 <- x$hubs$hubs2
      hubs2 <- hubs2[hubs2 %in% colnames(adja2)]
      
      border2[match(hubs2, rownames(adja2))] <- hubBorderCol
      labelFont2[match(hubs2, rownames(adja2))] <- hubLabelFont
      borderWidth2[match(hubs2, rownames(adja2))] <- hubBorderWidth
      
      if (nodeTransp != hubTransp) {
        hubcol2 <- colToTransp(nodecol2, hubTransp)
        nodecol2[rownames(adja2) %in% hubs2] <- 
          hubcol2[rownames(adja2) %in% hubs2]
      }
    }
  }
  
  #==========================================================================
  # Title
  
  if (showTitle) {
    if (twoNets) {
      if (is.null(title1) || is.null(title2)) {
        if (is.null(groupNames)) {
          main1 = paste0("group '" , xgroups[1], "'" )
          main2 = paste0("group '" , xgroups[2], "'" )
        } else {
          stopifnot(length(groupNames) == 2)
          
          main1 = as.character(groupNames[1])
          main2 = as.character(groupNames[2])
        }
      } else {
        main1 <- title1
        main2 <- title2
      }
    } else {
      if (is.null(title1)) {
        main1 <- ""
        warning('Argument "showTitle" is TRUE, but no title was defined.')
      } else {
        main1 <- title1
      }
    }
    
  }
  
  #==========================================================================
  # Define cut parameter for qgraph()
  
  if (!is.null(cut) & (length(cut) == 1)) {
    stopifnot(is.numeric(cut))
    cut1 <- cut2 <- cut
  } else if (!is.null(cut) & (length(cut) == 2)) {
    stopifnot(is.numeric(cut))
    cut1 <- cut[1]
    cut2 <- cut[2]
  } else {
    weights1 <- adja1[upper.tri(adja1)]
    weights1 <- weights1[abs(weights1) > 0]
    nNodes1 <- ncol(adja1)
    
    if (twoNets) {
      weights2 <- adja2[upper.tri(adja2)]
      weights2 <- weights2[abs(weights2) > 0]
      nNodes2 <- ncol(adja2)
    } else {
      weights2 <- 0
      nNodes2 <- 1
    }
    
    if (length(weights1) > (3 * nNodes1) || length(weights2) > (3 * nNodes2)) {
      cut1 <- max(sort(abs(weights1), decreasing = TRUE)[2 * nNodes1],
                  quantile(abs(weights1), 0.75))
      cut2 <- ifelse(twoNets, 
                     max(sort(abs(weights2), decreasing = TRUE)[2 * nNodes2],
                         quantile(abs(weights2), 0.75)), 0)
      
      
    } else if (length(weights1) > 1 || length(weights2) > 1) {
      cut1 <- quantile(abs(weights1), 0.75)
      cut2 <- ifelse(twoNets, quantile(abs(weights2), 0.75), 0)
      
    } else {
      cut1 <- 0
      cut2 <- 0
    }
    
    if (twoNets) {
      cut1 <- cut2 <- mean(c(cut1, cut2))
    }
    
  }
  
  #==========================================================================
  # Layout
  
  if (twoNets & sameLayout & layoutGroup == "union") {
    graph1 <- igraph::graph_from_adjacency_matrix(adja1, weighted = TRUE)
    graph2 <- igraph::graph_from_adjacency_matrix(adja2, weighted = TRUE)
    
    graph_u <- igraph::union(graph1, graph2)
    
    # element-wise minimum of edge weights
    E(graph_u)$weight <- pmin(E(graph_u)$weight_1, 
                              E(graph_u)$weight_2, 
                              na.rm = TRUE) 
    
    graph_u <- igraph::delete_edge_attr(graph_u, "weight_1")
    graph_u <- igraph::delete_edge_attr(graph_u, "weight_2")
    
    if (is.function(layout)) {
      lay1 <- layout(graph_u)
      
    } else {
      adja_u <- as.matrix(as_adjacency_matrix(graph_u, attr = "weight",
                                              sparse = T))
      
      lay1 <- qgraph(adja_u, color = nodecol1, layout = layout, 
                     vsize = nodeSize1, labels = colnames(adja1),
                     label.scale = labelScale, 
                     border.color = border1, border.width = borderWidth1,
                     label.font = labelFont1, label.cex = cexLabels,
                     edge.width = edgeWidth, repulsion = repulsion, cut = cut1,
                     shape = nodeShape1, DoNotPlot = TRUE, ...)$layout
    }
    
    rownames(lay1) <- rownames(adja1)
    lay2 <- lay1
    
  } else {
    # Group 1
    if (is.matrix(layout)) {
      lay1 <- layout
      
    } else if (is.list(layout)) {
      lay1 <- layout[[1]]
      
    } else if (is.function(layout)) {
      gr <- graph_from_adjacency_matrix(adja1, mode = "undirected",
                                        weighted = TRUE, diag = FALSE)
      lay1 <- layout(gr)
      rownames(lay1) <- colnames(adja1)
    } else {
      lay1 <- qgraph(adja1, color = nodecol1, layout = layout, vsize = nodeSize1,
                     labels = colnames(adja1),label.scale = labelScale,
                     border.color = border1, border.width = borderWidth1,
                     label.font = labelFont1, label.cex = cexLabels,
                     edge.width = edgeWidth, repulsion = repulsion, cut = cut1,
                     DoNotPlot = TRUE, shape = nodeShape1, ...)$layout
      rownames(lay1) <- rownames(adja1)
    }
    lay2 <- NULL
    
    #--------------------------------------------
    # Group 2
    
    if (twoNets) {
      if (is.matrix(layout)) {
        lay2 <- layout
        
      } else if (is.list(layout)) {
        lay2 <- layout[[2]]
        
      } else if (is.function(layout)) {
        gr <- igraph::graph_from_adjacency_matrix(adja2, mode = "undirected",
                                                  weighted = TRUE, diag = FALSE)
        
        lay2 <- layout(gr)
        rownames(lay2) <- colnames(adja2)
        
      } else {
        lay2 <- qgraph(adja2,
                       color = nodecol2,
                       layout = layout,
                       vsize = nodeSize2,
                       labels = colnames(adja2),
                       label.scale = labelScale,
                       border.color = border2,
                       border.width = borderWidth2,
                       label.font = labelFont2,
                       label.cex = cexLabels,
                       edge.width = edgeWidth,
                       repulsion = repulsion,
                       cut = cut2,
                       DoNotPlot = TRUE,
                       shape = nodeShape2,
                       ... )$layout
        
        rownames(lay2) <- rownames(adja2)
        
      }
      
      if (sameLayout) {
        if (layoutGroup == 1) {
          lay2 <- lay1
        } else {
          lay1 <- lay2
        }
      }
    }
  }
  
  #==========================================================================
  # Edge colors
  
  if (is.null(posCol)) {
    poscol1 <- "#009900"
    poscol2 <- "darkgreen"
  } else {
    if (length(posCol) == 2) {
      poscol1 <- posCol[1]
      poscol2 <- posCol[2]
    } else {
      poscol1 <- poscol2 <- posCol
    }
    
  }
  
  if (is.null(negCol)) {
    negcol1 <- "red"
    negcol2 <- "#BF0000"
  } else if (negCol == FALSE) {
    negcol1 <- poscol1
    negcol2 <- poscol2
  } else {
    if (length(negCol) == 2) {
      negcol1 <- negCol[1]
      negcol2 <- negCol[2]
    } else {
      negcol1 <- negcol2 <- negCol
    }
  }
  
  if (edgeTranspLow > 0) {  # transparency for values below cut
    poscol1 <- colToTransp(poscol1, edgeTranspLow)
    negcol1 <- colToTransp(negcol1, edgeTranspLow)
  }
  
  if (edgeTranspHigh > 0) {  # transparency for values above cut
    poscol2 <- colToTransp(poscol2, edgeTranspHigh)
    negcol2 <- colToTransp(negcol2, edgeTranspHigh)
  }
  
  if (!is.null(x$input$assoEst1)) {
    assoMat1 <- x$input$assoEst1
  } else {
    assoMat1 <- x$input$dissScale1
  }
  
  #--------------------------------------------
  # Edge colors (group 1)
  
  if (x$input$assoType == "dissimilarity") {
    assoMat1 <- x$input$dissEst1
  } else {
    assoMat1 <- x$input$assoEst1
  }
  
  if (negDiffCol) {
    colmat1 <- matrix(NA, ncol(assoMat1), ncol(assoMat1), 
                      dimnames = dimnames(assoMat1))
    colmat1[assoMat1 < 0] <- negcol1
    colmat1[assoMat1 >= 0] <- poscol1
    
    colmat1 <- colmat1[kept1, kept1]
    
    colmat1[abs(adja1) >= cut1 & colmat1 == poscol1] <- poscol2
    colmat1[abs(adja1) >= cut1 & colmat1 == negcol1] <- negcol2
    
  } else {
    colmat1 <- matrix(poscol1, ncol(assoMat1), ncol(assoMat1), 
                      dimnames = dimnames(assoMat1))
    colmat1 <- colmat1[kept1, kept1]
    colmat1[abs(adja1) >= cut1 & colmat1 == poscol1] <- poscol2
  }
  
  #--------------------------------------------
  
  # # plot dissimilarity values instead of similarities
  # if (plotdiss) {
  #   adja1 <- 1-adja1
  #   adja1[adja1 == 1] <- 0
  #   
  #   adja2 <- 1-adja2
  #   adja2[adja2 == 1] <- 0
  # }
  
  #--------------------------------------------
  # Labels (group 1)
  
  labels1.orig <- rownames(adja1)
  
  if (is.null(labels) || (is.logical(labels) && labels == TRUE)) {
    
    labels1 <- editLabels(rownames(adja1), 
                          shortenLabels = shortenLabels,
                          labelLength = labelLength, 
                          labelPattern = labelPattern,
                          addBrack = TRUE,
                          charToRm = charToRm,
                          verbose = FALSE)
    
  } else if (is.logical(labels) && labels == FALSE) {
    labels1 <- labels
    
  } else if (is.list(labels)) {
    stopifnot(all(colnames(adja1) %in% names(labels[[1]])))
    labels1 <- labels[[1]][colnames(adja1)]
    
  } else if (is.vector(labels)) {
    stopifnot(all(colnames(adja1) %in% names(labels)))
    labels1 <- labels[colnames(adja1)]
    
  } else {
    stop("Argument 'labels' must be either a named vector with a label for each 
         node, or a list of length 2 naming the nodes in each network.")
  }
  
  #--------------------------------------------
  # Label size (group 1)
  
  if (!is.null(labels) && is.logical(labels) && labels == FALSE) {
    cexLabels1 <- 0
    
  } else {
    if (highlightHubs) {
      if (is.null(cexHubLabels)) {
        cexHubLabels <- cexLabels
      }
      
      cexLabels1 <- rep(cexLabels, length(labels1))
      cexLabels1[match(hubs1, rownames(adja1))] <- cexHubLabels
      
    } else {
      cexLabels1 <- cexLabels
    }
    
    if (any(cexLabels1 == 0)) {
      labels1[cexLabels1 == 0] <- ""
    }
  }
  
  #--------------------------------------------
  # Filter edges without influencing the layout (group 1)
  
  if (edgeInvisFilter != "none") {
    adja1 <- .filterEdges(adja1, 
                          assoEst = x$input$assoMat1,
                          dissEst = dissMat1,
                          edgeFilter = edgeInvisFilter,
                          edgeFilterPar = edgeInvisPar)
  }
  
  #==========================================================================
  if (twoNets) {
    # Edge colors (group 2)
    
    if (x$input$assoType == "dissimilarity") {
      assoMat2 <- x$input$dissEst2
    } else {
      assoMat2 <- x$input$assoEst2
    }
    
    
    if (negDiffCol) {
      colmat2 <- matrix(NA, ncol(assoMat2), ncol(assoMat2),
                        dimnames = dimnames(assoMat2))
      
      colmat2[assoMat2 < 0] <- negcol1
      colmat2[assoMat2 >= 0] <- poscol1
      
      colmat2 <- colmat2[kept2, kept2]
      
      colmat2[abs(adja2) >= cut2 & colmat2 == poscol1] <- poscol2
      colmat2[abs(adja2) >= cut2 & colmat2 == negcol1] <- negcol2
      
    } else {
      colmat2 <- matrix(poscol1, ncol(assoMat2), ncol(assoMat2),
                        dimnames = dimnames(assoMat2))
      
      colmat2 <- colmat2[kept2, kept2]
      colmat2[abs(adja2) >= cut2 & colmat2 == poscol1] <- poscol2
    }
    
    #--------------------------------------------
    # Labels (group 2)
    labels2.orig <- rownames(adja2)
    
    if (is.null(labels) || (is.logical(labels) && labels == TRUE)) {
      labels2 <- editLabels(rownames(adja2), 
                            shortenLabels = shortenLabels,
                            labelLength = labelLength, 
                            labelPattern = labelPattern,
                            addBrack = TRUE,
                            charToRm = charToRm,
                            verbose = FALSE)
      
    } else if (is.logical(labels) && labels == FALSE) {
      labels2 <- labels
      
    } else if (is.list(labels)) {
      stopifnot(all(colnames(adja2) %in% names(labels[[2]])))
      labels2 <- labels[[2]][colnames(adja2)]
      
    } else if (is.vector(labels)) {
      stopifnot(all(colnames(adja2) %in% names(labels)))
      labels2 <- labels[colnames(adja2)]
    }
    
    # store label names to file
    if (!is.null(labelFile) && !is.logical(labels)) {
      stopifnot(is.character(labelFile))
      
      labelMat <- cbind(labels1, labels1.orig)
      labelMat <- rbind(c("Renamed", "Original"), labelMat)
      labelMat <- rbind(c(" ", " "), labelMat)
      labelMat <- rbind(c("Network 1:", " "), labelMat)
      
      dframe <- data.frame(labelMat, stringsAsFactors=FALSE)
      
      # apply format over each column for alignment
      if (shortenLabels != "none") {
        formWidth <- max(suppressWarnings(sum(as.numeric(labelPattern), 
                                              na.rm = TRUE)),
                         labelLength) + 4
      } else {
        formWidth <- 12
      }
      
      dframe <- apply(dframe, 2, format, width = formWidth)
      
      utils::write.table(dframe, labelFile, quote=FALSE, row.names=FALSE, 
                         col.names = FALSE, append = FALSE)
      
      labelMat <- cbind(labels2, labels2.orig)
      labelMat <- rbind(c("Renamed", "Original"), labelMat)
      labelMat <- rbind(c(" ", " "), labelMat)
      labelMat <- rbind(c("Network 2:", " "), labelMat)
      labelMat <- rbind(c(" ", " "), labelMat)
      
      dframe <- data.frame(labelMat, stringsAsFactors=FALSE)
      
      # apply format over each column for alignment
      dframe <- apply(dframe, 2, format, width = formWidth)
      
      utils::write.table(dframe, labelFile, quote=FALSE, row.names=FALSE,
                         col.names = FALSE, append = TRUE)
    }
    
    #--------------------------------------------
    # Label size (group 2)
    
    if (!is.null(labels) && is.logical(labels) && labels == FALSE) {
      cexLabels2 <- 0
      
    } else {
      if (highlightHubs) {
        cexLabels2 <- rep(cexLabels, length(labels2))
        cexLabels2[match(hubs2, rownames(adja2))] <- cexHubLabels
        
      } else {
        cexLabels2 <- cexLabels
      }
      
      if (any(cexLabels2 == 0)) {
        labels2[cexLabels2 == 0] <- ""
      }
    }
    
    #--------------------------------------------
    # Filter edges without influencing the layout (group 2)
    
    if (edgeInvisFilter != "none") {
      adja2 <- .filterEdges(adja2, 
                            assoEst = x$input$assoMat2,
                            dissEst = dissMat2,
                            edgeFilter = edgeInvisFilter,
                            edgeFilterPar = edgeInvisPar)
    }
    
  } else { #single network
    
    # store label names to file
    if (!is.null(labelFile) && !is.logical(labels)) {
      stopifnot(is.character(labelFile))
      
      labelMat <- cbind(labels1, labels1.orig)
      labelMat <- rbind(c("Renamed", "Original"), labelMat)
      
      dframe <- data.frame(labelMat, stringsAsFactors=FALSE)
      
      # apply format over each column for alignment
      if (shortenLabels != "none") {
        formWidth <- max(suppressWarnings(sum(as.numeric(labelPattern), 
                                              na.rm = TRUE)),
                         labelLength) + 4
      } else {
        formWidth <- 12
      }
      dframe <- apply(dframe, 2, format, width = formWidth)
      
      utils::write.table(dframe, labelFile, quote=FALSE, row.names=FALSE, 
                         col.names = FALSE, append = FALSE)
    }
  }
  
  #=============================================================================
  # Plot network(s)
  
  if (twoNets) {
    par(mfrow = c(1,2), xpd=NA)
    
    if (groupsChanged) {
      q2 <- qgraph(adja2,
                   color = nodecol2,
                   layout = lay2,
                   vsize = nodeSize2,
                   labels = labels2,
                   label.scale = labelScale,
                   border.color = border2,
                   border.width = borderWidth2,
                   label.font = labelFont2,
                   label.cex = cexLabels2,
                   edge.width = edgeWidth,
                   edge.color = colmat2,
                   repulsion = repulsion,
                   cut = cut2,
                   mar = mar,
                   shape = nodeShape2,
                   DoNotPlot = !doPlot,
                   ...)
      
      if (showTitle & doPlot) title(main = list(main2, cex = cexTitle))
    }
    
    q1 <- qgraph(adja1,
                 color = nodecol1,
                 layout = lay1,
                 vsize = nodeSize1,
                 labels = labels1,
                 label.scale = labelScale,
                 border.color = border1,
                 border.width = borderWidth1,
                 label.font = labelFont1,
                 label.cex = cexLabels1,
                 edge.width = edgeWidth,
                 edge.color = colmat1,
                 repulsion = repulsion,
                 cut = cut1,
                 mar = mar,
                 shape = nodeShape1,
                 DoNotPlot = !doPlot,
                 ...)
    
    if (showTitle & doPlot) title(main = list(main1, cex = cexTitle))
    
    if (!groupsChanged) {
      q2 <- qgraph(adja2,
                   color = nodecol2,
                   layout = lay2,
                   vsize = nodeSize2,
                   labels = labels2,
                   label.scale = labelScale,
                   border.color = border2,
                   border.width = borderWidth2,
                   label.font = labelFont2,
                   label.cex = cexLabels2,
                   edge.width = edgeWidth,
                   edge.color = colmat2,
                   repulsion = repulsion,
                   cut = cut2,
                   mar = mar,
                   shape = nodeShape2,
                   DoNotPlot = !doPlot,
                   ...)
      
      if (showTitle & doPlot) title(main = list(main2, cex = cexTitle))
    }
    
    #---------------------------------------------
  } else {
    par(xpd=NA)
    
    q1 <- qgraph(adja1,
                 color = nodecol1,
                 layout = lay1,
                 vsize = nodeSize1,
                 labels = labels1,
                 label.scale = labelScale,
                 border.color = border1,
                 border.width = borderWidth1,
                 label.font = labelFont1,
                 label.cex = cexLabels1,
                 edge.width = edgeWidth,
                 edge.color = colmat1,
                 repulsion = repulsion,
                 cut = cut1,
                 mar = mar,
                 shape = nodeShape1,
                 DoNotPlot = !doPlot,
                 ...)
    
    if (showTitle & doPlot) title(main = list(main1, cex = cexTitle))
    
    q2 <- NULL
    
    labels2 <- NULL
  }
  
  layout.out <- list(layout1 = lay1, layout2 = lay2)
  nodecol.out <- list(nodecol1 = nodecol1, nodecol2 = nodecol2)
  labels.out <- list(labels1 = labels1, labels2 = labels2)
  
  
  output <- list(q1 = q1, q2 = q2, layout = layout.out,
                 nodecolor = nodecol.out, labels = labels.out)
  
  #------------------------------------------------------------------
  #------------------------------------------------------------------
  
  par(mfrow = opar$mfrow, mar = opar$mar, xpd = opar$xpd)
  
  invisible(output)
}

