diff --git a/inst/scripts/pqtl-genome-scan.R b/inst/scripts/pqtl-genome-scan.R index 4e2b6539974ebb9614dac1c408a3e052f01a7822..0b297add8e816494ff7af3a6fefec8cdd28ccb04 100644 --- a/inst/scripts/pqtl-genome-scan.R +++ b/inst/scripts/pqtl-genome-scan.R @@ -8,11 +8,9 @@ library(qtl2) args <- commandArgs(trailingOnly = TRUE) #prev args -qtldata <- args[1] #path to yaml file required to create crossdata object (is that still needed ?) -crossdata_path <- args[2] #path to crossdata object -clean_pr_path <- args[3] #path to genoprobs objects -kloco_path <- args[4] #path to kinship (loco) objects -cores <- args[5] #nb of cores +crossdata_path <- args[1] #path to crossdata object +clean_pr_path <- args[2] #path to genoprobs objects +kloco_path <- args[3] #path to kinship (loco) objects #as of now (not using args) # qtldata <- "TEST.yaml" @@ -22,11 +20,14 @@ cores <- args[5] #nb of cores # cores <- 1 -#novel args re: covar -isSexCovar <- args[6] #sexcovar: drop-down menu perhaps ? two options, either yes or no. if yes, extraction of sexcovar and put in model. if no, nope. -whichAddCovars <- args[7] #addcovars: at most 2 non-sex covar -#rémi, re: addcovars (else see issue). This can be a drop down menu (or two drop-down menus? but whatif only one covar is selected?) listing all possible covars existing in covar file (BUT sex covar). Only two of them at most can be selected, perhaps a +# novel args re: covar +# sexcovar: drop-down menu perhaps ? two options, either yes or no. if yes, extraction of sexcovar and put in model. if no, nope. +isSexCovar <- args[4] +whichAddCovars <- args[5] #addcovars: at most 2 non-sex covar +#rémi, re: addcovars (else see issue). This can be a drop down menu (or two drop-down menus? but whatif only one covar is selected?) listing all possible covars existing in covar file (BUT sex covar). Only two of them at most can be selected, perhaps a ##remi: il faudrait pouvoir proposer un champ "none" (pour aucun ajout de covar), dans ce cas whichAddCovars <- NULL +cores <- args[6] #nb of cores + #as of now # isSexCovar <- 'yes' @@ -44,17 +45,21 @@ dir.create("outfiles/Rdata") ## ----def covars-------------------------------------------- ##---------extract Xcovar (under the condition that the user spec he wants to apply a sex covar)------ -if (isSexCovar=="yes") { +if (isSexCovar == "yes") { Xcovar <- get_x_covar(crossdata) -} +} ###-------- as for addcovars: if whichAddCovars is not null (exists at least 1 covar besides sex), generates the covar formula for model matrix and all what follows. -if (!is.null(whichAddCovars)==TRUE) { - covar.formula <- as.formula(sprintf("~ 0 + %s", gsub(",", "+", whichAddCovars))) #mk the covar formula for model matrix (tested, works with 1 & 2 addcovars) - covar.modmat <- model.matrix(covar.formula,data = crossdata$covar) #get the model matrix out of it - covar.modmat <- covar.modmat[,-1, drop=FALSE] # omit the first column bkos k-1 levels; drop=FALSE ensures it stays as a matrix - rownames(covar.modmat) <- rownames(crossdata$covar) #rownames should have remained identical but still, better safe than sorry -} +if (!is.null(whichAddCovars) == TRUE) { + covar.formula <- + as.formula(sprintf("~ 0 + %s", gsub(",", "+", whichAddCovars))) #mk the covar formula for model matrix (tested, works with 1 & 2 addcovars) + covar.modmat <- + model.matrix(covar.formula, data = crossdata$covar) #get the model matrix out of it + covar.modmat <- + covar.modmat[,-1, drop = FALSE] # omit the first column bkos k-1 levels; drop=FALSE ensures it stays as a matrix + rownames(covar.modmat) <- + rownames(crossdata$covar) #rownames should have remained identical but still, better safe than sorry +} #there could be a simpler way,actually ?: addcovar <- model.matrix(~Study, data = pheno)[,-1] (but I prefer the longer one...) @@ -69,13 +74,45 @@ if (!is.null(whichAddCovars)==TRUE) { if (exists("Xcovar") & exists("covar.modmat")) { - out <- scan1(genoprobs = clean_pr, pheno = crossdata$pheno, kinship = kloco, addcovar = covar.modmat, Xcovar = Xcovar, cores = cores, max_batch = 1) ##sex covar AND addcovar -} else if(exists("Xcovar") & !exists("covar.modmat")) { - out <- scan1(genoprobs = clean_pr, pheno = crossdata$pheno, kinship = kloco, Xcovar = Xcovar, cores = cores, max_batch = 1) ##sexcovar but no addcovar -} else if(!exists("Xcovar") & exists("covar.modmat")) { - out <- scan1(genoprobs = clean_pr, pheno = crossdata$pheno, kinship = kloco, addcovar = covar.modmat, cores = cores, max_batch = 1) ##addcovar but no sexcovar + out <- + scan1( + genoprobs = clean_pr, + pheno = crossdata$pheno, + kinship = kloco, + addcovar = covar.modmat, + Xcovar = Xcovar, + cores = cores, + max_batch = 1 + ) ##sex covar AND addcovar +} else if (exists("Xcovar") & !exists("covar.modmat")) { + out <- + scan1( + genoprobs = clean_pr, + pheno = crossdata$pheno, + kinship = kloco, + Xcovar = Xcovar, + cores = cores, + max_batch = 1 + ) ##sexcovar but no addcovar +} else if (!exists("Xcovar") & exists("covar.modmat")) { + out <- + scan1( + genoprobs = clean_pr, + pheno = crossdata$pheno, + kinship = kloco, + addcovar = covar.modmat, + cores = cores, + max_batch = 1 + ) ##addcovar but no sexcovar } else { - out <- scan1(genoprobs = clean_pr, pheno = crossdata$pheno, kinship = kloco, cores = cores, max_batch = 1) ##neither addcovar nor sexcovar + out <- + scan1( + genoprobs = clean_pr, + pheno = crossdata$pheno, + kinship = kloco, + cores = cores, + max_batch = 1 + ) ##neither addcovar nor sexcovar } @@ -103,17 +140,17 @@ test_datplot <- data.frame( pos = unlist(crossdata$pmap), lod = out ) -# !!! Remi !! the lod csv file now contains SEVERAL phenotypes - perhaps drop down menu to viz them ?? +# !!! Remi !! the lod csv file now contains SEVERAL phenotypes - perhaps drop down menu to viz them ?? write.csv(test_datplot, "outfiles/LOD/lod.csv", row.names = FALSE) save(out, file = "outfiles/Rdata/out.dat") #if exists, save the covariate files if (exists("covar.modmat")) { - save(covar.modmat,file="outfiles/Rdata/covar.modmat.dat") + save(covar.modmat, file = "outfiles/Rdata/covar.modmat.dat") } if (exists("Xcovar")) { - save(Xcovar,file="outfiles/Rdata/Xcovar.dat") -} -#save(pr, file = "outfiles/Rdata/pr.dat") #not needed, already existing -#save(crossdata, file = "outfiles/Rdata/crossdata.dat") #not needed, already existing -#save(kinship, file = "outfiles/Rdata/kinship.dat") #not needed, already existing + save(Xcovar, file = "outfiles/Rdata/Xcovar.dat") +} +#save(pr, file = "outfiles/Rdata/pr.dat") #not needed, already existing +#save(crossdata, file = "outfiles/Rdata/crossdata.dat") #not needed, already existing +#save(kinship, file = "outfiles/Rdata/kinship.dat") #not needed, already existing