#' @name compFGA

#' @title Comparison of fraction genome altered

#' @description This function calculates Fraction Genome Altered (FGA), Fraction Genome Gained (FGG), and Fraction Genome Lost (FGL) seperately, and compares them among curent subtypes identified from multi-omics integrative clustering algorithms.

#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.

#' @param segment A data frame containing segmented copy number and columns must exactly include the following elements: c('sample','chrom','start','end','value'). Column of `value` should be segments value when \code{iscopynumber = FALSE} but copy-number value when \code{iscopynumber = TRUE}. Copy-number will be converted to segments by log2(copy-number/2).

#' @param iscopynumber A logical value to indicate if the fifth column of segment input is copy-number. If segment file derived from CNV calling provides copy number instead of segment_mean value, this argument must be switched to TRUE. FALSE by default.

#' @param cnathreshold A numeric value to indicate the cutoff for identifying copy-number gain or loss. 0.2 by default.

#' @param test.method A string value to indicate the method for statistical testing. Allowed values contain c('nonparametric', 'parametric'); nonparametric means two-sample wilcoxon rank sum test for two subtypes and Kruskal-Wallis rank sum test for multiple subtypes; parametric means two-sample t-test when only two subtypes are identified, and anova for multiple subtypes comparison; "nonparametric" by default.

#' @param barcolor A string vector to indicate the mapping color for bars of FGA, FGG and FGL.

#' @param clust.col A string vector storing colors for each subtype.

#' @param fig.path A string value to indicate the output path for storing the barplot.

#' @param fig.name A string value to indicate the name of the barplot.

#' @param width A numeric value to indicate the width of barplot.

#' @param height A numeric value to indicate the height of barplot.

#' @return A list contains the following components:

#'

#'         \code{summary}           a table summarizing the measurements of FGA, FGG, and FGL per sample

#'

#'         \code{FGA.p.value}       a nominal p value quantifying the difference of FGA among current subtypes

#'

#'         \code{pairwise.FGA.test} a pairwise BH adjusted p value matrix for multiple comparisons of FGA if more than 2 subtypes were identified

#'

#'         \code{FGG.p.value}       a nominal p value quantifying the difference of FGG among current subtypes

#'

#'         \code{pairwise.FGG.test} a pairwise BH adjusted p value matrix for multiple comparisons of FGG if more than 2 subtypes were identified

#'

#'         \code{FGL.p.value}       a nominal p value quantifying the difference of FGL among current subtypes

#'

#'         \code{pairwise.FGL.test} a pairwise BH adjusted p value matrix for multiple comparisons of FGL if more than 2 subtypes were identified

#'

#'         \code{test.method}       a string value indicating the statistical testing method to calculate p values

#'

#' @export

#' @import ggplot2

#' @import patchwork

#' @importFrom dplyr group_by summarize %>%

#' @importFrom ggpubr stat_compare_means

#' @examples # There is no example and please refer to vignette.

#' @references Cerami E, Gao J, Dogrusoz U, et al. (2012). The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data. Cancer Discov, 2(5):401-404.

#'

#' Gao J, Aksoy B A, Dogrusoz U, et al. (2013). Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal, 6(269):pl1-pl1.

