---
title: "Examining background gene lists using simulation - comparison of different enrichment techniques"
author: "Mark Ziemann & Anusuiya Bora"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    fig_width: 7
    fig_height: 5
theme: cosmo
---

Source: https://github.com/markziemann/background

```{r,libs}

library("parallel")
library("edgeR")
library("DESeq2")
library("limma")
library("stringi")
library("kableExtra")
library("fgsea")
library("clusterProfiler")
source("simpw_func.R")

```

## Intro

TODO


xxx object slots

1. expression counts

2. ground truth up genes

3. ground truth down genes

4. ground truth up gene sets

5. ground truth down gene sets

6. DE result (DESeq2)

7. DE genes up observed

8. DE genes down observed

9. clusterprofiler_default up gene sets

10. clusterprofiler_default down gene sets

11. clusterprofiler bg fix up gene sets

12. clusterprofiler bg fix down gene sets

13. clusterprofiler FDR fix up gene sets

14. clusterprofiler FDR fix down gene sets

15. clusterprofiler BG and FDR fix up gene sets

16. clusterprofiler BG and FDR fix down gene sets

17. fora up gene sets

18. fora down gene sets

19. fgsea up gene sets

20. fgsea down gene sets

# Get count data

```{r,getcounts,warning=FALSE}

a <- countData()

```

# Generate gene sets

```{r,genesetgeneration}

gsets <- randomGeneSets(a,setsize=30,nsets=200)

```

# run simulations over a range of parameters

```{r,setparameters}

SIMS=1000
#run 500 sims
FRAC_DE=0.05
FC=0.5
N_REPS=3
DGE_FUNC="deseq2"
SUM_COUNT=2e7
VARIANCE=c(0,0.1,0.2,0.3,0.4,0.5,0.6)

mygrid <- expand.grid(FRAC_DE,FC,N_REPS,DGE_FUNC,SUM_COUNT,VARIANCE)
colnames(mygrid) <- c("FRAC_DE","FC","N_REPS","DGE_FUNC","SUM_COUNT","VARIANCE")

mygrid

```

Now run the analysis.

```{r,sims}

res <- lapply(1:nrow(mygrid), function(i) {
  FRAC_DE=mygrid[i,"FRAC_DE"]
  FC=mygrid[i,"FC"]
  N_REPS=mygrid[i,"N_REPS"]
  DGE_FUNC=as.character(mygrid[i,"DGE_FUNC"])
  SUM_COUNT=mygrid[i,"SUM_COUNT"]
  VARIANCE=mygrid[i,"VARIANCE"]
  x <- agg_dge(a,N_REPS,SUM_COUNT,VARIANCE,FRAC_DE,FC,SIMS,DGE_FUNC,gsets,cores=32)
  as.data.frame(do.call(rbind, x))
})

res <- do.call(rbind,res)

saveRDS(res,file="res.Rds")

```

Now show the results.

```{r,results}

res %>% kbl(caption="simulation_results") %>% kable_paper("hover", full_width = F)

```

## Plot

fig1_ORAjac.png

