library(WGCNA); library(stringi); library(rjson) options(stringsAsFactors = FALSE); cat("Running the wgcna analysis script\n") # load expression data **assumes from json files row(traits)(columns info+samples) # pass the file_path as arg # pass the file path to read json data args = commandArgs(trailingOnly=TRUE) if (length(args)==0) { stop("Argument for the file location is required", call.=FALSE) } else { # default output file json_file_path = args[1] } inputData <- fromJSON(file = json_file_path) imgDir = inputData$TMPDIR trait_sample_data <- do.call(rbind, inputData$trait_sample_data) dataExpr <- data.frame(apply(trait_sample_data, 2, function(x) as.numeric(as.character(x)))) # transform expressionData dataExpr <- data.frame(t(dataExpr)) gsg = goodSamplesGenes(dataExpr, verbose = 3) if (!gsg$allOK) { if (sum(!gsg$goodGenes)>0) printFlush(paste("Removing genes:", paste(names(dataExpr)[!gsg$goodGenes], collapse = ", "))); if (sum(!gsg$goodSamples)>0) printFlush(paste("Removing samples:", paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ", "))); # Remove the offending genes and samples from the data: dataExpr <- dataExpr[gsg$goodSamples, gsg$goodGenes] } ## network constructions and modules names(dataExpr) = inputData$trait_names # Allow multi-threading within WGCNA enableWGCNAThreads() # choose softthreshhold (Calculate soft threshold) # xtodo allow users to pass args powers <- c(c(1:10), seq(from = 12, to=20, by=2)) sft <- pickSoftThreshold(dataExpr, powerVector = powers, verbose = 5) # check the power estimate if (is.na(sft$powerEstimate)){ powerEst = 1 }else{ powerEst = sft$powerEstimate } # pass user options network <- blockwiseModules(dataExpr, #similarity matrix options corType = inputData$corType, #adjacency matrix options power = powerEst, networkType = "unsigned", #TOM options TOMtype = inputData$TOMtype, #module indentification verbose = 3, minmodulesSize = inputData$minModuleSize, deepSplit = 3, PamRespectsDendro = FALSE ) cat("Generated network \n") network genImageRandStr <- function(prefix){ randStr <- paste(prefix,stri_rand_strings(1, 9, pattern = "[A-Za-z0-9]"),sep="_") return(paste(randStr,".png",sep="")) } mergedColors <- labels2colors(network$colors) imageLoc <- file.path(imgDir,genImageRandStr("WGCNAoutput")) png(imageLoc,width=1000,height=600,type='cairo-png') cat("Generating the CLuster dendrogram\n") plotDendroAndColors(network$dendrograms[[1]],mergedColors[network$blockGenes[[1]]], "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05) json_data <- list(input = inputData, output = list(ModEigens=network$MEs, soft_threshold=sft$fitIndices, blockGenes =network$blockGenes[[1]], net_colors =network$colors, power_used_for_analy=powerEst, net_unmerged=network$unmergedColors, imageLoc=imageLoc)) json_data <- toJSON(json_data) write(json_data,file= json_file_path)