R/pick_soft_threshold.R

In jeffreyLbrabec/tinker: Tissue Network Knowledge Transfer

#' Pick soft threshold for WGCNA
#'
#' @param expression_data The output from `extract_expr_data`; a table of just expression data.
#' @param type_of_network The `pickSoftThreshold` function of the `WGCNA` package can use either signed, unsigned, or 
#'     hybrid signed or unsigned networks for soft thresholding. The `WGCNA` documentation recommends "signed" be used
#'     and that is the default for this function.
#'
#' @return A figure containing two graphs from which to find the number to use for soft-thresholding.
#' @export
#'
#' @examples
#' gse <- download_gse_data("GSE108000")
#' 
#' num_data <- extract_expr_data(gse)
#' 
#' soft_thresholds <- pick_soft_threshold(num_data, type_of_network = "signed")
pick_soft_threshold <- function(expression_data = NULL, type_of_network = "signed") {
  
  #The commaned pickSoftThreshold requires rows correspond to samples and columns correspond to genes so we need to transpose our expr table.
  transposed <- t(expression_data)
  
  #Set the powers for the soft-threshold to test at
  powers = c(c(1:10), seq(from = 12, to=20, by=2))
  
  sft = WGCNA::pickSoftThreshold(transposed, powerVector = powers, networkType = type_of_network, verbose = 5)
  
  # Plot the results:
  WGCNA::sizeGrWindow(9, 5)
  
  graphics::par(mfrow = c(1,2));
  
  cex1 = 0.9;
  
  # Scale-free topology fit index as a function of the soft-thresholding power. Run plot() and text() command together.
  graphics::plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
       main = paste("Scale independence"));
 
   graphics::text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
       labels=powers,cex=cex1,col="red");
  
  # this line corresponds to using an R^2 cut-off of h
  graphics::abline(h=0.80,col="red")
  
  # Mean connectivity as a function of the soft-thresholding power
  graphics::plot(sft$fitIndices[,1], sft$fitIndices[,5],
       xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
       main = paste("Mean connectivity"))
  
  graphics::text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
  
}