#' @name runGSVA

#' @title Run gene set variation analysis

#' @description Use gene set variation analysis to calculate enrichment score of each sample in each subtype based on given gene set list of interest.

#' @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 norm.expr A matrix of normalized expression data with rows for genes and columns for samples; FPKM or TPM without log2 transformation is recommended.

#' @param gset.gmt.path A string value to indicate ABSOULUTE PATH/NAME of gene sets of interest stored as GMT format \url{https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats#GMT:_Gene_Matrix_Transposed_file_format_.28.2A.gmt.29}.

#' @param annCol A data.frame storing annotation information for samples.

#' @param annColors A list of string vectors for colors matched with annCol.

#' @param clust.col A string vector storing colors for annotating each subtype at the top of heatmap.

#' @param halfwidth A numeric value to assign marginal cutoff for truncating enrichment scores; 1 by default.

#' @param centerFlag A logical vector to indicate if enrichment scores should be centered; TRUE by default.

#' @param scaleFlag A logical vector to indicate if enrichment scores should be scaled; TRUE by default.

#' @param distance A string value of distance measurement for hierarchical clustering; 'euclidean' by default.

#' @param linkage A string value of clustering method for hierarchical clustering; 'ward.D' by default.

#' @param show_rownames A logic value to indicate if showing rownames (feature names) in heatmap; TRUE by default.

#' @param show_colnames A logic value to indicate if showing colnames (sample ID) in heatmap; FALSE by default.

#' @param color A string vector storing colors for heatmap.

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

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

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

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

#' @param ... Additional parameters pass to \link[ComplexHeatmap]{pheatmap}.

#'

#' @return A figure of enrichment heatmap (.pdf) and a list with the following components:

#'

#'         \code{gset.list}  a list storing gene sets information converted from GMT format by \link[clusterProfiler]{read.gmt}.

#'

#'         \code{raw.es}     a data.frame storing raw enrichment score based on given gene sets of interest by using specified \code{gsva.method}.

#'

#'         \code{scaled.es}  a data.frame storing z-scored enrichment score based on given gene sets of interest by using specified \code{gsva.method}.

#'

#' @export

#' @importFrom ClassDiscovery distanceMatrix

#' @importFrom clusterProfiler read.gmt

#' @importFrom GSVA gsva

#' @importFrom ComplexHeatmap pheatmap draw ht_opt

#' @importFrom grDevices pdf dev.off colorRampPalette

#' @references Barbie, D.A. et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature, 462(5):108-112.

#'

#' Hänzelmann, S., Castelo, R. and Guinney, J. (2013). GSVA: Gene set variation analysis for microarray and RNA-Seq data. BMC Bioinformatics, 14(1):7.

#'

#' Lee, E. et al. (2008). Inferring pathway activity toward precise disease classification. PLoS Comp Biol, 4(11):e1000217.

#'

#' Tomfohr, J. et al. (2005). Pathway level analysis of gene expression using singular value decomposition. BMC Bioinformatics, 6(1):1-11.

#'

#' Yu G, Wang L, Han Y, He Q (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS, 16(5):284-287.

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

