#' @name runMarker
#' @title Run identification of unique biomarkers
#' @description This function aims to identify uniquely and significantly expressed (overexpressed or downexpressed) biomarkers for each subtype identified by multi-omics integrative clustering algorithms. Top markers will be chosen to generate a template so as to run nearest template prediction for subtype verification.
#' @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 dea.method A string value to indicate the algorithm for differential expression analysis. Allowed value contains c('deseq2', 'edger', 'limma').
#' @param dat.path A string value to indicate the path for saving the files of differential expression analysis.
#' @param res.path A string value to indicate the path for saving the results for identifying subtype-specific markers.
#' @param p.cutoff A numeric value to indicate the nominal p value for identifying significant markers; pvalue < 0.05 by default.
#' @param p.adj.cutoff A numeric value to indicate the adjusted p value for identifying significant markers; padj < 0.05 by default.
#' @param dirct A string value to indicate the direction of identifying significant marker. Allowed values contain c('up', 'down'); `up` means up-regulated marker, and `down` means down-regulated marker.
#' @param n.marker A integer value to indicate how many top markers sorted by log2fc should be identified for each subtype; 200 by default.
#' @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 prefix A string value to indicate the prefix of differential expression file (use for searching files).
#' @param doplot A logic value to indicate if generating heatmap by using subtype-specific markers. TRUE by default.
#' @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 vector to assign marginal cutoff for truncating values in data; 3 by default.
#' @param centerFlag A logical vector to indicate if expression data should be centered; TRUE by default.
#' @param scaleFlag A logical vector to indicate if expression data should be scaled; TRUE 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 marker heatmap.
#' @param fig.name A string value to indicate the name of the marker 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 show_rownames A logic value to indicate if showing rownames (feature names) in heatmap; FALSE by default.
#' @param show_colnames A logic value to indicate if showing colnames (sample ID) in heatmap; FALSE by default.
#' @param ... Additional parameters pass to \link[ComplexHeatmap]{pheatmap}.
#' @return A figure of subtype-specific marker heatmap (.pdf) if \code{doPlot = TRUE} and a list with the following components:
#'
#'         \code{unqlist}   a string vector storing the unique marker across all subtypes.
#'
#'         \code{templates} a data.frame storing the the template information for nearest template prediction, which is used for verification in external cohort.
#'
#'         \code{dirct}     a string value indicating the direction for identifying subtype-specific markers.
#'
#'         \code{heatmap}   a complexheatmap object.
#' @export
#' @importFrom grDevices pdf dev.off colorRampPalette
#' @importFrom ComplexHeatmap pheatmap draw
#' @references Gu Z, Eils R, Schlesner M (2016). Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics, 32(18):2847-2849.
#' @examples # There is no example and please refer to vignette.
runMarker <- function(moic.res      = NULL,
                      dea.method    = c("deseq2", "edger", "limma"),
                      prefix        = NULL,
                      dat.path      = getwd(),
                      res.path      = getwd(),
                      p.cutoff      = 0.05,
                      p.adj.cutoff  = 0.05,
                      dirct         = "up",
                      n.marker      = 200,
                      doplot        = TRUE,
                      norm.expr     = NULL,
                      annCol        = NULL,
                      annColors     = NULL,
                      clust.col     = c("#2EC4B6","#E71D36","#FF9F1C","#BDD5EA","#FFA5AB","#011627","#023E8A","#9D4EDD"),
                      halfwidth     = 3,
                      centerFlag    = TRUE,
                      scaleFlag     = TRUE,
                      show_rownames = FALSE,
                      show_colnames = FALSE,
                      color         = c("#5bc0eb", "black", "#ECE700"),
                      fig.path      = getwd(),
                      fig.name      = NULL,
                      width         = 8,
                      height        = 8,
                      ...) {

  n.moic <- length(unique(moic.res$clust.res$clust))
  mo.method <- moic.res$mo.method
  DEpattern <- paste(mo.method, "_", ifelse(is.null(prefix),"",paste0(prefix,"_")), dea.method, ".*._vs_Others.txt$", sep="")
  DEfiles <- dir(dat.path,pattern = DEpattern)

  if(length(DEfiles) == 0) {stop("no DEfiles!")}
  if(length(DEfiles) != n.moic) {stop("not all the multi-omics clusters have DEfile!")}
  if (!is.element(dirct, c("up", "down"))) {stop( "dirct type error! Allowed value contains c('up', 'down').") }

  if(dirct == "up") { outlabel <- "unique_upexpr_marker.txt" }
  if(dirct == "down") { outlabel <- "unique_downexpr_marker.txt" }

  genelist <- c()
  for (filek in DEfiles) {
    DEres <- read.table(file.path(dat.path, filek), header=TRUE, row.names=NULL, sep="\t", quote="", stringsAsFactors=FALSE)
    DEres <- DEres[!duplicated(DEres[, 1]),]
    DEres <- DEres[!is.na(DEres[, 1]), ]
    rownames(DEres) <- DEres[, 1]
    DEres <- DEres[, -1]

    #rownames(DEres) <- toupper(rownames(DEres))
    if (dirct == "up") {
      genelist <- c( genelist, rownames(DEres[!is.na(DEres$padj) & DEres$pvalue < p.cutoff & DEres$padj < p.adj.cutoff & !is.na(DEres$log2fc) & DEres$log2fc > 0, ]) )
    }
    if (dirct == "down") {
      genelist <- c( genelist, rownames(DEres[!is.na(DEres$padj) & DEres$pvalue < p.cutoff & DEres$padj < p.adj.cutoff & !is.na(DEres$log2fc) & DEres$log2fc < 0, ]) )
    }
  }
  unqlist <- setdiff(genelist,genelist[duplicated(genelist)])

  marker <- list()
  for (filek in DEfiles) {
    DEres <- read.table(file.path(dat.path, filek), header=TRUE, row.names=NULL, sep="\t", quote="", stringsAsFactors=FALSE)
    DEres <- DEres[!duplicated(DEres[, 1]),]
    DEres <- DEres[!is.na(DEres[, 1]), ]
    rownames(DEres) <- DEres[, 1]
    DEres <- DEres[, -1]

    #rownames(DEres) <- toupper(rownames(DEres))
    if(dirct == "up") {
      outk <- intersect( unqlist, rownames(DEres[!is.na(DEres$padj) & DEres$pvalue < p.cutoff & DEres$padj < p.adj.cutoff & !is.na(DEres$log2fc) & DEres$log2fc > 0, ]) )
      outk <- DEres[outk,]
      outk <- outk[order(outk$log2fc, decreasing = TRUE),]

      if(nrow(outk) > n.marker) {
        marker[[filek]] <- outk[1:n.marker,]
      } else {
        marker[[filek]] <- outk
      }
      marker$dirct <- "up"
    }
    if(dirct == "down") {
      outk <- intersect( unqlist, rownames(DEres[!is.na(DEres$padj) & DEres$pvalue < p.cutoff & DEres$padj < p.adj.cutoff & !is.na(DEres$log2fc) & DEres$log2fc < 0, ]) )
      outk <- DEres[outk,]
      outk <- outk[order(outk$log2fc, decreasing = FALSE),]

      if(nrow(outk) > n.marker) {
        marker[[filek]] <- outk[1:n.marker,]
      } else {
        marker[[filek]] <- outk
      }
      marker$dirct <- "down"
    }
    # write file
    write.table(outk, file=file.path(res.path, paste(gsub("_vs_Others.txt","", filek, fixed = TRUE), outlabel, sep = "_")), row.names = TRUE, col.names = NA, sep = "\t", quote = FALSE)
  }

  # generate templates for nearest template prediction
  templates <- NULL
  for (filek in DEfiles) {
    tmp <- data.frame(probe = rownames(marker[[filek]]),
                      class = sub("_vs_Others.txt","",sub(".*.result.","",filek)),
                      dirct = marker$dirct,
                      stringsAsFactors = FALSE)
    templates <- rbind.data.frame(templates, tmp, stringsAsFactors = FALSE)
  }
  write.table(templates, file=file.path(res.path, paste0(mo.method,"_",dea.method,"_",dirct,"regulated_marker_templates.txt")), row.names = FALSE, sep = "\t", quote = FALSE)

  # generate heatmap with subtype-specific markers
  if(doplot) {
    if(is.null(norm.expr)) {
      stop("please provide a matrix or data.frame of normalized expression data with rows for genes and columns for samples; FPKM or TPM without log2 transformation is recommended.")
    }

    # 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]

    if(is.null(fig.name)) {
      outFig <- paste0("markerheatmap_using_",dirct,"regulated_genes.pdf")
    } else {
      outFig <- paste0(fig.name, "_using_", dirct,"regulated_genes.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(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.")
      gset <- norm.expr
    }
    if(max(norm.expr) >= 25 & min(norm.expr) >= 0){
      message("--log2 transformation done for expression data.")
      gset <- log2(norm.expr + 1)
    }

    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)
    }

    plotdata <- standarize.fun(gset[intersect(templates$probe, rownames(gset)), sam.order],
                               halfwidth = halfwidth,
                               centerFlag = centerFlag,
                               scaleFlag = scaleFlag)

    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
        }
      }
    }

    hm <- ComplexHeatmap::pheatmap(mat               = plotdata,
                                   border_color      = NA,
                                   cluster_cols      = FALSE,
                                   cluster_rows      = FALSE,
                                   annotation_col    = annCol,
                                   annotation_colors = annColors,
                                   legend_breaks     = pretty(c(-halfwidth,halfwidth)),
                                   legend_labels     = pretty(c(-halfwidth,halfwidth)),
                                   show_rownames     = show_rownames,
                                   show_colnames     = show_colnames,
                                   treeheight_col    = 0,
                                   treeheight_row    = 0,
                                   color             = grDevices::colorRampPalette(color)(64),
                                   ...)

    # save to pdf
    pdf(file.path(fig.path, outFig), width = width, height = height)
    draw(hm,annotation_legend_side = "left",heatmap_legend_side = "left")
    invisible(dev.off())

    # print to screen
    draw(hm,annotation_legend_side = "left",heatmap_legend_side = "left")
    return(list(unqlist = unqlist, templates = templates, dirct = dirct, heatmap = hm))
  } else {
    return(list(unqlist = unqlist, templates = templates, dirct = dirct))
  }
}