From e9110f840c8817092d0264990a4c7f6d16e966ba Mon Sep 17 00:00:00 2001
From: zsloan
Date: Thu, 22 Jul 2021 23:47:58 +0000
Subject: Added option for running pairscan to rqtl_wrapper.R

---
 scripts/rqtl_wrapper.R | 1 +
 1 file changed, 1 insertion(+)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index 7518175..db0b83c 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -11,6 +11,7 @@ option_list = list(
   make_option(c("-c", "--addcovar"), action="store_true", default=NULL, help="Use covariates (included as extra columns in the phenotype input file)"),
   make_option(c("--model"), type="character", default="normal", help="Mapping Model - Normal or Non-Parametric"),
   make_option(c("--method"), type="character", default="hk", help="Mapping Method - hk (Haley Knott), ehk (Extended Haley Knott), mr (Marker Regression), em (Expectation-Maximization), imp (Imputation)"),
+  make_option(c("--pairscan"), action="store_true", default=NULL, help="Run Pair Scan - the R/qtl function scantwo"),
   make_option(c("-i", "--interval"), action="store_true", default=NULL, help="Use interval mapping"),
   make_option(c("--nperm"), type="integer", default=0, help="Number of permutations"),
   make_option(c("--pstrata"), action="store_true", default=NULL, help="Use permutation strata (stored as final column/vector in phenotype input file)"),
-- 
cgit v1.2.3


From 14d76308a779f6e4da4c9f4c17f971f5a0c938b0 Mon Sep 17 00:00:00 2001
From: zsloan
Date: Fri, 23 Jul 2021 19:40:04 +0000
Subject: - Added scan_func function that determines whether scanone or scantwo
 (pairscan) is used - For pairscan default to using step 20 (subject to
 change, but some step is required during calc.genoprob to make it run fast
 enough) - Added some new verbose prints

---
 scripts/rqtl_wrapper.R | 41 +++++++++++++++++++++++++++++------------
 1 file changed, 29 insertions(+), 12 deletions(-)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index db0b83c..2641271 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -156,6 +156,9 @@ cross_object = geno_to_csvr(geno_file, trait_names, trait_vals, cross_file)
 if (!is.null(opt$interval)) {
   verbose_print('Calculating genotype probabilities with interval mapping\n')
   cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
+} else if (!is.null(opt$pairscan)) {
+  verbose_print('Calculating genotype probabilities with interval mapping\n')
+  cross_object <- calc.genoprob(cross_object, step=20)
 } else {
   verbose_print('Calculating genotype probabilities\n')
   cross_object <- calc.genoprob(cross_object)
@@ -164,6 +167,7 @@ if (!is.null(opt$interval)) {
 # Pull covariates out of cross object, if they exist
 covars = vector(mode = "list", length = length(trait_names) - 1)
 if (!is.null(opt$addcovar)) {
+  verbose_print('Pulling covariates out of cross object\n')
   #If perm strata are being used, it'll be included as the final column in the phenotype file
   if (!is.null(opt$pstrata)) {
     covar_names = trait_names[3:length(trait_names) - 1]
@@ -176,16 +180,28 @@ if (!is.null(opt$addcovar)) {
 # Pull permutation strata out of cross object, if it is being used
 perm_strata = vector()
 if (!is.null(opt$pstrata)) {
+  verbose_print('Pulling permutation strata out of cross object\n')
   strata_col = trait_names[length(trait_names)]
   perm_strata <- pull.pheno(cross_object, strata_col)
 }
 
 # If a marker name is supplied as covariate, get its vector of values and add them as a covariate
 if (!is.null(opt$control)) {
+  verbose_print('Creating marker covariates and binding them to covariates vector\n')
   marker_covars = create_marker_covars(cross_object, opt$control)
   covars <- cbind(covars, marker_covars)
 }
 
+if (!is.null(opt$pairscan)) {
+  scan_func <- function(...){
+    scantwo(...)
+  }
+} else {
+  scan_func <- function(...){
+    scanone(...)
+  }
+}
+
 # Calculate permutations
 if (opt$nperm > 0) {
   if (!is.null(opt$filename)){
@@ -196,19 +212,19 @@ if (opt$nperm > 0) {
 
   if (!is.null(opt$addcovar) || !is.null(opt$control)){
     if (!is.null(opt$pstrata)) {
-      verbose_print('Running ', opt$nperm, ' permutations with cofactors and strata\n')
-      perm_results = scanone(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
+      verbose_print('Running permutations with cofactors and strata\n')
+      perm_results = scan_func(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
     } else {
-      verbose_print('Running ', opt$nperm, ' permutations with cofactors\n')
-      perm_results = scanone(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, model=opt$model, method=opt$method)
+      verbose_print('Running permutations with cofactors\n')
+      perm_results = scan_func(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, model=opt$model, method=opt$method)
     }
   } else {
     if (!is.null(opt$pstrata)) {
-      verbose_print('Running ', opt$nperm, ' permutations with strata\n')
-      perm_results = scanone(cross_object, pheno.col=1, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
+      verbose_print('Running permutations with strata\n')
+      perm_results = scan_func(cross_object, pheno.col=1, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
     } else {
-      verbose_print('Running ', opt$nperm, ' permutations\n')
-      perm_results = scanone(cross_object, pheno.col=1, n.perm=opt$nperm, model=opt$model, method=opt$method)
+      verbose_print('Running permutations\n')
+      perm_results = scan_func(cross_object, pheno.col=1, n.perm=opt$nperm, model=opt$model, method=opt$method)
     }
   }
   write.csv(perm_results, perm_out_file)
@@ -221,10 +237,11 @@ if (!is.null(opt$filename)){
 }
 
 if (!is.null(opt$addcovar) || !is.null(opt$control)){
-  verbose_print('Running scanone with cofactors\n')
-  qtl_results = scanone(cross_object, pheno.col=1, addcovar=covars, model=opt$model, method=opt$method)
+  verbose_print('Running scan with cofactors\n')
+  qtl_results = scan_func(cross_object, pheno.col=1, addcovar=covars, model=opt$model, method=opt$method)
 } else {
-  verbose_print('Running scanone\n')
-  qtl_results = scanone(cross_object, pheno.col=1, model=opt$model, method=opt$method)
+  verbose_print('Running scan\n')
+  qtl_results = scan_func(cross_object, pheno.col=1, model=opt$model, method=opt$method)
 }
+
 write.csv(qtl_results, out_file)
-- 
cgit v1.2.3


From 91804b4d4f5f82ddad30b68691456885b39309c9 Mon Sep 17 00:00:00 2001
From: zsloan
Date: Fri, 23 Jul 2021 20:08:36 +0000
Subject: Added line priting pair-scan results to CSV and changed the default
 step-size to 10cM for pair-scan

---
 scripts/rqtl_wrapper.R | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index 2641271..ebd5a91 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -158,7 +158,7 @@ if (!is.null(opt$interval)) {
   cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
 } else if (!is.null(opt$pairscan)) {
   verbose_print('Calculating genotype probabilities with interval mapping\n')
-  cross_object <- calc.genoprob(cross_object, step=20)
+  cross_object <- calc.genoprob(cross_object, step=10)
 } else {
   verbose_print('Calculating genotype probabilities\n')
   cross_object <- calc.genoprob(cross_object)
@@ -244,4 +244,9 @@ if (!is.null(opt$addcovar) || !is.null(opt$control)){
   qtl_results = scan_func(cross_object, pheno.col=1, model=opt$model, method=opt$method)
 }
 
-write.csv(qtl_results, out_file)
+verbose_print('Writing results to CSV file\n')
+if (!is.null(opt$pairscan)) {
+  write.csv(qtl_results[1], out_file)
+} else {
+  write.csv(qtl_results, out_file)
+}
-- 
cgit v1.2.3


From 2cb18bb7916e029d3aa7779df710c1d0cbc910f8 Mon Sep 17 00:00:00 2001
From: zsloan
Date: Mon, 2 Aug 2021 21:07:52 +0000
Subject: Updated rqtl_wrapper to also return a map file when doing a pair-scan
 (since we need the list of markers/pseudomarkers and their positions)

---
 scripts/rqtl_wrapper.R | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index ebd5a91..ffff5b9 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -157,7 +157,7 @@ if (!is.null(opt$interval)) {
   verbose_print('Calculating genotype probabilities with interval mapping\n')
   cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
 } else if (!is.null(opt$pairscan)) {
-  verbose_print('Calculating genotype probabilities with interval mapping\n')
+  verbose_print('Calculating genotype probabilities for pair-scan\n')
   cross_object <- calc.genoprob(cross_object, step=10)
 } else {
   verbose_print('Calculating genotype probabilities\n')
@@ -246,7 +246,9 @@ if (!is.null(opt$addcovar) || !is.null(opt$control)){
 
 verbose_print('Writing results to CSV file\n')
 if (!is.null(opt$pairscan)) {
+  map_out_file = file.path(opt$outdir, paste("MAP_", opt$filename, sep = "" ))
   write.csv(qtl_results[1], out_file)
+  write.csv(qtl_results[2], map_out_file)
 } else {
   write.csv(qtl_results, out_file)
 }
-- 
cgit v1.2.3


From e13e10d0f0dc96518e1d5efa9865ac2821b1c3f9 Mon Sep 17 00:00:00 2001
From: zsloan
Date: Thu, 22 Jul 2021 23:47:58 +0000
Subject: Added option for running pairscan to rqtl_wrapper.R

---
 scripts/rqtl_wrapper.R | 1 +
 1 file changed, 1 insertion(+)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index 13c2684..4c96ff2 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -12,6 +12,7 @@ option_list = list(
   make_option(c("--covarstruct"), type="character", help="File detailing which covariates are categorical or numerical"),
   make_option(c("--model"), type="character", default="normal", help="Mapping Model - Normal or Non-Parametric"),
   make_option(c("--method"), type="character", default="hk", help="Mapping Method - hk (Haley Knott), ehk (Extended Haley Knott), mr (Marker Regression), em (Expectation-Maximization), imp (Imputation)"),
+  make_option(c("--pairscan"), action="store_true", default=NULL, help="Run Pair Scan - the R/qtl function scantwo"),
   make_option(c("-i", "--interval"), action="store_true", default=NULL, help="Use interval mapping"),
   make_option(c("--nperm"), type="integer", default=0, help="Number of permutations"),
   make_option(c("--pstrata"), action="store_true", default=NULL, help="Use permutation strata (stored as final column/vector in phenotype input file)"),
-- 
cgit v1.2.3


From 2849792ddc2e76750421c75d926d5a7f18c57e97 Mon Sep 17 00:00:00 2001
From: zsloan
Date: Fri, 23 Jul 2021 19:40:04 +0000
Subject: - Added scan_func function that determines whether scanone or scantwo
 (pairscan) is used - For pairscan default to using step 20 (subject to
 change, but some step is required during calc.genoprob to make it run fast
 enough) - Added some new verbose prints

---
 scripts/rqtl_wrapper.R | 40 ++++++++++++++++++++++++++++------------
 1 file changed, 28 insertions(+), 12 deletions(-)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index 4c96ff2..0d13ccb 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -175,6 +175,9 @@ cross_object = geno_to_csvr(geno_file, trait_names, trait_vals, cross_file, type
 if (!is.null(opt$interval)) {
   verbose_print('Calculating genotype probabilities with interval mapping\n')
   cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
+} else if (!is.null(opt$pairscan)) {
+  verbose_print('Calculating genotype probabilities with interval mapping\n')
+  cross_object <- calc.genoprob(cross_object, step=20)
 } else {
   verbose_print('Calculating genotype probabilities\n')
   cross_object <- calc.genoprob(cross_object)
@@ -189,6 +192,7 @@ if (type == "4-way") {
 # Pull covariates out of cross object, if they exist
 covars <- c() # Holds the covariates which should be passed to R/qtl
 if (!is.null(opt$addcovar)) {
+  verbose_print('Pulling covariates out of cross object\n')
   # If perm strata are being used, it'll be included as the final column in the phenotype file
   if (!is.null(opt$pstrata)) {
     covar_names = trait_names[2:(length(trait_names)-1)]
@@ -224,16 +228,28 @@ if (!is.null(opt$addcovar)) {
 # Pull permutation strata out of cross object, if it is being used
 perm_strata = vector()
 if (!is.null(opt$pstrata)) {
+  verbose_print('Pulling permutation strata out of cross object\n')
   strata_col = trait_names[length(trait_names)]
   perm_strata <- pull.pheno(cross_object, strata_col)
 }
 
 # If a marker name is supplied as covariate, get its vector of values and add them as a covariate
 if (!is.null(opt$control)) {
+  verbose_print('Creating marker covariates and binding them to covariates vector\n')
   marker_covars = create_marker_covars(cross_object, opt$control)
   covars <- cbind(covars, marker_covars)
 }
 
+if (!is.null(opt$pairscan)) {
+  scan_func <- function(...){
+    scantwo(...)
+  }
+} else {
+  scan_func <- function(...){
+    scanone(...)
+  }
+}
+
 # Calculate permutations
 if (opt$nperm > 0) {
   if (!is.null(opt$filename)){
@@ -244,19 +260,19 @@ if (opt$nperm > 0) {
 
   if (!is.null(opt$addcovar) || !is.null(opt$control)){
     if (!is.null(opt$pstrata)) {
-      verbose_print('Running ', opt$nperm, ' permutations with cofactors and strata\n')
-      perm_results = scanone(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
+      verbose_print('Running permutations with cofactors and strata\n')
+      perm_results = scan_func(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
     } else {
-      verbose_print('Running ', opt$nperm, ' permutations with cofactors\n')
-      perm_results = scanone(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, model=opt$model, method=opt$method)
+      verbose_print('Running permutations with cofactors\n')
+      perm_results = scan_func(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, model=opt$model, method=opt$method)
     }
   } else {
     if (!is.null(opt$pstrata)) {
-      verbose_print('Running ', opt$nperm, ' permutations with strata\n')
-      perm_results = scanone(cross_object, pheno.col=1, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
+      verbose_print('Running permutations with strata\n')
+      perm_results = scan_func(cross_object, pheno.col=1, n.perm=opt$nperm, perm.strata=perm_strata, model=opt$model, method=opt$method)
     } else {
-      verbose_print('Running ', opt$nperm, ' permutations\n')
-      perm_results = scanone(cross_object, pheno.col=1, n.perm=opt$nperm, model=opt$model, method=opt$method)
+      verbose_print('Running permutations\n')
+      perm_results = scan_func(cross_object, pheno.col=1, n.perm=opt$nperm, model=opt$model, method=opt$method)
     }
   }
   write.csv(perm_results, perm_out_file)
@@ -269,11 +285,11 @@ if (!is.null(opt$filename)){
 }
 
 if (!is.null(opt$addcovar) || !is.null(opt$control)){
-  verbose_print('Running scanone with cofactors\n')
-  qtl_results = scanone(cross_object, pheno.col=1, addcovar=covars, model=opt$model, method=opt$method)
+  verbose_print('Running scan with cofactors\n')
+  qtl_results = scan_func(cross_object, pheno.col=1, addcovar=covars, model=opt$model, method=opt$method)
 } else {
-  verbose_print('Running scanone\n')
-  qtl_results = scanone(cross_object, pheno.col=1, model=opt$model, method=opt$method)
+  verbose_print('Running scan\n')
+  qtl_results = scan_func(cross_object, pheno.col=1, model=opt$model, method=opt$method)
 }
 
 #QTL main effects on adjusted longevity
-- 
cgit v1.2.3


From acfea0fc8793372723be0e64316002ed9bc2403b Mon Sep 17 00:00:00 2001
From: zsloan
Date: Fri, 23 Jul 2021 20:08:36 +0000
Subject: Added line priting pair-scan results to CSV and changed the default
 step-size to 10cM for pair-scan

---
 scripts/rqtl_wrapper.R | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index 0d13ccb..26d0a9c 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -177,7 +177,7 @@ if (!is.null(opt$interval)) {
   cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
 } else if (!is.null(opt$pairscan)) {
   verbose_print('Calculating genotype probabilities with interval mapping\n')
-  cross_object <- calc.genoprob(cross_object, step=20)
+  cross_object <- calc.genoprob(cross_object, step=10)
 } else {
   verbose_print('Calculating genotype probabilities\n')
   cross_object <- calc.genoprob(cross_object)
@@ -338,4 +338,9 @@ if (type == "4-way") {
   colnames(qtl_results)[4:7] <- c("AC", "AD", "BC", "BD")
 }
 
-write.csv(qtl_results, out_file)
+verbose_print('Writing results to CSV file\n')
+if (!is.null(opt$pairscan)) {
+  write.csv(qtl_results[1], out_file)
+} else {
+  write.csv(qtl_results, out_file)
+}
-- 
cgit v1.2.3


From 4c3359082fbb1eef93a9e54c633f938606cba5f6 Mon Sep 17 00:00:00 2001
From: zsloan
Date: Mon, 2 Aug 2021 21:07:52 +0000
Subject: Updated rqtl_wrapper to also return a map file when doing a pair-scan
 (since we need the list of markers/pseudomarkers and their positions)

---
 scripts/rqtl_wrapper.R | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index 26d0a9c..fb12012 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -176,7 +176,7 @@ if (!is.null(opt$interval)) {
   verbose_print('Calculating genotype probabilities with interval mapping\n')
   cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
 } else if (!is.null(opt$pairscan)) {
-  verbose_print('Calculating genotype probabilities with interval mapping\n')
+  verbose_print('Calculating genotype probabilities for pair-scan\n')
   cross_object <- calc.genoprob(cross_object, step=10)
 } else {
   verbose_print('Calculating genotype probabilities\n')
@@ -340,7 +340,9 @@ if (type == "4-way") {
 
 verbose_print('Writing results to CSV file\n')
 if (!is.null(opt$pairscan)) {
+  map_out_file = file.path(opt$outdir, paste("MAP_", opt$filename, sep = "" ))
   write.csv(qtl_results[1], out_file)
+  write.csv(qtl_results[2], map_out_file)
 } else {
   write.csv(qtl_results, out_file)
 }
-- 
cgit v1.2.3


From 0e6990940e22eb0431d6a27e45d29bc04d8ad582 Mon Sep 17 00:00:00 2001
From: zsloan
Date: Tue, 22 Mar 2022 20:28:58 +0000
Subject: Change order of if statements for running genoprob command

Now it checks for pairscan first, just in case interval ends up being
passed (which is an irrelevant parameter for pairscan)

Also added a couple more verbose prints
---
 scripts/rqtl_wrapper.R | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

(limited to 'scripts/rqtl_wrapper.R')

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index fb12012..5ecd774 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -172,12 +172,12 @@ verbose_print('Generating cross object\n')
 cross_object = geno_to_csvr(geno_file, trait_names, trait_vals, cross_file, type)
 
 # Calculate genotype probabilities
-if (!is.null(opt$interval)) {
-  verbose_print('Calculating genotype probabilities with interval mapping\n')
-  cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
-} else if (!is.null(opt$pairscan)) {
+if (!is.null(opt$pairscan)) {
   verbose_print('Calculating genotype probabilities for pair-scan\n')
   cross_object <- calc.genoprob(cross_object, step=10)
+} else if (!is.null(opt$interval)) {
+  verbose_print('Calculating genotype probabilities with interval mapping\n')
+  cross_object <- calc.genoprob(cross_object, step=5, stepwidth="max")
 } else {
   verbose_print('Calculating genotype probabilities\n')
   cross_object <- calc.genoprob(cross_object)
@@ -241,10 +241,12 @@ if (!is.null(opt$control)) {
 }
 
 if (!is.null(opt$pairscan)) {
+  verbose_print("Running scantwo")
   scan_func <- function(...){
     scantwo(...)
   }
 } else {
+  verbose_print("Running scanone")
   scan_func <- function(...){
     scanone(...)
   }
-- 
cgit v1.2.3