R/MAPE_internal.R

In metaOmic/MetaPath: Perform the Meta-Analysis for Pathway Enrichment

## Internal functions of MAPE
##
## @title Internal functions
## @author Kui Shen, Xiangrui Zeng and George Tseng.

###########################
cor.func <- function (x, y)
{
  n <- length(y)
  xbar <- x %*% rep(1/n, n)
  sxx <- ((x - as.vector(xbar))^2) %*% rep(1, n)
  sxy <- (x - as.vector(xbar)) %*% (y - mean(y))
  syy <- sum((y - mean(y))^2)
  numer <- sxy/sxx
  sd <- sqrt((syy/sxx - numer^2)/(n - 2))
  tt <- numer/sd
  return(list(tt = tt, numer = numer, sd = sd))
}


###########################
coxfunc <- function (x, y, censoring.status)
{
  junk <- coxscor(x, y, censoring.status)
  scor <- junk$scor
  sd <- sqrt(coxvar(x, y, censoring.status, coxstuff.obj = junk$coxstuff.obj))
  tt <- scor/sd
  return(list(tt = tt, numer = scor, sd = sd))
}


###########################
cox.perm.sample <- function (expr, testgroup, censoring, nperm)
{
  obs = coxfunc(expr, testgroup, censoring)$tt
  perms = matrix(NA, nrow(expr), nperm)
  perms.mtx = matrix(NA, nrow = nperm, ncol = length(testgroup))
  for (i in 1:nperm) {
    perms.mtx[i, ] = sample(1:length(testgroup), size = length(testgroup))
  }
  for (t1 in 1:nperm) {
    perms[, t1] = coxfunc(expr, testgroup[perms.mtx[t1, ]],
                          censoring)$tt
  }
  rownames(perms) = rownames(expr)
  colnames(perms) = paste("B", 1:nperm, sep = "")
  obs = abs(obs)
  perms = abs(perms)
  return(list(obs = obs, perms = perms, perms.mtx = perms.mtx))
}


###########################
coxscor <- function (x, y, ic, offset = rep(0, length(y)))
{
  n <- length(y)
  nx <- nrow(x)
  yy <- y + (ic == 0) * (1e-05)
  otag <- order(yy)
  y <- y[otag]
  ic <- ic[otag]
  x <- x[, otag, drop = F]
  offset <- offset[otag]
  a <- coxstuff(x, y, ic, offset = offset)
  nf <- a$nf
  fail.times <- a$fail.times
  s <- a$s
  d <- a$d
  dd <- a$dd
  nn <- a$nn
  nno <- a$nno
  w <- rep(0, nx)
  for (i in (1:nf)) {
    w <- w + s[, i]
    oo <- (1:n)[y >= fail.times[i]]
    r <- rowSums(x[, oo, drop = F] * exp(offset[oo]))
    w <- w - (d[i]/nno[i]) * r
  }
  return(list(scor = w, coxstuff.obj = a))
}


#########################
coxstuff <- function (x, y, ic, offset = rep(0, length(y)))
{
  fail.times <- unique(y[ic == 1])
  nf <- length(fail.times)
  n <- length(y)
  nn <- rep(0, nf)
  nno <- rep(0, nf)
  for (i in 1:nf) {
    nn[i] <- sum(y >= fail.times[i])
    nno[i] <- sum(exp(offset)[y >= fail.times[i]])
  }
  s <- matrix(0, ncol = nf, nrow = nrow(x))
  d <- rep(0, nf)
  for (i in 1:nf) {
    o <- (1:n)[(y == fail.times[i]) & (ic == 1)]
    d[i] <- length(o)
  }
  oo <- match(y, fail.times)
  oo[ic == 0] <- NA
  oo[is.na(oo)] <- max(oo[!is.na(oo)]) + 1
  s <- t(rowsum(t(x), oo))
  if (ncol(s) > nf) {
    s <- s[, -ncol(s)]
  }
  dd <- rep(0, n)
  for (j in 1:nf) {
    dd[(y == fail.times[j]) & (ic == 1)] <- d[j]
  }
  return(list(fail.times = fail.times, s = s, d = d, dd = dd,
              nf = nf, nn = nn, nno = nno))
}


##########################
coxvar <- function (x, y, ic, offset = rep(0, length(y)), coxstuff.obj = NULL)
{
  nx <- nrow(x)
  n <- length(y)
  yy <- y + (ic == 0) * (1e-06)
  otag <- order(yy)
  y <- y[otag]
  ic <- ic[otag]
  x <- x[, otag, drop = F]
  offset <- offset[otag]
  if (is.null(coxstuff.obj)) {
    coxstuff.obj <- coxstuff(x, y, ic, offset = offset)
  }
  nf <- coxstuff.obj$nf
  fail.times <- coxstuff.obj$fail.times
  s <- coxstuff.obj$s
  d <- coxstuff.obj$d
  dd <- coxstuff.obj$dd
  nn <- coxstuff.obj$nn
  nno <- coxstuff.obj$nno
  x2 <- x^2
  oo <- (1:n)[y >= fail.times[1]]
  sx <- (1/nno[1]) * rowSums(x[, oo] * exp(offset[oo]))
  s <- (1/nno[1]) * rowSums(x2[, oo] * exp(offset[oo]))
  w <- d[1] * (s - sx * sx)
  for (i in 2:nf) {
    oo <- (1:n)[y >= fail.times[i - 1] & y < fail.times[i]]
    sx <- (1/nno[i]) * (nno[i - 1] * sx - rowSums(x[, oo,
                                                    drop = F] * exp(offset[oo])))
    s <- (1/nno[i]) * (nno[i - 1] * s - rowSums(x2[, oo,
                                                   drop = F] * exp(offset[oo])))
    w <- w + d[i] * (s - sx * sx)
  }
  return(w)
}


#######################
reg.perm.sample <- function (expr, testgroup, nperm)
{
  obs = cor.func(expr, testgroup)$tt[, 1]
  perms = matrix(NA, nrow(expr), nperm)
  perms.mtx = matrix(NA, nrow = nperm, ncol = length(testgroup))
  for (i in 1:nperm) {
    perms.mtx[i, ] = sample(1:length(testgroup), size = length(testgroup))
  }
  for (t1 in 1:nperm) {
    perms[, t1] = cor.func(expr, testgroup[perms.mtx[t1,
                                                     ]])$tt[, 1]
  }
  rownames(perms) = rownames(expr)
  colnames(perms) = paste("B", 1:nperm, sep = "")
  obs = abs(obs)
  perms = abs(perms)
  return(list(obs = obs, perms = perms, perms.mtx = perms.mtx))
}


######################
Tperm.sample <- function (x, fac, nperm)
{
  obs = genefilter::rowttests(x, fac, tstatOnly = T)$statistic
  names(obs) = rownames(x)
  perms = matrix(NA, nrow(x), nperm)
  for (t1 in 1:nperm) {
    perms[, t1] = genefilter::rowttests(x, sample(fac), tstatOnly = T)$statistic
  }
  rownames(perms) = rownames(x)
  colnames(perms) = paste("B", 1:nperm, sep = "")
  obs = abs(obs)
  perms = abs(perms)
  return(list(obs = obs, perms = perms))
}


#######################
F.perm.sample <- function(x, fac, nperm)
{
  obs = genefilter::rowFtests(x, fac, var.equal = TRUE)$statistic
  names(obs)=rownames(x)
  perms= matrix(NA,  nrow(x),  nperm)

  for(t1 in 1:nperm){
    perms[,t1]=genefilter::rowFtests(x, sample(fac), var.equal = TRUE)$statistic
  }
  rownames(perms)=rownames(x)
  colnames(perms)=paste('B',1:nperm,sep="")

  obs=abs(obs); perms=abs(perms)

  return(list(obs=obs, perms=perms))
}


######################
pqvalues.compute <- function (Stat.0, Stat.B, Stat.type, PI0 = 1)
{
  Stat.0 = as.matrix(abs(Stat.0))
  Stat.B = as.matrix(abs(Stat.B))
  if (nrow(Stat.0) != nrow(Stat.B))
    stop("# of rows of Stat.0 and Stat.B should be same")
  B = ncol(Stat.B)
  G = nrow(Stat.B)
  if (Stat.type == "Tstat") {
    Stat.all = cbind(Stat.0, Stat.B)
    count.0 = apply(Stat.0, 1, function(x) sum(x <= Stat.0,
                                               na.rm = T))
    Stat.rank.all = rank(-as.vector(Stat.all), ties.method = "max")
    pvalue.0 = as.matrix(Stat.rank.all[1:G] - count.0)/(B *
                                                          G)
    if (is.null(PI0)) {
      PI0 = sum(pvalue.0 >= 0.5)/(0.5 * G)
    }
    else {
      PI0 = 1
    }
    qvalue.0 = PI0 * pvalue.0 * G/(apply(Stat.0, 1, function(x) sum(x <=
                                                                      Stat.0, na.rm = T)))
    qvalue.0 = ifelse(qvalue.0 <= 1, qvalue.0, 1)
    Stat.rank = rank(-as.vector(Stat.B), ties.method = "max")
    pvalue.B = matrix(Stat.rank/(B * G), G, B)
    rownames(pvalue.B) = rownames(Stat.B)
    colnames(pvalue.B) = colnames(Stat.B)
  }
  else if (Stat.type == "Pvalue") {
    Stat.all = cbind(Stat.0, Stat.B)
    count.0 = apply(Stat.0, 1, function(x) sum(x >= Stat.0))
    Stat.rank.all = rank(as.vector(Stat.all), ties.method = "max")
    pvalue.0 = as.matrix(Stat.rank.all[1:G] - count.0)/(B *
                                                          G)
    if (is.null(PI0)) {
      PI0 = sum(pvalue.0 >= 0.5)/(0.5 * G)
    }
    else {
      PI0 = 1
    }
    qvalue.0 = PI0 * pvalue.0 * G/(apply(Stat.0, 1, function(x) sum(x >=
                                                                      Stat.0, na.rm = T)))
    qvalue.0 = ifelse(qvalue.0 <= 1, qvalue.0, 1)
    Stat.rank = rank(as.vector(Stat.B), ties.method = "max")
    pvalue.B = matrix(Stat.rank/(B * G), G, B)
    rownames(pvalue.B) = rownames(Stat.B)
    colnames(pvalue.B) = colnames(Stat.B)
  }
  else {
    stop("wrong stat.type")
  }
  return(list(pvalue.0 = pvalue.0, pvalue.B = pvalue.B, qvalue.0 = qvalue.0,
              PI0 = PI0))
}


