---
title: "Examining FDR correction issue in clusterProfiler - RNA-seq example"
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: 7
theme: cosmo
---

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

## Introduction

This guide is a Rmarkdown script that conducts differential expression and enrichment analysis,
which are very popular workflows for transcriptome data.

The goal of this work is to understand the behaviour of ClusterProfiler's p-value adjustment method
as we have observed some weird behaviour.
In particular we see that clusterprofiler's criteria for inclusion/exclusion of gene sets is irregular.
For example clusterprofiler seems to include gene sets as detected if there is just 1 member detected
in the foreground list.
If a gene set has 10 in the background and  0 in the foreground it looks like it gets excluded.
This is a bit strange and contrary to best practice, where the detection threshold should be set by
the overlap size in the background list, rather than the foreground list.

Theoretically, this behaviour could have a negative impact on performance, as FDR is applied to the
wrong pathways.

```{r,libs}

suppressPackageStartupMessages({
  library("getDEE2")
  library("DESeq2")
  library("kableExtra")
  library("clusterProfiler")
  library("fgsea")
  library("eulerr")
  library("gplots")
  library("parallel")
})

```

For this guide I will be using bulk RNA-seq data from several previous studies.

## Load data

Gene sets were obtained from MSigDB.

```{r,load_data}

myfiles <- list.files("../dataprep",pattern="*Rds",full.names=TRUE)
d <- lapply(myfiles,readRDS)
names(d) <- basename(myfiles)

mygmts <- list.files("../gmt",pattern="*gmt",full.names=TRUE)
gs <- lapply(mygmts,gmtPathways)
names(gs) <- basename(mygmts)

names(gs)
names(gs) <- c("KEGG", "Reactome", "Wikipathways", "miR targets", "TFT GTRD",
  "GO", "HPO", "Cellmarkers", "Hallmark" )


```

## Basic geneset analysis

```{r,gs1}

# number of sets
lapply(gs,length)

# set size summary
lapply(gs,function(x) {
  y <- unlist(lapply(x, length ))
  summary(y)
} )

# gene coverage
lapply(gs,function(x) {
  y <- unique(unlist(x))
  length(y)
} )

# total number of annotations
lapply(gs, function(x) {
  sum(unlist(lapply(x,length)))
} )

```

## Enrichment with Clusterprofiler

Here is a function that conducts over-representation analysis with ClusterProfiler.
Genes that are not annotated as having any gene set category are discarded.

```{r,orabad}

# df = deseq output
# df must have "ensID geneSymbol" as rownames

# gs = gene sets
# gs format is list of symbols
orabad <- function( df, gs, ngenes=2000 ) {

  def <- df

  gs1 <- stack(gs)
  colnames(gs1) <- c("gene","term")
  gs1 <- gs1[,c(2,1)]

  defup <- head(rownames(subset(def,log2FoldChange>0)),ngenes)
  defup <- defup[lapply(strsplit(defup," "),length)==2]
  defup <- unique(sapply(strsplit(defup," "),"[[",2))

  defdn <- head(rownames(subset(def,log2FoldChange<0)),ngenes)
  defdn <- defdn[lapply(strsplit(defdn," "),length)==2]
  defdn <- unique(sapply(strsplit(defdn," "),"[[",2))

  bg <- rownames(def)
  bg <- bg[lapply(strsplit(bg," "),length)==2]
  bg <- unique(sapply(strsplit(bg," "),"[[",2))

  ora_up <- as.data.frame(enricher(gene = defup ,
    universe = bg, minGSSize = 5, maxGSSize = 500000, TERM2GENE = gs1,
    pAdjustMethod="fdr",  pvalueCutoff = 1, qvalueCutoff = 1  ))

  ora_up$geneID <- NULL
  ora_ups <- rownames(ora_up)

  gr <- as.numeric(sapply(strsplit(ora_up$GeneRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_up$GeneRatio,"/"),"[[",2))

  br <- as.numeric(sapply(strsplit(ora_up$BgRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_up$BgRatio,"/"),"[[",2))

  ora_up$es <- gr/br
  ora_up <- ora_up[order(-ora_up$es),]
  ora_up$Description=NULL

  ora_dn <- as.data.frame(enricher(gene = defdn ,
    universe = bg, minGSSize = 5, maxGSSize = 500000, TERM2GENE = gs1,
    pAdjustMethod="fdr",  pvalueCutoff = 1, qvalueCutoff = 1  ))

  ora_dn$geneID <- NULL
  ora_dns <- rownames(ora_dn)

  gr <- as.numeric(sapply(strsplit(ora_dn$GeneRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_dn$GeneRatio,"/"),"[[",2))

  br <- as.numeric(sapply(strsplit(ora_dn$BgRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_dn$BgRatio,"/"),"[[",2))

  ora_dn$es <- gr/br
  ora_dn <- ora_dn[order(-ora_dn$es),]
  ora_dn$Description=NULL

  result <- list("orabad_up"=ora_up,"orabad_dn"=ora_dn)
  return(result)
}

```

