R/mixed_models_sommer.R

   1 ## File was edited in ESS mode in emacs and some indents changed. sorry.
   2 ##First read in the arguments listed at the command line
   3 args=commandArgs(TRUE)
   4
   5 ##args is now a list of character vectors
   6 ## First check to see if arguments are passed.
   7 ## Then cycle through each element of the list and evaluate the expressions.
   8
   9 if(length(args)==0){
  10     print("No arguments supplied.")
  11     ##supply default values
  12     datafile=''
  13     paramfile=''
  14     fixed_factors<-c()
  15     random_factors<-c()
  16
  17 } else {
  18     for(i in 1:length(args)){
  19         print(paste("Processing arg ", args[[i]]));
  20         eval(parse(text=args[[i]]))
  21     }
  22 }
  23
  24
  25 workdir = dirname(datafile);
  26 setwd(workdir);
  27
  28 #library(lme4)
  29 #library(lmerTest)
  30 library(sommer)
  31 library(emmeans)
  32 #library(effects)
  33 library(stringr)
  34 library(dplyr)
  35
  36 pd = read.csv(datafile, sep="\t")
  37
  38 source(paramfile)  # should give us dependent_variable and the model
  39
  40 trait = dependent_variables
  41 pd$studyYear = as.factor(pd$studyYear)
  42 print(paste("FIXED MODEL :", fixed_model))
  43 print(paste("RANDOM MODEL: ", random_model))
  44 print(paste("FIXED FACTORS: ", fixed_factors))
  45 print(paste("RANDOM FACTORS: ", random_factors))
  46 #print(head(pd))
  47
  48 BLUE = as.data.frame(unique(pd$germplasmName))
  49 colnames(BLUE) = "germplasmName"
  50 adjustedBLUE = BLUE
  51
  52 BLUPS = BLUE
  53 adjusted_means = BLUE
  54
  55 for(i in 1:length(trait)){
  56
  57
  58     #print(paste("processing trait", trait[i]))
  59     dependent_variables = trait[i]
  60     dependent_variables = gsub(" ", "\\.", dependent_variables) # replace space with "." in variable name
  61     dependent_variables = gsub("\\|", "\\.", dependent_variables) # replace | with .
  62     dependent_variables = gsub("\\:", "\\.", dependent_variables)
  63     dependent_variables = gsub("\\-", "\\.", dependent_variables)
  64     dependent_variables = gsub("\\/", "\\.", dependent_variables)
  65
  66     print(paste("Dependent variables : ", dependent_variables))
  67
  68     pd <- pd[!(is.na(pd[c(dependent_variables)])), ]
  69
  70
  71
  72     genotypeEffectType = as.vector(str_match(random_model, 'germplasmName'))
  73     genotypeEffectType = ifelse(is.na(genotypeEffectType), 'fixed', 'random')
  74     print(paste('modeling genotypes as: ', genotypeEffectType))
  75
  76     if (genotypeEffectType=="random") {
  77
  78         mixmodel = mmer(as.formula(fixed_model), random = as.formula(random_model), rcov = ~ units, data=pd)
  79         print(paste("MIXED MODEL: ", mixmodel))
  80               varcomp<- summary(mixmodel)$varcomp
  81         print(varcomp)
  82
  83         print("---------")
  84
  85         print("---------")
  86
  87         ##BLUPS
  88         res <- randef(mixmodel) ##obtain the blups
  89         BLUP <- as.data.frame(res)<
  90         #print(BLUP)
  91
  92         ##ajusted means
  93         p0 = predict.mmer(mixmodel, D="germplasmName") ##runs the prediction
  94         summary(p0)
  95         adj = p0$pvals                                        ##obtains the predictions
  96         adjusted_means = as.data.frame(adj)
  97         #print(paste("adj", adj))
  98         print(adjusted_means)
  99
 100
 101     } else {
 102         if (random_model!="") {
 103         mixmodel = mmer(as.formula(fixed_model), random = as.formula(random_model), rcov = ~ units, data=pd)
 104         } else {
 105         mixmodel = mmer(as.formula(fixed_model),  rcov = ~ units, data=pd)
 106
 107         }
 108
 109         varcomp<-summary(mixmodel)$varcomp
 110         print(paste("MIXED MODEL: ", mixmodel))
 111
 112         #Computing fixed effects
 113         blue = summary(mixmodel)$beta
 114         BLUE<-as.data.frame(blue)
 115         #print(BLUE)
 116         #BLUE = merge(x = BLUE, y = fixedeff, by="germplasmName", all=TRUE)
 117
 118         # compute adjusted blues
 119         p0 = predict.mmer(mixmodel, D="germplasmName") ##runs the prediction
 120         summary(p0)
 121         adj = p0$pvals                                        ##obtains the predictions
 122         adjustedBLUE = as.data.frame(adj)
 123         #print(paste("adj", adj))
 124         print(adjustedBLUE)
 125
 126     }
 127 }
 128
 129 # save result files
 130 #
 131 # for random effects: file.BLUPs and file.adjustedBLUPs
 132 # for fixed effects: file.BLUEs and file.adjustedBLUEs
 133 #
 134 if (genotypeEffectType=="random") {
 135     outfile_blup = paste(datafile, ".BLUPs", sep="");
 136     sink(outfile_blup)
 137     write.table(BLUP, quote=F, sep='\t', row.names=FALSE)
 138     sink();
 139     outfile_adjmeans = paste(datafile, ".adjustedBLUPs", sep="")
 140     sink(outfile_adjmeans)
 141     write.table(adjusted_means, quote=F , sep='\t', row.names=FALSE)
 142     sink();
 143
 144 } else {   #fixed
 145     outfile_blue = paste(datafile, ".BLUEs", sep="")
 146     sink(outfile_blue)
 147     write.table(BLUE, quote=F , sep='\t', row.names=FALSE)
 148     sink();
 149     outfile_adjblues = paste(datafile, ".adjustedBLUEs", sep="")
 150     sink(outfile_adjblues)
 151     write.table(adjustedBLUE, quote=F , sep='\t', row.names=FALSE)
 152     sink();
 153 }