about summary refs log tree commit diff
path: root/scripts/rqtl_wrapper.R
diff options
context:
space:
mode:
Diffstat (limited to 'scripts/rqtl_wrapper.R')
-rw-r--r--scripts/rqtl_wrapper.R20
1 files changed, 15 insertions, 5 deletions
diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index eb660b5..ea2c345 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -172,7 +172,10 @@ 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)) {
+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 if (!is.null(opt$pairscan)) {
@@ -207,14 +210,19 @@ if (!is.null(opt$addcovar)) {
     name <- paste0("T_", covarDescr[x, 1]) # The covar description file doesn't have T_ in trait names (the cross object does)
     type <- covarDescr[x, 2]
     if(type == "categorical"){
-      if(length(table(covars[,name])) > 2){ # More then 2 levels create the model matrix for the factor
-        mdata <- data.frame(toExpand = as.factor(covars[, name]))
+      verbose_print('Binding covars to covars\n')
+      if (nrow(covarDescr) < 2){
+        this_covar = covars
+      } else {
+        this_covar = covars[,name]
+      }
+      if(length(table(this_covar)) > 2){ # More then 2 levels create the model matrix for the factor
+        mdata <- data.frame(toExpand = as.factor(this_covar))
         options(na.action='na.pass')
         modelmatrix <- model.matrix(~ toExpand + 0, mdata)[,-1]
         covars <- cbind(covars, modelmatrix)
       }else{ # 2 levels? just bind the trait as covar
-        verbose_print('Binding covars to covars\n')
-        covars <- cbind(covars, covars[,name])
+        covars <- cbind(covars, this_covar)
       }
     }
   }
@@ -236,10 +244,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(...)
   }