This function fixes the problem with not inlcuding sets for FDR correction.

```{r,oragood}

oragood <- function( df, gs, ngenes=2000 ) {

  def <- df

  gs1 <- stack(gs)
  colnames(gs1) <- c("gene","term")
  gs1 <- gs1[,c(2,1)]

  defup <- head(rownames(subset(def,log2FoldChange>0)),ngenes)
  defup <- defup[lapply(strsplit(defup," "),length)==2]
  defup <- unique(sapply(strsplit(defup," "),"[[",2))

  defdn <- head(rownames(subset(def,log2FoldChange<0)),ngenes)
  defdn <- defdn[lapply(strsplit(defdn," "),length)==2]
  defdn <- unique(sapply(strsplit(defdn," "),"[[",2))

  bg <- rownames(def)
  bg <- bg[lapply(strsplit(bg," "),length)==2]
  bg <- unique(sapply(strsplit(bg," "),"[[",2))

  gs1 <- gs1[which(gs1$gene %in% bg),]
  terms <- names(which(table(gs1$term)>5))
  gs1 <- gs1[gs1$term %in% terms,]

  ora_up <- as.data.frame(enricher(gene = defup ,
    universe = bg, minGSSize = 5, maxGSSize = 500000, TERM2GENE = gs1,
    pAdjustMethod="fdr",  pvalueCutoff = 1, qvalueCutoff = 1  ))

  ora_up$geneID <- NULL
  ora_ups <- rownames(ora_up)

  gr <- as.numeric(sapply(strsplit(ora_up$GeneRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_up$GeneRatio,"/"),"[[",2))

  br <- as.numeric(sapply(strsplit(ora_up$BgRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_up$BgRatio,"/"),"[[",2))

  ora_up$es <- gr/br
  ora_up <- ora_up[order(-ora_up$es),]
  ora_up$Description=NULL

  nsets <- length(which(table(gs1$term)>5))
  nres <- nrow(ora_up)
  diff <- nsets - nres
  pvals <- c(ora_up$pvalue,rep(1,diff))
  ora_up$p.adjust <- p.adjust(pvals,method="fdr")[1:nrow(ora_up)]

  ora_dn <- as.data.frame(enricher(gene = defdn ,
    universe = bg, minGSSize = 5, maxGSSize = 500000, TERM2GENE = gs1,
    pAdjustMethod="fdr",  pvalueCutoff = 1, qvalueCutoff = 1  ))

  ora_dn$geneID <- NULL
  ora_dns <- rownames(ora_dn)

  gr <- as.numeric(sapply(strsplit(ora_dn$GeneRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_dn$GeneRatio,"/"),"[[",2))

  br <- as.numeric(sapply(strsplit(ora_dn$BgRatio,"/"),"[[",1)) /
    as.numeric(sapply(strsplit(ora_dn$BgRatio,"/"),"[[",2))

  ora_dn$es <- gr/br
  ora_dn <- ora_dn[order(-ora_dn$es),]
  ora_dn$Description=NULL

  nsets <- length(which(table(gs1$term)>5))
  nres <- nrow(ora_dn)
  diff <- nsets - nres
  pvals <- c(ora_dn$pvalue,rep(1,diff))
  ora_dn$p.adjust <- p.adjust(pvals,method="fdr")[1:nrow(ora_dn)]

  result <- list("oragood_up"=ora_up,"oragood_dn"=ora_dn)
  return(result)
}


```