##########################
EnrichmentC_Exact<- function (madata, label, censoring = NULL, DB.matrix, DEgene.number = DEgene.number,
                              size.min = 15,size.max = 500, nperm = 500, gene.pvalues = NULL,
                              resp.type = NULL)
{
  if (!is.null(madata)) {
    genes.in.study = featureNames(madata)
    set2allgenes.mtx = DB.matrix
    gene.common = intersect(featureNames(madata), colnames(set2allgenes.mtx))
    if (nrow(set2allgenes.mtx) > 1) {
      set2allgenes.mtx = as.matrix(set2allgenes.mtx[, gene.common])
    }
    else {
      set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[,
                                                      gene.common]))
    }
    madata = madata[gene.common, ]
    if (resp.type == "twoclass") {
      tstat = genefilter::rowttests(exprs(madata), as.factor(label),
                                    tstatOnly = F)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "multiclass") {
      tstat = genefilter::rowFtests(exprs(madata), as.factor(label),
                                    var.equal = TRUE)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "survival") {
      if (is.null(censoring)) {
        stop("Error: censoring status is null")
      }
      cox.out <- coxfunc(exprs(madata), label, censoring)
      scores <- abs(cox.out$tt)
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else if (resp.type == "continuous") {
      cor.out <- cor.func(exprs(madata), label)
      scores <- abs(cor.out$tt[, 1])
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
  }
  DEgene = gene.name.sort[1:DEgene.number]
  pvalue.set.0 = matrix(NA,nrow = nrow(DB.matrix),ncol = 1)
  rownames(pvalue.set.0) = rownames(DB.matrix)
  for (i in 1:nrow(DB.matrix)){
    count_table<-matrix(0,2,2)
    p_value <-NA
    ####in the gene list and in the pathway
    count_table[1,1]<-sum(DEgene %in% colnames(DB.matrix[,DB.matrix[i,]==1]))
    ####in the gene list but not in the pathway
    count_table[1,2]<-length(DEgene)-count_table[1,1]
    ####not in the gene list but in the pathway
    count_table[2,1]<-sum(genes.in.study%in% colnames(DB.matrix[,DB.matrix[i,]==1]))
    ####not in the gene list and not in the pathway
    count_table[2,2]<-length(genes.in.study)-count_table[2,1]
    if(length(count_table)==4){
      pvalue.set.0[i,1] <- fisher.test(count_table, alternative="greater")$p}
  }
  pvalue.set.0 = pvalue.set.0[pvalue.set.0[,1]!=1,,drop=FALSE]
  qvalue.set.0 = matrix(p.adjust(pvalue.set.0[,1], "BH"),ncol = 1)
  return(list(pvalue.set.0 = pvalue.set.0,qvalue.set.0 = qvalue.set.0,DEgene = DEgene))
}


########################
CPI_Exact <- function (study, label, censoring.status, DB.matrix, size.min = 15,
                       DEgene.number = DEgene.number, size.max = 500, nperm = 500,
                       stat = NULL, rth.value = NULL,resp.type)
{
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  out = list()
  for (t1 in 1:length(study)) {
    madata = study[[t1]]
    testlabel = madata[[label]]
    out[[t1]] = list()
    if (resp.type == "survival") {
      censoring = madata[[censoring.status]]
    }
    out[[t1]] = EnrichmentC_Exact(madata = madata, label = testlabel,DEgene.number = DEgene.number,
                                  censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                                  size.max = size.max, nperm = nperm, resp.type = resp.type)
  }
  set.common = rownames(out[[1]]$pvalue.set.0)
  for (t1 in 2:length(study)) {
    set.common = intersect(set.common, rownames(out[[t1]]$pvalue.set.0))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  pvalue.0.mtx = matrix(NA, length(set.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[set.common,
                                                ]
  }
  genelist = list()
  for (t1 in 1:length(study)) {
    genelist[[t1]] = out[[t1]]$DEgene
  }
  rownames(pvalue.0.mtx) = set.common
  rm(out)

  AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
  p_value_meta = AW.fisher$pvalues
  weights_meta = AW.fisher$weights
  q_value_meta = p.adjust(p_value_meta, "BH")
  summary<-data.frame(q_value_meta = q_value_meta,p_value_meta = p_value_meta,
                      p_data = pvalue.0.mtx)
  return(list(summary = summary,genelist = genelist,AW.weights = weights_meta))
}


####################
CPI_gene_KS <- function (study, label, censoring.status, DB.matrix, size.min = 15,
                         size.max = 500, nperm = 500, stat = NULL, rth.value = NULL,
                         resp.type)
{if (is.null(names(study)))
  names(study) = paste("study.", 1:length(study), sep = "")
out = list()
for (t1 in 1:length(study)) {
  madata = study[[t1]]
  testlabel = madata[[label]]
  out[[t1]] = list()
  if (resp.type == "survival") {
    censoring = madata[[censoring.status]]
  }
  out[[t1]] = Enrichment_KS_gene(madata = madata, label = testlabel,
                                 censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                                 size.max = size.max, nperm = nperm, resp.type = resp.type)
}
set.common = rownames(out[[1]]$pvalue.set.0)
for (t1 in 2:length(study)) {
  set.common = intersect(set.common, rownames(out[[t1]]$pvalue.set.0))
}
if (is.null(names(study)))
  names(study) = paste("study.", 1:length(study), sep = "")
pvalue.B.array = array(data = NA, dim = c(length(set.common),
                                          nperm, length(study)))
dimnames(pvalue.B.array) = list(set.common, paste("perm",
                                                  1:nperm, sep = ""), names(study))
pvalue.0.mtx = matrix(NA, length(set.common), length(study))
qvalue.0.mtx = matrix(NA, length(set.common), length(study))
for (t1 in 1:length(study)) {
  pvalue.B.array[, , t1] = out[[t1]]$pvalue.set.B[set.common,]
  pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[set.common,]
  qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[set.common,]
}
rownames(qvalue.0.mtx) = set.common
rownames(pvalue.0.mtx) = set.common
rm(out)

AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
p_value_meta = AW.fisher$pvalues
weights_meta = AW.fisher$weights
q_value_meta = p.adjust(p_value_meta, "BH")
summary<-data.frame(q_value_meta = q_value_meta,p_value_meta = p_value_meta,
                    p_data = pvalue.0.mtx)
return(list(summary = summary,AW.weights = weights_meta))
}


#########################
CPI_sample_KS <- function (study, label, censoring.status = NULL, DB.matrix, size.min = 15,
                           size.max = 500, nperm = 500, stat, rth.value = NULL, resp.type)
{
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  out = list()
  for (t1 in 1:length(study)) {
    madata = study[[t1]]
    testlabel = madata[[label]]
    out[[t1]] = list()
    if (resp.type == "survival") {
      censoring = madata[[censoring.status]]
    }
    out[[t1]] = Enrichment_KS_sample(madata = madata, label = testlabel,
                                     censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                                     size.max = size.max, nperm = nperm, resp.type = resp.type)
  }
  common.pathway = rownames(out[[1]]$pvalue.set.0)
  for (t1 in 1:length(study)) {
    common.pathway = intersect(common.pathway, rownames(out[[t1]]$pvalue.set.0))
  }
  pvalue.B.array = array(data = NA, dim = c(length(common.pathway),
                                            nperm, length(study)))
  pvalue.0.mtx = matrix(NA, length(common.pathway), length(study))
  qvalue.0.mtx = matrix(NA, length(common.pathway), length(study))
  rownames(pvalue.0.mtx) = common.pathway
  colnames(pvalue.0.mtx) = names(study)
  rownames(qvalue.0.mtx) = common.pathway
  colnames(qvalue.0.mtx) = names(study)
  dimnames(pvalue.B.array) = list(common.pathway, paste("perm",
                                                        1:nperm, sep = ""), names(study))
  for (t1 in 1:length(study)) {
    pvalue.B.array[, , t1] = out[[t1]]$pvalue.set.B[common.pathway,
                                                    ]
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[common.pathway,
                                                ]
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[common.pathway,
                                                ]
  }

  AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
  p_value_meta = AW.fisher$pvalues
  weights_meta = AW.fisher$weights
  q_value_meta = p.adjust(p_value_meta, "BH")

  summary<-data.frame(q_value_meta = q_value_meta,p_value_meta = p_value_meta,
                      p_data = pvalue.0.mtx)

  return(list(summary = summary,AW.weights = weights_meta))
}


############################
EnrichmentM_Exact<- function (madata, label, censoring = NULL, DB.matrix = DB.matrix,
                              DEgene.number = DEgene.number,
                              size.min = 15,size.max = 500, nperm = 500, gene.pvalues = NULL,
                              resp.type = NULL,gene.name.sort,genes.in.study)
{
  if (!is.null(madata)) {
    genes.in.study = featureNames(madata)
    set2allgenes.mtx = DB.matrix
    gene.common = intersect(featureNames(madata), colnames(set2allgenes.mtx))
    if (nrow(set2allgenes.mtx) > 1) {
      set2allgenes.mtx = as.matrix(set2allgenes.mtx[, gene.common])
    }
    else {
      set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[,
                                                      gene.common]))
    }
    madata = madata[gene.common, ]
    if (resp.type == "twoclass") {
      tstat = genefilter::rowttests(exprs(madata), as.factor(label),
                                    tstatOnly = F)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "multiclass") {
      tstat = genefilter::rowFtests(exprs(madata), as.factor(label),
                                    var.equal = TRUE)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "survival") {
      if (is.null(censoring)) {
        stop("Error: censoring status is null")
      }
      cox.out <- coxfunc(exprs(madata), label, censoring)
      scores <- abs(cox.out$tt)
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else if (resp.type == "continuous") {
      cor.out <- cor.func(exprs(madata), label)
      scores <- abs(cor.out$tt[, 1])
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
  }
  DEgene = gene.name.sort[1:DEgene.number]
  pvalue.set.0 = matrix(NA,nrow = nrow(DB.matrix),ncol = 1)
  rownames(pvalue.set.0) = rownames(DB.matrix)
  for (i in 1:nrow(DB.matrix)){
    count_table<-matrix(0,2,2)
    p_value <-NA
    ####in the gene list and in the pathway
    count_table[1,1]<-sum(DEgene %in% colnames(DB.matrix[,DB.matrix[i,]==1]))
    ####in the gene list but not in the pathway
    count_table[1,2]<-length(DEgene)-count_table[1,1]
    ####not in the gene list but in the pathway
    count_table[2,1]<-sum(genes.in.study%in% colnames(DB.matrix[,DB.matrix[i,]==1]))
    ####not in the gene list and not in the pathway
    count_table[2,2]<-length(genes.in.study)-count_table[2,1]
    if(length(count_table)==4){
      pvalue.set.0[i,1] <- fisher.test(count_table, alternative="greater")$p}
  }
  pvalue.set.0 = pvalue.set.0[pvalue.set.0[,1]!=1,,drop=FALSE]
  qvalue.set.0 = matrix(p.adjust(pvalue.set.0[,1], "BH"),ncol = 1)
  rownames(qvalue.set.0) = rownames(pvalue.set.0)
  return(list(pvalue.set.0 = pvalue.set.0,qvalue.set.0 = qvalue.set.0,DEgene = DEgene))
  #      qvalue.set.0 = matrix(p.adjust(pvalue.set.0[,1], "BH"),ncol = 1)
  #      rownames(qvalue.set.0) = rownames(DB.matrix)
  #      return(list(pvalue.set.0 = pvalue.set.0,qvalue.set.0 = qvalue.set.0,
  #                  pvalue.set.B = pvalue.set.B))
}


######################
MAPE_P_Exact<-function (study, label, censoring.status, DB.matrix, size.min = 15,
                        DEgene.number = DEgene.number, size.max = 500, nperm = 500,
                        stat = NULL, rth.value = NULL,resp.type)
{
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  out = list()
  for (t1 in 1:length(study)) {
    madata = study[[t1]]
    testlabel = madata[[label]]
    out[[t1]] = list()
    if (resp.type == "survival") {
      censoring = madata[[censoring.status]]
    }
    out[[t1]] = EnrichmentM_Exact(madata = madata, label = testlabel, DEgene.number = DEgene.number,
                                  censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                                  size.max = size.max, nperm = nperm, resp.type = resp.type)
  }
  set.common = rownames(out[[1]]$pvalue.set.0)
  for (t1 in 2:length(study)) {
    set.common = intersect(set.common, rownames(out[[t1]]$pvalue.set.0))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  pvalue.0.mtx = matrix(NA, length(set.common), length(study))
  qvalue.0.mtx = matrix(NA, length(set.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[set.common, ]
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[set.common,]
  }
  rownames(qvalue.0.mtx) = set.common
  rownames(pvalue.0.mtx) = set.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,n,1)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,1,n)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = set.common
    qvalue.meta = p.adjust(pvalue.meta, "BH")
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  qvalue.meta = matrix(qvalue.meta,ncol = 1)
  rownames(qvalue.meta) = rownames(pvalue.meta)
  if (stat == "AW Fisher"){
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx,AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx))
  }
}


