R/solGS/genetic_gain.r

   1  #SNOPSIS
   2
   3  #visualizes genetic gain using boxplot.
   4
   5  #AUTHOR
   6  # Isaak Y Tecle (iyt2@cornell.edu)
   7
   8
   9 options(echo = FALSE)
  10
  11 library(data.table)
  12 #library(phenoAnalysis)
  13 library(dplyr)
  14 library(methods)
  15 library(ggplot2)
  16 #library(plotly)
  17
  18
  19
  20 allArgs <- commandArgs()
  21
  22 outputFiles <- scan(grep("output_files", allArgs, value = TRUE),
  23                     what = "character")
  24
  25 inputFiles  <- scan(grep("input_files", allArgs, value = TRUE),
  26                    what = "character")
  27
  28 trGebvFiles <- grep("training", inputFiles, value = TRUE)
  29
  30 message('training file: ', trGebvFiles)
  31
  32 slGebvFiles <- grep("selection", inputFiles, value = TRUE)
  33
  34 message('selection file: ', slGebvFiles)
  35
  36 boxplotFile <-  grep("genetic_gain_plot", outputFiles, value = TRUE)
  37 message('boxplot file: ', boxplotFile)
  38
  39 plotDataFile <-  grep("genetic_gain_data", outputFiles, value = TRUE)
  40 message('boxplot file: ', plotDataFile)
  41
  42 trGebv   <- c()
  43 gebvsCol <- 'GEBVs'
  44
  45 for (trGebvFile in trGebvFiles) {
  46     gebv <- data.frame(fread(trGebvFile, header = TRUE))
  47     trait <- names(gebv)[2]
  48     colnames(gebv)[2] <- gebvsCol
  49     gebv$trait <- trait
  50     trGebv <- bind_rows(trGebv, gebv)
  51 }
  52
  53 trGebv$pop <- 'training'
  54 slGebv <- c()
  55
  56 for (slGebvFile in slGebvFiles) {
  57     gebv <- data.frame(fread(slGebvFile, header = TRUE))
  58     trait <- names(gebv)[2]
  59     colnames(gebv)[2] <- gebvsCol
  60     gebv$trait <- trait
  61     slGebv <- bind_rows(slGebv, gebv)
  62 }
  63
  64 slGebv$pop <- 'selection'
  65
  66 boxplotData <- bind_rows(trGebv, slGebv)
  67 boxplotData$pop <- as.factor(boxplotData$pop)
  68 boxplotData$trait <- as.factor(boxplotData$trait)
  69
  70 boxplotData$pop <- with(boxplotData, relevel(pop, "training"))
  71
  72 ## pop       <- 'pop'
  73 ## training  <- 'training'
  74 ## selection <- 'selection'
  75 ## trait     <- 'trait'
  76 ## Gebvs     <- 'GEBVs'
  77
  78 bp <- ggplot(boxplotData,
  79              aes(y=GEBVs, x=pop, fill=pop)) +
  80              geom_boxplot(width=0.4) +
  81              stat_summary(geom="text", fun.y=quantile, size=5,
  82              aes(label=sprintf("%1.3f", ..y..), color=pop),
  83              position=position_nudge(x=0.35)) +
  84              theme_bw() +
  85              facet_wrap(~trait, ncol=2, scales="free") +
  86              theme(legend.position="none",
  87                axis.text=element_text(color='blue', size=12,  face='bold'),
  88                axis.title.y=element_text(color='blue'),
  89                axis.title.x=element_blank(),
  90                strip.text.x=element_text(color='blue', face='bold'))
  91
  92 wid <- 480 * if(length(trGebvFiles) > 1) 2 else 1;
  93
  94 png(boxplotFile, width=wid)
  95 bp
  96 dev.off()
  97
  98 if (length(plotDataFile) != 0 ) {
  99     fwrite(boxplotData,
 100        file      = plotDataFile,
 101        sep       = "\t",
 102        row.names = FALSE,
 103        quote     = FALSE,
 104        )
 105
 106 }
 107
 108 q(save = "no", runLast = FALSE)