## Compare the results of bad and good analysis

```{r,compare}

# bad=ORA analysis, a list with results of up and down regulated pathways
# good=same as bad, except fixed the background problem
compare_ora_plots <- function(bad,good,libname=NULL,setsize=NULL){
  badup <- bad$orabad_up
  baddn <- bad$orabad_dn
  goodup <- good$oragood_up
  gooddn <- good$oragood_dn
  mup <- merge(badup,goodup,by=0)

  MAX=max(-log10(c(mup$p.adjust.x,mup$p.adjust.y)))
  UPHEADER=paste(libname,"up",setsize,"genes")
  DNHEADER=paste(libname,"dn",setsize,"genes")

  par(mfrow=c(1,2))

  plot(-log10(mup$p.adjust.x), -log10(mup$p.adjust.y),
    xlim=c(0,MAX), ylim=c(0,MAX),
    xlab="original -log10 FDR", ylab="corrected -log10 FDR", main=UPHEADER)
  grid()
  abline(a = 0, b = 1,lwd=2,lty=2)

  mdn <- merge(baddn,gooddn,by=0)

  MAX=max(-log10(c(mdn$p.adjust.x,mdn$p.adjust.y)))

  plot(-log10(mdn$p.adjust.x), -log10(mdn$p.adjust.y),
    xlim=c(0,MAX), ylim=c(0,MAX),
    xlab="original -log10 FDR", ylab="corrected -log10 FDR", main=DNHEADER)
  grid()
  abline(a = 0, b = 1,lwd=2,lty=2)
}

compare_ora <- function(bad,good){
  badup <- bad$orabad_up
  baddn <- bad$orabad_dn
  goodup <- good$oragood_up
  gooddn <- good$oragood_dn
  bup <- rownames(subset(badup,p.adjust<0.05))
  gup <- rownames(subset(goodup,p.adjust<0.05))
  bdn <- rownames(subset(baddn,p.adjust<0.05))
  gdn <- rownames(subset(gooddn,p.adjust<0.05))
  sig_bad <- unique(c(bup,bdn))
  sig_good <- unique(c(gup,gdn))
  nsig_bad <- length(sig_bad)
  nsig_good <- length(sig_good)
  jac <- length(intersect(sig_bad,sig_good)) / length(union(sig_bad,sig_good))
  result=c("nsig_bad"=nsig_bad,"nsig_good"=nsig_good,"Jaccard"=jac)
  return(result)
}

```

## Run some analyses

