Added exp2gcn_blockwise() + increased unit test coverage

NAMESPACE

@@ -14,6 +14,7 @@ export(dfs2one)
 export(enrichment_analysis)
 export(exp2cor)
 export(exp2gcn)
+export(exp2gcn_blockwise)
 export(exp2grn)
 export(exp_genes2orthogroups)
 export(exp_preprocess)
@@ -67,6 +68,7 @@ importFrom(WGCNA,TOMsimilarity)
 importFrom(WGCNA,adjacency)
 importFrom(WGCNA,adjacency.fromSimilarity)
 importFrom(WGCNA,bicor)
+importFrom(WGCNA,blockwiseModules)
 importFrom(WGCNA,checkSets)
 importFrom(WGCNA,consensusMEDissimilarity)
 importFrom(WGCNA,corPvalueFisher)
@@ -79,7 +81,6 @@ importFrom(WGCNA,mergeCloseModules)
 importFrom(WGCNA,moduleEigengenes)
 importFrom(WGCNA,multiSetMEs)
 importFrom(WGCNA,pickSoftThreshold)
-importFrom(WGCNA,plotEigengeneNetworks)
 importFrom(WGCNA,pquantile)
 importFrom(WGCNA,sampledBlockwiseModules)
 importFrom(WGCNA,scaleFreeFitIndex)

R/gcn_inference.R