runGSVA <- function(moic.res      = NULL,

                    norm.expr     = NULL,

                    gset.gmt.path = NULL,

                    gsva.method   = "gsva",

                    centerFlag    = TRUE,

                    scaleFlag     = TRUE,

                    halfwidth     = 1,

                    annCol        = NULL,

                    annColors     = NULL,

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

                    distance      = "euclidean",

                    linkage       = "ward.D",

                    show_rownames = TRUE,

                    show_colnames = FALSE,

                    color         = c("#366A9B", "#4E98DE", "#DDDDDD", "#FBCFA7", "#F79C4A"),

                    fig.path      = getwd(),

                    fig.name      = NULL,

                    width         = 8,

                    height        = 8,

                    ...) {

  # standardize function

  standarize.fun <- function(indata=NULL, halfwidth=NULL, centerFlag=TRUE, scaleFlag=TRUE) {

    outdata=t(scale(t(indata), center=centerFlag, scale=scaleFlag))

    if (!is.null(halfwidth)) {

      outdata[outdata>halfwidth]=halfwidth

      outdata[outdata<(-halfwidth)]= -halfwidth

    }

    return(outdata)

  }

  # check data

  comsam <- intersect(moic.res$clust.res$samID, colnames(norm.expr))

  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."))

  }

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

  norm.expr <- norm.expr[,comsam]

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

  # load gene set data and convert gmt to data.frame

  gset <- try(clusterProfiler::read.gmt(gset.gmt.path), silent = TRUE)

  if(class(gset) == "try-error") {stop("please provide correct ABSOLUTE PATH for gene sets of interest.")}

  # convert data.frame to list

  term <- unique(gset[,1])

  gset.list <- list()

  for (i in term) {

    gset.list[[i]] <- gset[which(gset[,1] == i),2]

  }

  # calculate gene set enrichment scores

  if(max(norm.expr) < 25 | (max(norm.expr) >= 25 & min(norm.expr) < 0)) {

    message("--expression profile seems to have been standardised (z-score or log transformation), no more action will be performed.")

  }

  if(max(norm.expr) >= 25 & min(norm.expr) >= 0){

    message("--log2 transformation done for expression data.")

    norm.expr <- log2(norm.expr + 1)

  }

  es <- GSVA::gsva(expr          = as.matrix(norm.expr),

                   gset.idx.list = gset.list,

                   method        = gsva.method,

                   parallel.sz   = 1)

  es.backup <- es

  es <- standarize.fun(es, halfwidth = halfwidth, centerFlag = centerFlag, scaleFlag = scaleFlag)

  message(gsva.method," done...")

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

    outFig <- paste0("enrichment_heatmap_using_", gsva.method, ".pdf")

  } else {

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

  }

  sam.order <- moic.res$clust.res[order(moic.res$clust.res$clust, decreasing = FALSE), "samID"]

  colvec <- clust.col[1:n.moic]

  names(colvec) <- paste0("CS",1:n.moic)

  if(!is.null(annCol) & !is.null(annColors)) {

    annCol <- annCol[sam.order, , drop = FALSE]

    annCol$Subtype <- paste0("CS",moic.res$clust.res[sam.order,"clust"])

    annColors[["Subtype"]] <- colvec

  } else {

    annCol <- data.frame("Subtype" = paste0("CS",moic.res$clust.res[sam.order,"clust"]),

                         row.names = sam.order,

                         stringsAsFactors = FALSE)

    annColors <- list("Subtype" = colvec)

  }

  if(!is.null(annCol) & !is.null(annColors)) {

    for (i in names(annColors)) {

      if(is.function(annColors[[i]])) {

        annColors[[i]] <- annColors[[i]](pretty(range(annCol[,i]),n = 64)) # transformat colorRamp2 function to color vector

      }

    }

  }

  ht_opt$message = FALSE

  if(is.null(distance) | is.null(linkage)) {

    hcg <- FALSE

  } else {

    hcg <- hclust(ClassDiscovery::distanceMatrix(t(as.matrix(es[,sam.order])), distance), linkage)

  }

  hm <- ComplexHeatmap::pheatmap(mat               = es[,sam.order],

                                 border_color      = NA,

                                 cluster_cols      = FALSE,

                                 cluster_rows      = hcg,

                                 annotation_col    = annCol,

                                 annotation_colors = annColors,

                                 show_rownames     = show_rownames,

                                 show_colnames     = show_colnames,

                                 color             = grDevices::colorRampPalette(color)(64),

                                 ...)

  # save to pdf

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

  draw(hm)

  invisible(dev.off())

  # print to screen

  draw(hm)

  return(list(gset.list = gset.list, raw.es = es.backup, scaled.es = es))

}