```{r,run1}

params <- expand.grid(1:length(d),1:length(gs))
params2 <- expand.grid(names(d),names(gs))

colnames(params) <- c("d","gs")

ora_dat <- mclapply(1:nrow(params), function(i) {
  j=params[i,1]
  k=params[i,2]
  bad <- orabad(d[[j]],gs[[k]])
  good <- oragood(d[[j]],gs[[k]])
  res <- compare_ora(bad,good)
  return(res)
} , mc.cores=8)

ora_dat <- do.call(rbind,ora_dat)

res1 <- cbind(params2,ora_dat)

res1 %>%
  kbl(row.names = TRUE, caption="all results") %>%
  kable_paper("hover", full_width = F)

summary(res1$nsig_bad)
summary(res1$nsig_good)
summary(res1$Jaccard)

# NUMBER OF SIGNIFICANT SETS
ns1 <- matrix( res1$nsig_bad , ncol=9 )
rownames(ns1) <- names(d)
colnames(ns1) <- names(gs)

ns1 %>%
  kbl(row.names = TRUE, caption="significant genesets in original analysis") %>%
  kable_paper("hover", full_width = F)

ns1 <- apply(ns1,2,function(x) { x <- x[!is.na(x)] ; x <- x[is.finite(x)] ;  mean(x) } )

ns2 <- matrix( res1$nsig_good , ncol=9 )
rownames(ns2) <- names(d)
colnames(ns2) <- names(gs)

ns2 %>%
  kbl(row.names = TRUE, caption="significant genesets in corrected analysis") %>%
  kable_paper("hover", full_width = F)

ns2 <- apply(ns2,2,function(x) { x <- x[!is.na(x)] ; x <- x[is.finite(x)] ;  mean(x) } )
ns2

df <- data.frame(ns1,ns2)
df <- df[order(df$ns1),]

par(mar=c(c(5.1, 7.1, 2.1, 2.1) ))
barplot(t(df),beside=TRUE,horiz=TRUE,las=1,xlim=c(0,330),
  xlab="no. significant sets",
  col=c("gray30","gray80"))
text(df$ns1+20,((1:nrow(df))*3)-1.5,labels=signif(df[,1],3))
text(df$ns2+20,((1:nrow(df))*3)-0.5,labels=signif(df[,2],3))
legend("bottomright", inset=.02, legend=c("corrected","original"),
  fill=c("gray80","gray30"), horiz=FALSE, cex=1)

png("fig2_ORAfdr.png")
par(mar=c(c(5.1, 7.1, 2.1, 2.1) ))
barplot(t(df),beside=TRUE,horiz=TRUE,las=1,xlim=c(0,330),
  xlab="no. significant sets",
  col=c("gray30","gray80"))
text(df$ns1+20,((1:nrow(df))*3)-1.5,labels=signif(df[,1],3))
text(df$ns2+20,((1:nrow(df))*3)-0.5,labels=signif(df[,2],3))
legend("bottomright", inset=.02, legend=c("corrected","original"),
  fill=c("gray80","gray30"), horiz=FALSE, cex=1)
dev.off()

pdf("fig2_ORAfdr.pdf",width=7,height=7)
par(mar=c(c(5.1, 7.1, 2.1, 2.1) ))
barplot(t(df),beside=TRUE,horiz=TRUE,las=1,xlim=c(0,330),
  xlab="no. significant sets",
  col=c("gray30","gray80"))
text(df$ns1+20,((1:nrow(df))*3)-1.5,labels=signif(df[,1],3))
text(df$ns2+20,((1:nrow(df))*3)-0.5,labels=signif(df[,2],3))
legend("bottomright", inset=.02, legend=c("corrected","original"),
  fill=c("gray80","gray30"), horiz=FALSE, cex=1)
dev.off()

# RATIO
rat <- matrix( (res1$nsig_good+0) / (res1$nsig_bad+0) , ncol=9 )
rownames(rat) <- names(d)
colnames(rat) <- names(gs)

signif(rat,3) %>%
  kbl(row.names = TRUE, caption="Ratio of significant genesets Corrected:Original") %>%
  kable_paper("hover", full_width = F)

rat2 <- apply(rat,2,function(x) { x <- x[!is.na(x)] ; x <- x[is.finite(x)] ;  mean(x) } )
rat2

# JACCARD
jac <- matrix( res1$Jaccard , ncol=9 )
rownames(jac) <- names(d)
colnames(jac) <- names(gs)

signif(jac,3) %>%
  kbl(row.names = TRUE, caption="Jaccard similarity Corrected:Original") %>%
  kable_paper("hover", full_width = F)

jac2 <- apply(jac,2,function(x) { x <- x[!is.na(x)] ; x <- x[is.finite(x)] ;  mean(x) } )
jac2

barplot(sort(jac2),horiz=TRUE,las=1,xlim=c(0,1),
  main="Jaccard similarity Corrected:Original",xlab="Jaccard index")

# overall
oa <- matrix(paste(res1$nsig_bad,res1$nsig_good,signif(res1$Jaccard,2)),ncol=9)
rownames(oa) <- names(d)
colnames(oa) <- names(gs)

oa %>%
  kbl(row.names = TRUE, caption="Original:Corrected:Jaccard") %>%
  kable_paper("hover", full_width = F)

```

