1 files changed, 358 insertions, 0 deletions
diff --git a/scripts/rqtl2_wrapper.R b/scripts/rqtl2_wrapper.R
new file mode 100644
index 0000000..af13efa
--- /dev/null
+++ b/scripts/rqtl2_wrapper.R
@@ -0,0 +1,358 @@
+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 = "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.")
+)
+
+# Parse command-line arguments
+opt_parser <- OptionParser(option_list = option_list)
+opt <- parse_args(opt_parser)
+
+# Assign parsed arguments to variables
+NO_OF_CORES <- opt$cores
+SCAN_METHOD <- opt$method
+NO_OF_PERMUTATION <- opt$nperm
+
+NO_OF_CORES <- 20
+
+# 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))
+}
+
+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
+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")
+# Read and parse the input file
+cat("Reading and parsing the input file.\n")
+json_data <- fromJSON(file = INPUT_FILE_PATH)
+
+# 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)
+}
+
+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,
+    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,
+    founder_geno_transposed = json_data$founder_geno_transposed,
+    geno_transposed = json_data$geno_transposed
+  )
+}
+
+# Generate the cross object
+cat("Generating the cross object at", control_file_path, "\n")
+generate_cross_object(control_file_path, json_data)
+
+# 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)) {
+  cat("Cross meets required specifications for a cross.\n")
+} else {
+  cat("Cross does not meet required specifications.\n")
+}
+
+# Print cross summary
+cat("A summary about the cross you provided:\n")
+summary(cross)
+
+# 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)
+
+# Calculate genetic probabilities
+cat("Calculating the genetic probabilities.\n")
+Pr <- perform_genetic_pr(cross)
+
+# 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)
+}
+
+# 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)
+error_lod <- do.call("cbind", error_lod)
+
+# Get phenotypes and covariates
+cat("Getting the phenotypes and covariates.\n")
+pheno <- cross$pheno
+# covar <- match(cross$covar$sex, c("f", "m")) # make numeric
+# TODO rework on this
+covar <- NULL
+if (!is.null(covar)) {
+  names(covar) <- rownames(cross$covar)
+}
+
+Xcovar <- get_x_covar(cross)
+cat("The covariates are:\n")
+print(covar)
+cat("The Xcovar are:\n")
+print(Xcovar)
+
+# 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)
+}
+
+# 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")
+}
+
+# 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) {
+  if (method == "LMM") {
+    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)
+}
+
+# Perform the genome scan for the cross object
+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)
+}
+
+# Save scan results
+scan_file <- file.path(opt$directory, "scan_results.csv")
+write.csv(scan_results, scan_file)
+
+# 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) {
+  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 permutation strata is needed
+if (!is.null(opt$pstrata) && !is.null(Xcovar)) {
+  perm_strata <- mat2strata(Xcovar)
+} else {
+  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")
+
+if (NO_OF_PERMUTATION > 0) {
+  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)
+  
+  # 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)
+}
+
+
+
+# Function to get QTL effects
+get_qtl_effect <- function(chromosome, geno_prob, pheno, covar = NULL, LOCO = NULL) {
+  cat("Finding the QTL effect for chromosome", chromosome, "\n")
+  chr_Pr <- geno_prob[, chromosome]
+  if (!is.null(chr_Pr)) {
+    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)
+}
+
+# Get QTL effects for each chromosome
+# TODO
+
+# Prepare output data
+gmap_file <- file.path(opt$directory, json_data$geno_map_file)
+pmap_file <- file.path(opt$directory, json_data$physical_map_file)
+
+
+
+
+
+# Construct the Map object from cross with columns (Marker, chr, cM, Mb)
+gmap <- cross$gmap  # Genetic map in cM
+pmap <- cross$pmap  # Physical map in Mb
+# Convert lists to data frames
+gmap_df <- data.frame(
+  marker = unlist(lapply(gmap, names)), 
+  chr = rep(names(gmap), sapply(gmap, length)),  # Add chromosome info
+  CM = unlist(gmap), 
+  stringsAsFactors = FALSE
+)
+
+pmap_df <- data.frame(
+  marker = unlist(lapply(pmap, names)), 
+  chr = rep(names(pmap), sapply(pmap, length)),  # Add chromosome info
+  MB = unlist(pmap), 
+  stringsAsFactors = FALSE
+)
+# Merge using full outer join (by marker and chromosome)
+merged_map <- merge(gmap_df, pmap_df, by = c("marker", "chr"), all = TRUE)
+map_file <- file.path(opt$directory, "map.csv")
+write.csv(merged_map, map_file, row.names = FALSE)
+
+output <- list(
+  permutation_file = permutation_results_file,
+  significance_file = significance_results_file,
+  scan_file = scan_file,
+  gmap_file = gmap_file,
+  pmap_file = pmap_file,
+  map_file  = map_file,
+  permutations = NO_OF_PERMUTATION,
+  scan_method = SCAN_METHOD
+)
+
+# Write output to JSON file
+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)
\ No newline at end of file