#' 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")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.