aboutsummaryrefslogtreecommitdiff
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(...)
}