|
14
|
1
|
|
|
2 ########################################################
|
|
|
3 #
|
|
|
4 # creation date : 05/01/16
|
|
|
5 # last modification : 02/09/16
|
|
|
6 # author : Dr Nicolas Beaume
|
|
|
7 #
|
|
|
8 ########################################################
|
|
|
9
|
|
|
10
|
|
|
11 #log <- file(paste(getwd(), "log_rrBLUP.txt", sep="/"), open = "wt")
|
|
|
12 #sink(file = log, type="message")
|
|
|
13 library(rrBLUP)
|
|
|
14 ############################ helper functions #######################
|
|
|
15
|
|
|
16
|
|
|
17 ################################## main function ###########################
|
|
|
18
|
|
|
19 rrBLUP <- function(genotype, phenotype, outFile, evaluation=F, folds) {
|
|
|
20 if(evaluation) {
|
|
|
21 prediction <- NULL
|
|
|
22 for(i in 1:length(folds)) {
|
|
|
23 train <- genotype[-folds[[i]],]
|
|
|
24 test <- genotype[folds[[i]],]
|
|
|
25 phenoTrain <- phenotype[-folds[[i]]]
|
|
|
26 phenoTest <- phenotype[folds[[i]]]
|
|
|
27 model <-mixed.solve(phenoTrain, Z=genoTrain,K=NULL, SE=F,return.Hinv = F)
|
|
|
28 pred <- as.matrix(genoTest) %*% as.matrix(model$u)
|
|
|
29 pred <- pred[,1]+model$beta
|
|
|
30 prediction <- c(prediction, list(pred))
|
|
|
31 }
|
|
|
32 saveRDS(prediction, file=paste(outFile,".rds", sep=""))
|
|
|
33 # just create a model
|
|
|
34 } else {
|
|
|
35 model <-mixed.solve(phenotype, Z=genotype,K=NULL, SE=F,return.Hinv = F)
|
|
|
36 saveRDS(model, file = paste(outFile, ".rds", sep = ""))
|
|
|
37 }
|
|
|
38 }
|
|
|
39
|
|
|
40
|
|
|
41 ############################ main #############################
|
|
|
42 # running from terminal (supposing the OghmaGalaxy/bin directory is in your path) :
|
|
|
43 # rrBLUP.sh -i path_to_genotype -p phenotype_file -f folds -e -f fold_file -o out_file
|
|
|
44 ## -i : path to the file that contains the genotypes, must be a .rda file (as outputed by loadGenotype).
|
|
|
45 # please note that the table must be called "encoded" when your datafile is saved into .rda (automatic if encoded.R was used)
|
|
|
46
|
|
|
47 ## -p : file that contains the phenotype must be a .rda file (as outputed by loadGenotype.R).
|
|
|
48 # please note that the table must be called "phenotype" when your datafile is saved into .rda (automatic if loadGenotype.R was used)
|
|
|
49
|
|
|
50 ## -e : do evaluation instead of producing a model
|
|
|
51
|
|
|
52 ## -f : index of the folds to which belong each individual
|
|
|
53 # please note that the list must be called "folds" (automatic if folds.R was used)
|
|
|
54
|
|
|
55 ## -o : path to the output folder and generic name of the analysis. The extension related to each type of
|
|
|
56 # files are automatically added
|
|
|
57
|
|
|
58 cmd <- commandArgs(T)
|
|
|
59 source(cmd[1])
|
|
|
60 # load genotype and phenotype
|
|
|
61 con = file(genotype)
|
|
|
62 genotype <- readLines(con = con, n = 1, ok=T)
|
|
|
63 close(con)
|
|
|
64 genotype <- read.table(genotype, sep="\t", h=T)
|
|
|
65 # phenotype is written as a table (in columns) but it must be sent as a vector for mixed.solve
|
|
|
66 phenotype <- read.table(phenotype, sep="\t", h=T)[,1]
|
|
|
67 # run !
|
|
|
68 rrBLUP(genotype = genotype, phenotype = phenotype, outFile = out)
|
|
|
69 cat(paste(paste(out, ".rds", sep = ""), "\n", sep="")) |
