R/spatial_modeling.R

   1 ################################################################################
   2 # Modeling Spatial Variation
   3 ################################################################################
   4
   5 # There are ten main steps to this protocol:
   6 # 1. Load the software needed.
   7 # 2. Declare user-supplied variables.
   8 # 3. Process the phenotypic data.
   9 # 4. Fit the mixed models in sommer.
  10 # 5. Format the information needed for output.
  11
  12
  13 ################################################################################
  14 # 1. Load software needed
  15 ################################################################################
  16
  17 library(sommer)
  18
  19
  20 ################################################################################
  21 # 2. Declare user-supplied variables.
  22 ################################################################################
  23 #Get Arguments
  24 args = commandArgs(trailingOnly=TRUE)
  25 if (length(args) < 2) {
  26   stop('Two or more arguments are required.')
  27 }
  28 phenotypeFile = args[1]
  29 traits = args[2]
  30
  31
  32 ################################################################################
  33 # 3. Process the phenotypic data.
  34 ################################################################################
  35 #read in the phenotypic data
  36 userPheno <- read.delim(phenotypeFile, header = TRUE, sep="\t", fill=TRUE)
  37
  38 #The user should be able to select their response variables from a drop-down menu
  39 #    of the column names of the userPheno object. Then, those strings should be passed
  40 #    to this vector, 'userResponse'.
  41
  42 userResponse <- unlist(strsplit(traits, split=",", fixed=T))
  43 userID <- "germplasmName"
  44 row <- "rowNumber"
  45 col <- "colNumber"
  46 userPheno$R <- as.factor(userPheno$rowNumber)
  47 userPheno$C <- as.factor(userPheno$colNumber)
  48
  49
  50 ################################################################################
  51 # 4. Fit the 2D Spline model in sommer
  52 ################################################################################
  53 # Make a list to save the models.
  54
  55 userModels <- list()
  56
  57 for(i in 1:length(userResponse)){
  58
  59    fixedArg <- paste(userResponse[i], " ~ ", "1 +", userID,")", sep = "")
  60
  61    randArg <- paste("~vs(", R, ")+vs(", C, ")+ spl2Da(",col,"," ,row ,")", sep = "")
  62
  63
  64    m2.sommer <- mmer(fixed = as.formula(fixedArg),
  65                      random = as.formula(randArg),
  66                      rcov= ~units,
  67                      data=userPheno, verbose = FALSE)
  68
  69
  70
  71
  72    userModels[[i]] <- m2.sommer
  73
  74 }
  75
  76
  77 ################################################################################
  78 # 5. Format the information needed for output.
  79 ################################################################################
  80
  81
  82 for(i in 1:length(userModels)){
  83
  84   m2.sommer <- userModels[[i]]
  85
  86   #Variance_Components <- summary(m2.sommer)$varcomp
  87
  88  # outputFile= paste(userID, " Spatial Variance Components", ".out", sep="")
  89
  90   write.csv(Variance_Components, outputFile)
  91
  92   res <- (randef(m2.sommer)$`u:variety`)
  93   BLUP <- as.data.frame(res)
  94
  95  # adj = coef(m2.sommer)$Trait
  96   outfile_blup = paste(phenotypeFile, ".BLUPs", sep="");
  97   write.table(BLUP, outfile_blup)
  98
  99 }