R/stability/ammi_script.R

   1
   2
   3
   4 library("methods")
   5 library("dplyr")
   6 library("gridExtra")
   7 library("agricolae")
   8 library("gge")
   9
  10 ##### Get data #####
  11 args = commandArgs(trailingOnly = TRUE)
  12
  13
  14 pheno <- read.table(args[1], sep = "\t", header = TRUE)
  15 study_trait <- args[2]
  16 cat("Study trait is ", study_trait[1])
  17 figure1_file_name <- args[3]
  18 figure2_file_name <- args[4]
  19 AMMIFile <- args[5]
  20 method <- args[6]
  21
  22 #Making names standard
  23 names <- colnames(pheno)
  24 new_names <- gsub(".CO.*","", names)
  25 colnames(pheno) <- new_names
  26 cat(colnames(pheno),"\n")
  27
  28 #Finding which column is the study trait
  29 for (i in 1:ncol(pheno)){
  30         a = noquote(colnames(pheno[i]))
  31         b = study_trait
  32         cat("COL: ", a, " vs study_trait: ", b,"\n")
  33         if (a==b){
  34                 print(a)
  35                 col = i
  36                 i = ncol(pheno)
  37                 break
  38         }else{
  39                 cat("working ",i,"\n")
  40                 i=i+1
  41         }
  42 }
  43
  44 env <-as.factor(pheno$locationName)
  45 gen <-as.factor(pheno$germplasmName)
  46 rep <-as.factor(pheno$replicate)
  47 # trait <-as.numeric(pheno[,col])
  48
  49 message<-"The analysis could not be completed. Please set your dataset with more than 1 location."
  50 message2<-""
  51 locations <- unique(pheno$locationDbId)
  52 acc <- unique(pheno$germplasmDbId)
  53 subGen <- unique(subset(pheno, select=c(germplasmDbId, germplasmName)))
  54
  55
  56 if (! length(locations)>1){
  57
  58         png(AMMIFile, height = 130, width=800)
  59         z<-tableGrob(message)
  60         grid.arrange(z)
  61         dev.off()
  62
  63         sub1 <- unique(subset(pheno, locationDbId == locations[1], select=c(locationDbId, locationName)))
  64         png(figure1_file_name, height = 100, width=800)
  65         p<-tableGrob(sub1)
  66         grid.arrange(p)
  67         dev.off()
  68
  69         acc <- unique(pheno$germplasmDbId)
  70         png(figure2_file_name, height = (21*length(acc)), width = 800)
  71         sub2 <- unique(subset(pheno, locationDbId == locations[1], select=c(germplasmDbId, germplasmName)))
  72         q<-tableGrob(sub2)
  73         grid.arrange(q)
  74         dev.off()
  75 } else {
  76         cat("Starting stability analysis...","\n")
  77 }
  78
  79 cat("create model", "\n")
  80 summary(pheno)
  81 summary(env)
  82 summary(gen)
  83 summary(rep)
  84 model<- with(pheno,AMMI(env, gen, rep, pheno[,col], console=FALSE))
  85
  86 cat("marker 1", "\n")
  87 anova <-format(round(model$ANOVA, 3))
  88 cat("marker 2", "\n")
  89 analysis <- model$analysis
  90 cat("marker 3", "\n")
  91 anova
  92 cat("marker 4", "\n")
  93 analysis
  94
  95 cat("DONE WITH MODEL", "\n")
  96
  97 png(AMMIFile, height=130, width=800)
  98 p<-tableGrob(anova)
  99 grid.arrange(p)
 100 dev.off()
 101
 102
 103 if(method=="ammi"){
 104
 105         cat("running ammi", "\n")
 106         # Biplot and Triplot
 107         png(figure1_file_name,height=400)
 108         plot(model, first=0,second=1, number=TRUE, xlab = study_trait)
 109         dev.off()
 110
 111         tam <- nrow(model$analysis)
 112         if (tam >2){
 113           png(figure2_file_name, height= 300)
 114           plot(model, type=2, number=TRUE, xlab = study_trait)
 115           dev.off()
 116         }else{
 117           png(figure2_file_name, height=5, width=5)
 118           y <- tableGrob(message2)
 119           grid.arrange(y)
 120           dev.off()
 121         }
 122
 123         }else if( method=="gge"){
 124                 dat1 <- pheno %>% select(germplasmName, locationName, trait=study_trait)
 125                 head(dat1)
 126
 127                 model1 <- gge(dat1, trait~germplasmName*locationName, scale=FALSE)
 128
 129                 name1 = paste(study_trait,"- GGE Biplot", sep=" ")
 130
 131                 png(figure1_file_name, height=400, width=500)
 132                 biplot(model1, main=name1, flip=c(1,0), origin=0, hull=TRUE)
 133                 dev.off()
 134
 135                 model2 <- gge(dat1, trait~germplasmName*locationName, scale=TRUE)
 136                 png(figure2_file_name,height=400, width=500)
 137                 biplot(model2, main=name1, flip=c(1,1), origin=0)
 138                 dev.off()
 139         }