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         if (a==b){
  33                 print(a)
  34                 col = i
  35                 i = ncol(pheno)
  36         }else{
  37                 cat("working ",i,"\n")
  38                 i=i+1
  39         }
  40 }
  41
  42
  43 env <-as.factor(pheno$locationName)
  44 gen <-as.factor(pheno$germplasmName)
  45 rep <-as.factor(pheno$replicate)
  46 # trait <-as.numeric(pheno[,col])
  47
  48 message<-"The analysis could not be completed. Please set your dataset with more than 1 location."
  49 message2<-""
  50 locations <- unique(pheno$locationDbId)
  51 acc <- unique(pheno$germplasmDbId)
  52 subGen <- unique(subset(pheno, select=c(germplasmDbId, germplasmName)))
  53
  54 if (! length(locations)>1){
  55
  56         png(AMMIFile, height = 130, width=800)
  57         z<-tableGrob(message)
  58         grid.arrange(z)
  59         dev.off()
  60
  61         sub1 <- unique(subset(pheno, locationDbId == locations[1], select=c(locationDbId, locationName)))
  62         png(figure1_file_name, height = 100, width=800)
  63         p<-tableGrob(sub1)
  64         grid.arrange(p)
  65         dev.off()
  66
  67         acc <- unique(pheno$germplasmDbId)
  68         png(figure2_file_name, height = (21*length(acc)), width = 800)
  69         sub2 <- unique(subset(pheno, locationDbId == locations[1], select=c(germplasmDbId, germplasmName)))
  70         q<-tableGrob(sub2)
  71         grid.arrange(q)
  72         dev.off()
  73 } else {
  74         cat("Starting stability analysis...","\n")
  75 }
  76
  77 model<- with(pheno,AMMI(env, gen, rep, pheno[,col], console=FALSE))
  78
  79 anova <-format(round(model$ANOVA, 3))
  80 analysis <- model$analysis
  81 anova
  82 analysis
  83
  84 png(AMMIFile, height=130, width=800)
  85 p<-tableGrob(anova)
  86 grid.arrange(p)
  87 dev.off()
  88
  89
  90 if(method=="ammi"){
  91
  92         # Biplot and Triplot
  93         png(figure1_file_name,height=400)
  94         plot(model, first=0,second=1, number=TRUE, xlab = study_trait)
  95         dev.off()
  96
  97         tam <- nrow(model$analysis)
  98         if (tam >2){
  99           png(figure2_file_name, height= 300)
 100           plot(model, type=2, number=TRUE, xlab = study_trait)
 101           dev.off()
 102         }else{
 103           png(figure2_file_name, height=5, width=5)
 104           y <- tableGrob(message2)
 105           grid.arrange(y)
 106           dev.off()
 107         }
 108
 109         }else if( method=="gge"){
 110                 dat1 <- pheno %>% select(germplasmName, locationName, trait=study_trait)
 111                 head(dat1)
 112
 113                 model1 <- gge(dat1, trait~germplasmName*locationName, scale=FALSE)
 114
 115                 name1 = paste(study_trait,"- GGE Biplot", sep=" ")
 116
 117                 png(figure1_file_name, height=400, width=500)
 118                 biplot(model1, main=name1, flip=c(1,0), origin=0, hull=TRUE)
 119                 dev.off()
 120
 121                 model2 <- gge(dat1, trait~germplasmName*locationName, scale=TRUE)
 122                 png(figure2_file_name,height=400, width=500)
 123                 biplot(model2, main=name1, flip=c(1,1), origin=0)
 124                 dev.off()
 125         }