wrapper.GSEA=function(i, logicalVectors, data, cls.vector, GENESETS, datatype=NULL, clsname="", WD=getwd(), OUTDIR=getwd()){

  if(is.null(datatype)){stop("specify data type (B, N)")}

  # filter matrix and cls.vector

  data=data[,logicalVectors[[i]]&!is.na(logicalVectors)]

  cls.vector=cls.vector[logicalVectors[[i]]&!is.na(logicalVectors)]

  dataname=names(logicalVectors)[i]

  command=run.GSEA(data, cls.vector, GENESETS, datatype, dataname, clsname, WD, OUTDIR)

}

# function to go through

run.GSEA=function(data, cls.vector, GENESETS, datatype=NULL, dataname="", clsname="", WD=getwd(), OUTDIR=getwd(), GSEA_home="~/gsea-3.0.jar"){

  clsname=gsub("^_|_$", "", gsub("[[:punct:]]", "_", clsname)) # no special allowed here

  if(is.null(datatype)){stop("specify data type (B, N)")}

  # PATHS, create if needed

  if(!dir.exists(WD)){dir.create(WD, recursive = T); cat("Working directory made", sep="\n\n")}

  if(!dir.exists(OUTDIR)){dir.create(OUTDIR, recursive = T); cat("Output directory made", sep="\n\n")}

  setwd(WD)

  if(datatype=="N"){

    data=data[,order(as.numeric(cls.vector), decreasing=T)]

    cls.vector=cls.vector[order(as.numeric(cls.vector), decreasing=T)]

    CLUSTER_COMPARISON="continuous_phenotype"

  }

  if(datatype=="B"){

    data=data[,order(as.numeric(cls.vector), decreasing=T)]

    cls.vector=cls.vector[order(as.numeric(cls.vector), decreasing=T)]

    CLUSTER_COMPARISON=paste(unique(cls.vector), collapse="_versus_")

  }

  print("data sorted based on cls vector")

  #*********************** Need gct? ***************************

  NAME_OUT2=paste(dataname, gsub("-", "_", clsname), sep="_")

  GEXP=paste(WD, "/", NAME_OUT2, ".gct", sep="")

  print("making gct file")

  gsea.write.gct(data, GEXP)

  print("gct ready")

  #****************************** genesets ***********************************

  if(is.list(GENESETS)){

    writeGMT(GENESETS, "geneset_temp.gmt")

    print("gmt made from named list")

    GENESETS="geneset_temp.gmt"

  }

  # no need to do anything

  #*************************************************************************************

  if(datatype=="N"){

    print("this is numeric feature class, making class file")

    A1=paste("#numeric", sep = "\t")

    A2=paste("#feat", sep = "\t")

    A3=t(data.frame(as.character(cls.vector)))

    CLS=paste0(OUTDIR,"/",dataname, "_", CLUSTER_COMPARISON, "_", gsub("-", "_", clsname), ".cls")

    write(A1, file = CLS, append = F) # All .cls files require this dummy header line.

    write(A2, file = CLS, append = T)

    write.table(A3, file = CLS, append = T, quote = F, sep = "\t",

                na = "", row.names = F, col.names = F)

    PARAMS=" -collapse false -mode Max_probe -norm meandiv -nperm 1000 -permute sample -rnd_type no_balance -scoring_scheme weighted -metric Pearson -sort real -order descending -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 25 -rnd_seed timestamp -save_rnd_lists false -set_max 500 -set_min 5 -zip_report false -gui false"

  }

  if(datatype=="B"){

    print("this is binary/two group feature class")

    # make vector with labels

    cls.vector=cls.vector[!is.na(cls.vector), drop=F]

    A1=paste(length(cls.vector), 2, 1, sep = "\t")

    A2=paste("#", "1", "0", sep = "\t")

    A3=t(data.frame(as.character(cls.vector)))

    CLS=paste0(OUTDIR,"/", dataname, "_", CLUSTER_COMPARISON, "_", gsub("-", "_", clsname), "_temp.cls")

    write(A1, file = CLS, append = F) # All .cls files require this dummy header line.

    write(A2, file = CLS, append = T)

    write.table(A3, file = CLS, append = T, quote = F, sep = "\t",

                na = "", row.names = F, col.names = F)

    PARAMS=" -collapse false -mode Max_probe -norm meandiv -nperm 1000 -permute sample -rnd_type no_balance -scoring_scheme weighted -metric Signal2Noise -sort real -order descending -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 25 -rnd_seed timestamp -save_rnd_lists false -zip_report false -gui false"

  }

  RUN=paste0("java -Xmx5000m -cp ", GSEA_home, " xtools.gsea.Gsea -res ", GEXP,

             " -cls ", CLS, 

             " -gmx ", GENESETS,

             " -rpt_label ", NAME_OUT2, "_", CLUSTER_COMPARISON,

             " -out ", OUTDIR, " ", PARAMS)

}

