diff options
Diffstat (limited to 'scripts/ctl_analysis.R')
| -rw-r--r-- | scripts/ctl_analysis.R | 217 |
1 files changed, 217 insertions, 0 deletions
diff --git a/scripts/ctl_analysis.R b/scripts/ctl_analysis.R new file mode 100644 index 0000000..e06ac25 --- /dev/null +++ b/scripts/ctl_analysis.R @@ -0,0 +1,217 @@ + +library(ctl) +library(stringi); +library(rjson) + +options(stringsAsFactors = FALSE); + +args = commandArgs(trailingOnly=TRUE) + +imgDir = Sys.getenv("GENERATED_IMAGE_DIR") + +if (length(args)==0) { + stop("Argument for the data file", call.=FALSE) +} else { + + json_file_path = args[1] +} + +json_file_path +# add validation for the files +input <- fromJSON(file = json_file_path) + + +cat("The input data is \n") + + +genoData <- input$genoData +phenoData <- input$phenoData + +parametric <- switch( + input$parametric, + "True" = TRUE, + "False" = FALSE + + ) +# create the matixes + +# genotypes Matrix of genotypes. (individuals x markers) +# phenotypes Matrix of phenotypes. (individuals x phenotypes) + +geno_matrix = t(matrix(unlist(genoData$genotypes), + nrow=length(genoData$markernames), ncol=length(genoData$individuals), + dimnames=list(genoData$markernames, genoData$individuals), byrow=TRUE)) + + +pheno_matrix = t(matrix(as.numeric(unlist(phenoData$traits)), nrow=length(phenoData$trait_db_list), ncol=length( + phenoData$individuals), dimnames= list(phenoData$trait_db_list, phenoData$individuals), byrow=TRUE)) + +ctls <- CTLscan(geno_matrix,pheno_matrix,nperm=input$nperm,strategy=input$strategy,parametric=parametric,nthreads=3,verbose=TRUE) + +genImageRandStr <- function(prefix){ + + randStr <- paste(prefix,stri_rand_strings(1, 9, pattern = "[A-Za-z0-9]"),sep="_") + + return(paste(randStr,".png",sep="")) +} + + +genRandomFileName <- function(prefix,file_ext=".png"){ + + randStr = paste(prefix,stri_rand_strings(1, 9, pattern = "[A-Za-z0-9]"),sep="_") + + return(paste(randStr,file_ext,sep="")) +} + + +# #output matrix significant CTL interactions with 4 columns: trait, marker, trait, lod +ctl_significant <- CTLsignificant(ctls,significance = input$significance) + +colnames(ctl_significant) = c("trait","marker","trait_2","LOD","dcor") + + + +imageLoc = file.path(input$imgDir,genRandomFileName("CTLline")) + +png(imageLoc,width=1000,height=600,type='cairo-png') + +# Create the lineplot +ctl.lineplot(ctls,significance = input$significance, gap = 50, +col = "orange", bg.col = "lightgray", cex = 1, verbose = FALSE) + +dev.off() + + +n = 2 +ctl_plots = c() + +for (trait in phenoData$trait_db_list) +{ + image_loc = file.path(input$imgDir,genRandomFileName(paste("CTL",n,sep=""))) + png(image_loc,width=1000, height=600, type='cairo-png') + plot.CTLobject(ctls,n-1,significance= input$significance, main=paste("Phenotype",trait,sep="")) + + ctl_plots = append(ctl_plots,image_loc) + + dev.off() + n = n + 1 + +} + + + +network_file_name = paste("ctlnet",stri_rand_strings(1, 9, pattern = "[A-Za-z0-9]"),sep="_") +netfile = file.path(input$imgDir,paste(network_file_name,".sif",sep="")) + +nodefile = file.path(input$imgDir,paste(network_file_name,".nodes",sep="")) + +# fn overrides ctlnetwork function to target gn2 use case + +CTLnetworkGn<- function(CTLobject, mapinfo, significance = 0.05, LODdrop = 2, what = c("names","ids"), short = FALSE, add.qtls = FALSE,verbose = TRUE){ + if(missing(CTLobject) || is.null(CTLobject)) stop("argument 'CTLobject' is missing, with no default") + if("CTLscan" %in% class(CTLobject)) CTLobject = list(CTLobject) + if(length(what) > 1) what = what[1] + + results <- NULL + significant <- CTLsignificant(CTLobject, significance, what = "ids") + if(!is.null(significant)){ + all_m <- NULL; all_p <- NULL; + + + cat("",file=netfile); cat("",file=nodefile); + if(verbose) cat("NETWORK.SIF\n") + edges <- NULL + for(x in 1:nrow(significant)){ + data <- as.numeric(significant[x,]) + CTLscan <- CTLobject[[data[1]]] + markern <- rownames(CTLscan$dcor) + traitsn <- colnames(CTLscan$dcor) + name <- ctl.name(CTLscan) + if(what=="ids"){ + tid <- which(traitsn %in% ctl.name(CTLobject[[data[1]]])) + name <- paste("P",tid,sep="") + markern <- paste("M",1:nrow(CTLobject[[data[1]]]$dcor), sep="") + traitsn <- paste("P", 1:ncol(CTLobject[[data[1]]]$dcor), sep="") + } + if(add.qtls){ # Add QTL to the output SIF + bfc <- length(CTLscan$qtl) + above <- which(CTLscan$qtl > -log10(significance)) + qtlnms <- names(above); qtlmid <- 1 + for(m in above){ + cat(name,"\tQTL\t",markern[m],"\tQTL\t",CTLscan$qtl[m],"\n",sep="",file=netfile,append=TRUE) + all_m <- CTLnetwork.addmarker(all_m, mapinfo, markern[data[2]], qtlnms[qtlmid]) + qtlmid <- qtlmid+1 + } + } + lod <- CTLscan$ctl[data[2],data[3]] + qlod1 <- CTLscan$qtl[data[2]] + qlod2 <- qlod1 + edgetype <- NA + if(length(CTLobject) >= data[3]){ # Edge type based on QTL LOD scores + qlod2 <- CTLobject[[data[3]]]$qtl[data[2]] + if((qlod1-qlod2) > LODdrop){ + edgetype <- 1 + }else if((qlod1-qlod2) < -LODdrop){ + edgetype <- -1 + }else{ edgetype <- 0; } + } else { + cat("Warning: Phenotype", data[3], "from", data[1], "no CTL/QTL information\n") + qlod2 <- NA; + } + #Store the results + results <- rbind(results, c(data[1], data[2], data[3], lod, edgetype, qlod1, qlod2)) + + if(nodefile == "" && !verbose){ }else{ + if(short){ + edge <- paste(name,traitsn[data[3]]) + edgeI <- paste(traitsn[data[3]],name) + if(!edge %in% edges && !edgeI %in% edges){ + cat(name, "\t", markern[data[2]],"\t", traitsn[data[3]],"\n",file=netfile, append=TRUE,sep="") + edges <- c(edges,edge) + } + }else{ + cat(name, "\t", "CTL_", data[1],"_",data[3], "\t", markern[data[2]],file=netfile, append=TRUE,sep="") + cat("\tCTL\t", lod, "\n", file=netfile, append=TRUE,sep="") + cat(markern[data[2]], "\t", "CTL_", data[1],"_",data[3], "\t",file=netfile, append=TRUE,sep="") + cat(traitsn[data[3]],"\tCTL\t", lod, "\n", file=netfile,append=TRUE,sep="") + } + } + all_m <- CTLnetwork.addmarker(all_m, mapinfo, markern[data[2]], rownames(CTLscan$dcor)[data[2]]) + all_p <- unique(c(all_p, name, traitsn[data[3]])) + } + colnames(results) <- c("TRAIT1","MARKER","TRAIT2","LOD_C","CAUSAL","LOD_T1","LOD_T2") + if(verbose) cat("NODE.DESCRIPTION\n") + if(nodefile == "" && !verbose){ }else{ + for(m in all_m){ cat(m,"\n", sep="", file=nodefile, append=TRUE); } + for(p in all_p){ cat(p,"\tPHENOTYPE\n", sep="", file=nodefile, append=TRUE); } + } + } + if(!is.null(results)){ + class(results) <- c(class(results),"CTLnetwork") + } + invisible(results) +} + +CTLnetwork.addmarker <- function(markers, mapinfo, name, realname){ + if(!missing(mapinfo)){ + id <- which(rownames(mapinfo) %in% realname) + fname <- paste(name,"\tMARKER\t",mapinfo[id,1],"\t",mapinfo[id,2],sep="") + markers <- unique(c(markers, fname)) + } + return(markers) +} + + +# generating network + + +ctl_network = CTLnetworkGn(ctls, significance = input$significance, LODdrop = 2,short = FALSE, add.qtls = FALSE, verbose = TRUE) + + + +json_data <- list(phenotypes = input$phenoData$trait_db_list,significance_data = ctl_significant,image_loc = imageLoc,ctl_plots=ctl_plots,network_file_name = network_file_name) + +json_data <- toJSON(json_data) + +write(json_data,file= json_file_path) + |