--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/prediction.R	Fri Oct 28 08:47:48 2016 -0400
@@ -0,0 +1,52 @@
+########################################################
+#
+# creation date : 26/01/16
+# last modification : 02/06/16
+# author : Dr Nicolas Beaume
+# owner : IRRI
+#
+########################################################
+
+library(rrBLUP)
+suppressWarnings(suppressMessages(library(randomForest)))
+library(e1071)
+suppressWarnings(suppressMessages(library(glmnet)))
+library(methods)
+
+
+############################ main #############################
+classifierNames <- c("list", "randomForest", "svm", "glmnet")
+# load argument
+cmd <- commandArgs(trailingOnly = T)
+source(cmd[1])
+# load data
+con = file(genotype)
+genotype <- readLines(con = con, n = 1, ok=T)
+close(con)
+genotype <- read.table(genotype, sep="\t", h=T)
+con = file(model)
+model <- readLines(con = con, n = 1, ok=T)
+close(con)
+model <- readRDS(model)
+# check if the classifier name is valid
+if(all(is.na(match(class(model), classifierNames)))) {
+  stop(paste(class(model), "is not recognized as a valid model. Valid models are : ", classifierNames))
+}
+# run prediction according to the classifier
+# rrBLUP prediction
+if(any(class(model) == "list")) {
+  predictions <- as.matrix(genotype) %*% as.matrix(model$u)
+  predictions <- predictions[,1]+model$beta
+  predictions <- data.frame(lines=rownames(genotype), predictions=predictions)
+# LASSO prediction  
+} else if(any(class(model) == "glmnet")) {
+  predictions <- predict(model, as.matrix(genotype), type = "response")
+  predictions <- data.frame(lines=rownames(genotype), predictions=predictions)
+# SVM or RandomForest prediction (predict is a wrapper for many machine learning function)
+} else {
+  predictions <- predict(model, genotype)
+  predictions <- data.frame(lines=names(predictions), predictions=predictions)
+}
+# save results
+write.table(predictions, file=out, sep="\t", row.names = F)
+