|
87
|
1
|
|
|
2 ########################################################
|
|
|
3 #
|
|
|
4 # creation date : 05/01/16
|
|
|
5 # last modification : 25/10/16
|
|
|
6 # author : Dr Nicolas Beaume
|
|
|
7 #
|
|
|
8 ########################################################
|
|
|
9
|
|
|
10
|
|
|
11 library(rrBLUP)
|
|
|
12 ############################ helper functions #######################
|
|
|
13
|
|
|
14
|
|
|
15 ################################## main function ###########################
|
|
|
16 # do rrBLUP evaluation of classification.
|
|
|
17 # optimization of paramaters is included in rrBLUP package
|
|
|
18 rrBLUP <- function(genotype, phenotype, outFile, evaluation=F, folds) {
|
|
|
19 # Evaluation mode
|
|
|
20 if(evaluation) {
|
|
|
21 prediction <- NULL
|
|
|
22 # run over folds
|
|
|
23 for(i in 1:length(folds)) {
|
|
|
24 # create training and test set for this fold
|
|
|
25 train <- genotype[-folds[[i]],]
|
|
|
26 test <- genotype[folds[[i]],]
|
|
|
27 phenoTrain <- phenotype[-folds[[i]]]
|
|
|
28 phenoTest <- phenotype[folds[[i]]]
|
|
|
29 # create model
|
|
|
30 model <-mixed.solve(phenoTrain, Z=train,K=NULL, SE=F,return.Hinv = F)
|
|
|
31 # predict current test set
|
|
|
32 pred <- as.matrix(test) %*% as.matrix(model$u)
|
|
|
33 pred <- pred[,1]+model$beta
|
|
|
34 prediction <- c(prediction, list(pred))
|
|
|
35 }
|
|
|
36 # save results
|
|
|
37 saveRDS(prediction, file=paste(outFile,".rds", sep=""))
|
|
|
38 # just create a model
|
|
|
39 } else {
|
|
|
40 # create and save modle
|
|
|
41 model <-mixed.solve(phenotype, Z=genotype,K=NULL, SE=F,return.Hinv = F)
|
|
|
42 saveRDS(model, file = paste(outFile, ".rds", sep = ""))
|
|
|
43 }
|
|
|
44 }
|
|
|
45
|
|
|
46
|
|
|
47 ############################ main #############################
|
|
|
48 # get argument from xml file
|
|
|
49 cmd <- commandArgs(T)
|
|
|
50 source(cmd[1])
|
|
|
51 # for evaluation mode : set evaluation as True and load fold file
|
|
|
52 if(as.integer(evaluation) == 1) {
|
|
|
53 evaluation <- T
|
|
|
54 con = file(folds)
|
|
|
55 folds <- readLines(con = con, n = 1, ok=T)
|
|
|
56 close(con)
|
|
|
57 folds <- readRDS(folds)
|
|
|
58 } else{
|
|
|
59 evaluation <- F
|
|
|
60 }
|
|
|
61 # load genotype and phenotype
|
|
|
62 con = file(genotype)
|
|
|
63 genotype <- readLines(con = con, n = 1, ok=T)
|
|
|
64 close(con)
|
|
|
65 genotype <- read.table(genotype, sep="\t", h=T)
|
|
|
66 phenotype <- read.table(phenotype, sep="\t", h=T)[,1]
|
|
|
67 # run !
|
|
|
68 rrBLUP(genotype = genotype, phenotype = phenotype, outFile = out,
|
|
|
69 evaluation = evaluation, folds = folds)
|
|
|
70 # return path of the result to galaxy
|
|
|
71 cat(paste(paste(out, ".rds", sep = ""), "\n", sep="")) |
