| Previous changeset 78:5a56e1e3b235 (2016-10-28) Next changeset 80:4ce7b1bdb320 (2016-10-28) |
|
Commit message:
Uploaded |
|
added:
lasso.R |
| b |
| diff -r 5a56e1e3b235 -r 14b976f46889 lasso.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/lasso.R Fri Oct 28 08:46:49 2016 -0400 |
| [ |
| @@ -0,0 +1,122 @@ +######################################################## +# +# creation date : 08/01/16 +# last modification : 01/09/16 +# author : Dr Nicolas Beaume +# owner : IRRI +# +######################################################## + +suppressWarnings(suppressMessages(library(glmnet))) +library(methods) +############################ helper functions ####################### + + +# optimize alpha parameter +optimize <- function(genotype, phenotype, alpha=seq(0,1,0.1), repet=7) { + acc <- NULL + indexAlpha <- 1 + for(a in alpha) { + curAcc <- NULL + # repeat nfolds time each analysis + for(i in 1:repet) { + # draw at random 1/3 of the training set for testing and thus choose alpha + # note it is not a cross-validation + n <- ceiling(nrow(genotype)/3) + indexTest <- sample(1:nrow(genotype), size=n) + # create training set and test set + train <- genotype[-indexTest,] + test <- genotype[indexTest,] + phenoTrain <- phenotype[-indexTest] + phenoTest <- phenotype[indexTest] + # cv.glmnet allow to compute lambda at the current alpha + cv <- cv.glmnet(x=as.matrix(train), y=phenoTrain, alpha=a) + # create model + model <- glmnet(x=as.matrix(train), y=phenoTrain, alpha=a, lambda = cv$lambda.1se) + # predict test set + pred <- predict(model, test, type = "response") + # compute r2 for choosing the best alpha + curAcc <- c(curAcc, r2(phenoTest, pred)) + } + # add mean r2 for this value of lambda to the accuracy vector + acc <- c(acc, mean(curAcc)) + } + # choose best alpha + names(acc) <- alpha + return(as.numeric(names(acc)[which.max(acc)])) +} + +# compute r2 by computing the classic formula +# compare the sum of square difference from target to prediciton +# to the sum of square difference from target to the mean of the target +r2 <- function(target, prediction) { + sst <- sum((target-mean(target))^2) + ssr <- sum((target-prediction)^2) + return(1-ssr/sst) +} +################################## main function ########################### + +lasso <- function(genotype, phenotype, evaluation = T, outFile, folds, alpha=NULL) { + # go for optimization + if(is.null(alpha)) { + alpha <- seq(0,1,0.1) + alpha <- optimize(genotype=genotype, phenotype=phenotype, alpha = alpha) + } + # evaluation + if(evaluation) { + prediction <- NULL + # do cross-validation + for(i in 1:length(folds)) { + # create training and test set + train <- genotype[-folds[[i]],] + test <- genotype[folds[[i]],] + phenoTrain <- phenotype[-folds[[i]]] + phenoTest <- phenotype[folds[[i]]] + # cv.glmnet helps to compute the right lambda for a chosen alpha + cv <- cv.glmnet(x=as.matrix(train), y=phenoTrain, alpha=alpha) + # create model + lasso.fit <- glmnet(x=as.matrix(train), y=phenoTrain, alpha=alpha, lambda = cv$lambda.1se) + # predict value of the test set for further evaluation + prediction <- c(prediction, list(predict(lasso.fit, test, type = "response")[,1])) + } + # save predicted value for test set of each fold for further evaluation + saveRDS(prediction, file=paste(outFile,".rds", sep="")) + # just create a model + } else { + # cv.glmnet helps to compute the right lambda for a chosen alpha + cv <- cv.glmnet(x=genotype, y=phenotype, alpha=alpha) + # create model + model <- glmnet(x=genotype, y=phenotype, alpha=alpha, lambda=cv$lambda.1se) + # save model + saveRDS(model, file = paste(outFile, ".rds", sep = "")) + } +} + +############################ main ############################# +# load argument +cmd <- commandArgs(T) +source(cmd[1]) +# check if evaluation is required +evaluation <- F +if(as.integer(doEvaluation) == 1) { + evaluation <- T + con = file(folds) + folds <- readLines(con = con, n = 1, ok=T) + close(con) + folds <- readRDS(folds) +} +# load classifier parameters +alpha <- as.numeric(alpha) +if(alpha < 0 | alpha > 1) {alpha <- NULL} +# 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 is written as a table (in columns) but it must be sent as a vector for mixed.solve +phenotype <- read.table(phenotype, sep="\t", h=T)[,1] +# run ! +lasso(genotype = data.matrix(genotype), phenotype = phenotype, + evaluation = evaluation, outFile = out, folds = folds, alpha = alpha) +# return path of the result file to galaxy +cat(paste(paste(out, ".rds", sep = ""), "\n", sep="")) |