######################
MAPE_G_Exact = function (study, label, censoring.status, DB.matrix, size.min = 15,
                         DEgene.number = DEgene.number,size.max = 500, nperm = 500,
                         stat = NULL,rth.value = NULL,resp.type)
{
  gene.common = featureNames(study[[1]])
  for (t1 in 2:length(study)) {
    gene.common = intersect(gene.common, featureNames(study[[t1]]))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  madata = list()
  for (t1 in 1:length(study)) {
    madata[[t1]] = study[[t1]][gene.common, ]
  }
  Tstat.p = list()
  out = list()
  for (t1 in 1:length(study)) {
    Tstat.p[[t1]] = list()
    x = exprs(madata[[t1]])
    testlabel = madata[[t1]][[label]]
    if (resp.type == "twoclass") {
      tstat = genefilter::rowttests(x, as.factor(testlabel),
                                    tstatOnly = F)
      Tstat.p[[t1]] = (tstat$p.value)
    }
    else if (resp.type == "multiclass") {
      tstat = genefilter::rowFtests(x, as.factor(testlabel),
                                    var.equal = TRUE)
      Tstat.p[[t1]] = (tstat$p.value)
    }
    else if (resp.type == "survival") {
      if (is.null(censoring.status)) {
        stop("Error: censoring.status is null")
      }
      censoring = madata[[t1]][[censoring.status]]
      Tstat.p[[t1]] = cox.perm.sample(expr = x, testgroup = testlabel,
                                      censoring = censoring, nperm = 500)
    }
    else if (resp.type == "continuous") {
      Tstat.p[[t1]] = reg.perm.sample(expr = x, testgroup = testlabel,
                                      nperm = 500)
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
    out[[t1]] = list()
    if(resp.type %in% c("survival" , "continous")){
      out[[t1]] = pqvalues.compute(Tstat.p[[t1]]$obs, Tstat.p[[t1]]$perms,
                                   Stat.type = "Tstat")
    }
    else {
      out[[t1]]$pvalue.0 = Tstat.p[[t1]]
      out[[t1]]$qvalue.0 = p.adjust(Tstat.p[[t1]],"BH")
    }
  }
  pvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  qvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.0
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.0
  }
  rownames(qvalue.0.mtx) = gene.common
  rownames(pvalue.0.mtx) = gene.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,n,1)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,1,n)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = gene.common
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  gene.pvalues = as.vector(pvalue.meta)
  names(gene.pvalues) = gene.common
  gene.qvalues = p.adjust(gene.pvalues,"BH")
  gene.name.sort = names(sort(gene.pvalues, decreasing = F))
  meta.set = EnrichmentM_Exact(madata = NULL, label = NULL, DEgene.number = DEgene.number,
                               DB.matrix = DB.matrix, size.min = size.min, size.max = size.max,
                               nperm = nperm, gene.pvalues = gene.pvalues,
                               gene.name.sort = gene.name.sort ,genes.in.study = gene.common)
  if(stat == "AW Fisher"){
    return(list(pvalue.meta = meta.set$pvalue.set.0, qvalue.meta = meta.set$qvalue.set.0,
                gene.meta.qvalues = gene.qvalues,
                gene.meta.pvalues = gene.pvalues, gene.study.qvalues = qvalue.0.mtx,
                gene.study.pvalues = pvalue.0.mtx,AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = meta.set$pvalue.set.0, qvalue.meta = meta.set$qvalue.set.0,
                gene.meta.qvalues = gene.qvalues,
                gene.meta.pvalues = gene.pvalues, gene.study.qvalues = qvalue.0.mtx,
                gene.study.pvalues = pvalue.0.mtx))
  }
}


#########################
MAPE_P_gene_KS = function (study, label, censoring.status, DB.matrix, size.min = 15,
                           size.max = 500, nperm = 500, stat = NULL, rth.value = NULL,
                           resp.type)
{
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  out = list()
  for (t1 in 1:length(study)) {
    madata = study[[t1]]
    testlabel = madata[[label]]
    out[[t1]] = list()
    if (resp.type == "survival") {
      censoring = madata[[censoring.status]]
    }
    out[[t1]] = Enrichment_KS_gene(madata = madata, label = testlabel,
                                   censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                                   size.max = size.max, nperm = nperm, resp.type = resp.type)
  }
  set.common = rownames(out[[1]]$pvalue.set.0)
  for (t1 in 2:length(study)) {
    set.common = intersect(set.common, rownames(out[[t1]]$pvalue.set.0))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  pvalue.B.array = array(data = NA, dim = c(length(set.common),
                                            nperm, length(study)))
  dimnames(pvalue.B.array) = list(set.common, paste("perm",
                                                    1:nperm, sep = ""), names(study))
  pvalue.0.mtx = matrix(NA, length(set.common), length(study))
  qvalue.0.mtx = matrix(NA, length(set.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.B.array[, , t1] = out[[t1]]$pvalue.set.B[set.common,]
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[set.common, ]
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[set.common,]
  }
  rownames(qvalue.0.mtx) = set.common
  rownames(pvalue.0.mtx) = set.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), max)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), min)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) sort(x)[rth.value])
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) pchisq(-2 *
                                                                sum(log(x)), DF, lower.tail = T)))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) pchisq(-2 *
                                                              sum(log(x)), DF, lower.tail = T))
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  meta.out = pqvalues.compute(P.0, P.B, Stat.type = "Pvalue")
  return(list(pvalue.meta = meta.out$pvalue.0, qvalue.meta = meta.out$qvalue.0,
              pvalue.meta.B = meta.out$pvalue.B, qvalue.set.study = qvalue.0.mtx,
              pvalue.set.study = pvalue.0.mtx))
}


