R/gs_plotsHCL.R

#' Create a comparison plot for phenotypes
#' Function that creates boxplots comparing HCL color components phenotypes (e.g. Heritability or Prediction_ability).
#'
#' @param HCLTable Table with information from cross validation: GEBVs for each fold and each repetition.
#' @return None. Plot is saved to a PDF file
#' @import ggplot2
#' @import dplyr 
#' @export
gs_plotsHCL <- function (HCLFilename, outputDir, columnName, TITLE) {
	HCLTable = read.csv (HCLFilename)
	message (">>> Creating output plots and tables for best model ...")

	# Summary plot boxplot best models for traits
	meansTable = HCLTable %>% dplyr::group_by (Traits, Models) %>% 
		summarize (Means=mean(get(columnName)), .groups="keep") 

	# Get best model for each trait
	bestTable = meansTable %>% group_by (Traits) %>% slice_max (Means) %>% slice_head

	# Create table with correlations from best models
	getBestCorr <- function (trait, model) 
		return (dplyr::filter (HCLTable, Traits==trait, Models==model))
	bestCorrs = do.call (rbind, mapply (getBestCorr, bestTable$Traits, bestTable$Models, SIMPLIFY=F))

	traitColors = unlist (lapply(1:8, function(x) rep(x,3))) 
	nTraits     = length (unique (bestCorrs$Traits))
	traitColors = traitColors [1: nTraits] 
	outPlotname   = sprintf ("%s/out-%s-HCL-Comparison-PLOT.pdf", outputDir, columnName)
	ggplot (bestCorrs, aes(x=HCL_component, y=get(columnName))) + ylim (0,1) +
		    geom_boxplot (alpha=0.3, fill=traitColors) + 
			theme(axis.text.x=element_text(angle=0, hjust=1)) + 
			labs (title=TITLE, x="HCL_component", y=columnName) + 
			stat_summary(geom='text', label=round (bestTable$Means,2) , fun=max, vjust = -1, size=3, col="blue") +
			stat_summary(geom='text', label=bestTable$Models, fun=min, vjust = 1.5, angle=0,size=3, col="blue") + 
			facet_wrap (~Prefix, ncol=8) 
	ggsave (outPlotname, width=11)
}

#-------------------------------------------------------------
#-------------------------------------------------------------
main <- function () {
	library (dplyr)
	library (ggplot2)
	args = commandArgs (trailingOnly=T)
	resultsDir = args[1]
	gs_outputs (resultsDir)
}
luisgarreta/ppGS documentation built on Jan. 17, 2022, 9:23 p.m.