Introduction

In this report we establish a new method for generating simulated data with known ground truth. This will be used to test different gene methylation enrichment approaches systematically.

The general steps are:

Import GSE158422 data corresponding to control (non-tumour tissue).
From the 37 samples, create two groups of 18 samples. One of these will be considered “control” and the other “case”.
Create random gene sets that have similar sized to Reactome pathways.
Some gene sets will be selected to be differentially methylated. Half of these will be hypermethylated and the others will be hypomethylated.
The changes will be incorporated into the “case” samples.
The enrichment analysis will be conducted.
The accuracy will be calculated.

suppressPackageStartupMessages({
  library("stringi")
  library("limma")
  library("missMethyl")
  library("IlluminaHumanMethylation450kmanifest")
  library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
  library("HGNChelper")
  library('org.Hs.eg.db')
  library("psych")
  library("mitch")
  library("kableExtra")
})

# optimised for 128 GB sever with 32 threads
CORES=6

Load data

annotations
probe sets
gene sets
design matrix
mval matrix

# get probe-gene mapping
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","UCSC_RefGene_Group","Islands_Name","Relation_to_Island")])
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
gp2 <- lapply(gp2,unique)
gt <- stack(gp2)
colnames(gt) <- c("gene","probe")
dim(gt)

## [1] 684970      2

str(gt)

## 'data.frame':    684970 obs. of  2 variables:
##  $ gene : chr  "YTHDF1" "EIF2S3" "PKN3" "CCDC57" ...
##  $ probe: Factor w/ 865859 levels "cg18478105","cg09835024",..: 1 2 3 4 5 6 7 8 8 9 ...

head(gt)

##     gene      probe
## 1 YTHDF1 cg18478105
## 2 EIF2S3 cg09835024
## 3   PKN3 cg14361672
## 4 CCDC57 cg01763666
## 5   INF2 cg12950382
## 6  CDC16 cg02115394

#new.hgnc.table <- getCurrentHumanMap()
#new.hgnc.table <- readRDS("new.hgnc.table.rds")
#fix <- checkGeneSymbols(gt$gene,map=new.hgnc.table)
#fix2 <- fix[which(fix$x != fix$Suggested.Symbol),]
#length(unique(fix2$x))
#gt$gene <- fix$Suggested.Symbol

if (!file.exists("GSE158422_design.rds")) {
  download.file("https://ziemann-lab.net/public/gmea_prototype/GSE158422_design.rds", "GSE158422_design.rds")
}
design <- readRDS("GSE158422_design.rds")

if (!file.exists("GSE158422_design.rds")) {
  options(timeout=10000000000000)
  download.file("https://ziemann-lab.net/public/gmea_prototype/GSE158422_mx.rds","GSE158422_mx.rds")
}
mval <- readRDS("GSE158422_mx.rds")

Gene set database

We could use Reactome pathways, however these have a lot of overlapping sets, which could cause inflated false positives. After some testing, I found that the gene set size strongly impacts the accuracy.

For reference, Reactome sets have a mean of 48 and median of 15, while MSigDB has a median of 47 and mean of 116.

Therefore, a reasonable analysis would include small (20), medium (50) and large (100) sets.

To reflect coverage of Reactome and other pathway databases, only half the genes will be included in sets.

gene_catalog <- unique(gt$gene)

# set gene set size here
lengths <- rep(20,1000)

randomGeneSets <- function(gene_catalog, lengths, seed){
  set.seed(seed) ; gene_catalog_half <- sample(gene_catalog,length(gene_catalog)/2)
  num_gsets <- length(lengths)
  set.seed(seed) ; seeds <- sample(1:1e6, num_gsets)
  gsets <- lapply(1:num_gsets,function(i) {
    set.seed(seeds[i]) ; gs <- sample(gene_catalog_half,lengths[i])
    return(gs)
  } )
  names(gsets)<-stri_rand_strings(length(gsets), 15, pattern = "[A-Za-z]")
  return(gsets)
}

gsets <- randomGeneSets(gene_catalog,lengths,seed=100)

Incorporate changes

TODO: to incorporate changes to case samples.

Need to figure out what magnitude to change. Will refer to the cancer/normal comparison.

Select genes and probes to alter.

seed=100
frac_genes=0.5
frac_probes=0.5
delta=1
nsamples=10
normal_mval <- mval[,(1:(ncol(mval)/2)*2)]

incorp_dm <- function(genesets,myann,mval,seed, frac_genes,frac_probes,
  groupsize,delta=1,gene_catalog) {

  # divide gene sets between hyper and hypomethylated
  nset <- floor(length(genesets)/2)
  set.seed(seed) ; gtup <-sample(genesets,nset)
  set.seed(seed) ; gtdn <- sample(setdiff(genesets,gtup),nset)
  gup <- unname(unlist(gtup))
  gdn <- unname(unlist(gtdn))
  # add extra 10% DM genes
  gnon <- setdiff(gene_catalog,c(gup,gdn))
  gextra <- round(length(gnon)*0.1)
  set.seed(seed) ; gup <- c(gup,sample(gnon,gextra))
  gnon <- setdiff(gene_catalog,c(gup,gdn))
  set.seed(seed) ; gdn <- c(gdn,sample(gnon,gextra))
  # make probe-gene vector
  probe2gene <- strsplit(myann$UCSC_RefGene_Name,";")
  names(probe2gene) <- rownames(myann)
  probe2gene <- unlist(probe2gene)
  # select probes hypermethylated
  set.seed(seed) ; gup2 <- sample(gup,floor(length(gup)*frac_genes))
  pup <- names(probe2gene[which(probe2gene %in% gup2)])
  set.seed(seed) ; pup2 <- sample(pup,floor(length(pup)*frac_probes))
  # select probes hypomethylated
  set.seed(seed) ; gdn2 <- sample(gdn,floor(length(gdn)*frac_genes))
  pdn <- names(probe2gene[which(probe2gene %in% gdn2)])
  set.seed(seed) ; pdn2 <- sample(pdn,floor(length(pdn)*frac_probes))
  # add 10% DM probes as well
  probes <- rownames(myann)
  pnon <- setdiff(probes,c(pup,pdn))
  pextra <- round(length(pnon)*0.1)
  set.seed(seed) ; pup <- c(pup,sample(pnon,pextra))
  pnon <- setdiff(probes,c(pup,pdn))
  set.seed(seed) ; pdn <- c(pdn,sample(pnon,pextra))
  # divide samples between ctrl and case
  ncols <- ncol(mval)
  maxgroupsize=floor(ncols/2)
  if ( groupsize > maxgroupsize ) { stop("groupsize cannot be larger than half the ncols of mval") }
  set.seed(seed) ; ctrl <- sample(1:ncols,groupsize)
  set.seed(seed) ; case <- sample(setdiff(1:ncols,ctrl),groupsize)
  mval_ctrl <- mval[,ctrl]
  mval_case <- mval[,case]
  # incorporate altered signals - change by +1 or -1
  mval_case[rownames(mval_case) %in% pup2,] <-  mval_case[rownames(mval_case) %in% pup2,] + delta
  mval_case[rownames(mval_case) %in% pdn2,] <-  mval_case[rownames(mval_case) %in% pdn2,] - delta
  mval2 <- cbind(mval_ctrl,mval_case)
  result <- list("mval"=mval2,"probes up"=pup2,"probes down"=pdn2,
    "genes up"=gup2,"genes down"=gdn2,
    "genesets up"=gtup,"genesets down"=gtdn)
  return(result)
}

GSAMETH function

# limma
runlimma <- function(mval,design,myann) {
  fit.reduced <- lmFit(mval,design)
  fit.reduced <- eBayes(fit.reduced)
  dm <- topTable(fit.reduced,coef=ncol(design), number = Inf)
  dm <- merge(myann,dm,by=0)
  dm <- dm[order(dm$P.Value),]
  rownames(dm) <- dm$Row.names
  dm$Row.names=NULL
  return(dm)
}

This is how to use the function

This could be complicated as it requires translation of symbols to entrez IDs, but could be simplified if the entrez translation is done at the gsameth step.

