wrapper.wilcoxtest=function(genelist, data, logicalVectors, logicalVector_normals=NULL, ALTERNATIVE="greater", adj.method="BH", CORES=1, prettynum=T){

  library(parallel)

  data=do.call(rbind, mclapply(genelist, TestGeneWilcox, data, logicalVectors, logicalVector_normals=logicalVector_normals, ALTERNATIVE, prettynum, mc.cores = CORES))

  data$adj.p=p.adjust(data$p, method = adj.method)

  data$adj.method=adj.method

  # signifcode

  data$signifCode=""

  data$signifCode[as.numeric(data$adj.p)<0.1]="."

  data$signifCode[as.numeric(data$adj.p)<0.05]="*"

  data$signifCode[as.numeric(data$adj.p)<0.01]="**"

  data$signifCode[as.numeric(data$adj.p)<0.001]="***"

  data$signifCode[as.numeric(data$adj.p)<0.0001]="****"

  if(prettynum){

    data$adj.p=prettyNum(signif(data$adj.p,2))

  }

  if(adj.method=="BH"){

    colnames(data)[colnames(data)%in%"adj.p"]="FDR"

  }

  data=data[,c(1:5, 7,9, 6,8)]

  return(data)

}

FUN_GOTHROUGH=function(i, genelist, list_cancers, logicalVector_normals, TEST_FAILS=F, HG_PVAL=0.001, MW_PVAL=0.001, ALTERNATIVE="greater"){

  name=names(list_cancers)[i]

  logicalVector=as.logical(unlist(list_cancers[[i]]))

  logicalVector_normals=logicalVector_normals[!names(logicalVector_normals)%in%name]

  result=FUN_HGT_MW_TEST(genelist, logicalVector,logicalVector_normals, name, TEST_FAILS=F, data = data, profile = profile,HG_PVAL = HG_PVAL, MW_PVAL=MW_PVAL, ALTERNATIVE = ALTERNATIVE, CORES = 10)

}

fisher.wrapper=function(lv.list1, lv.list2, alternative="two.sided", method="bonferroni",log10=F, prettyNumbers=F,orderdata=F, cores=1){

  n1=names(lv.list1)

  n2=names(lv.list2)

  res=do.call(rbind,mclapply(seq(lv.list1), function(i){

    do.call(rbind,lapply(seq(lv.list2), function(j){

      data.frame("featureA"=n1[i], "featureB"=n2[j],fisher.2x2(lv.list1[[i]], lv.list2[[j]]))

    }))

  }, mc.cores=cores))

  res$adj.p=p.adjust(res$p, method=method)

  res$signifCode=""

  res$signifCode[as.numeric(res$adj.p)<0.1]="."

  res$signifCode[as.numeric(res$adj.p)<0.05]="*"

  res$signifCode[as.numeric(res$adj.p)<0.01]="**"

  res$signifCode[as.numeric(res$adj.p)<0.001]="***"

  res$signifCode[as.numeric(res$adj.p)<0.0001]="****"

  if(log10){

    res$p=abs(log10(res$p))

    res$adj.p=abs(log10(res$p))

  }

  if(prettyNumbers){

    res$log.odds=prettyNum(signif(res$log.odds,2))

    res$p=prettyNum(signif(res$p,2))

    res$adj.p=prettyNum(signif(res$adj.p,2))

  }

  if(orderdata)res=res[order(res$featureA, as.numeric(res$adj.p), decreasing = F),]

  return(res)

}

