diff options
| -rw-r--r-- | scripts/rqtl2_wrapper.R | 628 |
1 files changed, 195 insertions, 433 deletions
diff --git a/scripts/rqtl2_wrapper.R b/scripts/rqtl2_wrapper.R index 275370b..e9655d0 100644 --- a/scripts/rqtl2_wrapper.R +++ b/scripts/rqtl2_wrapper.R @@ -1,507 +1,269 @@ -# This file contains an implementation of qtl mapping using r-qtl2 -# For r-qtl1 implementation see: ./rqtl_wrapper.R - library(qtl2) library(rjson) library(stringi) library(optparse) - +# Define command-line options option_list <- list( - make_option(c("-d", "--directory"), action = "store", default = NULL, type = "character", help="Temporary working directory: should also host the input file ."), -make_option(c("-i", "--input_file"), action="store", default=NULL, type='character', help="a yaml or json file with required data to create the cross file"), - make_option(c("-o", "--output_file"), action="store", default=NULL, type='character', help="a file path of where to write the output json results"), - make_option(c("-c", "--cores"), type="integer", default=1, help="No of cores to use while making - computation"), - make_option(c("-p", "--nperm"), type="integer", default=0, action="store_true", help="No of permutations "), - make_option(c("-m", "--method"), action = "store", default = "HK", type = "character", help="Scan Mapping Method - HK (Haley Knott), LMM( Linear Mixed Model ), LOCO (Leave one Chromosome Out)"), - make_option(c("--pstrata"), action="store_true", default=NULL, help="Use permutation strata"), - make_option(c("-t", "--threshold"), type="integer", default= 1, action="store_true", help="Minimum LOD score for a Peak") + make_option(c("-d", "--directory"), action = "store", default = NULL, type = "character", + help = "Temporary working directory: should also host the input file."), + make_option(c("-i", "--input_file"), action = "store", default = NULL, type = 'character', + help = "A YAML or JSON file with required data to create the cross file."), + make_option(c("-o", "--output_file"), action = "store", default = NULL, type = 'character', + help = "A file path of where to write the output JSON results."), + make_option(c("-c", "--cores"), type = "integer", default = 1, + help = "Number of cores to use while making computation."), + make_option(c("-p", "--nperm"), type = "integer", default = 0, + help = "Number of permutations."), + make_option(c("-m", "--method"), action = "store", default = "HK", type = "character", + help = "Scan Mapping Method - HK (Haley Knott), LMM (Linear Mixed Model), LOCO (Leave One Chromosome Out)."), + make_option(c("--pstrata"), action = "store_true", default = NULL, + help = "Use permutation strata."), + make_option(c("-t", "--threshold"), type = "integer", default = 1, + help = "Minimum LOD score for a Peak.") ) -opt_parser = OptionParser(option_list=option_list); +# Parse command-line arguments +opt_parser <- OptionParser(option_list = option_list) opt <- parse_args(opt_parser) -NO_OF_CORES = opt$cores -SCAN_METHOD = opt$method -NO_OF_PERMUTATION = opt$nperm -# Step: check for mandatory file paths -# NOTE this is the working dir where the cross file will be generated -# NOTE this is where the cross file is generated +# Assign parsed arguments to variables +NO_OF_CORES <- opt$cores +SCAN_METHOD <- opt$method +NO_OF_PERMUTATION <- opt$nperm -if(is.null(opt$directory) || !(dir.exists(opt$directory))){ -# check if directory exists - print_help(opt_parser) - stop("The working directory does not exists or is NULL\n") +# Validate input and output file paths +validate_file_paths <- function(opt) { + if (is.null(opt$directory) || !dir.exists(opt$directory)) { + print_help(opt_parser) + stop("The working directory does not exist or is NULL.\n") + } + + INPUT_FILE_PATH <- opt$input_file + OUTPUT_FILE_PATH <- opt$output_file + + if (!file.exists(INPUT_FILE_PATH)) { + print_help(opt_parser) + stop("The input file ", INPUT_FILE_PATH, " you provided does not exist.\n") + } else { + cat("Input file exists. Reading the input file...\n") + } + + if (!file.exists(OUTPUT_FILE_PATH)) { + print_help(opt_parser) + stop("The output file ", OUTPUT_FILE_PATH, " you provided does not exist.\n") + } else { + cat("Output file exists...", OUTPUT_FILE_PATH, "\n") + } + + return(list(input = INPUT_FILE_PATH, output = OUTPUT_FILE_PATH)) } -INPUT_FILE_PATH = opt$input_file -OUTPUT_FILE_PATH = opt$output_file -if (!(file.exists(INPUT_FILE_PATH))) { - print_help(opt_parser) - stop("The input file", INPUT_FILE_PATH, " you provided does not exists\n") -} else { - cat("Input file exists Reading the input file .... \n") -} -if (!(file.exists(OUTPUT_FILE_PATH))) { - print_help(opt_parser) - stop("The output file ",OUTPUT_FILE_PATH, " you provided does not exists\n") -} else { - cat("Output file exists ...", OUTPUT_FILE_PATH, "\n") -} +file_paths <- validate_file_paths(opt) +INPUT_FILE_PATH <- file_paths$input +OUTPUT_FILE_PATH <- file_paths$output -# Utility function to generate random file names of size n: -genRandomFileName <- function(prefix, string_size = 9 , file_ext = ".txt") { - randStr = paste(prefix, stri_rand_strings(1, string_size, pattern = "[A-Za-z0-9]"), sep = - "_") +# Utility function to generate random file names +genRandomFileName <- function(prefix, string_size = 9, file_ext = ".txt") { + randStr <- paste(prefix, stri_rand_strings(1, string_size, pattern = "[A-Za-z0-9]"), sep = "_") return(paste(randStr, file_ext, sep = "")) } +# Generate control file path +control_file_path <- file.path(opt$directory, genRandomFileName(prefix = "control", file_ext = ".json")) +cat("Generated the control file path at", control_file_path, "\n") -# Step: Generate the control file name - -control_file_path <- file.path(opt$directory, - genRandomFileName(prefix = "control", file_ext = ".json")) -cat("Generated the control file path at ", control_file_path, "\n") +# Read and parse the input file +cat("Reading and parsing the input file.\n") +json_data <- fromJSON(file = INPUT_FILE_PATH) - -# Step: Reading and Parsing the input file -cat("Reading and parsing the input file \n") -json_data = fromJSON(file = INPUT_FILE_PATH) -if (is.null(json_data$sep)){ -cat("Using ',' as a default sep for cross file\n") -json_data$sep = "," -} -if (is.null(json_data$na.strings)){ -cat("Using '-' and 'NA' as the default na.strings\n") - json_data$na.strings = c("-" , "NA") -} -default_keys = c( - "geno_transposed", "founder_geno_transposed", - "pheno_transposed" , "covar_transposed", - "phenocovar_transposed") - -for (item in default_keys) { -if (!(item %in% names(json_data))){ - cat("Using FALSE as default parameter for ", item, "\n") - json_data[item] = FALSE -} +# Set default values for JSON data +set_default_values <- function(json_data) { + if (is.null(json_data$sep)) { + cat("Using ',' as a default separator for cross file.\n") + json_data$sep <- "," + } + if (is.null(json_data$na.strings)) { + cat("Using '-' and 'NA' as the default na.strings.\n") + json_data$na.strings <- c("-", "NA") + } + + default_keys <- c("geno_transposed", "founder_geno_transposed", "pheno_transposed", + "covar_transposed", "phenocovar_transposed") + + for (item in default_keys) { + if (!(item %in% names(json_data))) { + cat("Using FALSE as default parameter for", item, "\n") + json_data[item] <- FALSE + } + } + + return(json_data) } - - -# Note the files below should be in the same directory as the location for the crosss file -generate_cross_object <- function(control_file_path, json_data) { - # function to write the cross object from a json data object - 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, - pmap_file = json_data$physical_map_file, - phenocovar_file = json_data$phenocovar_file, - geno_codes = json_data$geno_codes, - alleles = json_data$alleles, - na.strings = json_data$na.strings, - geno_transposed = json_data$geno_transposed, - sex_file = json_data$sex_file, - founder_geno_file = json_data$founder_geno_file, - covar_file = json_data$covar_file, - sex_covar = json_data$sex_covar, - sex_codes = json_data$sex_codes, - crossinfo_file = json_data$crossinfo_file, - crossinfo_covar = json_data$crossinfo_covar, - crossinfo_codes = json_data$crossinfo_codes, - xchr = json_data$xchr, - overwrite = TRUE - ) +json_data <- set_default_values(json_data) + +# Function to generate the cross object +generate_cross_object <- function(control_file_path, json_data) { + 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, + pmap_file = json_data$physical_map_file, + phenocovar_file = json_data$phenocovar_file, + geno_codes = json_data$geno_codes, + alleles = json_data$alleles, + na.strings = json_data$na.strings, + geno_transposed = json_data$geno_transposed, + sex_file = json_data$sex_file, + founder_geno_file = json_data$founder_geno_file, + covar_file = json_data$covar_file, + sex_covar = json_data$sex_covar, + sex_codes = json_data$sex_codes, + crossinfo_file = json_data$crossinfo_file, + crossinfo_covar = json_data$crossinfo_covar, + crossinfo_codes = json_data$crossinfo_codes, + xchr = json_data$xchr, + overwrite = TRUE ) } -# Step: generate the cross file -cat("Generating the cross object at ", control_file_path, "\n") +# Generate the cross object +cat("Generating the cross object at", control_file_path, "\n") generate_cross_object(control_file_path, json_data) -cat("reading the cross object from", control_file_path, "\n") -cross <- read_cross2(control_file_path, quiet = FALSE) # replace this with a dynamic path -# check integrity of the cross -cat("Check the integrity of the cross object") -check_cross2(cross) +# Read the cross object +cat("Reading the cross object from", control_file_path, "\n") +cross <- read_cross2(control_file_path, quiet = FALSE) + +# Check the integrity of the cross object +cat("Checking the integrity of the cross object.\n") if (check_cross2(cross)) { - print("Cross meets required specifications for a cross") + cat("Cross meets required specifications for a cross.\n") } else { - print("Cross does not meet required specifications") + cat("Cross does not meet required specifications.\n") } - -# Cross Summarys -cat("A Summary about the Cross You Provided\n") +# Print cross summary +cat("A summary about the cross you provided:\n") summary(cross) -n_ind(cross) -n_chr(cross) -cat("names of markers in the object\n") -marker_names(cross) -cat("names of phenotypes in a the object") -pheno_names(cross) -cat("IDs for all individuals in the cross cross object that have genotype data\n") -ind_ids_geno(cross) -cat(" IDs for all individuals in the cross object that have phenotype data") -ind_ids_pheno(cross) -cat("Name of the founder Strains/n") -founders(cross) - -# Function for computing the genetic probabilities -perform_genetic_pr <- function(cross, - cores = NO_OF_CORES, - step=1, - map=NULL, - map_function=c("haldane", "kosambi", "c-f", "morgan"), - error_prob = 0.002 - ) { - #' Function to calculate the genetic probabilities - #' @description function to perform genetic probabilities - #' @param cores number no of cores to use Defaults to "1" - #' @param map Genetic map of markers. defaults to "NONE" - #' @param use_pseudomarkers option to insert pseudo markers in the gmap default "FALSE" - #' @param error_prob - #' @param map_function Character string indicating the map function to use to convert genetic - #' @param step for default "1" - #' @return a list of three-dimensional arrays of probabilities, individuals x genotypes x pst - - return(calc_genoprob(cross, map=map, - error_prob=error_prob, map_function=map_function, - quiet=FALSE, cores=cores)) +# Function to compute genetic probabilities +perform_genetic_pr <- function(cross, cores = NO_OF_CORES, step = 1, map = NULL, + map_function = c("haldane", "kosambi", "c-f", "morgan"), + error_prob = 0.002) { + calc_genoprob(cross, map = map, error_prob = error_prob, map_function = map_function, + quiet = FALSE, cores = cores) } +# Insert pseudomarkers to the genetic map +cat("Inserting pseudomarkers to the genetic map with step 1 and stepwidth fixed.\n") +MAP <- insert_pseudomarkers(cross$gmap, step = 1, stepwidth = "fixed", cores = NO_OF_CORES) -#Step: insert pseudomarkers to genetic map -# TODO need to review this to match rqtl1 -cat("Inserting pseudomarkers to the genetic map with steps", 1 , "and stepwidth" , "fixed\n") -MAP <- insert_pseudomarkers(cross$gmap, step= 1, stepwidth = "fixed", cores = NO_OF_CORES) - +# Calculate genetic probabilities +cat("Calculating the genetic probabilities.\n") +Pr <- perform_genetic_pr(cross) -# Step: calculate the genetic probabilities -cat("Calculating the genetic probabilities\n") -Pr = perform_genetic_pr(cross) - - - - -# Step: perform allele probabilites if cross ways -if (cross$crosstype == "4way"){ - cat("Calculating Allele Genetic probability for 4way cross\n") - aPr <- genoprob_to_alleleprob(pr) +# Calculate allele probabilities for 4-way cross +if (cross$crosstype == "4way") { + cat("Calculating allele genetic probability for 4-way cross.\n") + aPr <- genoprob_to_alleleprob(Pr) } -#Function to Calculate genotyping error LOD scores -cat("Calculating the genotype error LOD scores\n") +# Calculate genotyping error LOD scores +cat("Calculating the genotype error LOD scores.\n") error_lod <- calc_errorlod(cross, Pr, quiet = FALSE, cores = NO_OF_CORES) -# combine into one matrix error_lod <- do.call("cbind", error_lod) - - - -## grab phenotypes and covariates; ensure that covariates have names attribute -# TODO rework on this -cat("Getting the phenotypes and covariates\n") +# Get phenotypes and covariates +cat("Getting the phenotypes and covariates.\n") pheno <- cross$pheno covar <- match(cross$covar$sex, c("f", "m")) # make numeric -covar - -if (!is.null(covar)){ - names(covar) <- rownames(cross$covar) - +if (!is.null(covar)) { + names(covar) <- rownames(cross$covar) } -print("The covariates are") -covar Xcovar <- get_x_covar(cross) -print("The Xcovar are ") +cat("The covariates are:\n") +print(covar) +cat("The Xcovar are:\n") print(Xcovar) -# Function to calculate the kinship -get_kinship <- function(probability, method="LMM"){ -if (opt$method == "LMM"){ - kinship = calc_kinship(probability) -} else if (opt$method == "LOCO"){ - kinship = calc_kinship(probability, "loco") -}else { - kinship = NULL -} +# Function to calculate kinship +get_kinship <- function(probability, method = "LMM") { + if (method == "LMM") { + kinship <- calc_kinship(probability) + } else if (method == "LOCO") { + kinship <- calc_kinship(probability, "loco") + } else { + kinship <- NULL + } + return(kinship) } - -cat("Calculating the kinship for the genetic probability\n") -if (cross$crosstype == "4way"){ +# Calculate kinship for the genetic probability +cat("Calculating the kinship for the genetic probability.\n") +if (cross$crosstype == "4way") { kinship <- get_kinship(aPr, opt$method) } else { - kinship <- get_kinship(Pr, "loco") + kinship <- get_kinship(Pr, "loco") } - - # Function to perform genome scan -perform_genome_scan <- function(cross, - genome_prob, - method, - addcovar = NULL, - intcovar = NULL, - kinship = NULL, - model = c("normal","binary"), - Xcovar = NULL) { - #' perform genome scan - #' @description perform scan1 using haley-knott regression, perform scan1 using haley-knott #' or linear model, or LOCO linear model - #' the cross object required to pull the pheno - #' @param method to method to perform scan1 either by haley-knott regression(HL), - #' linear mixed model(LMM) or , for the LOCO method(LOCO) - #' @param intcovar A numeric optional matrix of interactive covariates. - #' @param addcovar An optional numeric matrix of additive covariates. - #' @param Xcovar An optional numeric matrix with additional additive covariates used for null #' used for null hypothesis when scanning the X chromosome. - #' @param model Indicates whether to use a normal model (least squares) or binary model - #' @return An object of class "scan1" - - +perform_genome_scan <- function(cross, genome_prob, method, addcovar = NULL, intcovar = NULL, + kinship = NULL, model = c("normal", "binary"), Xcovar = NULL) { if (method == "LMM") { - # provide parameters for this - cat("Performing scan1 using Linear mixed model\n") - - out <- scan1( - genome_prob, - cross$pheno, - kinship = kinship, - # addcovar = covar, - # Xcovar = Xcovar, - # intcovar = intcovar, - model = model, - cores = NO_OF_CORES - ) - } else if (method == "LOCO") { - cat("Performing scan1 using Leave one chromosome out\n") - out <- scan1( - genome_prob, - cross$pheno, - kinship = kinship, - # addcovar = covar, - # intcovar = intcovar, - model = model, - # Xcovar = Xcovar, - cores = NO_OF_CORES - ) - } - - else if (method == "HK"){ - cat("Performing scan1 using Haley Knott\n") - out <- scan1(genome_prob, - cross$pheno, - addcovar = NULL, - intcovar = intcovar, - model = model, - Xcovar = Xcovar, - cores = NO_OF_CORES - ) + cat("Performing scan1 using Linear Mixed Model.\n") + out <- scan1(genome_prob, cross$pheno, kinship = kinship, model = model, cores = NO_OF_CORES) + } else if (method == "LOCO") { + cat("Performing scan1 using Leave One Chromosome Out.\n") + out <- scan1(genome_prob, cross$pheno, kinship = kinship, model = model, cores = NO_OF_CORES) + } else if (method == "HK") { + cat("Performing scan1 using Haley Knott.\n") + out <- scan1(genome_prob, cross$pheno, addcovar = addcovar, intcovar = intcovar, + model = model, Xcovar = Xcovar, cores = NO_OF_CORES) } - - - return (out) + return(out) } # Perform the genome scan for the cross object -if (cross$crosstype == "4way"){ - sex <- (cross$covar$Sex == "male")*1 - names(sex) <- rownames(cross$covar) - sex <- setNames( (cross$covar$Sex == "male")*1, rownames(cross$covar)) - scan_results <- perform_genome_scan(aPr, cross, kinship=kinship, method = "LOCO", addcovar = sex) +if (cross$crosstype == "4way") { + sex <- setNames((cross$covar$Sex == "male") * 1, rownames(cross$covar)) + scan_results <- perform_genome_scan(aPr, cross, kinship = kinship, method = "LOCO", addcovar = sex) } else { - scan_results <- perform_genome_scan(cross = cross, - genome_prob = Pr, - kinship = kinship, - method = SCAN_METHOD) + scan_results <- perform_genome_scan(cross = cross, genome_prob = Pr, kinship = kinship, + method = SCAN_METHOD) } - - +# Save scan results scan_file <- file.path(opt$directory, "scan_results.csv") write.csv(scan_results, scan_file) -# function: perform permutation tests for single-QTL method -perform_permutation_test <- function(cross, - genome_prob, - n_perm, - method = opt$method, - covar = NULL, - Xcovar = NULL, - addcovar = NULL, - intcovar = NULL, - perm_Xsp = FALSE, - kinship = NULL, - model = c("normal", "binary"), - chr_lengths = NULL, +# Function to perform permutation tests +perform_permutation_test <- function(cross, genome_prob, n_perm, method = opt$method, + covar = NULL, Xcovar = NULL, addcovar = NULL, + intcovar = NULL, perm_Xsp = FALSE, kinship = NULL, + model = c("normal", "binary"), chr_lengths = NULL, perm_strata = NULL) { - - #' Function to peform permutation tests for single QTL method - #' @description The scan1perm() function takes the - #' same arguments as scan1(), plus additional a #rguments to control the permutations: - #' @param cross the cross object required to fetch the pheno - #' @param genome_prob the genomic probability matrix - #' @param method to computation method used to perform the genomic scan - #' @param intcovar - #' @param addcovar - #' @param Xcovar - #' @param perm_Xsp If TRUE, do separate permutations for the autosomes and the X chromosome. - #' @param perm_strata Vector of strata, for a stratified permutation test. - #' @param n_perm Number of permutation replicates. - #' @param chr_lengths engths of the chromosomes; - #' @return object of class "scan1perm". - return (scan1perm( - genome_prob, - cross$pheno, - kinship = kinship, - Xcovar = Xcovar, - intcovar = intcovar, - addcovar = addcovar, - n_perm = n_perm, - perm_Xsp = perm_Xsp, - model = model, - chr_lengths = chr_lengths, - cores = NO_OF_CORES - )) + scan1perm(genome_prob, cross$pheno, kinship = kinship, Xcovar = Xcovar, intcovar = intcovar, + addcovar = addcovar, n_perm = n_perm, perm_Xsp = perm_Xsp, model = model, + chr_lengths = chr_lengths, cores = NO_OF_CORES) } - - - -# check if pstrata -if (!(is.null(opt$pstrata)) && (!is.null(Xcovar))){ -perm_strata <- mat2strata(Xcovar) +# Check if permutation strata is needed +if (!is.null(opt$pstrata) && !is.null(Xcovar)) { + perm_strata <- mat2strata(Xcovar) } else { -perm_strata <- NULL + perm_strata <- NULL } +# Perform permutation test if requested +permutation_results_file <- file.path(opt$directory, "permutation.csv") +significance_results_file <- file.path(opt$directory, "significance.csv") -permutation_results_file = file.path(opt$directory, "permutation.csv") -significance_results_file = file.path(opt$directory, "significance.csv") - -# perform permutation if (NO_OF_PERMUTATION > 0) { - cat("performing permutation test for the cross object with permutations ", NO_OF_PERMUTATION, "\n") - perm <- perform_permutation_test( - cross, Pr, - n_perm = NO_OF_PERMUTATION, - perm_strata = perm_strata, - method = opt$method - ) - # Function to get LOD significance thresholds - get_lod_significance <- function(perm, thresholds = c(0.01, 0.05, 0.63)) { - cat("Getting the permutation summary with significance thresholds:", thresholds, "\n") - summary(perm, alpha = thresholds) - } - - # Compute LOD significance - lod_significance <- get_lod_significance(perm) - - # Save results - write.csv(lod_significance, significance_results_file) - write.csv(perm, permutation_results_file) -} - - - - -# step: get the lod peaks -# TODO fix the threshold here -cat("Fetching the lod peaks with threshold", opt$threshold, "\n") -lod_peaks = find_peaks( - scan_results, - threshold =opt$threshold, - map = cross$gmap, - cores = NO_OF_CORES -) - - -# step: get the estimated qtl effects -get_qtl_effect <- function(chromosome,geno_prob,pheno,covar=NULL,LOCO= NULL){ - cat("Finding the qtl effect\n") - chr_Pr <- geno_prob[,chromosome] - if (!is.null(chr_Pr)){ - cat("Finding qtl effect for chromosome ", chromosome, "\n") - if (!is.null(LOCO)) { - cat("Finding qtl effect for chromosome ", chromosome, "with LOCO \n") - kinship <- calc_kinship(chr_Pr, "loco")[[chromosome]] - return(scan1coef(chr_Pr, pheno, kinship, addcovar=covar)) - } - else { - return(scan1coef(chr_Pr, pheno, addcovar=covar)) - } - } - return(NULL) -} - - - -# take the first phenotype in the cross -# grab phenotypes and covariates; ensure that covariates have names attribute - -pheno <- cross$pheno[,1] -if (!is.null(cross$covar) && !is.null(cross$covar$sex)){ - covar <- match(cross$covar$sex, c("f", "m")) # make numeric - names(covar) <- rownames(cross$covar) -} else { -covar <- NULL -} - -meffects <- c() -meffects_plots <- c() -# TODO add plots for meffects -for (chr in chr_names(cross)){ - cat("Getting the qtl effect for chromosome", chr, "\n") - if (cross$crosstype == "4way"){ - coeff_results <- get_qtl_effect(chr, aPr, pheno, LOCO="LOCO", covar = sex) - cat("Generating the qtl effects plots\n") - file_name = genRandomFileName(prefix = "RQTL_EFFECT_", file_ext = ".png") - image_loc = file.path(base_dir , file_name) - par(mar=c(4.1, 4.1, 0.6, 0.6)) - png(image_loc, - width = 1000, - height = 600, - type = 'cairo-png') - plot( - coeff_results, - cross$gmap[chr], - bgcolor="gray95", - legend="bottomleft" - ) - meffects <- append(meffects_plots, image_loc) - } else { - coeff_results <- get_qtl_effect(chr, Pr, pheno) - } - meffects <- append(meffects, coeff_results) -} - - - -gmap_file <- file.path(opt$directory, json_data$geno_map_file) -pmap_file <- file.path(opt$directory, json_data$physical_map_file) -output = list(lod_peaks = lod_peaks, - scan_results =scan_results, - genetic_probabilities = Pr, - lod_peaks = lod_peaks, - permutation_file = permutation_results_file, - significance_file = significance_results_file, - scan_file = scan_file, - chromosomes = chr_names(cross), - error_lod = error_lod, - gmap_file = gmap_file, - pmap_file = pmap_file, - meffects_plots = meffects_plots, - permutations = NO_OF_PERMUTATION, - scan_method = SCAN_METHOD - ) -output_json_data <-toJSON(output) -cat("The output file path generated is", OUTPUT_FILE_PATH, "\n") -cat("Writing to the output file\n") -write(output_json_data, file=OUTPUT_FILE_PATH) - + cat("Performing permutation test for the cross object with", NO_OF_PERMUTATION, "permutations.\n") + perm <- perform_permutation_test(cross, Pr, n_perm = NO_OF_PERMUTATION, perm_strata = perm_strata, + method = opt$method) \ No newline at end of file |