Simplified exp2gcn() by writing wrappers + added the option to not re…

…turn correlation matrix in exp2gcn()
almeidasilvaf · Mar 22, 2024 · ad6be31 · ad6be31
1 parent 1a3c1ca
commit ad6be31
Show file tree

Hide file tree

Showing 11 changed files with 394 additions and 158 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -65,7 +65,7 @@ Imports:
     minet, 
     GENIE3,
     SummarizedExperiment
-RoxygenNote: 7.3.0
+RoxygenNote: 7.3.1
 Suggests: 
     knitr,
     rmarkdown,

diff --git a/NAMESPACE b/NAMESPACE
@@ -7,10 +7,12 @@ export(check_SFT)
 export(consensus_SFT_fit)
 export(consensus_modules)
 export(consensus_trait_cor)
+export(cor2adj)
 export(cormat_to_edgelist)
 export(detect_communities)
 export(dfs2one)
 export(enrichment_analysis)
+export(exp2cor)
 export(exp2gcn)
 export(exp2grn)
 export(exp_genes2orthogroups)
@@ -77,7 +79,6 @@ importFrom(WGCNA,mergeCloseModules)
 importFrom(WGCNA,moduleEigengenes)
 importFrom(WGCNA,multiSetMEs)
 importFrom(WGCNA,pickSoftThreshold)
-importFrom(WGCNA,plotDendroAndColors)
 importFrom(WGCNA,plotEigengeneNetworks)
 importFrom(WGCNA,pquantile)
 importFrom(WGCNA,sampledBlockwiseModules)

diff --git a/R/BioNERO-package.R b/R/BioNERO-package.R
@@ -0,0 +1,6 @@
+#' @keywords internal
+"_PACKAGE"
+
+## usethis namespace: start
+## usethis namespace: end
+NULL
diff --git a/R/gcn_inference.R b/R/gcn_inference.R
@@ -93,18 +93,26 @@ SFT_fit <- function(exp, net_type = "signed", rsquared = 0.8,
     return(result)
 }
 
-
-#' Reconstruct gene coexpression network from gene expression
+#' Infer gene coexpression network from gene expression
 #'
-#' @param exp A gene expression data frame with genes in row names and
-#' samples in column names or a `SummarizedExperiment` object.
-#' @param net_type Network type. One of 'signed', 'signed hybrid' or
-#' 'unsigned'. Default: 'signed'.
-#' @param module_merging_threshold Correlation threshold to merge similar
-#' modules into a single one. Default: 0.8.
-#' @param SFTpower SFT power generated by the function \code{SFT_fit}.
-#' @param cor_method Correlation method. One of "pearson", "biweight" or
-#' "spearman". Default is "spearman".
+#' @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", "biweight" or "spearman". Default: "spearman".
+#' @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 return_cormat Logical indicating whether the correlation matrix
+#' should be returned. If TRUE (default), an element named `correlation_matrix`
+#' containing the correlation matrix will be included in the result list.
 #' @param verbose Logical indicating whether to display progress
 #' messages or not. Default: FALSE.
 #'
