| Previous changeset 86:2212133e6a36 (2016-10-28) Next changeset 88:aef3240b58ac (2016-10-28) |
|
Commit message:
Uploaded |
|
added:
rrBLUP.R |
| b |
| diff -r 2212133e6a36 -r 2f423d8656ae rrBLUP.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rrBLUP.R Fri Oct 28 08:49:07 2016 -0400 |
| [ |
| @@ -0,0 +1,71 @@ + +######################################################## +# +# creation date : 05/01/16 +# last modification : 25/10/16 +# author : Dr Nicolas Beaume +# +######################################################## + + +library(rrBLUP) +############################ helper functions ####################### + + +################################## main function ########################### +# do rrBLUP evaluation of classification. +# optimization of paramaters is included in rrBLUP package +rrBLUP <- function(genotype, phenotype, outFile, evaluation=F, folds) { + # Evaluation mode + if(evaluation) { + prediction <- NULL + # run over folds + for(i in 1:length(folds)) { + # create training and test set for this fold + train <- genotype[-folds[[i]],] + test <- genotype[folds[[i]],] + phenoTrain <- phenotype[-folds[[i]]] + phenoTest <- phenotype[folds[[i]]] + # create model + model <-mixed.solve(phenoTrain, Z=train,K=NULL, SE=F,return.Hinv = F) + # predict current test set + pred <- as.matrix(test) %*% as.matrix(model$u) + pred <- pred[,1]+model$beta + prediction <- c(prediction, list(pred)) + } + # save results + saveRDS(prediction, file=paste(outFile,".rds", sep="")) + # just create a model + } else { + # create and save modle + model <-mixed.solve(phenotype, Z=genotype,K=NULL, SE=F,return.Hinv = F) + saveRDS(model, file = paste(outFile, ".rds", sep = "")) + } +} + + +############################ main ############################# +# get argument from xml file +cmd <- commandArgs(T) +source(cmd[1]) +# for evaluation mode : set evaluation as True and load fold file +if(as.integer(evaluation) == 1) { + evaluation <- T + con = file(folds) + folds <- readLines(con = con, n = 1, ok=T) + close(con) + folds <- readRDS(folds) +} else{ + evaluation <- F +} +# load genotype and phenotype +con = file(genotype) +genotype <- readLines(con = con, n = 1, ok=T) +close(con) +genotype <- read.table(genotype, sep="\t", h=T) +phenotype <- read.table(phenotype, sep="\t", h=T)[,1] +# run ! +rrBLUP(genotype = genotype, phenotype = phenotype, outFile = out, + evaluation = evaluation, folds = folds) +# return path of the result to galaxy +cat(paste(paste(out, ".rds", sep = ""), "\n", sep="")) \ No newline at end of file |