R/heritability/h2_blup_rscript.R

   1
   2
   3 library("rjson")
   4 library("dplyr")
   5 library("methods")
   6
   7
   8 ##### Get data #####
   9 args = commandArgs(trailingOnly = TRUE)
  10
  11 pheno <- read.table(args[1], sep = "\t", header = TRUE)
  12
  13 study_trait <- args[2]
  14 h2File <- args[3]
  15 h2CsvFile <- args[4]
  16 errorFile <- args[5]
  17
  18 errorMessages <- c()
  19
  20 study_trait <- gsub("\\."," ",study_trait)
  21 cat("study trait is ", study_trait,"\n")
  22
  23
  24 names <- colnames(pheno)
  25 allTraits <- grep("CO_", names, value = TRUE)
  26 library(dplyr)
  27 pheno <- pheno[,colnames(pheno) %in% c("locationName", "germplasmName", "studyYear", "blockNumber","replicate", allTraits)]
  28 names <- colnames(pheno)
  29 new_names <- gsub(".CO_.*","", names)
  30 new_names <- gsub("\\."," ",new_names)
  31 colnames(pheno) <- new_names
  32
  33 #Calculating missing data
  34 missingData <- data.frame(missingData = tapply(pheno[,6], pheno$germplasmName, function(x){
  35   return(length(which(is.na(x))))
  36 }))
  37
  38 missingData <- tibble::rownames_to_column(missingData, "germplasmName")
  39 md <- sum(missingData$missingData)/nrow(pheno)
  40
  41 #Removing trait with more than 60% of missing data
  42 if(md[1] > 0.6){
  43   message1 <- c("Please, check your dataset! There are more than 60% of missing data for selected trait.")
  44 }else{
  45   message1 <- c()
  46 }
  47
  48
  49 #Removing non numeric data
  50 isNumeric <- is.numeric(pheno[,6])
  51 if(isNumeric == "FALSE"){
  52   message2 <- c("Please, check your dataset! The selected trait is not numeric.")
  53 }else{
  54   message2 <- c()
  55 }
  56
  57 #checkning number of locations
  58 szloc <- length(unique(pheno$locationName))
  59 szreps <- length(unique(pheno$replicate))
  60 szyr <- length(unique(pheno$studyYear))
  61
  62 library(lmerTest)
  63 # Still need check temp data to ensure wright dimension
  64 an.error.occured <- FALSE
  65
  66 H2 = c()
  67 her = c()
  68 Vg = c()
  69 Ve = c()
  70 Vres = c()
  71 resp_var = c()
  72 j = 1
  73 tryCatch({ for(i in 6:ncol(pheno)){
  74     outcome = colnames(pheno)[i]
  75     if (szreps > 1){
  76       if (szloc == 1){
  77         if (szyr == 1){
  78           model <- lmer(get(outcome)~(1|germplasmName)+replicate,
  79                         na.action = na.exclude,
  80                         data=pheno)
  81           variance = as.data.frame(VarCorr(model))
  82           gvar = variance [1,"vcov"]
  83           envar = 0
  84           resvar = variance [2, "vcov"]
  85         }else{
  86           model <- lmer(get(outcome) ~ (1|germplasmName) + replicate + studyYear,
  87                         na.action = na.exclude,
  88                         data=pheno)
  89           variance = as.data.frame(VarCorr(model))
  90           gvar = variance [1,"vcov"]
  91           envar = 0
  92           resvar = variance [2, "vcov"]
  93         }
  94       }else if (szloc > 1) {
  95         if (szyr == 1){
  96           model <- lmer(get(outcome) ~ (1|germplasmName) + replicate + (1|locationName),
  97                         na.action = na.exclude,
  98                         data=pheno)
  99           variance = as.data.frame(VarCorr(model))
 100           gvar = variance [1,"vcov"]
 101           envar = variance [2, "vcov"]
 102           resvar = variance [3, "vcov"]
 103         }else{
 104           model <- lmer(get(outcome) ~ (1|germplasmName) + replicate + (1|locationName) + studyYear,
 105                         na.action = na.exclude,
 106                         data=pheno)
 107           variance = as.data.frame(VarCorr(model))
 108           gvar = variance [1,"vcov"]
 109           envar = variance [2, "vcov"]
 110           resvar = variance [3, "vcov"]
 111         }
 112       }
 113     }else if (szreps == 1){
 114       if (szloc ==1){
 115         if (szyr == 1){
 116           model <- lmer(get(outcome)~(1|germplasmName) + blockNumber,
 117                         na.action = na.exclude,
 118                         data=pheno)
 119           variance = as.data.frame(VarCorr(model))
 120           gvar = variance [1,"vcov"]
 121           envar = 0
 122           resvar = variance [2, "vcov"]
 123         }else{
 124           model <- lmer(get(outcome) ~ (1|germplasmName) + studyYear + blockNumber,
 125                         na.action = na.exclude,
 126                         data=pheno)
 127           variance = as.data.frame(VarCorr(model))
 128           gvar = variance [1,"vcov"]
 129           envar = 0
 130           resvar = variance [2, "vcov"]
 131         }
 132       }else if (szloc > 1){
 133         if (szyr ==1){
 134           model <- lmer(get(outcome)~(1|germplasmName)+ (1|locationName) +  blockNumber,
 135                         na.action = na.exclude,
 136                         data=pheno)
 137           variance = as.data.frame(VarCorr(model))
 138           gvar = variance [1,"vcov"]
 139           envar = variance [2, "vcov"]
 140           resvar = variance [3, "vcov"]
 141         }else{
 142           model <- lmer(get(outcome) ~ (1|germplasmName) + studyYear + (1|locationName) + blockNumber,
 143                         na.action = na.exclude,
 144                         data=pheno)
 145           variance = as.data.frame(VarCorr(model))
 146           gvar = variance [1,"vcov"]
 147           envar = variance [2, "vcov"]
 148           resvar = variance [3, "vcov"]
 149         }
 150       }
 151     }
 152
 153     H2 = append(H2, gvar/ (gvar + (envar) + (resvar)))
 154     H2nw = format(round(H2[j], 4), nsmall = 4)
 155     her = append(her, round(as.numeric(H2nw), digits =3))
 156     Vg = append(Vg, round(as.numeric(gvar), digits = 3))
 157     Ve = append(Ve, round(as.numeric(envar), digits = 2))
 158     Vres = append(Vres, round(as.numeric(resvar), digits = 3))
 159     resp_var = append(resp_var, colnames(pheno)[i])
 160     j=j+1
 161   }
 162
 163 }, error = function(e) {
 164   an.error.occured <<- TRUE
 165   errorMessages <<- c(errorMessages, as.character(e))
 166 })
 167
 168
 169
 170 #Prepare information to export data
 171 tryCatch({
 172
 173   Heritability = data.frame(trait = resp_var,
 174                           Hert = her,
 175                           Vg = Vg,
 176                           Ve = Ve,
 177                           Vres = Vres
 178                           )
 179
 180   Heritability = na.omit(Heritability)
 181
 182   h2_json <- jsonlite::toJSON(Heritability)
 183   jsonlite::write_json(h2_json, h2File)
 184   write.csv(Heritability, file = h2CsvFile)
 185
 186 }, error = function(e) {
 187   an.error.occured <<- TRUE
 188   errorMessages <<- c(errorMessages, as.character(e))
 189 })
 190
 191 errorMessages <- c(errorMessages, as.character(message1))
 192 errorMessages <- c(errorMessages, as.character(message2))
 193
 194 cat("Was there an error? ", an.error.occured,"\n")
 195 if ( length(errorMessages) > 0 ) {
 196   print(sprintf("Writing Error Messages to file: %s", errorFile))
 197   print(errorMessages)
 198   write(errorMessages, errorFile)
 199 }
 200
 201
 202 #-------------------------------------------------------------------------