about summary refs log tree commit diff
diff options
context:
space:
mode:
authorAlexander_Kabui2024-10-25 16:38:09 +0300
committerAlexander_Kabui2024-10-25 16:42:03 +0300
commit436163a37a2c2b0843b122508d1caee47049e7b9 (patch)
tree9978b398bc21ca2c2472f9b90e80c313f8a6c904
parent7c18d1e49db9ac565f9a46be7233e4b8e82c498d (diff)
downloadgenenetwork3-436163a37a2c2b0843b122508d1caee47049e7b9.tar.gz
Refactor: apply code formatting.
-rw-r--r--scripts/rqtl2_wrapper.R279
1 files changed, 174 insertions, 105 deletions
diff --git a/scripts/rqtl2_wrapper.R b/scripts/rqtl2_wrapper.R
index 1fe754e..c6d1911 100644
--- a/scripts/rqtl2_wrapper.R
+++ b/scripts/rqtl2_wrapper.R
@@ -7,16 +7,16 @@ library(rjson)
 library(stringi)
 library(stringr)
 
-options(stringsAsFactors = FALSE);
-args = commandArgs(trailingOnly=TRUE)
+options(stringsAsFactors = FALSE)
+
+args = commandArgs(trailingOnly = TRUE)
 
 # get the json file path with pre metadata required to create the cross
 
 
