aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--scripts/rqtl2_wrapper.R57
1 files changed, 56 insertions, 1 deletions
diff --git a/scripts/rqtl2_wrapper.R b/scripts/rqtl2_wrapper.R
index e9655d0..9d288af 100644
--- a/scripts/rqtl2_wrapper.R
+++ b/scripts/rqtl2_wrapper.R
@@ -266,4 +266,59 @@ 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) \ No newline at end of file
+ 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)
+
+output <- list(
+ permutation_file = permutation_results_file,
+ significance_file = significance_results_file,
+ scan_file = scan_file,
+ gmap_file = gmap_file,
+ pmap_file = pmap_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