aboutsummaryrefslogtreecommitdiff
path: root/scripts
diff options
context:
space:
mode:
authorzsloan2021-05-26 15:09:00 +0000
committerzsloan2021-06-18 22:08:04 +0000
commit9b8fb0f63114119e58663a3e7fee352ce7375fb4 (patch)
tree9251666b0b5ad60ce8cb2cf8797909a17f531554 /scripts
parent0bfe2990ce31fbd808b680a883619575d864aef5 (diff)
downloadgenenetwork3-9b8fb0f63114119e58663a3e7fee352ce7375fb4.tar.gz
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
Diffstat (limited to 'scripts')
-rw-r--r--scripts/rqtl_wrapper.R51
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)