#' 'compute N-of-1-kMEn'
#'
#' @param abs.logFC #abslogFC is calculated as |log2(X1+1)-log2(X2+1)|, where X1 is the expression value of a gene under condition 1, and X2 is the expression values of a gene under condition 2
#' @param gene.symbol gene symbol is a vector of gene names corresponding to abs.logFC
#' @param GeneSet.list GeneSet.list is a list of genesets. Each element of this geneset contains the gene names of the genes belonging to the gene set.
#' @param FDR.method method for multiplicity correction
#' @param alternative indicates the alternative hypothesis in Fisher's exact test.
#' 
#' @return A list contains two elements. The first element is the pathways and their corresponding p.values. The second element is all the genes and their corresponding DEG status
#' If either input vectors have no variance, the MD distance is chosen to the simple difference.


kMEn = function(abs.logFC, gene.symbol, GeneSet.list, FDR.method = 'BY', alternative = c('two.sided','greater','less')){
  set.seed(369)
  alternative = match.arg(alternative)
  res.kMeans = kmeans(abs.logFC, centers = 2, iter.max = 1000, nstart = 500)
  ind.DEG = res.kMeans$cluster==which.max(res.kMeans$centers) 
  DE.df = data.frame(symbol = gene.symbol, DEstatus = as.integer(ind.DEG))
  p_val = sapply(GeneSet.list, function(x) fisher.test(makeContingencyTable(x,DE.df),alternative = alternative)$p.value)
  or = sapply(GeneSet.list, function(x) cal.OR(makeContingencyTable(x,DE.df)))
  res.pathDE = data.frame(nom_pval = p_val, pval_adj = p.adjust(p_val,method = FDR.method), OR = or, check.rows = T)
  res.ls = list()
  res.ls[['Pathway']] = res.pathDE[order(res.pathDE$nom_pval),]
  DEG.status = as.matrix(ind.DEG)
  colnames(DEG.status) = 'DEG_status'
  res.ls[['DEG']] = DEG.status
  return(res.ls)
}

# define a function to calculate odds ratio when the contingeency table is given
cal.OR = function(cont.table){
  a = cont.table[1,1]
  b = cont.table[1,2]
  c = cont.table[2,1]
  d = cont.table[2,2]
  OR = (a*d)/(b*c)
  return(OR)
}

# define a function to create a contingency table for the enrichment analysis
makeContingencyTable = function(GeneSet, DE.df){ # GeneSet lists the gene symbols of this geneset; DE.df is the data frame which shows the DE status of every gene in the genome(first column is symbol; second column is DEstatus)
  ind.geneSet = DE.df$symbol %in% GeneSet
  DEGinSet = sum(DE.df$DEstatus[ind.geneSet])
  #NonDEGinSet = length(GeneSet) - DEGinSet
  NonDEGinSet = sum(ind.geneSet) - DEGinSet
  DEGoutSet = sum(DE.df$DEstatus) - DEGinSet
  NonDEGoutSet = dim(DE.df)[1] - sum(DE.df$DEstatus) - NonDEGinSet
  if(DEGinSet + NonDEGinSet + DEGoutSet + NonDEGoutSet != dim(DE.df)[1]){
    warning('check contigency table')
  }
  contingency.table = matrix(c(DEGinSet,DEGoutSet,NonDEGinSet,NonDEGoutSet),ncol=2)
  return(contingency.table)
}


# define a function to calculate absolute log Fold change
f.absLogFC = function(x,y){  # x and y are two vectors expression values
  return(abs(log2(x+1) - log2(y+1)))
}