fisher.matrix.pairwise=function(E, plot=F, alternative="two.sided"){

  # perform a fisher exact test for each pairs of gene and saves the oddsRatio and p-value

  res <- NULL

  for(i in seq(dim(E)[2])){

    for(j in seq(dim(E)[2])){

      if(i!=j){

        f <- fisher.test(E[,i], E[,j], alternative = alternative) 

        res <- rbind(res,cbind(geneA=colnames(E)[i],geneB=colnames(E)[j],oddsRatio=f$estimate,pvalue=f$p.value))

      }   

    }

  }

  # some formatting

  res <- as.data.frame(res)

  res$geneA <- factor(res$geneA,levels=unique(res$geneA))

  res$geneB <- factor(res$geneB,levels=unique(res$geneB))

  res$oddsRatio <- as.numeric(as.character(res$oddsRatio))

  # avoid inf values

  res$oddsRatio[res$oddsRatio>10]=10

  res$oddsRatio[res$oddsRatio<=0]=0.01

  res$pvalue <- as.numeric(as.character(res$pvalue))

  # use p.adjust to correct for multi testing using a FDR

  res <- cbind(res,fdr=p.adjust(res$pvalue,"fdr"))

  # change the FDR in labels for plotting

  res$stars <- cut(res$fdr, breaks=c(-Inf, 0.001, 0.01, 0.05, Inf), label=c("***", "**", "*", ""))

  if(plot){

    # plot with ggplot 2

    require(ggplot2)

    require(cowplot) # not necessary but the plot is nicer

    p <- ggplot(res, aes(geneA, geneB)) + geom_tile(aes(fill = log2(oddsRatio+0.01)),colour = "white") + scale_fill_gradient2(midpoint=1) + geom_text(aes(label=stars), color="black", size=5)

    p

    return(p)

  }

  return(res)

}

fun.get.cox=function(i, logicalv, DATA, univariate=T, pretty=T, TIME, STATUS){

  logicalvector=logicalv[[i]]

  n=names(logicalv)[i]

  # remove if all are NA

  DATA=DATA[,!apply(DATA[logicalvector,,drop=F], 2, function(v)all(is.na(v)))]

  if(univariate){

    univ=do.call(rbind, lapply(seq(dim(DATA)[2]), function(j){

      y=Surv(TIME[logicalvector], STATUS[logicalvector])

      fit=coxph(formula = as.formula(paste("y ~ ", paste(colnames(DATA)[j], collapse="+"))), DATA[logicalvector,,drop=F])

      b=cox.zph(fit)[[1]][3] # p.value

      a=summary(fit)

      pval=a$coefficients[1,5]

      coef=a$conf.int[c(1,3,4)]

      v=data.frame(colnames(DATA)[j], t(coef), pval,"","", a$concordance[1], b<0.05, stringsAsFactors = F)

      rownames(v)=NULL

      return(v)

    }))

    # adjust P-value depending on number of genes tested:

    univ[,6]=p.adjust(univ[,5], method="BH")

  }else{

    y=Surv(TIME[logicalvector], STATUS[logicalvector])

    fit=coxph(formula = as.formula(paste("y ~ ", paste(colnames(DATA), collapse="+"))), DATA[logicalvector,])

    b=cox.zph(fit)$table[,3] # p.value

    b=b[!names(b)%in%"GLOBAL"]

    a=summary(fit)

    pval=a$coefficients[,5]

    coef=a$conf.int[,c(1,3,4)]

    univ=data.frame(rownames(a$coefficients), coef, pval, "", "", a$concordance[1], b<0.05, stringsAsFactors = F)

    rownames(univ)=NULL

    # no need to adjust P-value here:

    univ[,6]=univ[,5]

  }

  univ=univ[order(univ[,5]),]

  univ[,7][as.numeric(univ[,6 ])<0.1]="."

  univ[,7][as.numeric(univ[,6 ])<0.05]="*"

  univ[,7][as.numeric(univ[,6 ])<0.01]="**"

  univ[,7][as.numeric(univ[,6 ])<0.001]="***"

  univ[,7][as.numeric(univ[,6 ])<0.0001]="****"

  univ[,7][is.na(univ[,7])]=""

  val=data.frame(univ, n, stringsAsFactors = F)

  colnames(val)=c("Feature", "exp(coef)", "lower .95", "upper .95", "P", "Adj.P", "Signif", "concordance", "zph P < 0.05", "Name")

  if(pretty){

    for(i in c(2:6,8)){

      val[,i]=prettyNum(signif(val[,i],2))

    }

  }

  return(val)

}