```{r,plot1}

cp <- subset(res,PWAY_FUNC == "clusterProfiler default")
cpbg <- subset(res,PWAY_FUNC == "clusterProfiler BG fix")
cpfdr <- subset(res,PWAY_FUNC == "clusterProfiler FDR fix")
cpbgfdr <- subset(res,PWAY_FUNC == "clusterProfiler BG and FDR fix")
fo <- subset(res,PWAY_FUNC == "fora")
fg <- subset(res,PWAY_FUNC == "fgsea")

cp %>% kbl(caption="clusterProfiler") %>% kable_paper("hover", full_width = F)
cpbg %>% kbl(caption="clusterProfiler with background bug fix") %>% kable_paper("hover", full_width = F)
cpfdr %>% kbl(caption="clusterProfiler with FDR bug fix") %>% kable_paper("hover", full_width = F)
cpbgfdr %>% kbl(caption="clusterProfiler with background and FDR fixes") %>% kable_paper("hover", full_width = F)
fo %>% kbl(caption="fora") %>% kable_paper("hover", full_width = F)
fg %>% kbl(caption="fg") %>% kable_paper("hover", full_width = F)

par(mfrow=c(1,3))
par(mar=c(c(5.1, 5.1, 2.1, 2.1) ))

plot(cp$VARIANCE,cp$p,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="precision")
points(cpbg$VARIANCE,cpbg$p,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$p,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$p,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$p,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$p,type="b",pch=19,col="blue")

legend("bottomleft", inset=.02, title="tool",
   c("CP default","CP BG fix","CP FDR fix","CP BG and FDR fix","fora","fgsea"),
   col=c("black","orange","darkgreen","purple","red","blue"),horiz=FALSE, cex=1.0, pch=19,lwd=2)

plot(cp$VARIANCE,cp$r,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="recall")
points(cpbg$VARIANCE,cpbg$r,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$r,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$r,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$r,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$r,type="b",pch=19,col="blue")

plot(cp$VARIANCE,cp$f,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="f1")
points(cpbg$VARIANCE,cpbg$f,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$f,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$f,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$f,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$f,type="b",pch=19,col="blue")


png("fig3_sim.png", width=7,height=5,units="in",res=150,pointsize=12)
par(mar=c(c(5.1, 5.1, 2.1, 2.1) ))
par(mfrow=c(1,3))

plot(cp$VARIANCE,cp$p,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="precision")
points(cpbg$VARIANCE,cpbg$p,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$p,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$p,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$p,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$p,type="b",pch=19,col="blue")
grid()
legend("bottomleft", inset=.02, title="tool",
   c("CP default","CP BG fix","CP FDR fix","CP BG and FDR fix","fora","fgsea"),
   col=c("black","orange","darkgreen","purple","red","blue"),
   horiz=FALSE, cex=1.0, pch=19,lwd=2)

plot(cp$VARIANCE,cp$r,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="recall")
points(cpbg$VARIANCE,cpbg$r,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$r,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$r,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$r,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$r,type="b",pch=19,col="blue")
grid()

plot(cp$VARIANCE,cp$f,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="f1")
points(cpbg$VARIANCE,cpbg$f,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$f,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$f,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$f,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$f,type="b",pch=19,col="blue")
grid()
dev.off()

#################################

pdf("fig3_sim.pdf", width=7,height=5,pointsize=12)
par(mar=c(c(5.1, 5.1, 2.1, 2.1) ))
par(mfrow=c(1,3))

plot(cp$VARIANCE,cp$p,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="precision")
points(cpbg$VARIANCE,cpbg$p,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$p,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$p,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$p,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$p,type="b",pch=19,col="blue")
grid()
legend("bottomleft", inset=.02, title="tool",
   c("CP default","CP BG fix","CP FDR fix","CP BG and FDR fix","fora","fgsea"),
   col=c("black","orange","darkgreen","purple","red","blue"),
   horiz=FALSE, cex=1.0, pch=19,lwd=2)

plot(cp$VARIANCE,cp$r,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="recall")
points(cpbg$VARIANCE,cpbg$r,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$r,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$r,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$r,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$r,type="b",pch=19,col="blue")
grid()

plot(cp$VARIANCE,cp$f,ylim=c(0,1),type="b",pch=19,xlab="variance added",ylab="index",main="f1")
points(cpbg$VARIANCE,cpbg$f,type="b",pch=19,col="orange")
points(cpfdr$VARIANCE,cpfdr$f,type="b",pch=19,col="darkgreen")
points(cpbgfdr$VARIANCE,cpbgfdr$f,type="b",pch=19,col="purple")
points(fo$VARIANCE,fo$f,type="b",pch=19,col="red")
points(fg$VARIANCE,fg$f,type="b",pch=19,col="blue")
grid()
dev.off()

```

Barplot is simpler.

```{r,bars}

prec <- c("CP default"=mean(cp$p),"CP BG fix"=mean(cpbg$p),"CP FDR fix"=mean(cpfdr$p),"CP BG and FDR fix"=mean(cpbgfdr$p),"fora"=mean(fo$p),"fgsea"=mean(fg$p))
rec <- c("CP default"=mean(cp$r),"CP BG fix"=mean(cpbg$r),"CP FDR fix"=mean(cpfdr$r),"CP BG and FDR fix"=mean(cpbgfdr$r),"fora"=mean(fo$r),"fgsea"=mean(fg$r))
f1 <- sapply(1:length(prec), function(i) { 2/(1/rec[i]+1/prec[i]) })

par(mfrow=c(1,3))
par(mar=c(c(9.1, 3.5, 2.1, 1.1) ))
barplot(prec,ylim=c(0,1.1),las=2,main="precision",ylab="")
text((1:6*1.18)-0.45,prec+0.02,labels=signif(prec,3))
barplot(rec,ylim=c(0,1.1),las=2,main="recall",ylab="")
text((1:6*1.18)-0.45,rec+0.02,labels=signif(rec,3))
barplot(f1,ylim=c(0,1.1),las=2,main="F1",ylab="")
text((1:6*1.18)-0.45,f1+0.02,labels=signif(f1,3))

png("fig3_bars.png", width=7,height=5,units="in",res=150,pointsize=12)
par(mar=c(c(9.1, 3.5, 2.1, 1.1) ))
par(mfrow=c(1,3))
barplot(prec,ylim=c(0,1.1),las=2,main="precision",ylab="")
text((1:6*1.18)-0.45,prec+0.02,labels=signif(prec,3))
barplot(rec,ylim=c(0,1.1),las=2,main="recall",ylab="")
text((1:6*1.18)-0.45,rec+0.02,labels=signif(rec,3))
barplot(f1,ylim=c(0,1.1),las=2,main="F1",ylab="")
text((1:6*1.18)-0.45,f1+0.02,labels=signif(f1,3))
dev.off()

#################################
pdf("fig3_bars.pdf", width=7,height=5,pointsize=12)
par(mar=c(c(9.1, 3.5, 2.1, 1.1) ))
par(mfrow=c(1,3))
barplot(prec,ylim=c(0,1.1),las=2,main="precision",ylab="")
text((1:6*1.18)-0.45,prec+0.02,labels=signif(prec,3))
barplot(rec,ylim=c(0,1.1),las=2,main="recall",ylab="")
text((1:6*1.18)-0.45,rec+0.02,labels=signif(rec,3))
barplot(f1,ylim=c(0,1.1),las=2,main="F1",ylab="")
text((1:6*1.18)-0.45,f1+0.02,labels=signif(f1,3))
dev.off()



```

## Session information


```{r,session}

sessionInfo()

```
