--- a
+++ b/common_scripts/statistics/statistics_wrappers.R
@@ -0,0 +1,325 @@
+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)
+  }
+}