########################
MAPE_G_gene_KS = function (study, label, censoring.status, DB.matrix, size.min = 15,
                           size.max = 500, nperm = 500, stat = NULL, rth.value = NULL,
                           resp.type)
{
  gene.common = featureNames(study[[1]])
  for (t1 in 2:length(study)) {
    gene.common = intersect(gene.common, featureNames(study[[t1]]))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  madata = list()
  for (t1 in 1:length(study)) {
    madata[[t1]] = study[[t1]][gene.common, ]
  }
  Tstat.perm = list()
  out = list()
  for (t1 in 1:length(study)) {
    Tstat.perm[[t1]] = list()
    x = exprs(madata[[t1]])
    testlabel = madata[[t1]][[label]]
    if (resp.type == "twoclass") {
      Tstat.perm[[t1]] = Tperm.sample(x = x, fac = as.factor(testlabel),
                                      nperm = nperm)
    }
    else if (resp.type == "multiclass") {
      Tstat.perm[[t1]] = F.perm.sample(x = x, fac = as.factor(testlabel),
                                       nperm = nperm)
    }
    else if (resp.type == "survival") {
      if (is.null(censoring.status)) {
        stop("Error: censoring.status is null")
      }
      censoring = madata[[t1]][[censoring.status]]
      Tstat.perm[[t1]] = cox.perm.sample(expr = x, testgroup = testlabel,
                                         censoring = censoring, nperm = nperm)
    }
    else if (resp.type == "continuous") {
      Tstat.perm[[t1]] = reg.perm.sample(expr = x, testgroup = testlabel,
                                         nperm = nperm)
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
    out[[t1]] = list()
    out[[t1]] = pqvalues.compute(Tstat.perm[[t1]]$obs, Tstat.perm[[t1]]$perms,
                                 Stat.type = "Tstat")
  }
  pvalue.B.array = array(data = NA, dim = c(length(gene.common), nperm, length(study)))
  dimnames(pvalue.B.array) = list(gene.common, paste("perm",
                                                     1:nperm, sep = ""), names(study))
  pvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  qvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.B.array[, , t1] = out[[t1]]$pvalue.B
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.0
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.0
  }
  rownames(qvalue.0.mtx) = gene.common
  rownames(pvalue.0.mtx) = gene.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), max)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), min)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) sort(x)[rth.value])
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) pchisq(-2 *
                                                                sum(log(x)), DF, lower.tail = T)))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) pchisq(-2 *
                                                              sum(log(x)), DF, lower.tail = T))
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  meta.out = pqvalues.compute(P.0, P.B, Stat.type = "Pvalue")
  gene.qvalues = as.vector(meta.out$qvalue.0)
  names(gene.qvalues) = rownames(meta.out$qvalue.0)
  gene.pvalues = as.vector(meta.out$pvalue.0)
  names(gene.pvalues) = rownames(meta.out$pvalue.0)
  meta.set = Enrichment_KS_gene(madata = NULL, label = NULL,
                                DB.matrix = DB.matrix, size.min = size.min, size.max = size.max,
                                nperm = nperm, gene.pvalues = gene.qvalues)
  return(list(pvalue.meta = meta.set$pvalue.set.0, qvalue.meta = meta.set$qvalue.set.0,
              pvalue.meta.B = meta.set$pvalue.set.B, gene.meta.qvalues = gene.qvalues,
              gene.meta.pvalues = gene.pvalues, gene.study.qvalues = qvalue.0.mtx,
              gene.study.pvalues = pvalue.0.mtx))
}


######################
MAPE_G_sample_KS = function (study, label, censoring.status = NULL, DB.matrix, size.min = 15,
                             size.max = 500, nperm = 500, stat, rth.value = NULL, resp.type)
{
  gene.common = featureNames(study[[1]])
  for (t1 in 2:length(study)) {
    gene.common = intersect(gene.common, featureNames(study[[t1]]))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  madata = list()
  for (t1 in 1:length(study)) {
    madata[[t1]] = study[[t1]][gene.common, ]
  }
  Tstat.perm = list()
  out = list()
  for (t1 in 1:length(study)) {
    Tstat.perm[[t1]] = list()
    x = exprs(madata[[t1]])
    testlabel = madata[[t1]][[label]]
    if (resp.type == "twoclass") {
      Tstat.perm[[t1]] = Tperm.sample(x = x, fac = as.factor(testlabel),
                                      nperm = nperm)
    }
    else if (resp.type == "multiclass") {
      Tstat.perm[[t1]] = F.perm.sample(x = x, fac = as.factor(testlabel),
                                       nperm = nperm)
    }
    else if (resp.type == "survival") {
      if (is.null(censoring.status)) {
        stop("Error: censoring.aus is null")
      }
      censoring = madata[[t1]][[censoring.status]]
      Tstat.perm[[t1]] = cox.perm.sample(expr = x, testgroup = testlabel,
                                         censoring = censoring, nperm = nperm)
    }
    else if (resp.type == "continuous") {
      Tstat.perm[[t1]] = reg.perm.sample(expr = x, testgroup = testlabel,
                                         nperm = nperm)
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
    out[[t1]] = list()
    out[[t1]] = pqvalues.compute(Tstat.perm[[t1]]$obs, Tstat.perm[[t1]]$perms,
                                 Stat.type = "Tstat")
  }
  pvalue.B.array = array(data = NA, dim = c(length(gene.common),
                                            nperm, length(study)))
  dimnames(pvalue.B.array) = list(gene.common, paste("perm",
                                                     1:nperm, sep = ""), names(study))
  pvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  qvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.B.array[, , t1] = out[[t1]]$pvalue.B
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.0
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.0
  }
  rownames(qvalue.0.mtx) = gene.common
  rownames(pvalue.0.mtx) = gene.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), max)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), min)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) sort(x)[rth.value])
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) pchisq(-2 *
                                                                sum(log(x)), DF, lower.tail = T)))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) pchisq(-2 *
                                                              sum(log(x)), DF, lower.tail = T))
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  colnames(P.0) = "perm0"
  colnames(P.B) = paste("perm", 1:ncol(P.B), sep = "")
  meta.out = pqvalues.compute(P.0, P.B, Stat.type = "Pvalue")
  gene.qvalues = as.vector(meta.out$qvalue.0)
  names(gene.qvalues) = rownames(meta.out$qvalue.0)
  gene.pvalues = as.vector(meta.out$pvalue.0)
  names(gene.pvalues) = rownames(meta.out$pvalue.0)
  score.0 = P.0
  score.B = P.B
  genes.in.study = rownames(score.0)
  set2allgenes.mtx = DB.matrix
  gene.common = intersect(genes.in.study, colnames(set2allgenes.mtx))
  set2allgenes.mtx = set2allgenes.mtx[, gene.common, drop = F]
  score.0 = score.0[gene.common, , drop = F]
  score.B = score.B[gene.common, , drop = F]
  set.length = apply(set2allgenes.mtx, 1, sum)
  idx.1 = which(set.length >= size.min)
  idx.2 = which(set.length <= size.max)
  set.idx = intersect(idx.1, idx.2)
  if (length(set.idx) < 1) {
    stop("no gene sets satisfying size.min<=size<=size.max")
  }
  else {
    set2allgenes.mtx = set2allgenes.mtx[set.idx, , drop = F]
  }
  gene.name.sort = names(sort(abs(score.0[, 1]), decreasing = F))
  order.mtx.1 = set2allgenes.mtx[, gene.name.sort, drop = F]
  order.mtx.0 = (1 - order.mtx.1)
  order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
  order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
  order.mtx = order.mtx.0 + order.mtx.1
  ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                         max))
  ES.B = matrix(NA, nrow(ES.0), ncol(score.B))
  for (t1 in 1:ncol(score.B)) {
    gene.name.sort = names(sort(abs(score.B[, t1]), decreasing = F))
    order.mtx.1 = set2allgenes.mtx[, gene.name.sort, drop = F]
    order.mtx.0 = (1 - order.mtx.1)
    order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
    order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
    order.mtx = order.mtx.0 + order.mtx.1
    order.cumsum = t(apply(order.mtx, 1, cumsum))
    ES.B[, t1] = apply(order.cumsum, 1, max)
  }
  rownames(ES.B) = rownames(order.mtx)
  N.X = apply(set2allgenes.mtx, 1, sum)
  N.Y = ncol(set2allgenes.mtx) - N.X
  N = N.X * N.Y/(N.X + N.Y)
  ES.pval.0 = exp(-2 * N * ES.0^2)
  ES.pval.B = exp(-2 * N * ES.B^2)
  enrich.out = pqvalues.compute(ES.pval.0, ES.pval.B, Stat.type = "Pvalue")
  colnames(enrich.out$pvalue.0) = "MAPE_G_sample"
  colnames(enrich.out$qvalue.0) = "MAPE_G_sample"
  colnames(enrich.out$pvalue.B) = paste("perm", 1:ncol(enrich.out$pvalue.B),
                                        sep = "")
  return(list(pvalue.meta = enrich.out$pvalue.0, qvalue.meta = enrich.out$qvalue.0,
              pvalue.meta.B = enrich.out$pvalue.B))
}


#######################
MAPE_P_sample_KS <- function (study, label, censoring.status = NULL, DB.matrix, size.min = 15,
                              size.max = 500, nperm = 500, stat, rth.value = NULL, resp.type)
{
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  out = list()
  for (t1 in 1:length(study)) {
    madata = study[[t1]]
    testlabel = madata[[label]]
    out[[t1]] = list()
    if (resp.type == "survival") {
      censoring = madata[[censoring.status]]
    }
    out[[t1]] = Enrichment_KS_sample(madata = madata, label = testlabel,
                                     censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                                     size.max = size.max, nperm = nperm, resp.type = resp.type)
  }
  common.pathway = rownames(out[[1]]$pvalue.set.0)
  for (t1 in 1:length(study)) {
    common.pathway = intersect(common.pathway, rownames(out[[t1]]$pvalue.set.0))
  }
  pvalue.B.array = array(data = NA, dim = c(length(common.pathway),
                                            nperm, length(study)))
  pvalue.0.mtx = matrix(NA, length(common.pathway), length(study))
  qvalue.0.mtx = matrix(NA, length(common.pathway), length(study))
  rownames(pvalue.0.mtx) = common.pathway
  colnames(pvalue.0.mtx) = names(study)
  rownames(qvalue.0.mtx) = common.pathway
  colnames(qvalue.0.mtx) = names(study)
  dimnames(pvalue.B.array) = list(common.pathway, paste("perm",
                                                        1:nperm, sep = ""), names(study))
  for (t1 in 1:length(study)) {
    pvalue.B.array[, , t1] = out[[t1]]$pvalue.set.B[common.pathway,
                                                    ]
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[common.pathway,
                                                ]
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[common.pathway,
                                                ]
  }
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    rownames(P.0) = rownames(pvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), max)
    rownames(P.B) = rownames(pvalue.0.mtx)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    rownames(P.0) = rownames(pvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), min)
    rownames(P.B) = rownames(pvalue.0.mtx)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    rownames(P.0) = rownames(pvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) sort(x)[rth.value])
    rownames(P.B) = rownames(pvalue.0.mtx)
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) pchisq(-2 *
                                                                sum(log(x)), DF, lower.tail = T)))
    rownames(P.0) = rownames(pvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) pchisq(-2 *
                                                              sum(log(x)), DF, lower.tail = T))
    rownames(P.B) = rownames(pvalue.0.mtx)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  colnames(P.0) = "perm0"
  colnames(P.B) = paste("perm", 1:ncol(P.B), sep = "")
  meta.out = pqvalues.compute(P.0, P.B, Stat.type = "Pvalue")
  colnames(meta.out$pvalue.0) = "MAPE_P_sample"
  colnames(meta.out$qvalue.0) = "MAPE_P_sample"
  return(list(pvalue.meta = meta.out$pvalue.0, qvalue.meta = meta.out$qvalue.0,
              pvalue.meta.B = meta.out$pvalue.B, pvalue.set.study = pvalue.0.mtx,
              qvalue.set.study = qvalue.0.mtx))
}


