diff options
Diffstat (limited to 'scripts/rqtl_wrapper.R')
-rw-r--r-- | scripts/rqtl_wrapper.R | 20 |
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(...) } |