R/gl.report.factorloadings.r

In dartR.base: Analysing 'SNP' and 'Silicodart' Data - Basic Functions

#' @name gl.report.factorloadings
#' @title Reports factor loadings for a PCA or PCoA
#' @family matched reports

#' @description
#' Extracts the factor loadings from a glPCA object (generated by gl.pcoa) and plots
#' their distribuion.

#' @param pca Name of the glPCA object containing factor loadings [required].
#' @param axis Axis in the ordination used to display the factor loadings [default 1]
#' @param n.display Number of loci for which to display factorloadings [default 15]
#' @param plot.display If TRUE, resultant plots are displayed in the plot window
#' [default TRUE].
#' @param plot.theme Theme for the plot. See Details for options
#' [default theme_dartR()].
#' @param plot.colors List of two color names for the borders and fill of the
#'  plots [default c("#2171B5","#6BAED6")].
#' @param plot.dir Directory to save the plot RDS files [default as specified 
#' by the global working directory or tempdir()]
#' @param plot.file Name for the RDS binary file to save (base name only, exclude extension) [default NULL]
#' @param bins Number of bins to display in histograms [default 25].
#' @param verbose Verbosity: 0, silent or fatal errors; 1, begin and end; 2,
#' progress log; 3, progress and results summary; 5, full report
#'  [default NULL, unless specified using gl.set.verbosity]
#' @param ... Parameters passed to function \link[ggplot2]{ggsave}, 
#'  such as width and height, when the ggplot is to be saved.

#' @details 
#' The function extracts the factor loadings for a given axis from a PCA object 
#' generated by gl.pcoa and plots their magnitudes. Useful for identifying loci
#' that load high for a given axis. 
#' 
#'  A color vector can be obtained with gl.select.colors() and then passed to the function
#'  with the plot.colors parameter.
#'  
#' Themes can be obtained from in \itemize{
#'  \item \url{https://ggplot2.tidyverse.org/reference/ggtheme.html} and \item
#'  \url{https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/}
#'  }
#'   If a plot.file is given, the ggplot arising from this function is saved as an "RDS" 
#' binary file using saveRDS(); can be reloaded with readRDS(). A file name must be 
#' specified for the plot to be saved.

#'  If a plot directory (plot.dir) is specified, the ggplot binary is saved to that
#'  directory; otherwise to the tempdir(). 
#'  
#' @author Custodian: Arthur Georges -- Post to
#' \url{https://groups.google.com/d/forum/dartr}
#' 
#' @examples
#' pca <- gl.pcoa(testset.gl)
#' gl.report.factorloadings(pca = pca)
#' 
#' @export
#' @return The unchanged genlight object
#' 
gl.report.factorloadings <- function(pca,
                            axis=1,
                            n.display=15,
                            plot.display=TRUE,
                            plot.theme = theme_dartR(),
                            plot.colors = NULL,
                            plot.file=NULL,
                            plot.dir=NULL,
                            bins=25,
                            verbose = NULL,
                            ...) {

  # SET VERBOSITY
  verbose <- gl.check.verbosity(verbose)
  if(verbose==0){plot.display <- FALSE}
  
  # SET WORKING DIRECTORY
  plot.dir <- gl.check.wd(plot.dir,verbose=0)
  
  # SET COLOURS
  if(is.null(plot.colors)){
    plot.colors <- c("#2171B5", "#6BAED6")
  } else {
    if(length(plot.colors) > 2){
      if(verbose >= 2){cat(warn("  More than 2 colors specified, only the first 2 are used\n"))}
      plot.colors <- plot.colors[1:2]
    }
  }
  
  # FLAG SCRIPT START
  funname <- match.call()[[1]]
  utils.flag.start(func = funname,
                   build = "v.2023.3",
                   verbose = verbose)
  
  # CHECK DATATYPE
  datatype <- class(pca)
  if(datatype != "glPca"){
    cat(error("To report factor loadings, require a glPca object\n"))
    stop()
  } else {
    if(verbose >= 2){cat(report("  Reading a glPca object\n"))}
  }

  # DO THE JOB
  
  # Pull the factor loadings into a dataframe
  factor.loadings <- data.frame(pca$loadings[,axis])
  df <- cbind(rownames(factor.loadings),factor.loadings[,1])
  df <- data.frame(df)
  colnames(df) <- c("locus","loading")
  df$loading <- as.numeric(df$loading)
  
  # Create a second dataframe with the top n.display rows
  cat(report("  Reporting factor loadings for the top",n.display,"loci for axis",axis,"\n"))
  tmp <- df
  tmp$abs <- abs(as.numeric(df$loading))
  tmp <- tmp[order(tmp$abs,decreasing=TRUE),]
  tmp$abs <- NULL
  tmp$loading <- round(as.numeric(tmp$loading),4)
  
  # Display the top n.display rows
  print(tmp[1:n.display, ,drop=FALSE])
  
# Prepare the plots
# get title for plots
  title1 <- paste("Distribution of factor loadings for axis",axis)

# Calculate minimum and maximum graph cutoffs for callrate
min <- min(df$loading)
max <- max(df$loading)

loading <- NULL
# Boxplot
p1 <-
  ggplot(df, aes(y = loading)) + 
  geom_boxplot(color = plot.colors[1], fill = plot.colors[2]) +
  coord_flip() +
  plot.theme + 
  xlim(range = c(-1, 1)) +
  ylim(min, max) + 
  ylab(" ") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  ggtitle(title1)

# Histogram
p2 <-
  ggplot(df, aes(x = loading)) + 
  geom_histogram(bins = bins, color = plot.colors[1],fill = plot.colors[2]) +
  coord_cartesian(xlim = c(min, max)) + 
  xlab("Factor Loading") +
  ylab("Count") +
  plot.theme

# using package patchwork
p3 <- (p1 / p2) + plot_layout(heights = c(1, 4))
if (plot.display) {print(p3)}

  # Optionally save the plot
  
  if(!is.null(plot.file)){
    tmp <- utils.plot.save(p3,
                           dir=plot.dir,
                           file=plot.file,
                           verbose=verbose)
  }
  
  # FLAG SCRIPT END 
  
  if (verbose >= 1) {
    cat(report("Completed:", funname, "\n"))
  }

  # RETURN
  invisible(factor.loadings)
}