@@ -118,22 +126,22 @@ SFT_fit <- function(exp, net_type = "signed", rsquared = 0.8,
 #'   \item Objects to plot the dendrogram in \code{plot_dendro_and_colors}.
 #' }
 #' @author Fabricio Almeida-Silva
-#' @seealso
-#'  \code{\link[WGCNA]{adjacency.fromSimilarity}},\code{\link[WGCNA]{TOMsimilarity}},\code{\link[WGCNA]{standardColors}},\code{\link[WGCNA]{labels2colors}},\code{\link[WGCNA]{moduleEigengenes}},\code{\link[WGCNA]{plotEigengeneNetworks}},\code{\link[WGCNA]{mergeCloseModules}},\code{\link[WGCNA]{plotDendroAndColors}},\code{\link[WGCNA]{intramodularConnectivity}}
-#'  \code{\link[dynamicTreeCut]{cutreeDynamicTree}}
 #' @rdname exp2gcn
 #' @export
-#' @importFrom WGCNA TOMsimilarity standardColors labels2colors moduleEigengenes plotEigengeneNetworks mergeCloseModules plotDendroAndColors intramodularConnectivity
+#' @importFrom WGCNA TOMsimilarity standardColors labels2colors
+#' moduleEigengenes plotEigengeneNetworks mergeCloseModules
+#' intramodularConnectivity
 #' @importFrom dynamicTreeCut cutreeDynamicTree
 #' @importFrom stats as.dist median cor fisher.test hclust na.omit prcomp qnorm qqplot quantile sd var
 #' @importFrom grDevices colorRampPalette
 #' @examples
 #' data(filt.se)
 #' # The SFT fit was previously calculated and the optimal power was 16
-#' gcn <- exp2gcn(filt.se, SFTpower = 18, cor_method = "pearson")
+#' gcn <- exp2gcn(exp = filt.se, SFTpower = 18, cor_method = "pearson")
 exp2gcn <- function(
         exp, net_type = "signed", module_merging_threshold = 0.8,
-        SFTpower = NULL, cor_method = "spearman", verbose = FALSE
+        SFTpower = NULL, cor_method = "spearman", TOM_type = NULL,
+        return_cormat = TRUE, verbose = FALSE
 ) {
 
     params <- list(
@@ -142,108 +150,81 @@ exp2gcn <- function(
         SFTpower = SFTpower,
         cor_method = cor_method
     )
-    norm.exp <- handleSE(exp)
+    exp <- handleSE(exp)
 
     if(is.null(SFTpower)) { stop("Please, specify the SFT power.") }
 
+    # Calculate adjacency matrix
     if(verbose) { message("Calculating adjacency matrix...") }
-    cor_matrix <- calculate_cor_adj(cor_method, norm.exp, SFTpower, net_type)$cor
-    adj_matrix <- calculate_cor_adj(cor_method, norm.exp, SFTpower, net_type)$adj
-
-    #Convert to matrix
-    gene_ids <- rownames(adj_matrix)
-    adj_matrix <- matrix(adj_matrix, nrow=nrow(adj_matrix))
-    rownames(adj_matrix) <- gene_ids
-    colnames(adj_matrix) <- gene_ids
+    cor_matrix <- NULL
+    if(return_cormat) {
+        cor_matrix <- exp2cor(exp, cor_method = cor_method)
+        adj_matrix <- cor2adj(cor_matrix, beta = SFTpower, net_type = net_type)
+    } else {
+        adj_matrix <- exp2cor(exp, cor_method = cor_method) |>
+            cor2adj(beta = SFTpower, net_type = net_type)
+    }
 
-    #Calculate TOM from adjacency matrix
+    # Calculate TOM from adjacency matrix
     if(verbose) { message("Calculating topological overlap matrix (TOM)...") }
-    tomtype <- get_TOMtype(net_type)
-    TOM <- WGCNA::TOMsimilarity(adj_matrix, TOMType = tomtype)
+    if(is.null(TOM_type)) { TOM_type <- guess_tom(net_type) }
+    TOM <- WGCNA::TOMsimilarity(adj_matrix, TOMType = TOM_type)
 
-    #Hierarchically cluster genes
-    dissTOM <- 1-TOM #hclust takes a distance structure
-    geneTree <- hclust(as.dist(dissTOM), method="average")
-    geneTree$height <- round(geneTree$height, 7) # to address rounding issue with cutree
+    # Hierarchically cluster genes
+    geneTree <- hclust(as.dist(1 - TOM), method = "average")
+    geneTree$height <- round(geneTree$height, 7)
 
-    #Detecting coexpression modules
+    # Detect coexpression modules
     if(verbose) { message("Detecting coexpression modules...") }
-    old.module_labels <- dynamicTreeCut::cutreeDynamicTree(
+    original_mods <- dynamicTreeCut::cutreeDynamicTree(
         dendro = geneTree, minModuleSize = 30, deepSplit = TRUE
     )
 
-    nmod <- length(unique(old.module_labels))
+    nmod <- length(unique(original_mods))
     palette <- rev(WGCNA::standardColors(nmod))
-    old.module_colors <- WGCNA::labels2colors(
-        old.module_labels, colorSeq = palette
-    )
+    original_colors <- WGCNA::labels2colors(original_mods, colorSeq = palette)
 
-    #Calculate eigengenes and merge the ones who are highly correlated
+    # Calculate module eigengenes and pairwise distances between them
     if(verbose) { message("Calculating module eigengenes (MEs)...") }
-    old.MElist <- WGCNA::moduleEigengenes(
-        t(norm.exp), colors = old.module_colors, softPower = SFTpower
+    me_list <- WGCNA::moduleEigengenes(
+        t(exp), colors = original_colors, softPower = SFTpower
+    )
+    me <- me_list$eigengenes
+    original_metree <- hclust(
+        as.dist(1 - cor(me, method = "spearman")), method = "average"
     )
-    old.MEs <- old.MElist$eigengenes
-
-    #Calculate dissimilarity of module eigengenes
-    MEDiss1 <- 1-cor(old.MEs)
-
-    #Hierarchically cluster module eigengenes to see how they're related
-    old.METree <- hclust(as.dist(MEDiss1), method="average")
-
-    #Then, choose a height cut
-    MEDissThreshold <- 1-module_merging_threshold
 
-    #Merge the modules.
+    # Merge very similar modules
     if(verbose) { message("Merging similar modules...") }
-    if(cor_method == "pearson") {
-        merge1 <- WGCNA::mergeCloseModules(
-            t(norm.exp), old.module_colors, cutHeight = MEDissThreshold,
-            verbose = 0, colorSeq = palette
-        )
-    } else if(cor_method == "spearman") {
-        merge1 <- WGCNA::mergeCloseModules(
-            t(norm.exp), old.module_colors, cutHeight = MEDissThreshold,
-            verbose = 0, corOptions = list(use = "p", method = "spearman"),
-            colorSeq = palette
-        )
-    } else if(cor_method == "biweight") {
-        merge1 <- WGCNA::mergeCloseModules(
-            t(norm.exp), old.module_colors, cutHeight = MEDissThreshold,
-            verbose = 0, corFnc = bicor, colorSeq = palette
-        )
-    } else {
-        stop("Please, specify a correlation method. One of 'spearman', 'pearson' or 'biweight'.")
-    }
-    new.module_colors <- merge1$colors
-    new.MEs <- merge1$newMEs
-
-    #Plot dendrogram of merged modules eigengenes
-    new.METree <- hclust(as.dist(1-cor(new.MEs)), method="average")
+    merged <- merge_modules(
+        exp, original_colors, me, palette,
+        dissimilarity = 1 - module_merging_threshold, cor_method = cor_method
+    )
+    new_colors <- merged$colors
+    new_mes <- merged$newMEs
 
-    #Get data frame with genes and modules they belong to
-    genes_and_modules <- as.data.frame(cbind(gene_ids, new.module_colors))
-    colnames(genes_and_modules) <- c("Genes", "Modules")
+    # Get data frame with genes and modules they belong to
+    genes_modules <- data.frame(Genes = rownames(exp), Modules = new_colors)
 
     if(verbose) { message("Calculating intramodular connectivity...") }
-    kwithin <- WGCNA::intramodularConnectivity(adj_matrix, new.module_colors)
+    kwithin <- WGCNA::intramodularConnectivity(adj_matrix, new_colors)
 
-    result.list <- list(
+    result_list <- list(
         adjacency_matrix = adj_matrix,
-        MEs = new.MEs,
-        genes_and_modules = genes_and_modules,
+        MEs = new_mes,
+        genes_and_modules = genes_modules,
         kIN = kwithin,
         correlation_matrix = cor_matrix,
         params = params,
         dendro_plot_objects = list(
-            tree = geneTree,
-            Unmerged = old.module_colors,
-            Merged = new.module_colors
+            tree = geneTree, Unmerged = original_colors, Merged = new_colors
         )
     )
-    return(result.list)
+
+    return(result_list)
 }
 
+
 #' Plot eigengene network
 #'
 #' @param gcn List object returned by \code{exp2gcn}.
@@ -421,7 +402,7 @@ module_stability <- function(exp, net, nRuns = 20) {
         checkSoftPower = FALSE,
         corType = cor_method,
         networkType = net$params$net_type,
-        TOMType = get_TOMtype(net$params$net_type),
+        TOMType = guess_tom(net$params$net_type),
         TOMDenom = "mean",
         mergeCutHeight = 1 - net$params$module_merging_threshold,
         reassignThreshold = 0,