simgsa <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(unlist(strsplit(myann$UCSC_RefGene_Name,";")))
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  dm3 <- runlimma(mval=mval2,design=d,myann=myann)
  pup3 <- rownames(subset(dm3,adj.P.Val<0.05 & logFC>0))
  pdn3 <- rownames(subset(dm3,adj.P.Val<0.05 & logFC<0))
  if ( length(pup3) < 250 ) { pup3 <- head(rownames(subset(dm3, logFC > 0)), 250) }
  if ( length(pdn3) < 250 ) { pdn3 <- head(rownames(subset(dm3, logFC < 0)), 250) }
  # convert gene sets to entrez
  suppressWarnings(suppressMessages({ gene2entrez <- mapIds(org.Hs.eg.db, gene_catalog, 'ENTREZID', 'SYMBOL') }))
  gsets_entrez <- lapply(gsets,function(gs) {
    gs2 <- unique(gene2entrez[names(gene2entrez) %in% gs])
    gs2 <- gs2[!is.na(gs2)]
    return(gs2)
  })
  suppressWarnings(suppressMessages({
    gsaup3 <- gsameth(sig.cpg=pup3, all.cpg=rownames(dm3), collection=gsets_entrez, array.type="EPIC")
    gsadn3 <- gsameth(sig.cpg=pdn3, all.cpg=rownames(dm3), collection=gsets_entrez, array.type="EPIC")
  }))
  gsig_up3 <- rownames(subset(gsaup3,FDR<0.05))
  gsig_dn3 <- rownames(subset(gsadn3,FDR<0.05))
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(gsadn3)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

LA function

This process runs limma first and then aggregates the results before doing an enrichment test.

# enrich parametric
ttenrich <- function(m,genesets,cores=1,testtype="selfcontained") {
  res <- mclapply( 1:length(genesets), function(i) {
    scores <- m[,1]
    gs <- genesets[i]
    name <- names(gs)
    n_members <- length(which(rownames(m) %in% gs[[1]]))
    if ( n_members > 4 ) {
      tstats <- m[which(rownames(m) %in% gs[[1]]),]
      myn <- length(tstats)
      mymean <- mean(tstats)
      mymedian <- median(tstats)
      if ( testtype == "selfcontained" ) { wt <- t.test(tstats) }
      if ( testtype == "competitive" ) { wt <- t.test(tstats,scores) }
      res <- c(name,myn,mymean,mymedian,wt$p.value)
    }
  } , mc.cores = cores)
  res_df <- do.call(rbind, res)
  rownames(res_df) <- res_df[,1]
  res_df <- res_df[,-1]
  colnames(res_df) <- c("n_genes","t_mean","t_median","pval")
  tmp <- apply(res_df,2,as.numeric)
  rownames(tmp) <- rownames(res_df)
  res_df <- tmp
  res_df <- as.data.frame(res_df)
  res_df <- res_df[order(res_df$pval),]
  res_df$logp <- -log10(res_df$pval )
  res_df$fdr <- p.adjust(res_df$pval,method="fdr")
  res_df[order(abs(res_df$pval)),]
  return(res_df)
}

# enrich non-parametric
wtenrich <- function(m,genesets,cores=1,testtype="selfcontained") {
  res <- mclapply( 1:length(genesets), function(i) {
    scores <- m[,1]
    gs <- genesets[i]
    name <- names(gs)
    n_members <- length(which(rownames(m) %in% gs[[1]]))
    if ( n_members > 4 ) {
      tstats <- m[which(rownames(m) %in% gs[[1]]),]
      myn <- length(tstats)
      mymean <- mean(tstats)
      mymedian <- median(tstats)
      if ( testtype == "selfcontained" ) { wt <- wilcox.test(tstats) }
      if ( testtype == "competitive" ) { wt <- wilcox.test(tstats,scores) }
      res <- c(name,myn,mymean,mymedian,wt$p.value)
    }
  } , mc.cores = cores)
  res_df <- do.call(rbind, res)
  rownames(res_df) <- res_df[,1]
  res_df <- res_df[,-1]
  colnames(res_df) <- c("n_genes","t_mean","t_median","pval")
  tmp <- apply(res_df,2,as.numeric)
  rownames(tmp) <- rownames(res_df)
  res_df <- tmp
  res_df <- as.data.frame(res_df)
  res_df <- res_df[order(res_df$pval),]
  res_df$logp <- -log10(res_df$pval )
  res_df$fdr <- p.adjust(res_df$pval,method="fdr")
  res_df[order(abs(res_df$pval)),]
  return(res_df)
}

