--- a
+++ b/R/compFGA.R
@@ -0,0 +1,284 @@
+#' @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))
+  }
+}