######################
MAPE_I_KS <- function (MAP_GENE.obj, MAP_SET.obj)
{
  set.common = intersect(rownames(MAP_GENE.obj$pvalue.meta),
                         rownames(MAP_SET.obj$pvalue.meta))
  nperm = ncol(MAP_SET.obj$pvalue.meta.B)
  pvalue.set.B.array = array(data = NA, dim = c(length(set.common),
                                                nperm, 2))
  dimnames(pvalue.set.B.array) = list(set.common, paste("perm",
                                                        1:nperm, sep = ""), c(1,2))
  pvalue.set.0.mtx = matrix(NA, length(set.common), 2)
  pvalue.set.B.array[, , 1] = MAP_SET.obj$pvalue.meta.B[set.common,]
  pvalue.set.B.array[, , 2] = MAP_GENE.obj$pvalue.meta.B[set.common,]
  pvalue.set.0.mtx[, 1] = MAP_SET.obj$pvalue.meta[set.common,]
  pvalue.set.0.mtx[, 2] = MAP_GENE.obj$pvalue.meta[set.common,]
  rownames(pvalue.set.0.mtx) = set.common
  minP.0 = as.matrix(apply(pvalue.set.0.mtx, 1, min))
  rownames(minP.0) = set.common
  minP.B = apply(pvalue.set.B.array, c(1, 2), min)
  rownames(minP.B) = set.common
  meta.out = pqvalues.compute(minP.0, minP.B, Stat.type = "Pvalue")
  return(list(pvalue.meta = meta.out$pvalue.0, qvalue.meta = meta.out$qvalue.0))
}


########################
MAPE_I <- function (MAP_GENE.obj, MAP_SET.obj)
{
  set.common = intersect(rownames(MAP_GENE.obj$qvalue.meta),
                         rownames(MAP_SET.obj$qvalue.meta))

  pvalue.set.0.mtx = matrix(NA, length(set.common), 2)
  pvalue.set.0.mtx[, 1] = MAP_SET.obj$pvalue.meta[set.common,]
  pvalue.set.0.mtx[, 2] = MAP_GENE.obj$pvalue.meta[set.common,]
  rownames(pvalue.set.0.mtx) = set.common
  minP.0 = as.matrix(apply(pvalue.set.0.mtx, 1, min))
  rownames(minP.0) = set.common
  pvalue.meta = pbeta(minP.0,1,2)
  rownames(pvalue.meta) = rownames(minP.0)
  qvalue.meta = p.adjust(pvalue.meta,"BH")
  qvalue.meta = matrix(qvalue.meta,ncol = 1)
  rownames(qvalue.meta) = rownames(pvalue.meta)
  return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta))
}



########################
Enrichment_KS_gene <- function (madata, label, censoring = NULL, DB.matrix, size.min = 15,
                                size.max = 500, nperm = 500, gene.pvalues = NULL, resp.type = NULL)
{
  if (!is.null(madata)) {
    genes.in.study = featureNames(madata)
    set2allgenes.mtx = DB.matrix
    gene.common = intersect(featureNames(madata), colnames(set2allgenes.mtx))
    if (nrow(set2allgenes.mtx) > 1) {
      set2allgenes.mtx = as.matrix(set2allgenes.mtx[, gene.common])
    }
    else {
      set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[,
                                                      gene.common]))
    }
    madata = madata[gene.common, ]
    if (resp.type == "twoclass") {
      tstat = genefilter::rowttests(exprs(madata), as.factor(label),
                                    tstatOnly = F)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "multiclass") {
      tstat = genefilter::rowFtests(exprs(madata), as.factor(label),
                                    var.equal = TRUE)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "survival") {
      if (is.null(censoring)) {
        stop("Error: censoring status is null")
      }
      cox.out <- coxfunc(exprs(madata), label, censoring)
      scores <- abs(cox.out$tt)
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else if (resp.type == "continuous") {
      cor.out <- cor.func(exprs(madata), label)
      scores <- abs(cor.out$tt[, 1])
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
  }
  if (!is.null(gene.pvalues)) {
    if ((!is.vector(gene.pvalues)) | (is.null(names(gene.pvalues))))
      stop("gene.pvalues should be a vector with gene names")
    genes.in.study = names(gene.pvalues)
    set2allgenes.mtx = DB.matrix
    gene.common = intersect(genes.in.study, colnames(set2allgenes.mtx))
    if (nrow(set2allgenes.mtx) > 1) {
      set2allgenes.mtx = as.matrix(set2allgenes.mtx[, gene.common])
    }
    else {
      set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[,
                                                      gene.common]))
    }
    gene.pvalues = gene.pvalues[gene.common]
    gene.name.sort = names(sort(gene.pvalues, decreasing = F))
  }
  set.length = apply(set2allgenes.mtx, 1, sum)
  idx.1 = which(set.length >= size.min)
  idx.2 = which(set.length <= size.max)
  set.idx = intersect(idx.1, idx.2)
  if (length(set.idx) < 1)
    stop("no gene sets satisfying size.min<=size<=size.max")
  if (length(set.idx) > 1) {
    set2allgenes.mtx = set2allgenes.mtx[set.idx, ]
  }
  else {
    set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[set.idx,
                                                    ]))
  }

  if (nrow(set2allgenes.mtx) > 1) {
    order.mtx.1 = (set2allgenes.mtx[, gene.name.sort])
  }
  else {
    order.mtx.1 = t(as.matrix(set2allgenes.mtx[, gene.name.sort]))
  }
  order.mtx.0 = (1 - order.mtx.1)
  order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
  order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
  order.mtx = order.mtx.0 + order.mtx.1
  ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                         max))
  ES.B = matrix(NA, nrow(ES.0), nperm)
  for (t1 in 1:nperm) {
    if (nrow(order.mtx) > 1) {
      order.mtx.perm = order.mtx[, sample(ncol(order.mtx))]
    }
    else {
      order.mtx.perm = t(as.matrix(order.mtx[, sample(ncol(order.mtx))]))
    }
    order.cumsum = t(apply(order.mtx.perm, 1, cumsum))
    ES.B[, t1] = apply(order.cumsum, 1, max)
  }
  rownames(ES.B) = rownames(order.mtx)
  N.X = apply(set2allgenes.mtx, 1, sum)
  N.Y = ncol(set2allgenes.mtx) - N.X
  N = N.X * N.Y/(N.X + N.Y)
  enrich.out = pqvalues.compute(ES.0, ES.B, Stat.type = "Tstat")
  return(list(pvalue.set.0 = enrich.out$pvalue.0, pvalue.set.B = enrich.out$pvalue.B,
              qvalue.set.0 = enrich.out$qvalue.0))
}


######################
Enrichment_KS_sample <- function (madata, label, censoring = NULL, DB.matrix, size.min = 15,
                                  size.max = 500, nperm = 500, resp.type = NULL)
{
  genes.in.study = featureNames(madata)
  set2allgenes.mtx = DB.matrix
  gene.common = intersect(featureNames(madata), colnames(set2allgenes.mtx))
  set2allgenes.mtx = as.matrix(set2allgenes.mtx[, gene.common],
                               drop = F)
  madata = madata[gene.common, ]
  set.length = apply(set2allgenes.mtx, 1, sum)
  idx.1 = which(set.length >= size.min)
  idx.2 = which(set.length <= size.max)
  set.idx = intersect(idx.1, idx.2)
  if (length(set.idx) < 1) {
    stop("no gene sets satisfying size.min<=size<=size.max")
  }
  else {
    set2allgenes.mtx = set2allgenes.mtx[set.idx, , drop = F]
  }
  x = exprs(madata)
  testlabel = label
  if (resp.type == "twoclass") {
    Tstat.perm = Tperm.sample(x = x, fac = as.factor(testlabel),
                              nperm = nperm)
  }
  else if (resp.type == "multiclass") {
    Tstat.perm = F.perm.sample(x = x, fac = as.factor(testlabel),
                               nperm = nperm)
  }
  else if (resp.type == "survival") {
    if (is.null(censoring)) {
      stop("Error: censoring.status is null")
    }
    Tstat.perm = cox.perm.sample(expr = x, testgroup = testlabel,
                                 censoring = censoring, nperm = nperm)
  }
  else if (resp.type == "continuous") {
    Tstat.perm = reg.perm.sample(expr = x, testgroup = testlabel,
                                 nperm = nperm)
  }
  else {
    stop("Error: Wrong input augument for resp.type")
  }
  score.0 = Tstat.perm$obs
  score.B = Tstat.perm$perms
  gene.name.sort = names(sort(abs(score.0), decreasing = T))
  order.mtx.1 = set2allgenes.mtx[, gene.name.sort, drop = F]
  order.mtx.0 = (1 - order.mtx.1)
  order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
  order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
  order.mtx = order.mtx.0 + order.mtx.1
  ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                         max))
  ES.B = matrix(NA, nrow(ES.0), ncol(score.B))
  for (t1 in 1:ncol(score.B)) {
    gene.name.sort = names(sort(abs(score.B[, t1]), decreasing = T))
    order.mtx.1 = set2allgenes.mtx[, gene.name.sort, drop = F]
    order.mtx.0 = (1 - order.mtx.1)
    order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
    order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
    order.mtx = order.mtx.0 + order.mtx.1
    order.cumsum = t(apply(order.mtx, 1, cumsum))
    ES.B[, t1] = apply(order.cumsum, 1, max)
  }
  rownames(ES.B) = rownames(order.mtx)
  N.X = apply(set2allgenes.mtx, 1, sum)
  N.Y = ncol(set2allgenes.mtx) - N.X
  N = N.X * N.Y/(N.X + N.Y)
  enrich.out = pqvalues.compute(ES.0, ES.B, Stat.type = "Tstat")
  return(list(pvalue.set.0 = enrich.out$pvalue.0, pvalue.set.B = enrich.out$pvalue.B,
              qvalue.set.0 = enrich.out$qvalue.0))
}