compFGA <- function(moic.res     = NULL,

                    segment      = NULL,

                    iscopynumber = FALSE,

                    cnathreshold = 0.2,

                    test.method  = "nonparametric",

                    barcolor     = c("#008B8A", "#F2042C", "#21498D"),

                    clust.col    = c("#2EC4B6","#E71D36","#FF9F1C","#BDD5EA","#FFA5AB","#011627","#023E8A","#9D4EDD"),

                    fig.path     = getwd(),

                    fig.name     = NULL,

                    width        = 8,

                    height       = 4) {

  # check arguments

  if(!all(is.element(c("sample","chrom","start","end","value"), colnames(segment)))) {

    stop("segment data must have the following columns: sample, chrom, start, end, value.")

  }

  if(iscopynumber) {

    segment$value <- log2(segment$value/2) # convert copy-number to segment-mean value

  }

  comsam <- intersect(moic.res$clust.res$samID, unique(segment[,1]))

  # check data

  if(length(comsam) == nrow(moic.res$clust.res)) {

    message("--all samples matched.")

  } else {

    message(paste0("--",(nrow(moic.res$clust.res)-length(comsam))," samples mismatched from current subtypes."))

  }

  if(!is.element(test.method, c("nonparametric","parametric"))) {

    stop("test.method can be one of nonparametric or parametric.")

  }

  clust.res <- moic.res$clust.res[comsam, , drop = FALSE]

  segment <- segment[which(segment$sample %in% comsam),]

  n.moic <- length(unique(clust.res$clust))

  # data process

  segment$bases <- segment$end - segment$start

  # calculate FGA, FGG and FGL

  display.progress = function (index, totalN, breakN=20) {

    if ( index %% ceiling(totalN/breakN)  == 0  ) {

      cat(paste(round(index*100/totalN), "% ", sep = ""))

    }

  }

  std <- function(x, na.rm = TRUE) {

    if(na.rm) {

      x <- as.numeric(na.omit(x))

      sd(x)/sqrt(length(x))

    } else {sd(x)/sqrt(length(x))}

  }

  outTab <- data.frame()

  for (i in 1:length(unique(segment$sample))) {

    display.progress(index = i, totalN = length(unique(segment$sample)))

    tmp <- segment[segment$sample == names(table(segment$sample))[i],]

    if (length(tmp[abs(tmp$value) > cnathreshold,"bases"][6]) == 0) {

      FGA = 0

    } else {

      FGA = sum(tmp[abs(tmp$value) > cnathreshold,"bases"]) / sum(tmp[,"bases"])

    }

    if (length(tmp[tmp$value > cnathreshold,"bases"][6]) == 0) {

      FGG = 0

    } else {

      FGG = sum(tmp[tmp$value > cnathreshold,"bases"]) / sum(tmp[,"bases"])

    }

    if (length(tmp[tmp$value < (-cnathreshold),"bases"][6]) == 0) {

      FGL = 0

    } else {

      FGL = sum(tmp[tmp$value < (-cnathreshold),"bases"]) / sum(tmp[,"bases"])

    }

    tmp <- data.frame(samID            = names(table(segment$sample))[i],

                      FGA              = FGA,

                      FGG              = FGG,

                      FGL              = FGL,

                      stringsAsFactors = FALSE)

    outTab <- rbind.data.frame(outTab, tmp, stringsAsFactors = FALSE)

  }

  outTab$Subtype <- paste0("CS", clust.res[outTab$samID, "clust"])

  # calculate mean and se

  summaryFGA  <- outTab %>% group_by(Subtype) %>% dplyr::summarize(mean = mean(FGA, na.rm = TRUE), se = std(FGA, na.rm = TRUE))

  summaryFGG  <- outTab %>% group_by(Subtype) %>% dplyr::summarize(mean = mean(FGG, na.rm = TRUE), se = std(FGG, na.rm = TRUE))

  summaryFGL  <- outTab %>% group_by(Subtype) %>% dplyr::summarize(mean = mean(FGL, na.rm = TRUE), se = std(FGL, na.rm = TRUE))

  summaryFGGL <- data.frame(rbind.data.frame(summaryFGG,summaryFGL),class = rep(c("FGG","FGL"),c(nrow(summaryFGG),nrow(summaryFGL))),stringsAsFactors = FALSE)

  # statistical testing

  # generate boxviolin plot with statistical testing

  if(n.moic == 2 & test.method == "nonparametric") {

    statistic <- "wilcox.test"

    FGA.test  <- wilcox.test(outTab$FGA ~ outTab$Subtype)$p.value

    FGG.test  <- wilcox.test(outTab$FGG ~ outTab$Subtype)$p.value

    FGL.test  <- wilcox.test(outTab$FGL ~ outTab$Subtype)$p.value

  }

  if(n.moic == 2 & test.method == "parametric") {

    statistic <- "t.test"

    FGA.test  <- t.test(outTab$FGA ~ outTab$Subtype)$p.value

    FGG.test  <- t.test(outTab$FGG ~ outTab$Subtype)$p.value

    FGL.test  <- t.test(outTab$FGL ~ outTab$Subtype)$p.value

  }

  if(n.moic > 2 & test.method == "nonparametric") {

    statistic <- "kruskal.test"

    FGA.test  <- kruskal.test(outTab$FGA ~ outTab$Subtype)$p.value

    FGG.test  <- kruskal.test(outTab$FGG ~ outTab$Subtype)$p.value

    FGL.test  <- kruskal.test(outTab$FGL ~ outTab$Subtype)$p.value

    pairwise.FGA.test <- pairwise.wilcox.test(outTab$FGA,outTab$Subtype,p.adjust.method = "BH")

    pairwise.FGG.test <- pairwise.wilcox.test(outTab$FGG,outTab$Subtype,p.adjust.method = "BH")

    pairwise.FGL.test <- pairwise.wilcox.test(outTab$FGL,outTab$Subtype,p.adjust.method = "BH")

  }

  if(n.moic > 2 & test.method == "parametric") {

    statistic <- "anova"

    FGA.test  <- summary(aov(outTab$FGA ~ outTab$Subtype))[[1]][["Pr(>F)"]][1]

    FGG.test  <- summary(aov(outTab$FGG ~ outTab$Subtype))[[1]][["Pr(>F)"]][1]

    FGL.test  <- summary(aov(outTab$FGL ~ outTab$Subtype))[[1]][["Pr(>F)"]][1]

    pairwise.FGA.test <- pairwise.t.test(outTab$FGA,outTab$Subtype,p.adjust.method = "BH")

    pairwise.FGG.test <- pairwise.t.test(outTab$FGG,outTab$Subtype,p.adjust.method = "BH")

    pairwise.FGL.test <- pairwise.t.test(outTab$FGL,outTab$Subtype,p.adjust.method = "BH")

  }

  # generate barplot

  FGA.col <- barcolor[1]

  FGG.col <- barcolor[2]

  FGL.col <- barcolor[3]

  p1 <- ggplot(summaryFGA, aes(x = Subtype, y = mean,fill=rep("0",nrow(summaryFGA)))) +

    geom_bar(stat = 'identity') +

    geom_errorbar(aes(ymax = mean+se, ymin = mean-se),position = position_dodge(0.9), width = 0.15) +

    annotate(geom="text",

             #hjust = ifelse(n.moic %% 2 == 0, 0.5, 0),

             x = n.moic / 2 + 0.5,

             #y = as.numeric(summaryFGA[which.max(summaryFGA$mean),"mean"] + summaryFGA[which.max(summaryFGA$mean),"se"]),

             y = as.numeric(summaryFGA[which.max(summaryFGA$mean),"mean"]),

             size = 8, angle = 90, fontface = "bold",

             #label = paste0(statistic, " p = ", formatC(FGA.test, digits = 1, format = "e")),

             label = cut(FGA.test,c(0,0.001,0.01,0.05,0.1,1),labels = c("****","***","**","*","."))) +

    scale_x_discrete(name = "",position = "top") +

    theme_bw() +

    theme(axis.line.y = element_line(size = 0.8),

          axis.ticks.y = element_line(size = 0.2),

          axis.text.y = element_blank(),

          axis.title.x = element_text(vjust = -0.3,size = 12),

          axis.text.x = element_text(size = 10, color = "black"),

          plot.margin = unit(c(0.3, -1.7, 0.3, 0.3), "lines"),

          legend.title = element_blank()) +

    coord_flip() +

    scale_fill_manual(values  = FGA.col, breaks = c("0"), labels = c("Copy number-altered genome")) +

    scale_y_reverse(expand = c(0.01,0),

                    name = "FGA (Fraction of Genome Altered)", position = "left")

  p2 <- ggplot(summaryFGGL, aes(x = Subtype, y = ifelse(class == 'FGG',mean,-mean),fill=class)) +

    geom_bar(stat = 'identity') +

    geom_errorbar(data=summaryFGGL[summaryFGGL$class=='FGG',],aes(ymax = mean+se, ymin =mean-se),position = position_dodge(0.9), width = 0.15) +

    geom_errorbar(data=summaryFGGL[summaryFGGL$class=='FGL',],aes(ymax = -mean-se, ymin =-mean+se),position = position_dodge(0.9), width = 0.15) +

    annotate(geom="text",

             x = n.moic / 2 + 0.5,

             y = -as.numeric(summaryFGL[which.max(summaryFGL$mean),"mean"]),

             size = 8, angle = 90, fontface = "bold",

             label = cut(FGL.test,c(0,0.001,0.01,0.05,0.1,1),labels = c("****","***","**","*","."))) +

    annotate(geom="text",

             x = n.moic / 2 + 0.5,

             y = as.numeric(summaryFGG[which.max(summaryFGG$mean),"mean"]),

             size = 8, angle = 90, fontface = "bold",

             label = cut(FGG.test,c(0,0.001,0.01,0.05,0.1,1),labels = c("****","***","**","*","."))) +

    scale_x_discrete(name = "") +

    theme_bw() +

    theme(axis.line.y = element_line(size = 0.8),

          axis.ticks.y = element_line(size = 0.2),

          axis.text.y = element_blank(),

          axis.title.x = element_text(vjust = -0.3, size = 12),

          axis.text.x = element_text(size = 10, color = "black"),

          plot.margin = unit(c(0.3, 0.3, 0.3, -1), "lines"),

          legend.title = element_blank()) +

    coord_flip() +

    scale_fill_manual(values  = c(FGL.col, FGG.col), breaks = c("FGL","FGG"),

                      labels = c("Copy number-lost genome","Copy number-gained genome")) +

    scale_y_continuous(expand = c(0.01,0),

                       name = "FGL or FGG (Fraction of Genome Lost or Gained)")

  pp <- ggplot() +

    # geom_text(data = summaryFGGL,

    #           aes(label = Subtype, x=Subtype), y = 0.5,

    #           size = 0.8*11/.pt, # match font size to theme

    #           hjust = 0.5, vjust = 0.5) +

    geom_label(data = summaryFGGL,

               aes(label = Subtype, x = Subtype, fill = Subtype),

               y = 0.5,

               color = "white",

               size = 0.9*11/.pt, # match font size to theme

               hjust = 0.4, vjust = 0.5) +

    scale_fill_manual(values = clust.col) +

    theme_minimal()+

    theme(axis.line.y =element_blank(),

          axis.ticks.y =element_blank(),

          axis.text.y =element_blank(),

          axis.title.y =element_blank(),

          axis.title.x =element_blank(),

          plot.margin = unit(c(0.3, 0, 0.3, 0), "lines")

    ) +

    guides(fill = FALSE) +

    coord_flip() +

    scale_y_reverse()

  pal <- p1 + pp + p2 +

    plot_layout(widths = c(7,1,7), guides = 'collect') & theme(legend.position = 'top')

  # save to pdf

  if(is.null(fig.name)) {

    outFig <- "barplot of FGA.pdf"

  } else {

    outFig <- paste0(fig.name,".pdf")

  }

  ggsave(file.path(fig.path, outFig), width = width, height = height)

  # print to screen

  print(pal)

  if(n.moic > 2) {

    return(list(summary = outTab,

                FGA.p.value = FGA.test,

                pairwise.FGA.test = pairwise.FGA.test,

                FGG.p.value = FGG.test,

                pairwise.FGG.test = pairwise.FGG.test,

                FGL.p.value = FGL.test,

                pairwise.FGL.test = pairwise.FGL.test,

                test.method = statistic))

  } else {

    return(list(summary = outTab,

                FGA.p.value = FGA.test,

                FGG.p.value = FGG.test,

                FGL.p.value = FGL.test,

                test.method = statistic))

  }

}