# LA parametric competitive
simlac <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  dm3 <- runlimma(mval=mval2,design=d,myann=myann)
  dd <- merge(dm3,gt,by.x=0,by.y="probe")
  m1 <- aggregate(t ~ gene,dd,mean)
  rownames(m1) <- m1$gene
  m1$gene=NULL
  lares1 <- ttenrich(m=m1,genesets=gsets,cores=2,testtype="competitive")
  gsig_up3 <- rownames(subset(lares1, fdr < 0.05 & t_mean > 0))
  gsig_dn3 <- rownames(subset(lares1, fdr < 0.05 & t_mean < 0))
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(lares1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

# LA parametric competitive top
simlactop <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  dm3 <- runlimma(mval=mval2,design=d,myann=myann)
  dd <- merge(dm3,gt,by.x=0,by.y="probe")
  m1 <- aggregate(t ~ gene,dd, function(x) {
    if (abs(max(x)) > abs(min(x))) { max(x) } else { min(x) }
  })
  rownames(m1) <- m1$gene
  m1$gene=NULL
  lares1 <- ttenrich(m=m1,genesets=gsets,cores=2,testtype="competitive")
  gsig_up3 <- rownames(subset(lares1, fdr < 0.05 & t_mean > 0))
  gsig_dn3 <- rownames(subset(lares1, fdr < 0.05 & t_mean < 0))
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(lares1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

# LA nonparametric competitive
simnlac <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  dm3 <- runlimma(mval=mval2,design=d,myann=myann)
  dd <- merge(dm3,gt,by.x=0,by.y="probe")
  m1 <- aggregate(t ~ gene,dd,mean)
  rownames(m1) <- m1$gene
  m1$gene=NULL
  lares1 <- wtenrich(m=m1,genesets=gsets,cores=2,testtype="competitive")
  gsig_up3 <- rownames(subset(lares1, fdr < 0.05 & t_mean > 0))
  gsig_dn3 <- rownames(subset(lares1, fdr < 0.05 & t_mean < 0))
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(lares1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

# LA competitive mitch
simlacm <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  dm3 <- runlimma(mval=mval2,design=d,myann=myann)
  dd <- merge(dm3,gt,by.x=0,by.y="probe")
  m1 <- aggregate(t ~ gene,dd,mean)
  rownames(m1) <- m1$gene
  m1$gene=NULL
  lamres1 <- runmitch(m=m1,genesets=gsets,cores=2)
  gsig_up3 <- rownames( subset( lamres1, p.adjustANOVA < 0.05 & s.dist > 0 ) )
  gsig_dn3 <- rownames( subset( lamres1, p.adjustANOVA < 0.05 & s.dist < 0 ) )
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(lamres1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

# LA rank competitive mitch
simlrm <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  dm3 <- runlimma(mval=mval2,design=d,myann=myann)
  # rank probes first
  dm3$rank <-  rank(dm3$t) - nrow(subset(dm3,t<0))
  mm <- merge(dm3,gt,by.x=0,by.y="probe")
  head(mm)
  mma <-aggregate(mm$rank ~ gene,mm,mean)
  rownames(mma) <- mma$gene
  mma$gene = NULL
  colnames(mma) <- "meanrank"
  lrmres <- runmitch(m=mma,genesets=gsets,cores=2)
  gsig_up3 <- rownames( subset( lrmres, p.adjustANOVA < 0.05 & s.dist > 0 ) )
  gsig_dn3 <- rownames( subset( lrmres, p.adjustANOVA < 0.05 & s.dist < 0 ) )
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(lrmres)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

AL: Aggregate limma

Functions for aggregate-limma-enrich approach.

# AL approach parametric competitive
simalc <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  # al pipeline
  mm <- merge(mval2,gt,by.x=0,by.y="probe")
  mm$Row.names = NULL
  a <- aggregate(. ~ gene,mm,mean)
  rownames(a) <- a$gene
  a$gene=NULL
  fit.reduced <- lmFit(a,d)
  fit.reduced <- eBayes(fit.reduced)
  al <- topTable(fit.reduced,coef=ncol(d), number = Inf)
  m1 <- as.data.frame(al$t)
  rownames(m1) <- rownames(al)
  colnames(m1) <- "t"
  alres1 <- ttenrich(m=m1,genesets=gsets,cores=2,testtype="competitive")
  # summarise results
  gsig_up3 <- rownames(subset(alres1, fdr < 0.05 & t_mean > 0))
  gsig_dn3 <- rownames(subset(alres1, fdr < 0.05 & t_mean < 0))
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(alres1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

# AL approach nonparametric competitive
simnalc <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  # al pipeline
  mm <- merge(mval2,gt,by.x=0,by.y="probe")
  mm$Row.names = NULL
  a <- aggregate(. ~ gene,mm,mean)
  rownames(a) <- a$gene
  a$gene=NULL
  fit.reduced <- lmFit(a,d)
  fit.reduced <- eBayes(fit.reduced)
  al <- topTable(fit.reduced,coef=ncol(d), number = Inf)
  m1 <- as.data.frame(al$t)
  rownames(m1) <- rownames(al)
  colnames(m1) <- "t"
  alres1 <- wtenrich(m=m1,genesets=gsets,cores=2,testtype="competitive")
  # summarise results
  gsig_up3 <- rownames(subset(alres1, fdr < 0.05 & t_mean > 0))
  gsig_dn3 <- rownames(subset(alres1, fdr < 0.05 & t_mean < 0))
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(alres1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

Agg-limma-mitch function

This approach uses the aggregated mvals, limma and instead of a 1-sample t-test it uses mitch which is a competitive test and could give more interpretable results.

runmitch <- function(m,genesets,cores=1) {
  suppressMessages({ mres <- mitch_calc(m,genesets,minsetsize=5,cores=cores) })
  mres <- mres$enrichment_result
  rownames(mres) <- mres$set
  mres$set=NULL
  return(mres)
}

simalm <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  # alm pipeline
  mm <- merge(mval2,gt,by.x=0,by.y="probe")
  mm$Row.names = NULL
  a <- aggregate(. ~ gene,mm,mean)
  rownames(a) <- a$gene
  a$gene=NULL
  fit.reduced <- lmFit(a,d)
  fit.reduced <- eBayes(fit.reduced)
  al <- topTable(fit.reduced,coef=ncol(d), number = Inf)
  m1 <- as.data.frame(al$t)
  rownames(m1) <- rownames(al)
  colnames(m1) <- "t"
  almres1 <- runmitch(m=m1,genesets=gsets,cores=2)
  # summarise results
  gsig_up3 <- rownames( subset( almres1, p.adjustANOVA < 0.05 & s.dist > 0 ) )
  gsig_dn3 <- rownames( subset( almres1, p.adjustANOVA < 0.05 & s.dist < 0 ) )
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(almres1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

AA Aggregate-aggregate-limma functions

Use mean value works well here.

gsagg <- function(x,genesets,cores=1) {
  meds <- mclapply(1:length(genesets), function(i) {
    gs = genesets[[i]]
    xx <- x[rownames(x) %in% gs,]
    med <- apply(xx,2,mean)
  },mc.cores=cores)
  mymed <- do.call(rbind,meds)
  rownames(mymed) <- names(genesets)
  as.data.frame(mymed)
}

aalimma <- function(agag,design) {
  fit.reduced <- lmFit(agag,design)
  fit.reduced <- eBayes(fit.reduced)
  dmagg <- topTable(fit.reduced,coef=ncol(design), number = Inf)
  return(dmagg)
}

aal <- function(mval,myann,genesets,design,cores=1) {
  medf <- chragg(mval,myann,cores=cores)
  agag <- gsagg(x=medf,genesets=genesets,cores=cores)
  aalres <- aalimma(agag=agag,design=design)
  return(aalres)
}

simaa <- function(genesetdatabase, myann, mval, seed, frac_genes, frac_probes, groupsize, delta=1, num_dm_sets=50) {
  # generate gene sets
  gene_catalog <- unique(gt$gene)
  lengths <- unname(unlist(lapply(genesetdatabase,length)))
  gsets <- randomGeneSets(gene_catalog,lengths,seed=seed)
  # select gene sets to alter
  set.seed(seed) ; gset_mod <- sample(gsets,num_dm_sets)
  # incorporate select changes
  sim <- incorp_dm(genesets=gset_mod, myann=myann, mval=mval, seed=seed,
    frac_genes=0.5,frac_probes=0.5,groupsize=groupsize,delta=delta,
    gene_catalog=gene_catalog)
  # set up limma
  mval2 <- sim$mval
  ncols <- ncol(mval2)
  groupsize <- ncols/2
  ss <- data.frame(colnames(mval2))
  colnames(ss) <- "sample"
  ss$case <- c(rep(0,groupsize),rep(1,groupsize))
  d <- model.matrix(~ ss$case )
  mm <- merge(mval2,gt,by.x=0,by.y="probe")
  mm$Row.names = NULL
  a <- aggregate(. ~ gene,mm,mean)
  rownames(a) <- a$gene
  a$gene=NULL
  mystack <- stack(gsets)
  mmm <- merge(a,mystack,by.x=0,by.y="values")
  mmm$Row.names=NULL
  aa <- aggregate(. ~ ind,mmm,mean)
  rownames(aa) <- aa$ind
  aa$ind=NULL
  fit.reduced <- lmFit(aa,d)
  fit.reduced <- eBayes(fit.reduced)
  aares1 <- topTable(fit.reduced,coef=ncol(d), number = Inf)
  # summarise results
  gsig_up3 <- rownames(subset(aares1, adj.P.Val < 0.05 & logFC > 0))
  gsig_dn3 <- rownames(subset(aares1, adj.P.Val < 0.05 & logFC < 0))
  gtup <- names(sim[[6]])
  gtdn <- names(sim[[7]])
  UPTP=length(intersect(gsig_up3 ,gtup))
  UPFP=length(setdiff(gsig_up3 ,gtup))
  UPFN=length(setdiff(gtup,gsig_up3))
  DNTP=length(intersect(gsig_dn3 ,gtdn))
  DNFP=length(setdiff(gsig_dn3 ,gtdn))
  DNFN=length(setdiff(gtdn,gsig_dn3))
  TP=UPTP+DNTP
  FP=UPFP+DNFP
  FN=UPFN+DNFN
  TN=nrow(aares1)-DNTP-DNFP-DNFN-UPTP-UPFP-UPFN
  PREC=TP/(TP+FP)
  REC=TP/(TP+FN)
  F1=TP/(TP+(0.5*(FP+FN)))
  result <- c("TP"=TP,"FP"=FP,"FN"=FN,"TN"=TN,"PREC"=PREC,"REC"=REC)
  return(result)
}

F1 <- function(x,y) {
  ( 2 * x * y ) / ( x + y )
}

Run analyses

Set assumptions.

num_dm_sets=50
sims=20

groupsizes=c(3,5,8,12)
deltas=c(0.1,0.2,0.3,0.4,0.5)
#groupsizes=5
#deltas=0.4

params <- expand.grid("groupsizes"=groupsizes,"deltas"=deltas)

params %>% kbl(caption="Parameters to test") %>% kable_paper("hover", full_width = F)

Parameters to test
groupsizes	deltas
3	0.1
5	0.1
8	0.1
12	0.1
3	0.2
5	0.2
8	0.2
12	0.2
3	0.3
5	0.3
8	0.3
12	0.3
3	0.4
5	0.4
8	0.4
12	0.4
3	0.5
5	0.5
8	0.5
12	0.5

GSA meth sim

Cannot be run in multicore due to fragility of AnnotationDbi SQLite objects.

gres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- lapply(1:sims,function(i) {
    simgsa(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5,
      frac_probes=0.5, groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  })
  res <- do.call(rbind,res)
  return(res)
})

gres2 <- do.call(rbind,lapply(gres,colMeans))
gres2[,"PREC"] <- gres2[,"TP"] / ( gres2[,"TP"] + gres2[,"FP"] )

gres2 %>% kbl(caption="GSAmeth results") %>% kable_paper("hover", full_width = F)

GSAmeth results
TP	FP	FN	TN	PREC	REC
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.50	0.00	49.50	950.00	1.0000000	0.010
0.05	0.00	49.95	950.00	1.0000000	0.001
0.15	0.05	49.85	949.95	0.7500000	0.003
0.45	0.00	49.55	950.00	1.0000000	0.009
1.30	0.10	48.70	949.90	0.9285714	0.026
0.10	0.00	49.90	950.00	1.0000000	0.002
0.35	0.05	49.65	949.95	0.8750000	0.007
2.55	0.00	47.45	950.00	1.0000000	0.051
19.00	0.45	31.00	949.55	0.9768638	0.380
0.20	0.05	49.80	949.95	0.8000000	0.004
0.60	0.05	49.40	949.95	0.9230769	0.012
15.70	0.50	34.30	949.50	0.9691358	0.314
25.15	0.70	24.85	949.30	0.9729207	0.503

gres3p <- do.call(rbind,lapply(groupsizes, function (g) { gres2[params$groupsizes==g,"PREC"] }))
colnames(gres3p) <- deltas
rownames(gres3p) <- groupsizes
gres3p %>% kbl(caption="GSAmeth precision") %>% kable_paper("hover", full_width = F)

GSAmeth precision
	0.1	0.2	0.3	0.4	0.5
3	NaN	NaN	1.0000000	1.0000000	0.8000000
5	NaN	NaN	0.7500000	0.8750000	0.9230769
8	NaN	NaN	1.0000000	1.0000000	0.9691358
12	NaN	1	0.9285714	0.9768638	0.9729207

gres3r <- do.call(rbind,lapply(groupsizes, function (g) { gres2[params$groupsizes==g,"REC"] }))
colnames(gres3r) <- deltas
rownames(gres3r) <- groupsizes
gres3r %>% kbl(caption="GSAmeth recall") %>% kable_paper("hover", full_width = F)

GSAmeth recall
	0.2	0.3	0.4	0.5
3	0.00	0.001	0.002	0.004
5	0.00	0.003	0.007	0.012
8	0.00	0.009	0.051	0.314
12	0.01	0.026	0.380	0.503

F1(gres3p,gres3r) %>% kbl(caption="GSAmeth F1") %>% kable_paper("hover", full_width = F)

GSAmeth F1
	0.1	0.2	0.3	0.4	0.5
3	NaN	NaN	0.0019980	0.0039920	0.0079602
5	NaN	NaN	0.0059761	0.0138889	0.0236920
8	NaN	NaN	0.0178394	0.0970504	0.4743202
12	NaN	0.019802	0.0505837	0.5471562	0.6631510

LA sim

parametric competitive

lacres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simlac(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

lacres2 <- do.call(rbind,lapply(lacres,colMeans))
lacres2[,"PREC"] <- lacres2[,"TP"] / ( lacres2[,"TP"] + lacres2[,"FP"] )

lacres2 %>% kbl(caption="LA parametric results") %>% kable_paper("hover", full_width = F)

LA parametric results
TP	FP	FN	TN	PREC	REC
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.05	0.05	49.95	949.95	0.5	0.001
0.00	0.00	50.00	950.00	NaN	0.000
0.05	0.00	49.95	950.00	1.0	0.001
0.00	0.00	50.00	950.00	NaN	0.000
0.05	0.00	49.95	950.00	1.0	0.001
0.00	0.00	50.00	950.00	NaN	0.000
0.15	0.00	49.85	950.00	1.0	0.003
0.05	0.00	49.95	950.00	1.0	0.001
0.05	0.00	49.95	950.00	1.0	0.001
0.05	0.00	49.95	950.00	1.0	0.001
0.10	0.00	49.90	950.00	1.0	0.002
0.05	0.00	49.95	950.00	1.0	0.001
0.05	0.00	49.95	950.00	1.0	0.001
0.05	0.00	49.95	950.00	1.0	0.001
0.10	0.00	49.90	950.00	1.0	0.002
0.05	0.00	49.95	950.00	1.0	0.001
0.05	0.00	49.95	950.00	1.0	0.001

lacres3p <- do.call(rbind,lapply(groupsizes, function (g) { lacres2[params$groupsizes==g,"PREC"] }))
colnames(lacres3p) <- deltas
rownames(lacres3p) <- groupsizes
lacres3p %>% kbl(caption="LA parametric precision") %>% kable_paper("hover", full_width = F)

LA parametric precision
	0.1	0.2	0.3	0.4	0.5
3	NaN	NaN	NaN	1	1
5	NaN	1	1	1	1
8	NaN	NaN	1	1	1
12	0.5	1	1	1	1

lacres3r <- do.call(rbind,lapply(groupsizes, function (g) { lacres2[params$groupsizes==g,"REC"] }))
colnames(lacres3r) <- deltas
rownames(lacres3r) <- groupsizes
lacres3r %>% kbl(caption="LA parametric recall") %>% kable_paper("hover", full_width = F)

LA parametric recall
	0.1	0.2	0.3	0.4	0.5
3	0.000	0.000	0.000	0.001	0.001
5	0.000	0.001	0.003	0.002	0.002
8	0.000	0.000	0.001	0.001	0.001
12	0.001	0.001	0.001	0.001	0.001

F1(lacres3p,lacres3r) %>% kbl(caption="LA parametric F1") %>% kable_paper("hover", full_width = F)

LA parametric F1
	0.1	0.2	0.3	0.4	0.5
3	NaN	NaN	NaN	0.001998	0.001998
5	NaN	0.001998	0.0059821	0.003992	0.003992
8	NaN	NaN	0.0019980	0.001998	0.001998
12	0.001996	0.001998	0.0019980	0.001998	0.001998

parametric competitive with “top” aggregation

Too many false positives.

lactopres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simlactop(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

lactopres2 <- do.call(rbind,lapply(lactopres,colMeans))
lactopres2[,"PREC"] <- lactopres2[,"TP"] / ( lactopres2[,"TP"] + lactopres2[,"FP"] )

lactopres2 %>% kbl(caption="LA parametric results (top agg)") %>% kable_paper("hover", full_width = F)

LA parametric results (top agg)
TP	FP	FN	TN	PREC	REC
0.10	0.85	49.90	949.15	0.1052632	0.002
0.20	2.40	49.80	947.60	0.0769231	0.004
0.25	3.20	49.75	946.80	0.0724638	0.005
0.25	1.80	49.75	948.20	0.1219512	0.005
0.15	1.10	49.85	948.90	0.1200000	0.003
0.45	2.25	49.55	947.75	0.1666667	0.009
0.70	3.00	49.30	947.00	0.1891892	0.014
0.75	1.80	49.25	948.20	0.2941176	0.015
0.35	1.10	49.65	948.90	0.2413793	0.007
0.50	1.90	49.50	948.10	0.2083333	0.010
1.20	2.90	48.80	947.10	0.2926829	0.024
1.45	1.65	48.55	948.35	0.4677419	0.029
0.65	1.05	49.35	948.95	0.3823529	0.013
0.75	1.75	49.25	948.25	0.3000000	0.015
1.65	2.65	48.35	947.35	0.3837209	0.033
1.40	1.30	48.60	948.70	0.5185185	0.028
0.80	0.80	49.20	949.20	0.5000000	0.016
0.70	1.10	49.30	948.90	0.3888889	0.014
1.60	2.30	48.40	947.70	0.4102564	0.032
1.45	1.00	48.55	949.00	0.5918367	0.029

lactopres3p <- do.call(rbind,lapply(groupsizes, function (g) { lactopres2[params$groupsizes==g,"PREC"] }))
colnames(lactopres3p) <- deltas
rownames(lactopres3p) <- groupsizes
lactopres3p %>% kbl(caption="LA parametric precision (top agg)") %>% kable_paper("hover", full_width = F)

LA parametric precision (top agg)
	0.1	0.2	0.3	0.4	0.5
3	0.1052632	0.1200000	0.2413793	0.3823529	0.5000000
5	0.0769231	0.1666667	0.2083333	0.3000000	0.3888889
8	0.0724638	0.1891892	0.2926829	0.3837209	0.4102564
12	0.1219512	0.2941176	0.4677419	0.5185185	0.5918367

lactopres3r <- do.call(rbind,lapply(groupsizes, function (g) { lactopres2[params$groupsizes==g,"REC"] }))
colnames(lactopres3r) <- deltas
rownames(lactopres3r) <- groupsizes
lactopres3r %>% kbl(caption="LA parametric recall (top agg)") %>% kable_paper("hover", full_width = F)

LA parametric recall (top agg)
	0.1	0.2	0.3	0.4	0.5
3	0.002	0.003	0.007	0.013	0.016
5	0.004	0.009	0.010	0.015	0.014
8	0.005	0.014	0.024	0.033	0.032
12	0.005	0.015	0.029	0.028	0.029

F1(lactopres3p,lactopres3r) %>% kbl(caption="LA parametric F1 (top agg)") %>% kable_paper("hover", full_width = F)

LA parametric F1 (top agg)
	0.1	0.2	0.3	0.4	0.5
3	0.0039254	0.0058537	0.0136054	0.0251451	0.0310078
5	0.0076046	0.0170778	0.0190840	0.0285714	0.0270270
8	0.0093545	0.0260708	0.0443623	0.0607735	0.0593692
12	0.0096061	0.0285442	0.0546139	0.0531309	0.0552908

nonparametric competitive

nlacres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simnlac(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

nlacres2 <- do.call(rbind,lapply(nlacres,colMeans))
nlacres2[,"PREC"] <- nlacres2[,"TP"] / ( nlacres2[,"TP"] + nlacres2[,"FP"] )

nlacres2 %>% kbl(caption="LA nonparametric results") %>% kable_paper("hover", full_width = F)

LA nonparametric results
TP	FP	FN	TN	PREC	REC
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.05	0.00	49.95	950.00	1.0000000	0.001
0.25	0.00	49.75	950.00	1.0000000	0.005
0.35	0.05	49.65	949.95	0.8750000	0.007
0.60	0.05	49.40	949.95	0.9230769	0.012
0.15	0.00	49.85	950.00	1.0000000	0.003
0.50	0.00	49.50	950.00	1.0000000	0.010
1.15	0.05	48.85	949.95	0.9583333	0.023
1.70	0.05	48.30	949.95	0.9714286	0.034
0.30	0.00	49.70	950.00	1.0000000	0.006
0.65	0.00	49.35	950.00	1.0000000	0.013
2.05	0.10	47.95	949.90	0.9534884	0.041
2.60	0.10	47.40	949.90	0.9629630	0.052
0.90	0.00	49.10	950.00	1.0000000	0.018
1.30	0.00	48.70	950.00	1.0000000	0.026
3.05	0.15	46.95	949.85	0.9531250	0.061
3.45	0.15	46.55	949.85	0.9583333	0.069

nlacres3p <- do.call(rbind,lapply(groupsizes, function (g) { nlacres2[params$groupsizes==g,"PREC"] }))
colnames(nlacres3p) <- deltas
rownames(nlacres3p) <- groupsizes
nlacres3p %>% kbl(caption="LA nonparametric precision") %>% kable_paper("hover", full_width = F)

LA nonparametric precision
	0.1	0.2	0.3	0.4	0.5
3	NaN	1.0000000	1.0000000	1.0000000	1.0000000
5	NaN	1.0000000	1.0000000	1.0000000	1.0000000
8	0	0.8750000	0.9583333	0.9534884	0.9531250
12	0	0.9230769	0.9714286	0.9629630	0.9583333

nlacres3r <- do.call(rbind,lapply(groupsizes, function (g) { nlacres2[params$groupsizes==g,"REC"] }))
colnames(nlacres3r) <- deltas
rownames(nlacres3r) <- groupsizes
nlacres3r %>% kbl(caption="LA nonparametric recall") %>% kable_paper("hover", full_width = F)

LA nonparametric recall
	0.2	0.3	0.4	0.5
3	0.001	0.003	0.006	0.018
5	0.005	0.010	0.013	0.026
8	0.007	0.023	0.041	0.061
12	0.012	0.034	0.052	0.069

F1(nlacres3p,nlacres3r) %>% kbl(caption="LA nonparametric F1") %>% kable_paper("hover", full_width = F)

LA nonparametric F1
	0.1	0.2	0.3	0.4	0.5
3	NaN	0.0019980	0.0059821	0.0119284	0.0353635
5	NaN	0.0099502	0.0198020	0.0256663	0.0506823
8	NaN	0.0138889	0.0449219	0.0786194	0.1146617
12	NaN	0.0236920	0.0657005	0.0986717	0.1287313

LA competitive mitch

lacmres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simlacm(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

lacmres2 <- do.call(rbind,lapply(lacmres,colMeans))
lacmres2[,"PREC"] <- lacmres2[,"TP"] / ( lacmres2[,"TP"] + lacmres2[,"FP"] )

lacmres2 %>% kbl(caption="LA mitch results") %>% kable_paper("hover", full_width = F)

LA mitch results
TP	FP	FN	TN	PREC	REC
0.05	0.00	49.95	950.00	1.0000000	0.001
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.05	0.00	49.95	950.00	1.0000000	0.001
0.25	0.00	49.75	950.00	1.0000000	0.005
0.35	0.05	49.65	949.95	0.8750000	0.007
0.60	0.05	49.40	949.95	0.9230769	0.012
0.20	0.00	49.80	950.00	1.0000000	0.004
0.50	0.00	49.50	950.00	1.0000000	0.010
1.15	0.05	48.85	949.95	0.9583333	0.023
1.70	0.05	48.30	949.95	0.9714286	0.034
0.30	0.00	49.70	950.00	1.0000000	0.006
0.75	0.00	49.25	950.00	1.0000000	0.015
2.05	0.10	47.95	949.90	0.9534884	0.041
2.60	0.10	47.40	949.90	0.9629630	0.052
0.90	0.00	49.10	950.00	1.0000000	0.018
1.30	0.00	48.70	950.00	1.0000000	0.026
3.05	0.15	46.95	949.85	0.9531250	0.061
3.50	0.15	46.50	949.85	0.9589041	0.070

lacmres3p <- do.call(rbind,lapply(groupsizes, function (g) { lacmres2[params$groupsizes==g,"PREC"] }))
colnames(lacmres3p) <- deltas
rownames(lacmres3p) <- groupsizes
lacmres3p %>% kbl(caption="LA mitch precision") %>% kable_paper("hover", full_width = F)

LA mitch precision
	0.1	0.2	0.3	0.4	0.5
3	1	1.0000000	1.0000000	1.0000000	1.0000000
5	NaN	1.0000000	1.0000000	1.0000000	1.0000000
8	0	0.8750000	0.9583333	0.9534884	0.9531250
12	0	0.9230769	0.9714286	0.9629630	0.9589041

lacmres3r <- do.call(rbind,lapply(groupsizes, function (g) { lacmres2[params$groupsizes==g,"REC"] }))
colnames(lacmres3r) <- deltas
rownames(lacmres3r) <- groupsizes
lacmres3r %>% kbl(caption="LA mitch recall") %>% kable_paper("hover", full_width = F)

LA mitch recall
	0.1	0.2	0.3	0.4	0.5
3	0.001	0.001	0.004	0.006	0.018
5	0.000	0.005	0.010	0.015	0.026
8	0.000	0.007	0.023	0.041	0.061
12	0.000	0.012	0.034	0.052	0.070

F1(lacmres3p,lacmres3r) %>% kbl(caption="LA mitch F1") %>% kable_paper("hover", full_width = F)

LA mitch F1
	0.1	0.2	0.3	0.4	0.5
3	0.001998	0.0019980	0.0079681	0.0119284	0.0353635
5	NaN	0.0099502	0.0198020	0.0295567	0.0506823
8	NaN	0.0138889	0.0449219	0.0786194	0.1146617
12	NaN	0.0236920	0.0657005	0.0986717	0.1304753

LA rank probes

lrmres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simlrm(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

lrmres2 <- do.call(rbind,lapply(lrmres,colMeans))
lrmres2[,"PREC"] <- lrmres2[,"TP"] / ( lrmres2[,"TP"] + lrmres2[,"FP"] )

lrmres2 %>% kbl(caption="LA rank mitch results") %>% kable_paper("hover", full_width = F)

LA rank mitch results
TP	FP	FN	TN	PREC	REC
0.05	0.00	49.95	950.00	1.0000000	0.001
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.05	0.00	49.95	950.00	1.0000000	0.001
0.20	0.00	49.80	950.00	1.0000000	0.004
0.20	0.05	49.80	949.95	0.8000000	0.004
0.40	0.05	49.60	949.95	0.8888889	0.008
0.20	0.00	49.80	950.00	1.0000000	0.004
0.40	0.00	49.60	950.00	1.0000000	0.008
0.45	0.05	49.55	949.95	0.9000000	0.009
0.55	0.05	49.45	949.95	0.9166667	0.011
0.20	0.00	49.80	950.00	1.0000000	0.004
0.40	0.00	49.60	950.00	1.0000000	0.008
0.45	0.05	49.55	949.95	0.9000000	0.009
0.60	0.05	49.40	949.95	0.9230769	0.012
0.40	0.00	49.60	950.00	1.0000000	0.008
0.50	0.00	49.50	950.00	1.0000000	0.010
0.45	0.05	49.55	949.95	0.9000000	0.009
0.70	0.05	49.30	949.95	0.9333333	0.014

lrmres3p <- do.call(rbind,lapply(groupsizes, function (g) { lrmres2[params$groupsizes==g,"PREC"] }))
colnames(lrmres3p) <- deltas
rownames(lrmres3p) <- groupsizes
lrmres3p %>% kbl(caption="LA rank mitch precision") %>% kable_paper("hover", full_width = F)

LA rank mitch precision
	0.1	0.2	0.3	0.4	0.5
3	1	1.0000000	1.0000000	1.0000000	1.0000000
5	NaN	1.0000000	1.0000000	1.0000000	1.0000000
8	0	0.8000000	0.9000000	0.9000000	0.9000000
12	NaN	0.8888889	0.9166667	0.9230769	0.9333333

lrmres3r <- do.call(rbind,lapply(groupsizes, function (g) { lrmres2[params$groupsizes==g,"REC"] }))
colnames(lrmres3r) <- deltas
rownames(lrmres3r) <- groupsizes
lrmres3r %>% kbl(caption="LA rank mitch recall") %>% kable_paper("hover", full_width = F)

LA rank mitch recall
	0.1	0.2	0.3	0.4	0.5
3	0.001	0.001	0.004	0.004	0.008
5	0.000	0.004	0.008	0.008	0.010
8	0.000	0.004	0.009	0.009	0.009
12	0.000	0.008	0.011	0.012	0.014

F1(lrmres3p,lacmres3r) %>% kbl(caption="LA rank mitch F1") %>% kable_paper("hover", full_width = F)

LA rank mitch F1
	0.1	0.2	0.3	0.4	0.5
3	0.001998	0.0019980	0.0079681	0.0119284	0.0353635
5	NaN	0.0099502	0.0198020	0.0295567	0.0506823
8	NaN	0.0138786	0.0448537	0.0784272	0.1142560
12	NaN	0.0236803	0.0655680	0.0984538	0.1302326

AL sim

parametric competitive

alcres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simalc(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

alcres2 <- do.call(rbind,lapply(alcres,colMeans))
alcres2[,"PREC"] <- alcres2[,"TP"] / ( alcres2[,"TP"] + alcres2[,"FP"] )

alcres2 %>% kbl(caption="AL parametric results") %>% kable_paper("hover", full_width = F)

AL parametric results
TP	FP	FN	TN	PREC	REC
0.00	0.10	50.00	949.90	0.0000000	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.10	50.00	949.90	0.0000000	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.15	0.00	49.85	950.00	1.0000000	0.003
0.05	0.05	49.95	949.95	0.5000000	0.001
0.05	0.05	49.95	949.95	0.5000000	0.001
0.10	0.05	49.90	949.95	0.6666667	0.002
0.15	0.00	49.85	950.00	1.0000000	0.003
0.05	0.05	49.95	949.95	0.5000000	0.001
0.00	0.00	50.00	950.00	NaN	0.000
0.05	0.05	49.95	949.95	0.5000000	0.001
0.10	0.00	49.90	950.00	1.0000000	0.002
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.00	50.00	950.00	NaN	0.000

alcres3p <- do.call(rbind,lapply(groupsizes, function (g) { alcres2[params$groupsizes==g,"PREC"] }))
colnames(alcres3p) <- deltas
rownames(alcres3p) <- groupsizes
alcres3p %>% kbl(caption="AL parametric precision") %>% kable_paper("hover", full_width = F)

AL parametric precision
	0.1	0.2	0.3	0.4	0.5
3	0	0	0.0	0.6666667	0.5
5	NaN	NaN	1.0	1.0000000	1.0
8	0	0	0.5	0.5000000	NaN
12	0	0	0.5	NaN	NaN

alcres3r <- do.call(rbind,lapply(groupsizes, function (g) { alcres2[params$groupsizes==g,"REC"] }))
colnames(alcres3r) <- deltas
rownames(alcres3r) <- groupsizes
alcres3r %>% kbl(caption="AL parametric recall") %>% kable_paper("hover", full_width = F)

AL parametric recall
	0.3	0.4	0.5
3	0.000	0.002	0.001
5	0.003	0.003	0.002
8	0.001	0.001	0.000
12	0.001	0.000	0.000

F1(alcres3p,alcres3r) %>% kbl(caption="AL parametric F1") %>% kable_paper("hover", full_width = F)

AL parametric F1
	0.1	0.2	0.3	0.4	0.5
3	NaN	NaN	NaN	0.0039880	0.001996
5	NaN	NaN	0.0059821	0.0059821	0.003992
8	NaN	NaN	0.0019960	0.0019960	NaN
12	NaN	NaN	0.0019960	NaN	NaN

nonparametric competitive

nalcres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simnalc(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

nalcres2 <- do.call(rbind,lapply(nalcres,colMeans))
nalcres2[,"PREC"] <- nalcres2[,"TP"] / ( nalcres2[,"TP"] + nalcres2[,"FP"] )

nalcres2 %>% kbl(caption="AL nonparametric results") %>% kable_paper("hover", full_width = F)

AL nonparametric results
TP	FP	FN	TN	PREC	REC
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.10	50.00	949.90	0.0000000	0.000
0.05	0.00	49.95	950.00	1.0000000	0.001
0.25	0.00	49.75	950.00	1.0000000	0.005
0.15	0.05	49.85	949.95	0.7500000	0.003
0.60	0.05	49.40	949.95	0.9230769	0.012
0.05	0.00	49.95	950.00	1.0000000	0.001
0.40	0.00	49.60	950.00	1.0000000	0.008
0.75	0.15	49.25	949.85	0.8333333	0.015
1.65	0.05	48.35	949.95	0.9705882	0.033
0.30	0.10	49.70	949.90	0.7500000	0.006
0.85	0.00	49.15	950.00	1.0000000	0.017
1.20	0.15	48.80	949.85	0.8888889	0.024
3.05	0.10	46.95	949.90	0.9682540	0.061
0.85	0.15	49.15	949.85	0.8500000	0.017
1.30	0.05	48.70	949.95	0.9629630	0.026
2.20	0.15	47.80	949.85	0.9361702	0.044
3.30	0.10	46.70	949.90	0.9705882	0.066

nalcres3p <- do.call(rbind,lapply(groupsizes, function (g) { nalcres2[params$groupsizes==g,"PREC"] }))
colnames(nalcres3p) <- deltas
rownames(nalcres3p) <- groupsizes
nalcres3p %>% kbl(caption="AL nonparametric precision") %>% kable_paper("hover", full_width = F)

AL nonparametric precision
	0.1	0.2	0.3	0.4	0.5
3	0	1.0000000	1.0000000	0.7500000	0.8500000
5	NaN	1.0000000	1.0000000	1.0000000	0.9629630
8	0	0.7500000	0.8333333	0.8888889	0.9361702
12	0	0.9230769	0.9705882	0.9682540	0.9705882

nalcres3r <- do.call(rbind,lapply(groupsizes, function (g) { nalcres2[params$groupsizes==g,"REC"] }))
colnames(nalcres3r) <- deltas
rownames(nalcres3r) <- groupsizes
nalcres3r %>% kbl(caption="AL nonparametric recall") %>% kable_paper("hover", full_width = F)

AL nonparametric recall
	0.2	0.3	0.4	0.5
3	0.001	0.001	0.006	0.017
5	0.005	0.008	0.017	0.026
8	0.003	0.015	0.024	0.044
12	0.012	0.033	0.061	0.066

F1(nalcres3p,nalcres3r) %>% kbl(caption="AL nonparametric F1") %>% kable_paper("hover", full_width = F)

AL nonparametric F1
	0.1	0.2	0.3	0.4	0.5
3	NaN	0.0019980	0.0019980	0.0119048	0.0333333
5	NaN	0.0099502	0.0158730	0.0334317	0.0506329
8	NaN	0.0059761	0.0294695	0.0467381	0.0840497
12	NaN	0.0236920	0.0638298	0.1147695	0.1235955

ALM competitive test

almres <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simalm(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

almres2 <- do.call(rbind,lapply(almres,colMeans))
almres2[,"PREC"] <- almres2[,"TP"] / ( almres2[,"TP"] + almres2[,"FP"] )

almres2 %>% kbl(caption="ALM results") %>% kable_paper("hover", full_width = F)

ALM results
TP	FP	FN	TN	PREC	REC
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.00	50.00	950.00	NaN	0.000
0.00	0.05	50.00	949.95	0.0000000	0.000
0.00	0.10	50.00	949.90	0.0000000	0.000
0.05	0.00	49.95	950.00	1.0000000	0.001
0.30	0.00	49.70	950.00	1.0000000	0.006
0.20	0.10	49.80	949.90	0.6666667	0.004
0.60	0.05	49.40	949.95	0.9230769	0.012
0.05	0.00	49.95	950.00	1.0000000	0.001
0.40	0.00	49.60	950.00	1.0000000	0.008
0.75	0.15	49.25	949.85	0.8333333	0.015
1.65	0.05	48.35	949.95	0.9705882	0.033
0.35	0.05	49.65	949.95	0.8750000	0.007
0.85	0.00	49.15	950.00	1.0000000	0.017
1.40	0.15	48.60	949.85	0.9032258	0.028
3.05	0.10	46.95	949.90	0.9682540	0.061
0.90	0.10	49.10	949.90	0.9000000	0.018
1.30	0.05	48.70	949.95	0.9629630	0.026
2.20	0.15	47.80	949.85	0.9361702	0.044
3.30	0.10	46.70	949.90	0.9705882	0.066

almres3p <- do.call(rbind,lapply(groupsizes, function (g) { almres2[params$groupsizes==g,"PREC"] }))
colnames(almres3p) <- deltas
rownames(almres3p) <- groupsizes
almres3p %>% kbl(caption="ALM precision") %>% kable_paper("hover", full_width = F)

ALM precision
	0.1	0.2	0.3	0.4	0.5
3	0	1.0000000	1.0000000	0.8750000	0.9000000
5	NaN	1.0000000	1.0000000	1.0000000	0.9629630
8	0	0.6666667	0.8333333	0.9032258	0.9361702
12	0	0.9230769	0.9705882	0.9682540	0.9705882

almres3r <- do.call(rbind,lapply(groupsizes, function (g) { almres2[params$groupsizes==g,"REC"] }))
colnames(almres3r) <- deltas
rownames(almres3r) <- groupsizes
almres3r %>% kbl(caption="ALM recall") %>% kable_paper("hover", full_width = F)

ALM recall
	0.2	0.3	0.4	0.5
3	0.001	0.001	0.007	0.018
5	0.006	0.008	0.017	0.026
8	0.004	0.015	0.028	0.044
12	0.012	0.033	0.061	0.066

F1(almres3p,almres3r) %>% kbl(caption="ALM F1") %>% kable_paper("hover", full_width = F)

ALM F1
	0.1	0.2	0.3	0.4	0.5
3	NaN	0.0019980	0.0019980	0.0138889	0.0352941
5	NaN	0.0119284	0.0158730	0.0334317	0.0506329
8	NaN	0.0079523	0.0294695	0.0543162	0.0840497
12	NaN	0.0236920	0.0638298	0.1147695	0.1235955

AA sim

aares <- lapply(1:nrow(params) , function(j) {
  groupsize = params[j,1]
  delta = params[j,2]
  res <- mclapply(1:sims,function(i) {
    simaa(genesetdatabase=gsets, myann=myann, mval=normal_mval, seed=i*100, frac_genes=0.5, frac_probes=0.5,
      groupsize=groupsize, delta=delta, num_dm_sets=num_dm_sets)
  },mc.cores=6)
  res <- do.call(rbind,res)
  return(res)
})

aares2 <- do.call(rbind,lapply(aares,colMeans))
aares2[,"PREC"] <- aares2[,"TP"] / ( aares2[,"TP"] + aares2[,"FP"] )

aares2 %>% kbl(caption="AA results") %>% kable_paper("hover", full_width = F)

AA results
TP	FP	FN	TN	PREC	REC
2.35	82.90	47.65	867.10	0.0275660	0.047
1.05	38.70	48.95	911.30	0.0264151	0.021
0.35	10.50	49.65	939.50	0.0322581	0.007
0.00	0.10	50.00	949.90	0.0000000	0.000
2.45	82.60	47.55	867.40	0.0288066	0.049
1.05	37.80	48.95	912.20	0.0270270	0.021
0.50	10.20	49.50	939.80	0.0467290	0.010
0.00	0.10	50.00	949.90	0.0000000	0.000
2.45	81.15	47.55	868.85	0.0293062	0.049
1.15	36.85	48.85	913.15	0.0302632	0.023
0.80	11.35	49.20	938.65	0.0658436	0.016
0.00	0.10	50.00	949.90	0.0000000	0.000
2.50	79.80	47.50	870.20	0.0303767	0.050
1.85	36.75	48.15	913.25	0.0479275	0.037
1.60	12.30	48.40	937.70	0.1151079	0.032
0.30	0.10	49.70	949.90	0.7500000	0.006
3.00	79.00	47.00	871.00	0.0365854	0.060
2.30	36.85	47.70	913.15	0.0587484	0.046
2.50	12.95	47.50	937.05	0.1618123	0.050
2.65	0.10	47.35	949.90	0.9636364	0.053

aares3p <- do.call(rbind,lapply(groupsizes, function (g) { aares2[params$groupsizes==g,"PREC"] }))
colnames(aares3p) <- deltas
rownames(aares3p) <- groupsizes
aares3p %>% kbl(caption="AA precision") %>% kable_paper("hover", full_width = F)

AA precision
	0.1	0.2	0.3	0.4	0.5
3	0.0275660	0.0288066	0.0293062	0.0303767	0.0365854
5	0.0264151	0.0270270	0.0302632	0.0479275	0.0587484
8	0.0322581	0.0467290	0.0658436	0.1151079	0.1618123
12	0.0000000	0.0000000	0.0000000	0.7500000	0.9636364

aares3r <- do.call(rbind,lapply(groupsizes, function (g) { aares2[params$groupsizes==g,"REC"] }))
colnames(aares3r) <- deltas
rownames(aares3r) <- groupsizes
aares3r %>% kbl(caption="AA recall") %>% kable_paper("hover", full_width = F)

AA recall
	0.1	0.2	0.3	0.4	0.5
3	0.047	0.049	0.049	0.050	0.060
5	0.021	0.021	0.023	0.037	0.046
8	0.007	0.010	0.016	0.032	0.050
12	0.000	0.000	0.000	0.006	0.053

F1(aares3p,almres3r) %>% kbl(caption="AA F1") %>% kable_paper("hover", full_width = F)

AA F1
	0.1	0.2	0.3	0.4	0.5
3	0	0.0019329	0.0019340	0.0113780	0.0241287
5	0	0.0098200	0.0126547	0.0250978	0.0360469
8	0	0.0073692	0.0244337	0.0450432	0.0691867
12	NaN	0.0000000	0.0000000	0.1128237	0.1235388

Session information

save.image("GSE158422_simulate020.Rdata")

sessionInfo()

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.4                                   
##  [2] mitch_1.15.0                                       
##  [3] psych_2.3.9                                        
##  [4] org.Hs.eg.db_3.18.0                                
##  [5] AnnotationDbi_1.64.1                               
##  [6] HGNChelper_0.8.1                                   
##  [7] IlluminaHumanMethylation450kmanifest_0.4.0         
##  [8] missMethyl_1.36.0                                  
##  [9] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [10] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1 
## [11] minfi_1.48.0                                       
## [12] bumphunter_1.44.0                                  
## [13] locfit_1.5-9.8                                     
## [14] iterators_1.0.14                                   
## [15] foreach_1.5.2                                      
## [16] Biostrings_2.70.1                                  
## [17] XVector_0.42.0                                     
## [18] SummarizedExperiment_1.32.0                        
## [19] Biobase_2.62.0                                     
## [20] MatrixGenerics_1.14.0                              
## [21] matrixStats_1.2.0                                  
## [22] GenomicRanges_1.54.1                               
## [23] GenomeInfoDb_1.38.2                                
## [24] IRanges_2.36.0                                     
## [25] S4Vectors_0.40.2                                   
## [26] BiocGenerics_0.48.1                                
## [27] limma_3.58.1                                       
## [28] stringi_1.8.3                                      
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.2             later_1.3.2              
##   [3] BiocIO_1.12.0             bitops_1.0-7             
##   [5] filelock_1.0.3            BiasedUrn_2.0.11         
##   [7] tibble_3.2.1              preprocessCore_1.64.0    
##   [9] XML_3.99-0.16             lifecycle_1.0.4          
##  [11] lattice_0.22-5            MASS_7.3-60              
##  [13] base64_2.0.1              scrime_1.3.5             
##  [15] magrittr_2.0.3            sass_0.4.8               
##  [17] rmarkdown_2.25            jquerylib_0.1.4          
##  [19] yaml_2.3.8                httpuv_1.6.13            
##  [21] doRNG_1.8.6               askpass_1.2.0            
##  [23] DBI_1.1.3                 RColorBrewer_1.1-3       
##  [25] abind_1.4-5               zlibbioc_1.48.0          
##  [27] rvest_1.0.3               quadprog_1.5-8           
##  [29] purrr_1.0.2               RCurl_1.98-1.13          
##  [31] rappdirs_0.3.3            GenomeInfoDbData_1.2.11  
##  [33] genefilter_1.84.0         annotate_1.80.0          
##  [35] svglite_2.1.3             DelayedMatrixStats_1.24.0
##  [37] codetools_0.2-19          DelayedArray_0.28.0      
##  [39] xml2_1.3.6                tidyselect_1.2.0         
##  [41] beanplot_1.3.1            BiocFileCache_2.10.1     
##  [43] webshot_0.5.5             illuminaio_0.44.0        
##  [45] GenomicAlignments_1.38.0  jsonlite_1.8.8           
##  [47] multtest_2.58.0           ellipsis_0.3.2           
##  [49] survival_3.5-7            systemfonts_1.0.5        
##  [51] tools_4.3.2               progress_1.2.3           
##  [53] Rcpp_1.0.11               glue_1.6.2               
##  [55] gridExtra_2.3             mnormt_2.1.1             
##  [57] SparseArray_1.2.2         xfun_0.41                
##  [59] dplyr_1.1.4               HDF5Array_1.30.0         
##  [61] fastmap_1.1.1             GGally_2.2.0             
##  [63] rhdf5filters_1.14.1       fansi_1.0.6              
##  [65] openssl_2.1.1             caTools_1.18.2           
##  [67] digest_0.6.33             R6_2.5.1                 
##  [69] mime_0.12                 colorspace_2.1-0         
##  [71] gtools_3.9.5              biomaRt_2.58.0           
##  [73] RSQLite_2.3.4             utf8_1.2.4               
##  [75] tidyr_1.3.0               generics_0.1.3           
##  [77] data.table_1.14.10        rtracklayer_1.62.0       
##  [79] prettyunits_1.2.0         httr_1.4.7               
##  [81] htmlwidgets_1.6.4         S4Arrays_1.2.0           
##  [83] ggstats_0.5.1             pkgconfig_2.0.3          
##  [85] gtable_0.3.4              blob_1.2.4               
##  [87] siggenes_1.76.0           htmltools_0.5.7          
##  [89] echarts4r_0.4.5           scales_1.3.0             
##  [91] png_0.1-8                 rstudioapi_0.15.0        
##  [93] knitr_1.45                reshape2_1.4.4           
##  [95] tzdb_0.4.0                rjson_0.2.21             
##  [97] nlme_3.1-163              curl_5.2.0               
##  [99] cachem_1.0.8              rhdf5_2.46.1             
## [101] stringr_1.5.1             KernSmooth_2.23-22       
## [103] restfulr_0.0.15           GEOquery_2.70.0          
## [105] pillar_1.9.0              grid_4.3.2               
## [107] reshape_0.8.9             vctrs_0.6.5              
## [109] gplots_3.1.3              promises_1.2.1           
## [111] dbplyr_2.4.0              xtable_1.8-4             
## [113] beeswarm_0.4.0            evaluate_0.23            
## [115] readr_2.1.4               GenomicFeatures_1.54.1   
## [117] cli_3.6.2                 compiler_4.3.2           
## [119] Rsamtools_2.18.0          rlang_1.1.2              
## [121] crayon_1.5.2              rngtools_1.5.2           
## [123] nor1mix_1.3-2             mclust_6.0.1             
## [125] plyr_1.8.9                viridisLite_0.4.2        
## [127] BiocParallel_1.36.0       munsell_0.5.0            
## [129] Matrix_1.6-1.1            hms_1.1.3                
## [131] sparseMatrixStats_1.14.0  bit64_4.0.5              
## [133] ggplot2_3.4.4             Rhdf5lib_1.24.1          
## [135] KEGGREST_1.42.0           statmod_1.5.0            
## [137] shiny_1.8.0               highr_0.10               
## [139] memoise_2.0.1             bslib_0.6.1              
## [141] bit_4.0.5

Notes

simla: parametric self-contained
simlac: parametric competitive
simnla: nonparametric self-contained
simnlac: nonparametric competitive

simal: parametric self-contained
simalc: parametric competitive
simnal: nonparametric self-contained
simnalc: nonparametric competitive

Testing gene methylation enrichment analysis approaches using simulated data

The GMEA team

2023-12-28

Introduction

Load data

Gene set database

Incorporate changes

GSAMETH function

LA function

AL: Aggregate limma

Agg-limma-mitch function

AA Aggregate-aggregate-limma functions

Run analyses

GSA meth sim

LA sim

parametric competitive

parametric competitive with “top” aggregation

nonparametric competitive

LA competitive mitch

LA rank probes

AL sim

parametric competitive

nonparametric competitive

ALM competitive test

AA sim

Session information

Notes

groupsizes	deltas
3	0.1
5	0.1
8	0.1
12	0.1
3	0.2
5	0.2
8	0.2
12	0.2
3	0.3
5	0.3
8	0.3
12	0.3
3	0.4
5	0.4
8	0.4
12	0.4
3	0.5
5	0.5
8	0.5
12	0.5

groupsizes	deltas
3	0.1
5	0.1
8	0.1
12	0.1
3	0.2
5	0.2
8	0.2
12	0.2
3	0.3
5	0.3
8	0.3
12	0.3
3	0.4
5	0.4
8	0.4
12	0.4
3	0.5
5	0.5
8	0.5
12	0.5

groupsizes	deltas
3	0.1
5	0.1
8	0.1
12	0.1
3	0.2
5	0.2
8	0.2
12	0.2
3	0.3
5	0.3
8	0.3
12	0.3
3	0.4
5	0.4
8	0.4
12	0.4
3	0.5
5	0.5
8	0.5
12	0.5