fun.cox.elasticnet=function(DATA_ORG, time, status, cores=8, nfold=10,min.elnet=0.01,max.elnet=0.99, summary.km="85th_15th_percentile", percentage=0.25, REPEATS=100){

  library(glmnet)

  # run regularized regression

  a <- seq(min.elnet, max.elnet, 0.05)

  if(percentage==0){

    y=Surv(time, status)

    formdf1 <- as.formula(paste(" ~ ", paste(colnames(DATA_ORG),collapse="+")))

    x=model.matrix(formdf1, data=DATA_ORG)

    # go through alpha sequence with nfolds and repeat (100 times etc) cross validation on different sets of samples

    enet.rep=do.call(rbind, parallel::mclapply(seq(REPEATS), function(r){

      s <- do.call(rbind, lapply(a, function(i){

        cv <- cv.glmnet(x=x,y=y, family = "cox", nfold = nfold, type.measure = "deviance", alpha = i, standardize=F)

        data.frame(cvm = cv$cvm[cv$lambda == cv$lambda.min], lambda.min = cv$lambda.min, alpha = i)

      }))

      return(s)

    }, mc.cores=cores))

    # find mean alpha and lambda using all repeats

    search_m=do.call(rbind, lapply(a, function(alpha){

      cvm=mean(enet.rep$cvm[enet.rep$alpha==alpha])

      lambda=mean(enet.rep$lambda.min[enet.rep$alpha==alpha])

      data.frame(alpha, cvm, lambda)

    }))

    # minimum cvm error:

    cv3 <- search_m[search_m$cvm == min(search_m$cvm), ]

    print(cv3$alpha)

    print(cv3$lambda)

    print(cv3$cvm)

    # fit model with tuned alpha, use lambda sequence here, not single value

    md3 <- glmnet(x=x,y=y, family = "cox", alpha = cv3$alpha, standardize=F)

    coefficients <- coef(md3, s = cv3$lambda)[-1,,drop=F]

    coefficients=coefficients[order(abs(coefficients[,1]), decreasing = T),,drop=F]

    active_coefficients <- coefficients[,1] != 0

    coefficients_all=coefficients[active_coefficients,,drop=F]

    comprisk=DATA_ORG[,match(rownames(coefficients_all), colnames(DATA_ORG))]

    risk_patient=as.numeric(as.numeric(coefficients_all) %*% data.matrix(t(comprisk)))

    df=data.frame(x=risk_patient)

    ggplot(df, aes(x=x))+

      geom_density(color="lightblue", fill="lightblue")

    fun.kapplanMeier(time, status, CONTINUOUS=risk_patient, CONTINUOUS_SUMMARY = summary.km, MONTHS=F, PVAL=1, INDIVIDUAL_GROUPS=F, NAME = "predicted risk")

    a=list(summary(coxph(y ~ PI.train, data.frame("PI.train"=risk_patient))))

    c=list(coefficients_all)

    out=c(a,c)

    names(out)=c("PI.test", "coefficients")

    return(out)

  }else{

    print("Splitting to test and train data")

    # divide data so that there are equal proportions of outcomes with training and test set!

    nr.samples=floor(percentage*dim(DATA_ORG)[1])

    nr.events=floor(table(status)*percentage)

    set.seed(0)

    x <- which(status==1)

    vf <- logical(length(status==1))

    vf[x[sample.int(length(x), nr.events[2])]] <- TRUE

    set.seed(0)

    x <- which(status==0)

    vf2 <- logical(length(status==0))

    vf2[x[sample.int(length(x), nr.events[1])]] <- TRUE

    find=vf|vf2

    # define training and testing datasets

    TRAIN=DATA_ORG[!find,]

    TEST=DATA_ORG[find,]

    ytrain=Surv(time[!find], status[!find])

    ytest=Surv(time[find], status[find])

    # this must look the same:

    log=fun.kapplanMeier(time, status, GROUPS=(rownames(DATA_ORG)%in%rownames(TRAIN)*1), MONTHS=F, PVAL=1, INDIVIDUAL_GROUPS=T, NAME = "These must be similar\nand high pval: 1 TRAIN, 0 TEST")

    formdf1 <- as.formula(paste(" ~ ", paste(colnames(TRAIN),collapse="+")))

    x=model.matrix(formdf1, data=TRAIN)

    # go through alpha sequence with nfolds and repeat (100 times etc) cross validation on different sets of samples

    enet.rep=do.call(rbind, parallel::mclapply(seq(REPEATS), function(r){

      s <- do.call(rbind, lapply(a, function(i){

        cv <- cv.glmnet(x=x,y=ytrain, family = "cox", nfold = floor((dim(DATA_ORG)[1]/10)), type.measure = "deviance", alpha = i, standardize=F)

        data.frame(cvm = cv$cvm[cv$lambda == cv$lambda.min], lambda.min = cv$lambda.min, alpha = i)

      }))

      return(s)

    }, mc.cores=cores))

    # find mean alpha and lambda using all repeats

    search_m=do.call(rbind, lapply(a, function(alpha){

      cvm=mean(enet.rep$cvm[enet.rep$alpha==alpha])

      lambda=mean(enet.rep$lambda.min[enet.rep$alpha==alpha])

      data.frame(alpha, cvm, lambda)

    }))

    # minimum cvm error:

    cv3 <- search_m[search_m$cvm == min(search_m$cvm), ]

    print(cv3$alpha)

    print(cv3$lambda)

    print(cv3$cvm)

    # fit model with tuned alpha, use lambda sequence here, not single value

    md3 <- glmnet(x=x,y=ytrain, family = "cox", alpha = cv3$alpha, standardize=F)

    coefficients <- coef(md3, s = cv3$lambda)[-1,,drop=F]

    coefficients=coefficients[order(abs(coefficients[,1]), decreasing = T),,drop=F]

    active_coefficients <- coefficients[,1] != 0

    coefficients_all=coefficients[active_coefficients,,drop=F]

    comprisk=TRAIN[,match(rownames(coefficients_all), colnames(TRAIN))]

    risk_patient=as.numeric(as.numeric(coefficients_all) %*% data.matrix(t(comprisk)))

    comprisk=TEST[,match(rownames(coefficients_all), colnames(TEST))]

    risk_patient_test=as.numeric(as.numeric(coefficients_all) %*% data.matrix(t(comprisk)))

    df=data.frame(x=risk_patient)

    ggplot(df, aes(x=x))+

      geom_density(color="lightblue", fill="lightblue")

    df=data.frame(x=risk_patient_test)

    ggplot(df, aes(x=x))+

      geom_density(color="indianred", fill="indianred")

    fun.kapplanMeier(time[!find], status[!find], CONTINUOUS=risk_patient, CONTINUOUS_SUMMARY = summary.km, MONTHS=F, PVAL=1, INDIVIDUAL_GROUPS=F, NAME = "predicted risk")

    fun.kapplanMeier(time[find], status[find], CONTINUOUS=risk_patient_test, CONTINUOUS_SUMMARY = summary.km, MONTHS=F, PVAL=1, INDIVIDUAL_GROUPS=F, NAME = "validated risk")

    a=list(summary(coxph(ytrain ~ PI.train, data.frame("PI.train"=risk_patient))))

    b=list(summary(coxph(ytest ~ PI.test, data.frame("PI.test"=risk_patient_test))))

    c=list(coefficients_all)

    out=c(a,b,c)

    names(out)=c("PI.train", "PI.test", "coefficients")

    return(out)

  }

}