gsea.write.gct <- function(exprMat, gctFn){

  nGenes = nrow(exprMat)

  nConds = ncol(exprMat)

  write("#1.2", file = gctFn, append = F) # All .gct files require this dummy header line.

  write(paste(nGenes, nConds, sep = "\t"), file = gctFn, append = T)

  write(paste("Name", "Description", paste(colnames(exprMat), collapse = "\t"), sep = "\t"), file = gctFn, append = T)

  # The second column of the .gct file, "Description", is filled out with "na"'s. 

  rownames(exprMat) = paste(rownames(exprMat), "na", sep = "\t") # Append "\tna" to every gene clsname. 

  write.table(exprMat, file = gctFn, append = T, quote = F, sep = "\t",

              na = "", row.names = T, col.names = F)

}

writeGMT <- function #Create a gmt (gene matrix transposed) file

### Createss a gmt (gene matrix transposed) file such as those

### provided by mSigDB or geneSigDB, from an R list object.

### Function by Levi Waldron.

(object,

 ### R list object that will be converted to GMT file.  Each element

 ### should contain a vector of gene names, and the names of the

 ### elements will used for the gene set names

 fname

 ### Output file name for .gmt file

){

  if (class(object) != "list") stop("object should be of class 'list'")

  if(file.exists(fname)) unlink(fname)

  for (iElement in 1:length(object)){

    write.table(t(c(make.names(rep(names(object)[iElement],2)),object[[iElement]])),

                sep="\t",quote=FALSE,

                file=fname,append=TRUE,col.names=FALSE,row.names=FALSE)

  }

  ### Called for the effect of writing a .gmt file

}

GSEA = function(gene_list, genesets, pval) {

  set.seed(54321)

  library(dplyr)

  library(gage)

  library(fgsea)

  if ( any( duplicated(names(gene_list)) )  ) {

    warning("Duplicates in gene names")

    gene_list = gene_list[!duplicated(names(gene_list))]

  }

  if  ( !all( order(gene_list, decreasing = TRUE) == 1:length(gene_list)) ){

    warning("Gene list not sorted")

    gene_list = sort(gene_list, decreasing = TRUE)

  }

  # if not list:

  if(is.character(genesets)){

    genesets = fgsea::gmtPathways(genesets)

  }

  fgRes <- fgsea::fgsea(pathways = genesets, 

                        stats = gene_list,

                        minSize=15,

                        maxSize=600,

                        nperm=10000) %>% 

    as.data.frame() %>% 

    dplyr::filter(padj < !!pval)

  #print(dim(fgRes))

  ## Filter FGSEA by using gage results. Must be significant and in same direction to keep 

  gaRes = gage::gage(gene_list, gsets=genesets, same.dir=TRUE, set.size =c(15,600))

  ups = as.data.frame(gaRes$greater) %>% 

    tibble::rownames_to_column("Pathway") %>% 

    dplyr::filter(!is.na(p.geomean) & q.val < pval ) %>%

    dplyr::select("Pathway")

  downs = as.data.frame(gaRes$less) %>% 

    tibble::rownames_to_column("Pathway") %>% 

    dplyr::filter(!is.na(p.geomean) & q.val < pval ) %>%

    dplyr::select("Pathway")

  #print(dim(rbind(ups,downs)))

  keepups = fgRes[fgRes$NES > 0 & !is.na(match(fgRes$pathway, ups$Pathway)), ]

  keepdowns = fgRes[fgRes$NES < 0 & !is.na(match(fgRes$pathway, downs$Pathway)), ]

  ### Collapse redundant pathways

  Up = fgsea::collapsePathways(keepups, pathways = genesets, stats = gene_list,  nperm = 500, pval.threshold = 0.05)

  Down = fgsea::collapsePathways(keepdowns, genesets, gene_list,  nperm = 500, pval.threshold = 0.05) 

  fgRes = fgRes[ !is.na(match(fgRes$pathway, 

                              c( Up$mainPathways, Down$mainPathways))), ] %>% 

    arrange(desc(NES))

  fgRes$pathway = stringr::str_replace(fgRes$pathway, "GO_" , "")

  fgRes$Enrichment = ifelse(fgRes$NES > 0, "Up-regulated", "Down-regulated")

  filtRes = rbind(head(fgRes, n = 10),

                  tail(fgRes, n = 10 ))

  g = ggplot(filtRes, aes(reorder(pathway, NES), NES)) +

    geom_segment( aes(reorder(pathway, NES), xend=pathway, y=0, yend=NES)) +

    geom_point( size=5, aes( fill = Enrichment),

                shape=21, stroke=2) +

    scale_fill_manual(values = c("Down-regulated" = "dodgerblue",

                                 "Up-regulated" = "firebrick") ) +

    coord_flip() +

    labs(x="Pathway", y="Normalized Enrichment Score",

         title="GSEA - Biological Processes") + 

    theme_minimal()

  output = list("Results" = fgRes, "Plot" = g)

  return(output)

}