Added option for enabling permutation strata

Fixed issue where covariates with numerical names were being ignored

Fixed issue where the hash filename wasn't being used for the permutation output

1 files changed, 38 insertions, 13 deletions

diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index ae970d4..f7e0406 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -13,6 +13,7 @@ option_list = list(
   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("-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)"),
   make_option(c("-s", "--scale"), type="character", default="mb", help="Mapping scale - Megabases (Mb) or Centimorgans (cM)"),
   make_option(c("--control"), type="character", default=NULL, help="Name of marker (contained in genotype file) to be used as a control"),
   make_option(c("-o", "--outdir"), type="character", default=file.path(tmp_dir, "output"), help="Directory in which to write result file"),
@@ -135,15 +136,17 @@ df <- read.table(pheno_file, na.strings = "x", header=TRUE, check.names=FALSE)
 sample_names <- df$Sample
 trait_names <- colnames(df)[2:length(colnames(df))]
 
+# Since there will always only be one non-covar phenotype, its name will be in the first column
+pheno_name = unlist(trait_names)[1]
+
 trait_vals <- vector(mode = "list", length = length(trait_names))
 for (i in 1:length(trait_names)) {
   this_trait <- trait_names[i]
   this_vals <- df[this_trait]
   trait_vals[[i]] <- this_vals
-}
 
-# Since there will always only be one non-covar phenotype, its name will be in the first column
-pheno_name = unlist(trait_names)[1]
+  trait_names[i] = paste("T_", this_trait, sep = "")
+}
 
 verbose_print('Generating cross object\n')
 cross_object = geno_to_csvr(geno_file, trait_names, trait_vals, cross_file)
@@ -158,12 +161,24 @@ if (!is.null(opt$interval)) {
 }
 
 # Pull covariates out of cross object, if they exist
-covars = vector()
+covars = vector(mode = "list", length = length(trait_names) - 1)
 if (!is.null(opt$addcovar)) {
-  covar_names = trait_names[2:length(trait_names)]
+  #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]
+  } else {
+    covar_names = trait_names[3:length(trait_names)]
+  }
   covars <- pull.pheno(cross_object, covar_names)
 }
 
+# Pull permutation strata out of cross object, if it is being used
+perm_strata = vector()
+if (!is.null(opt$pstrata)) {
+  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)) {
   marker_covars = create_marker_covars(cross_object, opt$control)
@@ -173,17 +188,27 @@ if (!is.null(opt$control)) {
 # Calculate permutations
 if (opt$nperm > 0) {
   if (!is.null(opt$filename)){
-    perm_out_file = file.path(opt$outdir, paste("PERM_", opt$filename, sep = ""))
+    perm_out_file = file.path(opt$outdir, paste("PERM_", opt$filename, sep = "" ))
   } else {
-    perm_out_file = file.path(opt$outdir, paste(pheno_name, "_PERM_", stri_rand_strings(1, 8), ".csv", sep = ""))
+    perm_out_file = file.path(opt$outdir, paste(pheno_name, "_PERM_", stri_rand_strings(1, 8), sep = ""))
   }
 
   if (!is.null(opt$addcovar) || !is.null(opt$control)){
-    verbose_print('Running permutations with cofactors\n')
-    perm_results = scanone(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, model=opt$model, method=opt$method)
+    if (!is.null(opt$pstrata)) {
+      verbose_print('Running 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)
+    } else {
+      verbose_print('Running permutations with cofactors\n')
+      perm_results = scanone(cross_object, pheno.col=1, addcovar=covars, n.perm=opt$nperm, model=opt$model, method=opt$method)
+    }
   } else {
-    verbose_print('Running permutations\n')
-    perm_results = scanone(cross_object, pheno.col=1, n.perm=opt$nperm, model=opt$model, method=opt$method)
+    if (!is.null(opt$pstrata)) {
+      verbose_print('Running 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)
+    } else {
+      verbose_print('Running permutations\n')
+      perm_results = scanone(cross_object, pheno.col=1, n.perm=opt$nperm, model=opt$model, method=opt$method)
+    }
   }
   write.csv(perm_results, perm_out_file)
 }
@@ -191,7 +216,7 @@ if (opt$nperm > 0) {
 if (!is.null(opt$filename)){
   out_file = file.path(opt$outdir, opt$filename)
 } else {
-  out_file = file.path(opt$outdir, paste(pheno_name, "_", stri_rand_strings(1, 8), ".csv", sep = ""))
+  out_file = file.path(opt$outdir, paste(pheno_name, "_", stri_rand_strings(1, 8), sep = ""))
 }
 
 if (!is.null(opt$addcovar) || !is.null(opt$control)){
@@ -201,4 +226,4 @@ if (!is.null(opt$addcovar) || !is.null(opt$control)){
   verbose_print('Running scanone\n')
   qtl_results = scanone(cross_object, pheno.col=1, model=opt$model, method=opt$method)
 }
-write.csv(qtl_results, out_file)
\ No newline at end of file
+write.csv(qtl_results, out_file)