Now run an exhaustive analysis, changing the gene list length from 125 to 2000.

```{r,run2}

ngenes <- c(125,250,500,1000,2000)
params <- expand.grid(1:length(d),1:length(gs),ngenes)
params2 <- expand.grid(names(d),names(gs),ngenes)

colnames(params) <- c("d","gs","ngenes")

ora_dat <- mclapply(1:nrow(params), function(i) {
  j=params[i,1]
  k=params[i,2]
  l=params[i,3]
  bad <- orabad(d[[j]],gs[[k]],l)
  good <- oragood(d[[j]],gs[[k]],l)
  res <- compare_ora(bad,good)
  return(res)
} , mc.cores=8)

ora_dat <- do.call(rbind,ora_dat)

res1 <- cbind(params2,ora_dat)

res1 %>%
  kbl(row.names = TRUE, caption="all results") %>%
  kable_paper("hover", full_width = F)

summary(res1$nsig_bad)
summary(res1$nsig_good)
summary(res1$Jaccard)

lapply(ngenes, function(n) {
  res2 <- subset(res1,Var3==n) ; summary(res2$nsig_bad)
} )

lapply(ngenes, function(n) {
  res2 <- subset(res1,Var3==n) ; summary(res2$nsig_good)
} )

lapply(ngenes, function(n) {
  res2 <- subset(res1,Var3==n) ; summary(res2$Jaccard)
} )

# NUMBER OF SIGNIFICANT SETS

ns1 <- lapply(ngenes ,function(n) {
  res2 <- subset(res1,Var3==n)
  ns <- matrix( res2$nsig_bad , ncol=9 )
  rownames(ns) <- names(d)
  colnames(ns) <- names(gs)
  return(ns)
} )

names(ns1) <- ngenes

ns1

ns1 <- lapply(ns1,function(y) {
  res <- apply(y,2,function(x) {
    x <- x[!is.na(x)]
    x <- x[is.finite(x)]
    mean(x)
  } )
  return(res)
})

ns1

ns2 <- lapply(ngenes ,function(n) {
  res2 <- subset(res1,Var3==n)
  ns <- matrix( res2$nsig_good , ncol=9 )
  rownames(ns) <- names(d)
  colnames(ns) <- names(gs)
  return(ns)
} )

names(ns2) <- ngenes

ns2

ns2 <- lapply(ns2,function(y) {
  res <- apply(y,2,function(x) {
    x <- x[!is.na(x)]
    x <- x[is.finite(x)]
    mean(x)
  } )
  return(res)
})

ns2

# JACCARD
jac <- lapply(ngenes ,function(n) {
  res2 <- subset(res1,Var3==n)
  ns <- matrix( res2$Jaccard , ncol=9 )
  rownames(ns) <- names(d)
  colnames(ns) <- names(gs)
  return(ns)
} )
names(jac) <- ngenes

jac

jac <- lapply(jac,function(y) {
  res <- apply(y,2,function(x) {
    x <- x[!is.na(x)]
    x <- x[is.finite(x)]
    mean(x)
  } )
  return(res)
})

jac

jac <- do.call(rbind,jac)

plot(jac[,"GO"],type="b", pch=19,cex=1.5,lwd=2,
  xaxt = "n",ylim=c(0,1),ylab="Jaccard index",xlab="no. genes selected")
axis(1, at = seq(1, 5, by = 1),
     labels = ngenes)
points(jac[,"Reactome"],type="b",col="red",pch=19,cex=1.5,lwd=2)
points(jac[,"Cellmarkers"],type="b",col="blue",pch=19,cex=1.5,lwd=2)
points(jac[,"TFT GTRD"],type="b",col="gray50",pch=19,cex=1.5,lwd=2)

legend("bottomright", legend=c("GO","Reactome", "Cellmarkers", "TFT"),
       col=c("black", "red", "blue", "gray50"), lty=1, cex=1.2,lwd=3,pch=19)

png("fig2_ORAsz.png")
plot(jac[,"GO"],type="b", pch=19,cex=1.5,lwd=2,
  xaxt = "n",ylim=c(0,1),ylab="Jaccard index",xlab="no. genes selected")
axis(1, at = seq(1, 5, by = 1),
     labels = ngenes)
points(jac[,"Reactome"],type="b",col="red",pch=19,cex=1.5,lwd=2)
points(jac[,"Cellmarkers"],type="b",col="blue",pch=19,cex=1.5,lwd=2)
points(jac[,"TFT GTRD"],type="b",col="gray50",pch=19,cex=1.5,lwd=2)

legend("bottomright", legend=c("GO","Reactome", "Cellmarkers", "TFT"),
       col=c("black", "red", "blue", "gray50"), lty=1, cex=1.2,lwd=3,pch=19)
dev.off()

pdf("fig2_ORAsz.pdf",width=7, height=7)
plot(jac[,"GO"],type="b", pch=19,cex=1.5,lwd=2,
  xaxt = "n",ylim=c(0,1),ylab="Jaccard index",xlab="no. genes selected")
axis(1, at = seq(1, 5, by = 1),
     labels = ngenes)
points(jac[,"Reactome"],type="b",col="red",pch=19,cex=1.5,lwd=2)
points(jac[,"Cellmarkers"],type="b",col="blue",pch=19,cex=1.5,lwd=2)
points(jac[,"TFT GTRD"],type="b",col="gray50",pch=19,cex=1.5,lwd=2)

legend("bottomright", legend=c("GO","Reactome", "Cellmarkers", "TFT"),
       col=c("black", "red", "blue", "gray50"), lty=1, cex=1.2,lwd=3,pch=19)
dev.off()

```