######################
Enrichment_KS <- function (madata, label, censoring = NULL, DB.matrix, size.min = 15,
                           size.max = 500, gene.pvalues = NULL, resp.type = NULL)
{
  if (!is.null(madata)) {
    genes.in.study = featureNames(madata)
    set2allgenes.mtx = DB.matrix
    gene.common = intersect(featureNames(madata), colnames(set2allgenes.mtx))
    if (nrow(set2allgenes.mtx) > 1) {
      set2allgenes.mtx = as.matrix(set2allgenes.mtx[, gene.common])
    }
    else {
      set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[,
                                                      gene.common]))
    }
    madata = madata[gene.common, ]
    if (resp.type == "twoclass") {
      tstat = genefilter::rowttests(exprs(madata), as.factor(label),
                                    tstatOnly = F)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "multiclass") {
      tstat = genefilter::rowFtests(exprs(madata), as.factor(label),
                                    var.equal = TRUE)
      p.values = (tstat$p.value)
      names(p.values) = rownames(tstat)
      gene.name.sort = names(sort(p.values, decreasing = F))
    }
    else if (resp.type == "survival") {
      if (is.null(censoring)) {
        stop("Error: censoring status is null")
      }
      cox.out <- coxfunc(exprs(madata), label, censoring)
      scores <- abs(cox.out$tt)
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else if (resp.type == "continuous") {
      cor.out <- cor.func(exprs(madata), label)
      scores <- abs(cor.out$tt[, 1])
      gene.name.sort = names(sort(scores, decreasing = T))
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
  }
  if (!is.null(gene.pvalues)) {
    if ((!is.vector(gene.pvalues)) | (is.null(names(gene.pvalues))))
      stop("gene.pvalues should be a vector with gene names")
    genes.in.study = names(gene.pvalues)
    set2allgenes.mtx = DB.matrix
    gene.common = intersect(genes.in.study, colnames(set2allgenes.mtx))
    if (nrow(set2allgenes.mtx) > 1) {
      set2allgenes.mtx = as.matrix(set2allgenes.mtx[, gene.common])
    }
    else {
      set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[,
                                                      gene.common]))
    }
    gene.pvalues = gene.pvalues[gene.common]
    gene.name.sort = names(sort(gene.pvalues, decreasing = F))
  }
  set.length = apply(set2allgenes.mtx, 1, sum)
  idx.1 = which(set.length >= size.min)
  idx.2 = which(set.length <= size.max)
  set.idx = intersect(idx.1, idx.2)
  if (length(set.idx) < 1)
    stop("no gene sets satisfying size.min<=size<=size.max")
  if (length(set.idx) > 1) {
    set2allgenes.mtx = set2allgenes.mtx[set.idx, ]
  }
  else {
    set2allgenes.mtx = t(as.matrix(set2allgenes.mtx[set.idx,
                                                    ]))
  }
  if (nrow(set2allgenes.mtx) > 1) {
    order.mtx.1 = (set2allgenes.mtx[, gene.name.sort])
  }
  else {
    order.mtx.1 = t(as.matrix(set2allgenes.mtx[, gene.name.sort]))
  }
  order.mtx.0 = (1 - order.mtx.1)
  n_hit = rowSums(order.mtx.1)
  n_miss = rowSums(order.mtx.0)
  n_genes = ncol(order.mtx.1)
  nn = (n_hit*n_miss)/n_genes
  order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
  order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
  order.mtx = order.mtx.0 + order.mtx.1
  ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                         max))
  pvalue.0 = matrix(exp(-2*nn*(ES.0[,1]^2)), ncol = 1)
  qvalue.0 = matrix(p.adjust(pvalue.0[,1],"BH"), ncol = 1)
  rownames(pvalue.0) = rownames(ES.0)
  rownames(qvalue.0) = rownames(ES.0)
  return(list(pvalue.set.0 = pvalue.0, qvalue.set.0 = qvalue.0))
}



##########################
MAPE_P_KS = function (study, label, censoring.status, DB.matrix, size.min = 15,
                      size.max = 500, stat = NULL, rth.value = NULL,
                      resp.type)
{
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  out = list()
  for (t1 in 1:length(study)) {
    madata = study[[t1]]
    testlabel = madata[[label]]
    out[[t1]] = list()
    if (resp.type == "survival") {
      censoring = madata[[censoring.status]]
    }
    out[[t1]] = Enrichment_KS (madata = madata, label = testlabel,
                               censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                               size.max = size.max, resp.type = resp.type)
  }
  set.common = rownames(out[[1]]$pvalue.set.0)
  for (t1 in 2:length(study)) {
    set.common = intersect(set.common, rownames(out[[t1]]$pvalue.set.0))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  pvalue.0.mtx = matrix(NA, length(set.common), length(study))
  qvalue.0.mtx = matrix(NA, length(set.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[set.common, ]
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[set.common,]
  }
  rownames(qvalue.0.mtx) = set.common
  rownames(pvalue.0.mtx) = set.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,n,1)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,1,n)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = set.common
    qvalue.meta = p.adjust(pvalue.meta, "BH")
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  qvalue.meta = matrix(qvalue.meta,ncol = 1)
  rownames(qvalue.meta) = rownames(pvalue.meta)

  if(stat == "AW Fisher"){
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx,AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx))
  }
}



##########################
MAPE_G_KS = function (study, label, censoring.status, DB.matrix, size.min = 15,
                      size.max = 500, stat = NULL, rth.value = NULL,
                      resp.type)
{
  gene.common = featureNames(study[[1]])
  for (t1 in 2:length(study)) {
    gene.common = intersect(gene.common, featureNames(study[[t1]]))
  }
  if (is.null(names(study)))
    names(study) = paste("study.", 1:length(study), sep = "")
  madata = list()
  for (t1 in 1:length(study)) {
    madata[[t1]] = study[[t1]][gene.common, ]
  }
  Tstat.p = list()
  out = list()
  for (t1 in 1:length(study)) {
    Tstat.p[[t1]] = list()
    x = exprs(madata[[t1]])
    testlabel = madata[[t1]][[label]]
    if (resp.type == "twoclass") {
      tstat = genefilter::rowttests(x, as.factor(testlabel),
                                    tstatOnly = F)
      Tstat.p[[t1]] = (tstat$p.value)
    }
    else if (resp.type == "multiclass") {
      tstat = genefilter::rowFtests(x, as.factor(testlabel),
                                    var.equal = TRUE)
      Tstat.p[[t1]] = (tstat$p.value)
    }
    else if (resp.type == "survival") {
      if (is.null(censoring.status)) {
        stop("Error: censoring.status is null")
      }
      censoring = madata[[t1]][[censoring.status]]
      Tstat.p[[t1]] = cox.perm.sample(expr = x, testgroup = testlabel,
                                      censoring = censoring, nperm = 500)
    }
    else if (resp.type == "continuous") {
      Tstat.p[[t1]] = reg.perm.sample(expr = x, testgroup = testlabel,
                                      nperm = 500)
    }
    else {
      stop("Error: Wrong input augument for resp.type")
    }
    out[[t1]] = list()
    if(resp.type %in% c("survival" , "continous")){
      out[[t1]] = pqvalues.compute(Tstat.p[[t1]]$obs, Tstat.p[[t1]]$perms,
                                   Stat.type = "Tstat")
    }
    else {
      out[[t1]]$pvalue.0 = Tstat.p[[t1]]
      out[[t1]]$qvalue.0 = p.adjust(Tstat.p[[t1]],"BH")
    }
  }
  pvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  qvalue.0.mtx = matrix(NA, length(gene.common), length(study))
  for (t1 in 1:length(study)) {
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.0
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.0
  }
  rownames(qvalue.0.mtx) = gene.common
  rownames(pvalue.0.mtx) = gene.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,n,1)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,1,n)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
  }
  else if (stat == "Fisher") {
    DF = 2 * length(study)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = gene.common
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  gene.pvalues = as.vector(pvalue.meta)
  names(gene.pvalues) = gene.common
  gene.qvalues = p.adjust(gene.pvalues,"BH")
  meta.set = Enrichment_KS(madata = NULL, label = NULL,
                           DB.matrix = DB.matrix, size.min = size.min, size.max = size.max,
                           gene.pvalues = gene.pvalues)
  if(stat == "AW Fisher"){
    return(list(pvalue.meta = meta.set$pvalue.set.0, qvalue.meta = meta.set$qvalue.set.0,
                gene.meta.qvalues = gene.qvalues,
                gene.meta.pvalues = gene.pvalues, gene.study.qvalues = qvalue.0.mtx,
                gene.study.pvalues = pvalue.0.mtx, AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = meta.set$pvalue.set.0, qvalue.meta = meta.set$qvalue.set.0,
                gene.meta.qvalues = gene.qvalues,
                gene.meta.pvalues = gene.pvalues, gene.study.qvalues = qvalue.0.mtx,
                gene.study.pvalues = pvalue.0.mtx))
  }
}