-if (length(args)==0) {
-  stop("Argument for the metadata file is Missing ", call.=FALSE)
+if (length(args) == 0) {
+  stop("Argument for the metadata file is Missing ", call. = FALSE)
 } else {
-
   json_file_path = args[1]
 }
 
@@ -24,18 +24,19 @@ if (length(args)==0) {
 
 
 if (!(file.exists(json_file_path))) {
-   stop("The input file path does not exists")
+  stop("The input file path does not exists")
 } else {
-str_glue("The input path for the metadata >>>>>>> {json_file_path}")
-json_data  = fromJSON(file = json_file_path)
+  str_glue("The input path for the metadata >>>>>>> {json_file_path}")
+  json_data  = fromJSON(file = json_file_path)
 }
 
 
 
 # generate random string file path here
-genRandomFileName <- function(prefix,file_ext=".txt"){
-	randStr = paste(prefix,stri_rand_strings(1, 9, pattern = "[A-Za-z0-9]"),sep="_")
-	return(paste(randStr,file_ext,sep=""))
+genRandomFileName <- function(prefix, file_ext = ".txt") {
+  randStr = paste(prefix, stri_rand_strings(1, 9, pattern = "[A-Za-z0-9]"), sep =
+                    "_")
+  return(paste(randStr, file_ext, sep = ""))
 }
 
 # this should be read from the json file assigned to variables
@@ -43,11 +44,11 @@ genRandomFileName <- function(prefix,file_ext=".txt"){
 
 
 
-# should put constraints for items data required for this 
+# should put constraints for items data required for this
 crosstype <- json_data$crosstype
 geno_file <- json_data$geno_file
 pheno_file <- json_data$pheno_file
-geno_map_file <- json_data$geno_map_file 
+geno_map_file <- json_data$geno_map_file
 pheno_covar_file <- json_data$phenocovar_file
 alleles <- json_data$alleles
 founder_geno_file = json_data$founder_geno_file
@@ -62,53 +63,58 @@ gmap_file = json_data$gmap_file
 # parsing the required data for example the geno_codes
 
 
-# geno_codes  handle the geno codes here 
+# geno_codes  handle the geno codes here
 
 # make assertion for the geno_file and the geno_file exists
 # make assertion for the physical map file or the geno map file exists
 
-# create temp directory for this workspace 
-control_file_path  <- file.path("/home/kabui", genRandomFileName(prefix="control_", file_ext=".json"))
+# create temp directory for this workspace
+control_file_path  <- file.path("/home/kabui",
+                                genRandomFileName(prefix = "control_", file_ext = ".json"))
 
-str_glue(
-  "Generated control file path is  {control_file_path}"
-)
+str_glue("Generated control file path is  {control_file_path}")
 
 
 
 # create the cross file here from the arguments provided
-# todo add more validation checks here 
+# todo add more validation checks here
 
-# issue I can no define the different paths for files for example pheno_file 
+# issue I can no define the different paths for files for example pheno_file
 
 # think about the issue about geno codes ~~~~
 # function to generate a cross file from a json list
 
-generate_cross_object  <- function(json_data){
-return (write_control_file(control_file_path,
-    crosstype= json_data$crosstype,
-    geno_file= json_data$geno_file,
-    pheno_file= json_data$pheno_file,    
-    gmap_file= json_data$geno_map_file,
-    phenocovar_file= json_data$phenocovar_file,
-    geno_codes= json_data$geno_codes,
-    alleles= json_data$alleles,    
-    na.strings=json_data$na.strings,
-    overwrite = TRUE))
+generate_cross_object  <- function(json_data) {
+  return (
+    write_control_file(
+      control_file_path,
+      crosstype = json_data$crosstype,
+      geno_file = json_data$geno_file,
+      pheno_file = json_data$pheno_file,
+      gmap_file = json_data$geno_map_file,
+      phenocovar_file = json_data$phenocovar_file,
+      geno_codes = json_data$geno_codes,
+      alleles = json_data$alleles,
+      na.strings = json_data$na.strings,
+      overwrite = TRUE
+    )
+  )
 }
 
 
 # alternatively pass a  yaml file with
-dataset <- write_control_file(control_file_path,
-    crosstype= crosstype,
-    geno_file= geno_file,
-    pheno_file= pheno_file,    
-    gmap_file= geno_map_file,
-    phenocovar_file= pheno_covar_file,
-    geno_codes=c(L=1L, C=2L),
-    alleles= alleles,    
-    na.strings=c("-", "NA"),
-    overwrite = TRUE)
+dataset <- write_control_file(
+  control_file_path,
+  crosstype = crosstype,
+  geno_file = geno_file,
+  pheno_file = pheno_file,
+  gmap_file = geno_map_file,
+  phenocovar_file = pheno_covar_file,
+  geno_codes = c(L = 1L, C = 2L),
+  alleles = alleles,
+  na.strings = c("-", "NA"),
+  overwrite = TRUE
+)
 
 # make validation for the data
 dataset  <- read_cross2(control_file_path, quiet = FALSE) # replace this with a dynamic path
@@ -116,13 +122,13 @@ dataset  <- read_cross2(control_file_path, quiet = FALSE) # replace this with a
 # check integrity of the cross
 cat("Check the integrity of the cross object")
 check_cross2(dataset)
-if (check_cross2(dataset)){
+if (check_cross2(dataset)) {
   print("Dataset meets required specifications for a cross")
 } else {
   print("Dataset does not meet required specifications")
 }
 
-# Dataset Summarys 
+# Dataset Summarys
 cat("A Summary about the Dataset You Provided\n")
 summary(dataset)
 n_ind(dataset)
@@ -144,12 +150,18 @@ founders(dataset)
 analysis_type <- "single"
 
 
-perform_genetic_pr <- function(cross, cores= 1, error_prob=0.002, analysis_type="single"){
-# improve on this 
-if (analysis_type == "single"){
-    pr <- calc_genoprob(cross, error_prob= error_prob, quiet=FALSE, cores=cores)
+perform_genetic_pr <- function(cross,
+                               cores = 1,
+                               error_prob = 0.002,
+                               analysis_type = "single") {
+  # improve on this
+  if (analysis_type == "single") {
+    pr <- calc_genoprob(cross,
+                        error_prob = error_prob,
+                        quiet = FALSE,
+                        cores = cores)
     return (pr)
-}
+  }
 }
 
 # get the genetic probability
@@ -178,52 +190,89 @@ print(pheno)
 print(covar)
 print(Xcovar)
 
-# rework on fetching th Xcovar and getting the covar data 
+# rework on fetching th Xcovar and getting the covar data
 
 # perform kinship
 
 
-perform_genome_scan <- function(cross, genome_prob, method, covar =NULL, xCovar=NULL)
-# perform scan1 using haley-knott regression or linear model, or LOCO linear model
+perform_genome_scan <- function(cross,
+                                genome_prob,
+                                method,
+                                covar = NULL,
+                                xCovar = NULL)
+  # perform scan1 using haley-knott regression or linear model, or LOCO linear model
 {
-if (method == "LMM") {
-   # provide parameters for this
-   kinship = calc_kinship(genome_prob)
-   out  <- scan1(genome_prob, cross$pheno, kinship=kinship, addcovar=covar, Xcovar=Xcovar)    
+  if (method == "LMM") {
+    # provide parameters for this
+    kinship = calc_kinship(genome_prob)
+    out  <- scan1(
+      genome_prob,
+      cross$pheno,
+      kinship = kinship,
+      addcovar = covar,
+      Xcovar = Xcovar
+    )
+  }
+  
+  if (method == "LOCO") {
+    # perform kinship inside better option
+    kinship = calc_kinship(genome_prob, "loco")
+    out <- scan1(
+      genome_prob,
+      cross$pheno,
+      kinship = kinship,
+      addcovar = covar,
+      Xcovar = Xcovar
+    )
+  }
+  else {
+    # perform using Haley Knott
+    out <- scan1(genome_prob,
+                 cross$pheno,
+                 addcovar = NULL,
+                 Xcovar = Xcovar)
+  }
+  return (out)
 }
 
-if (method == "LOCO") {
-# perform kinship inside better option
-  kinship = calc_kinship(genome_prob, "loco")
-  out <- scan1(genome_prob, cross$pheno, kinship=kinship,addcovar=covar, Xcovar=Xcovar)
-}
-else {
-# perform using Haley Knott 
-out <- scan1(genome_prob, cross$pheno, addcovar=NULL, Xcovar=Xcovar)
-}
-return (out)
-}
-
-results <- perform_genome_scan(cross=dataset, genome_prob=Pr, method = "HMM")
+results <- perform_genome_scan(cross = dataset,
+                               genome_prob = Pr,
+                               method = "HMM")
 
 
-results # this should probably return the method use here 
+results # this should probably return the method use here
 
 # plot for the LOD scores from  performing the genome scan
 
-generate_lod_plot <- function(cross, scan_result, method, base_dir="."){
-# Plot LOD curves for a genome scan
-color <- c("slateblue", "violetred", "green3")
-par(mar=c(4.1, 4.1, 1.6, 1.1))
-ymx <- maxlod(scan_result)
-file_name = genRandomFileName(prefix="RQTL_LOD_SCORE_",file_ext=".png")
-image_loc = file.path(base_dir ,file_name)
-png(image_loc, width=1000, height=600, type='cairo-png')
-plot(scan_result, cross$gmap, lodcolumn=1, col=color[1], main=colnames(cross$pheno)[1],
-              ylim=c(0, ymx*1.02))
-legend("topleft", lwd=2, col=color[1], method, bg="gray90", lty=c(1,1,2))
-dev.off()
-return (image_loc)
+generate_lod_plot <- function(cross, scan_result, method, base_dir = ".") {
+  # Plot LOD curves for a genome scan
+  color <- c("slateblue", "violetred", "green3")
+  par(mar = c(4.1, 4.1, 1.6, 1.1))
+  ymx <- maxlod(scan_result)
+  file_name = genRandomFileName(prefix = "RQTL_LOD_SCORE_", file_ext = ".png")
+  image_loc = file.path(base_dir , file_name)
+  png(image_loc,
+      width = 1000,
+      height = 600,
+      type = 'cairo-png')
+  plot(
+    scan_result,
+    cross$gmap,
+    lodcolumn = 1,
+    col = color[1],
+    main = colnames(cross$pheno)[1],
+    ylim = c(0, ymx * 1.02)
+  )
+  legend(
+    "topleft",
+    lwd = 2,
+    col = color[1],
+    method,
+    bg = "gray90",
+    lty = c(1, 1, 2)
+  )
+  dev.off()
+  return (image_loc)
 }
 
 lod_file_path <- generate_lod_plot(dataset, results, "HK")
@@ -234,31 +283,51 @@ lod_file_path
 # Q  how do we perform geno scan with Genome scan with a single-QTL model ????
 
 # perform  permutation tests for single-QTL method
- 
 
-perform_permutation_test <- function(cross, genome_prob, n_perm, method="HKK", covar=NULL, Xcovar=NULL, perm_strata=NULL){
-# todo add chr_lengths and perm_Xsp
 
-if (method == "HKK") {
-    perm <- scan1perm(genome_prob, cross$pheno, Xcovar=Xcovar, n_perm= n_perm, perm_strata=perm_strata)
-}
-else if (method == "LMM") {
-   kinship = calc_kinship(genome_prob)
-   perm <- scan1perm(genome_prob, cross$pheno,
-                    kinship=kinship, Xcovar=Xcovar,
-		    n_perm=n_perm)
-}
-else if (method == "LOCO") {
- kinship = calc_kinship(genome_prob, "loco")
-operm3 <- scan1perm(genome_prob, cross$pheno,
-                   kinship=kinship , perm_strata=perm_strata,
-                   Xcovar= Xcovar, n_perm=n_perm)
-}
-return (perm) 
+perform_permutation_test <- function(cross,
+                                     genome_prob,
+                                     n_perm,
+                                     method = "HKK",
+                                     covar = NULL,
+                                     Xcovar = NULL,
+                                     perm_strata = NULL) {
+  # todo add chr_lengths and perm_Xsp
+  
+  if (method == "HKK") {
+    perm <- scan1perm(
+      genome_prob,
+      cross$pheno,
+      Xcovar = Xcovar,
+      n_perm = n_perm,
+      perm_strata = perm_strata
+    )
+  }
+  else if (method == "LMM") {
+    kinship = calc_kinship(genome_prob)
+    perm <- scan1perm(
+      genome_prob,
+      cross$pheno,
+      kinship = kinship,
+      Xcovar = Xcovar,
+      n_perm = n_perm
+    )
+  }
+  else if (method == "LOCO") {
+    kinship = calc_kinship(genome_prob, "loco")
+    operm3 <- scan1perm(
+      genome_prob,
+      cross$pheno,
+      kinship = kinship ,
+      perm_strata = perm_strata,
+      Xcovar = Xcovar,
+      n_perm = n_perm
+    )
+  }
+  return (perm)
 }
 
 # TODO ! get these parameters from argument from the user
-perm <- perform_permutation_test(dataset,Pr, n_perm=2, method="LMM")
+perm <- perform_permutation_test(dataset, Pr, n_perm = 2, method = "LMM")
 # get the permutation summary with a significance threshold
-summary(perm, alpha=c(0.2, 0.05))
-
+summary(perm, alpha = c(0.2, 0.05))
\ No newline at end of file