## Make charts showing changes in pvalues

For bulk RNA-seq data set #1, run enrichment with all the gene set libraries.

```{r,charts1,fig.height=4}

names(gs)

lapply(1:length(gs), function(i) {
  bad <- orabad(d[[1]],gs[[i]],ngenes=125)
  good <- oragood(d[[1]],gs[[i]],ngenes=125)
  compare_ora_plots(bad,good,libname=names(gs)[i],setsize=125)
} )

lapply(1:length(gs), function(i) {
  bad <- orabad(d[[1]],gs[[i]],ngenes=250)
  good <- oragood(d[[1]],gs[[i]],ngenes=250)
  compare_ora_plots(bad,good,libname=names(gs)[i],setsize=250)
} )

lapply(1:length(gs), function(i) {
  bad <- orabad(d[[1]],gs[[i]],ngenes=500)
  good <- oragood(d[[1]],gs[[i]],ngenes=500)
  compare_ora_plots(bad,good,libname=names(gs)[i],setsize=500)
} )

lapply(1:length(gs), function(i) {
  bad <- orabad(d[[1]],gs[[i]],ngenes=1000)
  good <- oragood(d[[1]],gs[[i]],ngenes=1000)
  compare_ora_plots(bad,good,libname=names(gs)[i],setsize=1000)
} )

lapply(1:length(gs), function(i) {
  bad <- orabad(d[[1]],gs[[i]],ngenes=2000)
  good <- oragood(d[[1]],gs[[i]],ngenes=2000)
  compare_ora_plots(bad,good,libname=names(gs)[i],setsize=2000)
} )

```

## Session information

For reproducibility

```{r,session}

sessionInfo()

```