######################
MAPE_P_KS_DE = function (ind.p = ind.p, DB.matrix, size.min = 15, gene.common,
                         size.max = 500, stat = NULL, rth.value = NULL)
{
  out = list()
  out2 = list()
  for (t1 in 1:ncol(ind.p)) {
    gene.name.sort = names(sort(ind.p[,t1],decreasing = F))
    gene.name.sort = gene.name.sort[gene.name.sort%in%gene.common]
    gene.name.sort = toupper(gene.name.sort)
    set2allgenes.mtx = DB.matrix
    order.mtx.1 = (set2allgenes.mtx[, gene.name.sort])
    order.mtx.0 = (1 - order.mtx.1)
    n_hit = rowSums(order.mtx.1)
    n_miss = rowSums(order.mtx.0)
    n_genes = ncol(order.mtx.1)
    nn = (n_hit*n_miss)/n_genes
    order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
    order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
    order.mtx = order.mtx.0 + order.mtx.1
    ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                           max))
    pvalue.0 = matrix(exp(-2*nn*(ES.0[,1]^2)), ncol = 1)
    rownames(pvalue.0) = rownames(ES.0)
    out[[t1]] = pvalue.0
    out2[[t1]] = as.matrix(p.adjust(pvalue.0,"BH"),ncol = 1)
    rownames(out2[[t1]]) = rownames(pvalue.0)
  }
  set.common = rownames(out[[1]])
  pvalue.0.mtx = matrix(NA, length(set.common), ncol(ind.p))
  qvalue.0.mtx = matrix(NA, length(set.common), ncol(ind.p))
  for (t1 in 1:ncol(ind.p)) {
    pvalue.0.mtx[, t1] = out[[t1]][set.common, ]
    qvalue.0.mtx[, t1] = out2[[t1]][set.common,]
  }
  rownames(qvalue.0.mtx) = set.common
  rownames(pvalue.0.mtx) = set.common
  qvalue.0.mtx = qvalue.0.mtx[complete.cases(qvalue.0.mtx),]
  pvalue.0.mtx = pvalue.0.mtx[complete.cases(pvalue.0.mtx),]
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,n,1)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,1,n)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "Fisher") {
    DF = 2 * ncol(ind.p)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = rownames(pvalue.0.mtx)
    qvalue.meta = p.adjust(pvalue.meta, "BH")
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  qvalue.meta = matrix(qvalue.meta,ncol = 1)
  rownames(qvalue.meta) = rownames(pvalue.meta)
  if(stat == "AW Fisher"){
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx,AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx))
  }
}


####################
CPI_KS <- function (study, label, censoring.status, DB.matrix, size.min = 15,
                    size.max = 500, stat = NULL, rth.value = NULL, resp.type)
{if (is.null(names(study)))
  names(study) = paste("study.", 1:length(study), sep = "")
out = list()
for (t1 in 1:length(study)) {
  madata = study[[t1]]
  testlabel = madata[[label]]
  out[[t1]] = list()
  if (resp.type == "survival") {
    censoring = madata[[censoring.status]]
  }
  out[[t1]] = Enrichment_KS (madata = madata, label = testlabel,
                             censoring = censoring, DB.matrix = DB.matrix, size.min = size.min,
                             size.max = size.max, resp.type = resp.type)
}
set.common = rownames(out[[1]]$pvalue.set.0)
for (t1 in 2:length(study)) {
  set.common = intersect(set.common, rownames(out[[t1]]$pvalue.set.0))
}
if (is.null(names(study)))
  names(study) = paste("study.", 1:length(study), sep = "")

pvalue.0.mtx = matrix(NA, length(set.common), length(study))
qvalue.0.mtx = matrix(NA, length(set.common), length(study))
for (t1 in 1:length(study)) {
  pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[set.common,]
  qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[set.common,]
}
rownames(qvalue.0.mtx) = set.common
rownames(pvalue.0.mtx) = set.common
rm(out)

AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
p_value_meta = AW.fisher$pvalues
weights.meta = AW.fisher$weights

q_value_meta = p.adjust(p_value_meta, "BH")
summary<-data.frame(q_value_meta = q_value_meta,p_value_meta = p_value_meta,
                    p_data = pvalue.0.mtx)
return(list(summary = summary,AW.weights = weights.meta))
}



#########################
MAPE_P_Exact_DE = function (ind.p = ind.p, DB.matrix, size.min = 15, gene.common,
                            size.max = 500, stat = NULL, rth.value = NULL, DEgene.number)
{
  out = list()
  out2 = list()
  out3 = list()
  for (t1 in 1:ncol(ind.p)) {
    genes.in.study = names(ind.p[,t1])
    gene.name.sort = names(sort(ind.p[,t1],decreasing = F))
    gene.name.sort = gene.name.sort[gene.name.sort%in%gene.common]
    gene.name.sort = toupper(gene.name.sort)
    DEgene = gene.name.sort[1:DEgene.number]
    pvalue.0 = matrix(NA,nrow = nrow(DB.matrix),ncol = 1)
    rownames(pvalue.0) = rownames(DB.matrix)
    for (i in 1:nrow(DB.matrix)){
      count_table<-matrix(0,2,2)
      p_value <-NA
      ####in the gene list and in the pathway
      count_table[1,1]<-sum(DEgene %in% colnames(DB.matrix[,DB.matrix[i,]==1]))
      ####in the gene list but not in the pathway
      count_table[1,2]<-length(DEgene)-count_table[1,1]
      ####not in the gene list but in the pathway
      count_table[2,1]<-sum(genes.in.study%in% colnames(DB.matrix[,DB.matrix[i,]==1]))
      ####not in the gene list and not in the pathway
      count_table[2,2]<-length(genes.in.study)-count_table[2,1]
      if(length(count_table)==4){
        pvalue.0[i,1] <- fisher.test(count_table, alternative="greater")$p}
    }
    #    pvalue.0 = pvalue.0[pvalue.0[,1]!=1,,drop=FALSE]
    qvalue.0 = matrix(p.adjust(pvalue.0[,1], "BH"),ncol = 1)
    rownames(qvalue.0) = rownames(pvalue.0)
    out[[t1]] = pvalue.0
    out2[[t1]] = as.matrix(p.adjust(pvalue.0,"BH"),ncol = 1)
    rownames(out2[[t1]]) = rownames(pvalue.0)
    out3[[t1]] = DEgene
  }
  set.common = rownames(out[[1]])
  pvalue.0.mtx = matrix(NA, length(set.common), ncol(ind.p))
  qvalue.0.mtx = matrix(NA, length(set.common), ncol(ind.p))
  for (t1 in 1:ncol(ind.p)) {
    pvalue.0.mtx[, t1] = out[[t1]][set.common, ]
    qvalue.0.mtx[, t1] = out2[[t1]][set.common,]
    pvalue.0.mtx = pvalue.0.mtx
  }
  rownames(qvalue.0.mtx) = set.common
  rownames(pvalue.0.mtx) = set.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,n,1)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,1,n)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "Fisher") {
    DF = 2 * ncol(ind.p)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(pvalue.0.mtx)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
    qvalue.meta = p.adjust(pvalue.meta,"BH")
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = set.common
    qvalue.meta = p.adjust(pvalue.meta, "BH")
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  qvalue.meta = matrix(qvalue.meta,ncol = 1)
  rownames(qvalue.meta) = rownames(pvalue.meta)
  genelist = list()
  for (t1 in 1:ncol(ind.p)) {
    genelist[[t1]] = out3[[t1]]
  }
  if(stat=="AW Fisher"){
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx,genelist = genelist, AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = pvalue.meta, qvalue.meta = qvalue.meta,
                qvalue.set.study = qvalue.0.mtx, pvalue.set.study = pvalue.0.mtx,genelist = genelist))
  }
}