@@ -116,23 +116,36 @@ SFT_fit <- function(exp, net_type = "signed", rsquared = 0.8,
 #' @param verbose Logical indicating whether to display progress
 #' messages or not. Default: FALSE.
 #'
-#' @return List containing: \itemize{
-#'   \item Adjacency matrix
-#'   \item Data frame of module eigengenes
-#'   \item Data frame of genes and their corresponding modules
-#'   \item Data frame of intramodular connectivity
-#'   \item Correlation matrix
-#'   \item Parameters used for network reconstruction
-#'   \item Objects to plot the dendrogram in \code{plot_dendro_and_colors}.
+#' @return A list containing the following elements: \itemize{
+#'   \item \emph{adjacency_matrix} Numeric matrix with network adjacencies.
+#'   \item \emph{MEs} Data frame of module eigengenes, with samples
+#'   in rows, and module eigengenes in columns.
+#'   \item \emph{genes_and_modules} Data frame with columns 'Genes' (character)
+#'   and 'Modules' (character) indicating the genes and the modules to
+#'   which they belong.
+#'   \item \emph{kIN} Data frame of degree centrality for each gene,
+#'   with columns 'kTotal' (total degree), 'kWithin' (intramodular degree),
+#'   'kOut' (extra-modular degree), and 'kDiff' (difference between
+#'   the intra- and extra-modular degree).
+#'   \item \emph{correlation_matrix} Numeric matrix with pairwise
+#'   correlation coefficients between genes.
+#'   If parameter \strong{return_cormat} is FALSE, this will be NULL.
+#'   \item \emph{params} List with network inference parameters passed
+#'   as input.
+#'   \item \emph{dendro_plot_objects} List with objects to plot the dendrogram
+#'   in \code{plot_dendro_and_colors}. Elements are named 'tree' (an hclust
+#'   object with gene dendrogram), 'Unmerged' (character with per-gene module
+#'   assignments before merging similar modules), and 'Merged' (character
+#'   with per-gene module assignments after merging similar modules).
 #' }
 #' @author Fabricio Almeida-Silva
 #' @rdname exp2gcn
 #' @export
 #' @importFrom WGCNA TOMsimilarity standardColors labels2colors
-#' moduleEigengenes plotEigengeneNetworks mergeCloseModules
-#' intramodularConnectivity
+#' moduleEigengenes mergeCloseModules intramodularConnectivity
 #' @importFrom dynamicTreeCut cutreeDynamicTree
-#' @importFrom stats as.dist median cor fisher.test hclust na.omit prcomp qnorm qqplot quantile sd var
+#' @importFrom stats as.dist median cor fisher.test hclust na.omit prcomp
+#' qnorm qqplot quantile sd var
 #' @importFrom grDevices colorRampPalette
 #' @examples
 #' data(filt.se)
@@ -224,6 +237,114 @@ exp2gcn <- function(
     return(result_list)
 }
 
+#' Infer gene coexpression network from gene expression in a blockwise manner
+#'
+#' @param exp Either a `SummarizedExperiment` object, or a gene expression
+#' matrix/data frame with genes in row names and samples in column names.
+#' @param net_type Character indicating the type of network to infer.
+#' One of 'signed', 'signed hybrid' or 'unsigned'. Default: 'signed'.
+#' @param module_merging_threshold Numeric indicating the minimum correlation
+#' threshold to merge similar modules into a single one. Default: 0.8.
+#' @param SFTpower Numeric scalar indicating the value of the \eqn{\beta}
+#' power to which correlation coefficients will be raised to ensure
+#' scale-free topology fit. This value can be obtained with
+#' the function \code{SFT_fit()}.
+#' @param cor_method Character with correlation method to use.
+#' One of "pearson" or "biweight". Default: "pearson".
+#' @param TOM_type Character specifying the method to use to calculate a
+#' topological overlap matrix (TOM). If NULL, TOM type will be automatically
+#' inferred from network type specified in \strong{net_type}. Default: NULL.
+#' @param max_block_size Numeric indicating the maximum block size for module
+#' detection.
+#' @param ... Additional arguments to \code{WGCNA::blockwiseModules()}.
+#'
+#' @return A list containing the following elements: \itemize{
+#'   \item \emph{MEs} Data frame of module eigengenes, with samples
+#'   in rows, and module eigengenes in columns.
+#'   \item \emph{genes_and_modules} Data frame with columns 'Genes' (character)
+#'   and 'Modules' (character) indicating the genes and the modules to
+#'   which they belong.
+#'   \item \emph{params} List with network inference parameters passed
+#'   as input.
+#'   \item \emph{dendro_plot_objects} List with objects to plot the dendrogram
+#'   in \code{plot_dendro_and_colors}. Elements are named 'tree' (an hclust
+#'   object with gene dendrogram), 'Unmerged' (character with per-gene module
+#'   assignments before merging similar modules), and 'Merged' (character
+#'   with per-gene module assignments after merging similar modules).
+#' }
+#'
+#' @author Fabricio Almeida-Silva
+#' @rdname exp2gcn_blockwise
+#' @export
+#' @importFrom WGCNA blockwiseModules
+#' @examples
+#' data(filt.se)
+#' # The SFT fit was previously calculated and the optimal power was 16
+#' cor <- WGCNA::cor
+#' gcn <- exp2gcn_blockwise(
+#'     exp = filt.se, SFTpower = 18, cor_method = "pearson"
+#' )
+exp2gcn_blockwise <- function(
+        exp, net_type = "signed", module_merging_threshold = 0.8,
+        SFTpower = NULL, cor_method = "pearson", TOM_type = NULL,
+        max_block_size = 5000, ...
+) {
+
+    params <- list(
+        net_type = net_type,
+        module_merging_threshold = module_merging_threshold,
+        SFTpower = SFTpower,
+        cor_method = cor_method
+    )
+    exp <- handleSE(exp)
+
+    if(is.null(SFTpower)) { stop("Please, specify the SFT power.") }
+
+    cortype <- "pearson"
+    if(cor_method == "biweight") {
+        cortype <- "bicor"
+    } else if(!cor_method %in% c("pearson", "biweight")) {
+        stop("Argument to `cor_method` must be one of 'pearson' or 'biweight'.")
+    }
+
+    if(is.null(TOM_type)) { TOM_type <- guess_tom(net_type) }
+
+    # Detect modules
+    cor <- WGCNA::cor
+    bmod <- WGCNA::blockwiseModules(
+        datExpr = as.matrix(t(exp)),
+        checkMissingData = FALSE,
+        maxBlockSize = max_block_size,
+        corType = cortype,
+        maxPOutliers = 0.1,
+        power = SFTpower,
+        networkType = net_type,
+        TOMType = TOM_type,
+        minModuleSize = 30,
+        mergeCutHeight = 1 - module_merging_threshold,
+        ...
+    )
+
+    # Create result list
+    genes_modules <- data.frame(
+        Genes = names(bmod$colors), Modules = bmod$colors, row.names = NULL
+    )
+
+    result_list <- list(
+        MEs = bmod$MEs,
+        genes_and_modules = genes_modules,
+        params = params,
+        dendro_plot_objects = list(
+            tree = bmod$dendrograms[[1]],
+            Unmerged = bmod$unmergedColors,
+            Merged = bmod$colors
+        )
+    )
+
+    return(result_list)
+}
+
+
 
 #' Plot eigengene network
 #'

R/network_comparison.R

@@ -135,16 +135,14 @@ modPres_WGCNA <- function(explist, ref_net, nPerm = 200) {
 
     # Calculate module preservation
     if(cor_method == "pearson") {
-      corOptions <- "use = 'p'"
-      corFnc <- "cor"
+        corOptions <- "use = 'p'"
+        corFnc <- "cor"
     } else if(cor_method == "spearman") {
-      corOptions <- list(use = 'p', method="spearman")
-      corFnc <- "cor"
+        corOptions <- "use = 'p', method = 'spearman'"
+        corFnc <- "cor"
     } else if(cor_method == "biweight") {
-      corOptions <- list(use = 'p', maxPOutliers = 0.05)
-      corFnc <- "bicor"
-    } else {
-      stop("Please, specify a valid correlation method.")
+        corOptions <- "use = 'p', maxPOutliers = 0.05"
+        corFnc <- "bicor"
     }
 
     pres <- WGCNA::modulePreservation(
@@ -261,12 +259,13 @@ modPres_netrep <- function(explist, ref_net = NULL, test_net = NULL,
     explist <- lapply(explist, function(x) return(t(x)))
 
     # Set data set names
+    data_names <- names(explist)
     if(is.null(names(explist))) {
         data_names <- c("cohort1", "cohort2")
-    } else {
-        data_names <- names(explist)
+        names(explist) <- data_names
     }
 
+
     # Create correlation list
     correlation_list <- list(
         as.matrix(ref_net$correlation_matrix),
@@ -281,28 +280,26 @@ modPres_netrep <- function(explist, ref_net = NULL, test_net = NULL,
     )
     names(network_list) <- data_names
 
-    # Set background label
-    if("grey" %in% ref_net$moduleColors) {
-        backgroundLabel <- "grey"
-    } else {
-        backgroundLabel <- 0
-    }
-
     # Set module assignments
     modAssignments <- ref_net$genes_and_modules$Modules
     names(modAssignments) <- ref_net$genes_and_modules$Genes
 
+    # Set background label
+    backgroundLabel <- "grey"
+    if(!"grey" %in% modAssignments) { backgroundLabel <- 0 }
+
     # Calculate preservation statistics
     pres <- NetRep::modulePreservation(
         network = network_list, data = explist, correlation = correlation_list,
         moduleAssignments = modAssignments,
         discovery = data_names[1], test = data_names[2],
-        nPerm = nPerm, nThreads = nThreads
+        nPerm = nPerm, nThreads = nThreads, verbose = FALSE
     )
 
     # Get preserved modules (p < 0.05 for all statistics)
-    max_pval <- apply(pres$p.value, 1, max)
+    max_pval <- apply(pres$p.value, 1, max, na.rm = TRUE)
     preservedmodules <- names(max_pval[max_pval < 0.05])
+
     if(length(preservedmodules) > 0) {
         message(length(preservedmodules), " modules in ", data_names[1],
                 " were preserved in ", data_names[2], ":", "\n",

tests/testthat/test_comparison.R

@@ -60,7 +60,9 @@ test_that("modPres_WGCNA() calculates module preservation", {
 # })
 
 test_that("module_preservation() calculates module preservation", {
+
     explist <- exp_ortho
+    names(explist) <- NULL
     ref_net <- gcn_osa
     test_net <- gcn_zma
     pres <- suppressWarnings(
@@ -70,6 +72,16 @@ test_that("module_preservation() calculates module preservation", {
         explist, ref_net, nPerm = 2, algorithm = "WGCNA"
     )
 
+    ref_net$params$cor_method <- "biweight"
+    pres_wgcna2 <- module_preservation(
+        explist, ref_net, nPerm = 2, algorithm = "WGCNA"
+    )
+
+    ref_net$params$cor_method <- "spearman"
+    pres_wgcna2 <- module_preservation(
+        explist, ref_net, nPerm = 2, algorithm = "WGCNA"
+    )
+
     expect_equal(class(pres), "list")
     expect_true("ggplot" %in% class(pres_wgcna))
     expect_error(module_preservation(