aboutsummaryrefslogtreecommitdiff
path: root/scripts
diff options
context:
space:
mode:
Diffstat (limited to 'scripts')
-rw-r--r--scripts/rqtl_wrapper.R43
1 files changed, 43 insertions, 0 deletions
diff --git a/scripts/rqtl_wrapper.R b/scripts/rqtl_wrapper.R
index 20c0e06..523888f 100644
--- a/scripts/rqtl_wrapper.R
+++ b/scripts/rqtl_wrapper.R
@@ -269,4 +269,47 @@ 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)
}
+
+#QTL main effects on adjusted longevity
+getEffects <- function(sdata, gtsprob, marker = "1_24042124", model = "longevity ~ sex + site + cohort + treatment", trait = "longevity"){
+ rownames(sdata) <- 1:nrow(sdata)
+ rownames(gtsprob) <- 1:nrow(gtsprob)
+ mp <- gtsprob[, grep(marker, colnames(gtsprob))]
+ gts <- unlist(lapply(lapply(lapply(apply(mp,1,function(x){which(x > 0.85)}),names), strsplit, ":"), function(x){
+ if(length(x) > 0){ return(x[[1]][2]); }else{ return(NA) }
+ }))
+
+ ismissing <- which(apply(sdata, 1, function(x){any(is.na(x))}))
+ if(length(ismissing) > 0){
+ sdata <- sdata[-ismissing, ]
+ gts <- gts[-ismissing]
+ }
+
+ mlm <- lm(as.formula(model), data = sdata)
+ pheAdj <- rep(NA, nrow(sdata))
+ adj <- residuals(mlm) + mean(sdata[, trait])
+ pheAdj[as.numeric(names(adj))] <- adj
+ means <- c(mean(pheAdj[which(gts == "AC")],na.rm=TRUE),mean(pheAdj[which(gts == "AD")],na.rm=TRUE),mean(pheAdj[which(gts == "BC")],na.rm=TRUE),mean(pheAdj[which(gts == "BD")],na.rm=TRUE))
+ std <- function(x) sd(x,na.rm=TRUE)/sqrt(length(x))
+ stderrs <- c(std(pheAdj[which(gts == "AC")]),std(pheAdj[which(gts == "AD")]),std(pheAdj[which(gts == "BC")]),std(pheAdj[which(gts == "BD")]))
+ paste0(round(means,0), " ± ", round(stderrs,2))
+}
+
+if (type == "4-way") {
+ verbose_print("Get phenotype name + genoprob + all phenotypes + models for 4-way crosses")
+ traitname <- colnames(pull.pheno(cross_object))[1]
+ gtsp <- pull.genoprob(cross_object)
+ allpheno <- pull.pheno(cross_object)
+ model <- paste0(traitname, " ~ ", paste0(covar_names, sep="", collapse=" + "))
+
+ meffects <- c()
+ verbose_print("Getting QTL main effects for 4-way crosses")
+ for(marker in rownames(qtl_results)){
+ meff <- getEffects(allpheno, gtsp, marker = marker, model, trait = traitname)
+ meffects <- rbind(meffects, meff)
+ }
+ qtl_results <- cbind(data.frame(qtl_results[,1:3]), meffects)
+ colnames(qtl_results)[4:7] <- c("AC", "AD", "BC", "BD")
+}
+
write.csv(qtl_results, out_file)