R/run_wgcna.R

In jeffreyLbrabec/tinker: Tissue Network Knowledge Transfer

#' Run WGCNA to group genes into related modules
#'
#' @param expression_data Your expression data from `extract_expr_data`.
#' @param sft_thresh_power The soft-thresholding power output from `pick_soft_threshold`. 
#'     This should be picked from the graph labeled Scale independence. The number closest
#'     to the "knee" of the number graph indicated by the red line.
#' @param type_of_network WGCNA has several network types you could work under. The `WGCNA`
#'     package recommends you use a signed network for your analysis though their are other
#'     options including unsigned and hybrid signed and unsigned.
#' @param correlation_type Type of correlation algorithm to use. The `WGCNA` package has 
#'     several options for `corrType` but bicor is the default.
#'
#' @return A figure with two graphs highlighting mean connectivity and scale free independence.
#' @export
#'
#' @examples
#' WGCNA::cor
#' gse <- download_gse_data("GSE108000")
#' num_data <- extract_expr_data(gse)
#' wgcna_out <- run_wgcna(num_data, sft_thresh_power = 12, type_of_network = "signed", 
#'                        correlation_type = "bicor")
#' 
run_wgcna <- function(expression_data = NULL, sft_thresh_power = NULL, type_of_network = "signed", correlation_type = "bicor") {
  
  WGCNA::cor()
  
  #WGCNA using the power we determined from pickSoftThreshold and a transposed expression table. 
  wgcna_out = WGCNA::blockwiseModules(t(expression_data), power = sft_thresh_power, networkType = type_of_network, corType = correlation_type)
  
  #Since WGCNA takes a bit to run we can save the out put to the disk and then call it back as a data object in our R environtment.
  #saveRDS(wgcnaout, file = "~/Documents/RotationProject/R_code/Drug Repositioning Pipeline/data/WGCNA_output2.RData") #format so this ends up being work_dir/project_dir/data_dir
  
  #Also probably don't need to save this to disk in the function? but might help.
  
  #Read the wgcna data back into R as an accesible data object.
  #wgcna_out = readRDS("~/Documents/RotationProject/R_code/Drug Repositioning Pipeline/data/WGCNA_output2.RData") #ditto to above
  #Moving forward we will only use the wgcna.out object for our analysis.
  
  return(wgcna_out)
  
  on.exit(cat("Genes successfully grouped into modules!\n"))
}