diff options
| -rw-r--r-- | scripts/rqtl_wrapper.R | 51 |
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) |