diff options
Diffstat (limited to 'scripts')
| -rw-r--r-- | scripts/rqtl2_wrapper.R | 279 |
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 |