##########################
MAPE_P_KS_gene = function (ind.p = ind.p, DB.matrix, size.min = 15, gene.common,
                           size.max = 500, stat = NULL, rth.value = NULL,nperm = nperm)
{
  out = list()
  for (t1 in 1:ncol(ind.p)) {
    gene.name.sort = names(sort(ind.p[,t1],decreasing = F))
    gene.name.sort = gene.name.sort[gene.name.sort%in%gene.common]
    gene.name.sort = toupper(gene.name.sort)
    set2allgenes.mtx = DB.matrix
    order.mtx.1 = (set2allgenes.mtx[, gene.name.sort])
    order.mtx.0 = (1 - order.mtx.1)
    order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
    order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
    order.mtx = order.mtx.0 + order.mtx.1
    ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                           max))
    ES.B = matrix(NA, nrow(ES.0), nperm)
    for (i in 1:nperm) {
      if (nrow(order.mtx) > 1) {
        order.mtx.perm = order.mtx[, sample(ncol(order.mtx))]
      }
      else {
        order.mtx.perm = t(as.matrix(order.mtx[, sample(ncol(order.mtx))]))
      }
      order.cumsum = t(apply(order.mtx.perm, 1, cumsum))
      ES.B[, i] = apply(order.cumsum, 1, max)
    }
    rownames(ES.B) = rownames(order.mtx)
    N.X = apply(set2allgenes.mtx, 1, sum)
    N.Y = ncol(set2allgenes.mtx) - N.X
    N = N.X * N.Y/(N.X + N.Y)
    enrich.out = pqvalues.compute(ES.0, ES.B, Stat.type = "Tstat")
    out[[t1]] = list()
    out[[t1]]$pvalue.set.0 = enrich.out$pvalue.0
    out[[t1]]$pvalue.set.B = enrich.out$pvalue.B
    out[[t1]]$qvalue.set.0 = enrich.out$qvalue.0
  }
  set.common = rownames(out[[1]]$pvalue.set.0)
  for (t1 in 2:ncol(ind.p)) {
    set.common = intersect(set.common, rownames(out[[t1]]$pvalue.set.0))
  }
  pvalue.B.array = array(data = NA, dim = c(length(set.common),
                                            nperm, ncol(ind.p)))
  dimnames(pvalue.B.array) = list(set.common, paste("perm",
                                                    1:nperm, sep = ""), colnames(ncol(ind.p)))
  pvalue.0.mtx = matrix(NA, length(set.common), ncol(ind.p))
  qvalue.0.mtx = matrix(NA, length(set.common), ncol(ind.p))
  for (t1 in 1:ncol(ind.p)) {
    pvalue.B.array[, , t1] = out[[t1]]$pvalue.set.B[set.common,]
    pvalue.0.mtx[, t1] = out[[t1]]$pvalue.set.0[set.common, ]
    qvalue.0.mtx[, t1] = out[[t1]]$qvalue.set.0[set.common,]
  }
  rownames(qvalue.0.mtx) = set.common
  rownames(pvalue.0.mtx) = set.common
  rm(out)
  if (stat == "maxP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, max))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), max)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, min))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), min)
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) sort(x)[rth.value]))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) sort(x)[rth.value])
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else if (stat == "Fisher") {
    DF = 2 * ncol(ind.p)
    P.0 = as.matrix(apply(pvalue.0.mtx, 1, function(x) pchisq(-2 *
                                                                sum(log(x)), DF, lower.tail = T)))
    rownames(P.0) = rownames(qvalue.0.mtx)
    P.B = apply(pvalue.B.array, c(1, 2), function(x) pchisq(-2 *
                                                              sum(log(x)), DF, lower.tail = T))
    rownames(P.B) = rownames(qvalue.0.mtx)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  meta.out = pqvalues.compute(P.0, P.B, Stat.type = "Pvalue")
  return(list(pvalue.meta = meta.out$pvalue.0, qvalue.meta = meta.out$qvalue.0,
              pvalue.meta.B = meta.out$pvalue.B, qvalue.set.study = qvalue.0.mtx,
              pvalue.set.study = pvalue.0.mtx))
}


########################
MAPE_G_Exact_DE = function (ind.p = ind.p,DB.matrix,gene.common,stat = NULL, rth.value = NULL,
                            size.min = 15, size.max = 500,DEgene.number)
{
  if (stat == "maxP") {
    P.0 = as.matrix(apply(ind.p, 1, max))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,n,1)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(ind.p, 1, min))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,1,n)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(ind.p, 1, function(x) sort(x)[rth.value]))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
  }
  else if (stat == "Fisher") {
    DF = 2 * ncol(ind.p)
    P.0 = as.matrix(apply(ind.p, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(ind.p)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = rownames(ind.p)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  gene.name.sort = names(sort(pvalue.meta[,1],decreasing = F))
  genes.in.study = gene.name.sort
  gene.name.sort = gene.name.sort[gene.name.sort%in%gene.common]
  gene.name.sort = toupper(gene.name.sort)
  DEgene = gene.name.sort[1:DEgene.number]
  pvalue.0 = matrix(NA,nrow = nrow(DB.matrix),ncol = 1)
  rownames(pvalue.0) = rownames(DB.matrix)
  for (i in 1:nrow(DB.matrix)){
    count_table<-matrix(0,2,2)
    p_value <-NA
    ####in the gene list and in the pathway
    count_table[1,1]<-sum(DEgene %in% colnames(DB.matrix[,DB.matrix[i,]==1]))
    ####in the gene list but not in the pathway
    count_table[1,2]<-length(DEgene)-count_table[1,1]
    ####not in the gene list but in the pathway
    count_table[2,1]<-sum(genes.in.study%in% colnames(DB.matrix[,DB.matrix[i,]==1]))
    ####not in the gene list and not in the pathway
    count_table[2,2]<-length(genes.in.study)-count_table[2,1]
    if(length(count_table)==4){
      pvalue.0[i,1] <- fisher.test(count_table, alternative="greater")$p}
  }
  #  pvalue.0 = pvalue.0[pvalue.0[,1]!=1,,drop=FALSE]
  qvalue.0 = matrix(p.adjust(pvalue.0[,1], "BH"),ncol = 1)
  rownames(qvalue.0) = rownames(pvalue.0)
  pvalue.set.0 = pvalue.0
  qvalue.set.0 = as.matrix(p.adjust(pvalue.0,"BH"),ncol = 1)
  rownames(qvalue.set.0) = rownames(pvalue.0)
  if(stat == "AW Fisher"){
    return(list(pvalue.meta = pvalue.set.0, qvalue.meta = qvalue.set.0, AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = pvalue.set.0, qvalue.meta = qvalue.set.0))
  }
}



#########################
MAPE_G_KS_gene = function (ind.p=ind.p,DB.matrix,gene.common,stat,rth.value,
                           size.min = 15, size.max = 500,nperm)
{
  if (stat == "maxP") {
    P.0 = as.matrix(apply(ind.p, 1, max))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,n,1)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(ind.p, 1, min))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,1,n)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(ind.p, 1, function(x) sort(x)[rth.value]))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
  }
  else if (stat == "Fisher") {
    DF = 2 * ncol(ind.p)
    P.0 = as.matrix(apply(ind.p, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(ind.p)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = rownames(ind.p)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  gene.name.sort = names(sort(pvalue.meta[,1],decreasing = F))
  gene.name.sort = gene.name.sort[gene.name.sort%in%gene.common]
  gene.name.sort = toupper(gene.name.sort)
  set2allgenes.mtx = DB.matrix
  order.mtx.1 = (set2allgenes.mtx[, gene.name.sort])
  order.mtx.0 = (1 - order.mtx.1)
  order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
  order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
  order.mtx = order.mtx.0 + order.mtx.1
  ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                         max))
  ES.B = matrix(NA, nrow(ES.0), nperm)
  for (t1 in 1:nperm) {
    if (nrow(order.mtx) > 1) {
      order.mtx.perm = order.mtx[, sample(ncol(order.mtx))]
    }
    else {
      order.mtx.perm = t(as.matrix(order.mtx[, sample(ncol(order.mtx))]))
    }
    order.cumsum = t(apply(order.mtx.perm, 1, cumsum))
    ES.B[, t1] = apply(order.cumsum, 1, max)
  }
  rownames(ES.B) = rownames(order.mtx)
  N.X = apply(set2allgenes.mtx, 1, sum)
  N.Y = ncol(set2allgenes.mtx) - N.X
  N = N.X * N.Y/(N.X + N.Y)
  enrich.out = pqvalues.compute(ES.0, ES.B, Stat.type = "Tstat")
  if(stat == "AW Fisher"){
    return(list(pvalue.meta = enrich.out$pvalue.0, pvalue.meta.B = enrich.out$pvalue.B,
                qvalue.meta = enrich.out$qvalue.0,AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = enrich.out$pvalue.0, pvalue.meta.B = enrich.out$pvalue.B,
                qvalue.meta = enrich.out$qvalue.0))
  }
}



#########################
MAPE_G_KS_DE = function (ind.p=ind.p,DB.matrix,gene.common,stat,rth.value,
                         size.min = 15, size.max = 500)
{
  if (stat == "maxP") {
    P.0 = as.matrix(apply(ind.p, 1, max))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,n,1)
  }
  else if (stat == "minP") {
    P.0 = as.matrix(apply(ind.p, 1, min))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,1,n)
  }
  else if (stat == "rth") {
    P.0 = as.matrix(apply(ind.p, 1, function(x) sort(x)[rth.value]))
    n = ncol(ind.p)
    pvalue.meta = pbeta(P.0,rth.value,(n - rth.value + 1))
  }
  else if (stat == "Fisher") {
    DF = 2 * ncol(ind.p)
    P.0 = as.matrix(apply(ind.p, 1, function(x) (-2 * sum(log(x)))))
    n = ncol(ind.p)
    pvalue.meta = pchisq(P.0,DF,lower.tail = F)
  }
  else if (stat == "AW Fisher") {
    AW.fisher = aw.fisher.pvalue(pvalue.0.mtx, method="original", weight.matrix=T)
    pvalue.meta = matrix(AW.fisher$pvalues,ncol = 1)
    weights.meta = AW.fisher$weights
    rownames(pvalue.meta) = rownames(ind.p)
  }
  else {
    stop("Please check: the selection of stat should be one of the following options: maxP,minP,rth and Fisher")
  }
  gene.name.sort = names(sort(pvalue.meta[,1],decreasing = F))
  gene.name.sort = gene.name.sort[gene.name.sort%in%gene.common]
  gene.name.sort = toupper(gene.name.sort)
  set2allgenes.mtx = DB.matrix
  order.mtx.1 = (set2allgenes.mtx[, gene.name.sort])
  order.mtx.0 = (1 - order.mtx.1)
  n_hit = rowSums(order.mtx.1)
  n_miss = rowSums(order.mtx.0)
  n_genes = ncol(order.mtx.1)
  nn = (n_hit*n_miss)/n_genes
  order.mtx.1 = t(apply(order.mtx.1, 1, function(x) x/sum(x)))
  order.mtx.0 = -t(apply(order.mtx.0, 1, function(x) x/sum(x)))
  order.mtx = order.mtx.0 + order.mtx.1
  ES.0 = as.matrix(apply(t(apply(order.mtx, 1, cumsum)), 1,
                         max))
  pvalue.0 = matrix(exp(-2*nn*(ES.0[,1]^2)), ncol = 1)
  rownames(pvalue.0) = rownames(ES.0)
  pvalue.set.0 = pvalue.0
  qvalue.set.0 = as.matrix(p.adjust(pvalue.0,"BH"),ncol = 1)
  rownames(qvalue.set.0) = rownames(pvalue.0)
  if(stat == "AW Fisher"){
    return(list(pvalue.meta = pvalue.set.0, qvalue.meta = qvalue.set.0, AW.weights = weights.meta))
  }else{
    return(list(pvalue.meta = pvalue.set.0, qvalue.meta = qvalue.set.0))
  }
}