#' @title Graph PC Components
#'
#' @description This function identifies the regions that contribute the most to an individual principal component. It then graphs
#' the signal in those regions in each individual input brain, separated by phenotype
#' @param PCAobject this is the output of your prcomp() function. It is a list that includes "rotation"
#' @param LabeledMatrix this is a matrix whose rows are brains and columns are brain regions, with descriptive rownames and colnames
#' @param PCnum this is which principal component you want to examine (e.g. PC1, PC2). Enter a number.
#' @param num2graph this is the number of contributing brain regions you want in your graph. For example, if you want to graph the signal in
#' individual brains of the top 5 contributing brain regions to PC3, you would enter PCnum = 3, num2graph = 5.
#' @export
#' @examples


graph_PC_components <- function(PCAobject, LabeledMatrix, PCnum, num2graph) {

  #PCAobject is the output from performing PCA, a list containing "rotation"
  # LabeledMatrix contains the signal values for brain regions plus any relevant variables (pheno, geno)
  # PC num is the particular PC you wish to graph, e.g. PC1, PC2, etc. input as a number
  # num2graph is the number of brain regions from the PC you wish to graph. They will be sorted according to their contribution,
  # so if you enter "10" you will get the top 10 regions

  PCnum = as.numeric(PCnum)
  num2graph = as.numeric(num2graph)

  loading_scores_PC <- PCAobject$rotation[,PCnum]
  loading_scores_PC
  loading_scores_PC_ranked <- sort(loading_scores_PC, decreasing=TRUE)
  top_regions_PC <- names(loading_scores_PC_ranked[1:num2graph])

  PC_Values <- LabeledMatrix %>%
  select(geno, pheno, all_of(top_regions_PC))

  PC_Values <- pivot_longer(PC_Values, cols = c(top_regions_PC), names_to= "Region", values_to = "Signal")

  PC_Values <- PC_Values %>%
    mutate("Signal" = as.numeric(Signal))

  position <- position_jitter(0.1, 0)

  title <- paste("Signal in top ", num2graph, " regions of PC", PCnum, sep = "")

  ggplot(PC_Values, aes(x = pheno, y = Signal, color = pheno)) +
    geom_point(position = position, size = 3, alpha = 0.5) +
    geom_violin(alpha = 0.1, scale = "count") +
    facet_wrap(~Region, ncol = 5) +
    scale_color_manual(values = c("red", "blue")) + #if using mid too, include "blueviolet" in the middle here
    scale_y_continuous(expand = c(0,0)) +
    ggtitle(title) +
    theme(
      panel.background = element_blank(),
      axis.line = element_line(color = "grey40"),
      panel.grid.major.y = element_line(color = "grey90"),
      strip.text.x = element_text(size = 8)
    )

}