comparison rrBLUP.R @ 87:2f423d8656ae draft

Uploaded

87:2f423d8656ae

########################################################
#
# 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=""))