Examining background gene lists in clusterProfiler

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

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 how ClusterProfiler manages the background list, as we have observed some weird behaviour. In particular we see that when we provide a custom background, those genes do not appear in the final analysis unless they are also included in the gene list.

In the code chunk below called libs, you can add and remove required R library dependancies. Check that the libraries listed here match the Dockerfile, otherwise you might get errors.

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

## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationDbi'

For this guide I will be using bulk RNA-seq data from a previous study, which is deposited at NCBI GEO and SRA under accession numbers: GSE55123 and SRP038101 (Lund et al, 2014). The experiment is designed to investigate the effect of Azacitidine treatment on AML3 cells.

The raw data have been processed by the DEE2 project, and the summary gene expression counts are available at the dee2.io website, and programmatically with the getDEE2 bioconductor package (Ziemann et al, 2019).

Alternatively, you could fetch data from another resource like NCBI GEO, Zenodo or from the host storage drive.

mdat <- getDEE2Metadata("hsapiens")

# get sample sheet
ss <- subset(mdat,SRP_accession=="SRP038101")

# fetch the whole set of RNA-seq data
x <- getDEE2("hsapiens",ss$SRR_accession , metadata=mdat, legacy=TRUE)

## For more information about DEE2 QC metrics, visit
##     https://github.com/markziemann/dee2/blob/master/qc/qc_metrics.md

mx <- x$GeneCounts
rownames(mx) <- paste(rownames(mx),x$GeneInfo$GeneSymbol)
dim(mx)

## [1] 58302     6

# aza no filtering
ss$trt <- grepl("Treated",ss$Experiment_title)

ss %>%
  kbl(caption="sample sheet for Aza treatment in AML3 cells") %>%
  kable_paper("hover", full_width = F)

sample sheet for Aza treatment in AML3 cells
	SRR_accession	QC_summary	SRX_accession	SRS_accession	SRP_accession	Experiment_title	trt
156267	SRR1171523	PASS	SRX472607	SRS559064	SRP038101	GSM1329859: Untreated.1; Homo sapiens; RNA-Seq	FALSE
156268	SRR1171524	WARN(3,4)	SRX472608	SRS559066	SRP038101	GSM1329860: Untreated.2; Homo sapiens; RNA-Seq	FALSE
156269	SRR1171525	WARN(3,4)	SRX472609	SRS559065	SRP038101	GSM1329861: Untreated.3; Homo sapiens; RNA-Seq	FALSE
156270	SRR1171526	WARN(3,4)	SRX472610	SRS559068	SRP038101	GSM1329862: Treated.1; Homo sapiens; RNA-Seq	TRUE
156271	SRR1171527	WARN(3,4)	SRX472611	SRS559067	SRP038101	GSM1329863: Treated.2; Homo sapiens; RNA-Seq	TRUE
156272	SRR1171528	WARN(3,4)	SRX472612	SRS559069	SRP038101	GSM1329864: Treated.3; Homo sapiens; RNA-Seq	TRUE

Data quality control

QC is important, even if you are using public transcriptome data. For RNA-seq it is a good idea to show the number of reads for each sample.

par(mar=c(5,7,5,1))
barplot(rev(colSums(mx)),horiz=TRUE,las=1,main="number of reads per sample in SRP038101")

Now make a MDS plot.

mds <- cmdscale(dist(t(mx)))

# expand plot area
XMIN=min(mds[,1])*1.3
XMAX=max(mds[,1])*1.3
YMIN=min(mds[,2])*1.3
YMAX=max(mds[,2])*1.3

cols <- as.character(grepl("Treated",ss$Experiment_title))
cols <- gsub("FALSE","lightblue",cols)
cols <- gsub("TRUE","pink",cols)
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  xlim=c(XMIN,XMAX),ylim=c(YMIN,YMAX),
  pch=19,cex=2,col=cols, main="MDS plot")
text(cmdscale(dist(t(mx))), labels=colnames(mx) )

Differential expression analysis

I will be using DESeq2 for DE analysis, the count matrix is prefiltered using a detection threshold of 10 reads per sample across all samples. All genes that meet the detection threshold will comprise the background list. The first 6 rows of the count matrix are shows.

mxf <- mx[which(rowMeans(mx)>=10),]
dim(mxf)

## [1] 13168     6

head(mxf,6) %>%
  kbl(caption="Count matrix format") %>%
  kable_paper("hover", full_width = F)

Count matrix format
	SRR1171523	SRR1171524	SRR1171525	SRR1171526	SRR1171527	SRR1171528
ENSG00000225630 MTND2P28	494	396	340	333	415	418
ENSG00000237973 MTCO1P12	52	39	40	30	37	29
ENSG00000248527 MTATP6P1	853	544	537	582	702	716
ENSG00000228327 AL669831.1	17	13	21	21	22	12
ENSG00000228794 LINC01128	42	27	30	32	40	23
ENSG00000230699 AL645608.3	20	11	13	10	15	22

dds <- DESeqDataSetFromMatrix(countData = mxf , colData = ss, design = ~ trt )
res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <-cbind(as.data.frame(z),mxf)
def <-as.data.frame(zz[order(zz$pvalue),])

head(def,10) %>%
  kbl(caption="Top DE genes for Aza treatment") %>%
  kable_paper("hover", full_width = F)

Top DE genes for Aza treatment
	baseMean	log2FoldChange	lfcSE	stat	SRR1171523	SRR1171524	SRR1171525	SRR1171526	SRR1171527	SRR1171528
ENSG00000165949 IFI27	1960.1970	-3.384492	0.0938869	-36.04861	4689	3583	2758	309	384	334
ENSG00000090382 LYZ	7596.0299	-1.650342	0.0561143	-29.41036	14212	10237	10789	3476	3764	3738
ENSG00000115461 IGFBP5	531.2217	-5.071157	0.1795239	-28.24781	1320	823	1025	23	26	43
ENSG00000157601 MX1	827.1511	-2.877795	0.1047823	-27.46450	1732	1501	1206	184	223	186
ENSG00000111331 OAS3	2127.2010	-2.661214	0.0972124	-27.37525	4204	3977	2972	562	614	560
ENSG00000070915 SLC12A3	424.5509	-3.374852	0.1298671	-25.98697	1012	721	653	63	85	76
ENSG00000234745 HLA-B	3197.0159	-1.431566	0.0604169	-23.69479	6085	4256	4023	1590	1872	1719
ENSG00000137965 IFI44	409.0957	-2.978581	0.1319352	-22.57608	829	740	635	76	111	89
ENSG00000204525 HLA-C	1631.6421	-1.461550	0.0660214	-22.13750	3112	2150	2106	791	923	891
ENSG00000110042 DTX4	524.1318	-2.470219	0.1173182	-21.05572	1166	883	688	166	168	145

Make a smear plot.

sigf <- subset(def,padj<=0.05)

DET=nrow(mxf)
NSIG=nrow(sigf)
NUP=nrow(subset(sigf,log2FoldChange>0))
NDN=nrow(subset(sigf,log2FoldChange<0))

HEADER=paste(DET,"detected genes;",NSIG,"w/FDR<0.05;",NUP,"up;",NDN,"down")

plot(log10(def$baseMean) ,def$log2FoldChange,
  cex=0.6,pch=19,col="darkgray",
  main="5-azacitidine treatment in AML3 cells",
  xlab="log10(basemean)",ylab="log2(fold change)")

points(log10(sigf$baseMean) ,sigf$log2FoldChange,
  cex=0.6,pch=19,col="red")

mtext(HEADER)

In the next sections I will run enrichment analysis with over-representation analysis (ORA) test and compare it to functional class scoring. I will also investigate some strange behviour of the ORA tool clusterprofiler.

ORA with Clusterprofiler custom background with otherwise default analysis

I’ve compiled a reporting checklist:

Reporting criteria	Method/resource used
Origin of gene sets	Reactome (2023-03-06)
Tool used	ClusterProfiler (check version at foot of report)
Statistical test used	hypergeometric test
P-value correction for multiple comparisons	FDR method
Background list	Genes with >=10 reads per sample on average across all samples
Gene Selection Criteria	DESeq2 FDR<0.05
ORA directionality	Separate tests for up- and down-regulation
Data availability	via DEE2 at accession SRP038101 (human)
Other parameters	Min gene set size of 10

Here I provide a custom background gene list.

Get the gene sets loaded in R.

# from https://reactome.org/download/current/ReactomePathways.gmt.zip
genesets2 <- read.gmt("ReactomePathways_2023-03-06.gmt")
gsets <- gmtPathways("ReactomePathways_2023-03-06.gmt")

message(paste("number of genes described in the annotation set:",length(unique(genesets2$gene))))

## number of genes described in the annotation set: 11457

Now filter the gene names into three lists, up-regulated, down-regulated and background. Background is simply all genes that were detected.

defup <- rownames(subset(def,padj<=0.05 & log2FoldChange>0))
defup <- unique(sapply(strsplit(defup," "),"[[",1))

defdn <- rownames(subset(def,padj<=0.05 & log2FoldChange<0))
defdn <- unique(sapply(strsplit(defdn," "),"[[",1))

bg <- rownames(def)
bg <- unique(sapply(strsplit(bg," "),"[[",1))

# write ensembl IDs for DAVID 
writeLines(text=bg,con="ebg.txt")
writeLines(text=defup,con="eup.txt")
writeLines(text=defdn,con="edn.txt")

defup <- rownames(subset(def,padj<=0.05 & log2FoldChange>0))
defup <- unique(sapply(strsplit(defup," "),"[[",2))

defdn <- rownames(subset(def,padj<=0.05 & log2FoldChange<0))
defdn <- unique(sapply(strsplit(defdn," "),"[[",2))

bg <- rownames(def)
bg <- unique(sapply(strsplit(bg," "),"[[",2))

message(paste("number of genes in background:",length(bg)))

## number of genes in background: 13164

Enrichment firstly with upregulated genes

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

ora_up$geneID <- NULL
ora_up <- subset(ora_up,p.adjust<0.05 & Count >=10)
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

head(ora_up) %>%
  kbl(row.names = FALSE, caption="Top upregulated pathways in Aza treatment") %>%
  kable_paper("hover", full_width = F)

Top upregulated pathways in Aza treatment
ID	GeneRatio	BgRatio	p.adjust	qvalue	Count	es
Interactions of Rev with host cellular proteins	24/1050	36/6944	0e+00	0e+00	24	4.408889
Nuclear import of Rev protein	22/1050	33/6944	0e+00	0e+00	22	4.408889
Postmitotic nuclear pore complex (NPC) reformation	18/1050	27/6944	1e-07	1e-07	18	4.408889
Rev-mediated nuclear export of HIV RNA	22/1050	34/6944	0e+00	0e+00	22	4.279216
Export of Viral Ribonucleoproteins from Nucleus	20/1050	31/6944	0e+00	0e+00	20	4.266667
NEP/NS2 Interacts with the Cellular Export Machinery	20/1050	31/6944	0e+00	0e+00	20	4.266667

topup2 <- rev(head(ora_up$es,10))
names(topup2) <- rev(head(ora_up$ID,10))

The reported background size does not match the dataset.

Now repeat with the downregulated genes

n_pw=10

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

ora_dn$geneID <- NULL
ora_dn <- subset(ora_dn,p.adjust<0.05 & Count >=10)
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

head(ora_dn,n_pw) %>%
  kbl(row.names = FALSE, caption="Top downregulated pathways in Aza treatment") %>%
  kable_paper("hover", full_width = F)

Top downregulated pathways in Aza treatment
ID	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	Count	es
Interleukin-20 family signaling	10/1178	12/6944	0.0000009	0.0000594	0.0000530	10	4.912281
Interferon alpha/beta signaling	41/1178	55/6944	0.0000000	0.0000000	0.0000000	41	4.394258
Interferon gamma signaling	40/1178	63/6944	0.0000000	0.0000000	0.0000000	40	3.742690
Interleukin-10 signaling	17/1178	27/6944	0.0000001	0.0000104	0.0000093	17	3.711501
Other interleukin signaling	10/1178	18/6944	0.0002230	0.0080676	0.0072011	10	3.274854
Growth hormone receptor signaling	10/1178	19/6944	0.0003996	0.0140422	0.0125340	10	3.102493
Smooth Muscle Contraction	13/1178	26/6944	0.0001058	0.0040676	0.0036307	13	2.947368
FGFR1 mutant receptor activation	10/1178	20/6944	0.0006785	0.0213979	0.0190997	10	2.947368
RHO GTPases Activate NADPH Oxidases	10/1178	20/6944	0.0006785	0.0213979	0.0190997	10	2.947368
Interleukin-4 and Interleukin-13 signaling	34/1178	70/6944	0.0000000	0.0000001	0.0000001	34	2.863158

topdn2 <- head(ora_dn$es,n_pw)
names(topdn2) <- head(ora_dn$ID,n_pw)

Make a barplot

par(mar=c(5,20,5,1))

cols <- c(rep("blue",n_pw),rep("red",n_pw))

barplot(c(topdn2,topup2),
  horiz=TRUE,las=1,cex.names=0.65,col=cols,
  main="top DE Reactomes",
  xlab="ES",
  cex.main=0.9)

mtext("ORA test")

ORA with Clusterprofiler with custom background and special gene set list to fix the potential bug

I’ve compiled a reporting checklist:

Reporting criteria	Method/resource used
Origin of gene sets	Reactome (2023-03-06)
Tool used	ClusterProfiler (check version at foot of report)
Statistical test used	hypergeometric test
P-value correction for multiple comparisons	FDR method
Background list	Genes with >=10 reads per sample on average across all samples
Gene Selection Criteria	DESeq2 FDR<0.05
ORA directionality	Separate tests for up- and down-regulation
Data availability	via DEE2 at accession SRP038101 (human)
Other parameters	Min gene set size of 10

Here I provide a background gene list for cluterprofiler to use, like a user should.

I have also done some modification to the Reactome gene sets, to en

Now filter the gene names into three lists, up-regulated, down-regulated and background. Background is simply all genes that were detected.

defup <- rownames(subset(def,padj<=0.05 & log2FoldChange>0))
defup <- unique(sapply(strsplit(defup," "),"[[",2))

defdn <- rownames(subset(def,padj<=0.05 & log2FoldChange<0))
defdn <- unique(sapply(strsplit(defdn," "),"[[",2))

bg <- rownames(def)
bg <- unique(sapply(strsplit(bg," "),"[[",2))
message(paste("number of genes in background:",length(bg)))

## number of genes in background: 13164

Adding all detected genes to the background appears to improve results!

bgdf <- data.frame("background",bg)
colnames(bgdf) <- c("term","gene")
genesets2 <- rbind(genesets2,bgdf)

Enrichment firstly with upregulated genes

orafix_up <- as.data.frame(enricher(gene = defup ,
  universe = bg,  maxGSSize = 5000, TERM2GENE = genesets2,
  pAdjustMethod="fdr",  pvalueCutoff = 1, qvalueCutoff = 1  ))

orafix_up$geneID <- NULL
orafix_up <- subset(orafix_up,p.adjust<0.05 & Count >=10)
orafix_ups <- rownames(orafix_up)

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

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

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

head(orafix_up) %>%
  kbl(row.names = FALSE, caption="Top upregulated pathways in Aza treatment") %>%
  kable_paper("hover", full_width = F)

Top upregulated pathways in Aza treatment
ID	GeneRatio	BgRatio	Count	es
Interactions of Rev with host cellular proteins	24/1672	36/13164	24	5.248804
Nuclear import of Rev protein	22/1672	33/13164	22	5.248804
Postmitotic nuclear pore complex (NPC) reformation	18/1672	27/13164	18	5.248804
Rev-mediated nuclear export of HIV RNA	22/1672	34/13164	22	5.094427
Export of Viral Ribonucleoproteins from Nucleus	20/1672	31/13164	20	5.079488
NEP/NS2 Interacts with the Cellular Export Machinery	20/1672	31/13164	20	5.079488

topup2 <- rev(head(orafix_up$es,10))
names(topup2) <- rev(head(orafix_up$ID,10))

Now with the downregulated genes

n_pw=10

orafix_dn <- as.data.frame(enricher(gene = defdn ,
  universe = bg,  maxGSSize = 5000, TERM2GENE = genesets2,
  pAdjustMethod="fdr",  pvalueCutoff = 1, qvalueCutoff = 1  ))

orafix_dn$geneID <- NULL
orafix_dn <- subset(orafix_dn,p.adjust<0.05 & Count >=10)
orafix_dns <- rownames(orafix_dn)

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

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

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

head(orafix_dn,n_pw) %>%
  kbl(row.names = FALSE, caption="Top downregulated pathways in Aza treatment") %>%
  kable_paper("hover", full_width = F)

Top downregulated pathways in Aza treatment
ID	GeneRatio	BgRatio	pvalue	p.adjust	qvalue	Count	es
Interleukin-20 family signaling	10/1926	12/13164	0.0000002	0.0000108	0.0000089	10	5.695742
Interferon alpha/beta signaling	41/1926	55/13164	0.0000000	0.0000000	0.0000000	41	5.095101
Interferon gamma signaling	40/1926	63/13164	0.0000000	0.0000000	0.0000000	40	4.339613
Interleukin-10 signaling	17/1926	27/13164	0.0000000	0.0000009	0.0000007	17	4.303450
Other interleukin signaling	10/1926	18/13164	0.0000622	0.0021864	0.0018093	10	3.797162
Growth hormone receptor signaling	10/1926	19/13164	0.0001141	0.0036942	0.0030571	10	3.597311
FGFR1 mutant receptor activation	10/1926	20/13164	0.0001984	0.0055459	0.0045896	10	3.417445
RHO GTPases Activate NADPH Oxidases	10/1926	20/13164	0.0001984	0.0055459	0.0045896	10	3.417445
Smooth Muscle Contraction	13/1926	26/13164	0.0000216	0.0008068	0.0006677	13	3.417445
Interleukin-4 and Interleukin-13 signaling	34/1926	70/13164	0.0000000	0.0000000	0.0000000	34	3.319804

topdn2 <- head(orafix_dn$es,n_pw)
names(topdn2) <- head(orafix_dn$ID,n_pw)

Make a barplot

par(mar=c(5,20,5,1))

cols <- c(rep("blue",n_pw),rep("red",n_pw))

barplot(c(topdn2,topup2),
  horiz=TRUE,las=1,cex.names=0.65,col=cols,
  main="top DE Reactomes",
  xlab="ES",
  cex.main=0.9)

mtext("ORA fix test")

Enrichment analysis with functional class scoring

Reporting criteria	Method/resource used
Origin of gene sets	Reactome (2023-03-06)
Tool used	FGSEA (check version at foot of report)
Statistical test used	Kolmogorov-Smirnov test
P-value correction for multiple comparisons	FDR method
Background list	Genes with >=10 reads per sample on average across all samples
Gene Selection Criteria	DESeq2 FDR<0.05
ORA directionality	Separate tests for up- and down-regulation
Data availability	via DEE2 at accession SRP038101 (human)
Other parameters	Min gene set size of 10

def2 <- def
def2$genename <- sapply(strsplit(rownames(def2)," "),"[[",2)
def2 <- def2[,c("stat","genename")]
def2 <- aggregate(. ~ genename, def2, sum)
stat <- def2$stat
names(stat) <- def2$genename
fgseaRes <- fgsea(pathways=gsets, stats=stat, minSize=10, nPermSimple = 10000)
fgseaRes <- fgseaRes[order(fgseaRes$pval),]

fgsea_up <- subset(fgseaRes,padj<0.05 & ES>0)
fgsea_ups <- fgsea_up$pathway
fgsea_dn <- subset(fgseaRes,padj<0.05 & ES<0)
fgsea_dns <- fgsea_dn$pathway

fgsea_up <-  data.frame(fgsea_up[order(-fgsea_up$ES),])

fgsea_dn <-  data.frame(fgsea_dn[order(fgsea_dn$ES),])

head(fgsea_up,n_pw) %>%
  kbl(row.names = FALSE, caption="FGSEA:Top upregulated pathways in Aza treatment") %>%
  kable_paper("hover", full_width = F)

FGSEA:Top upregulated pathways in Aza treatment
pathway	pval	padj	log2err	ES	NES	size	leadingEdge
Condensation of Prometaphase Chromosomes	2.49e-05	0.0002378	0.5756103	0.8415269	2.138218	11	SMC4 , NCAPG , CSNK2A2, CDK1 , CCNB1 , SMC2 , NCAPH , CCNB2 , NCAPD2
Postmitotic nuclear pore complex (NPC) reformation	0.00e+00	0.0000000	0.8012156	0.7968889	2.601757	27	KPNB1 , NDC1 , AHCTF1 , RCC1 , NUP155 , TNPO1 , NUP160 , POM121 , NUP205 , NUP54 , NUP93 , NUP188 , NUP107 , NUP98 , RAN , NUP85 , SEH1L , NUP43 , SUMO1 , NUP133 , NUP62 , RANGAP1
Interactions of Rev with host cellular proteins	0.00e+00	0.0000000	0.8870750	0.7869546	2.762050	36	KPNB1 , NDC1 , NUP153 , NPM1 , RCC1 , NUP155 , TPR , RANBP1 , RANBP2 , NUP160 , POM121 , NUP50 , XPO1 , NUP205 , NUP54 , NUP93 , NUP188 , NUP214 , NUP107 , NUP98 , RAN , NUP85 , SEH1L , NUP43 , NUP133 , POM121C, NUP62 , RANGAP1, NUP88 , RAE1
Nuclear import of Rev protein	0.00e+00	0.0000000	0.8390889	0.7824653	2.704171	33	KPNB1 , NDC1 , NUP153 , NPM1 , RCC1 , NUP155 , TPR , RANBP2 , NUP160 , POM121 , NUP50 , NUP205 , NUP54 , NUP93 , NUP188 , NUP214 , NUP107 , NUP98 , RAN , NUP85 , SEH1L , NUP43 , NUP133 , POM121C, NUP62 , NUP88 , RAE1
Rev-mediated nuclear export of HIV RNA	0.00e+00	0.0000000	0.8513391	0.7821972	2.723299	34	NDC1 , NUP153 , RCC1 , NUP155 , TPR , RANBP1 , RANBP2 , NUP160 , POM121 , NUP50 , XPO1 , NUP205 , NUP54 , NUP93 , NUP188 , NUP214 , NUP107 , NUP98 , RAN , NUP85 , SEH1L , NUP43 , NUP133 , POM121C, NUP62 , RANGAP1, NUP88 , RAE1
SLBP Dependent Processing of Replication-Dependent Histone Pre-mRNAs	4.34e-04	0.0028733	0.4984931	0.7809511	1.984302	11	NCBP1 , LSM11 , SNRPD3, SNRPF , ZNF473, NCBP2 , SLBP , SNRPG , SNRPE
SLBP independent Processing of Histone Pre-mRNAs	5.12e-04	0.0032766	0.4772708	0.7801965	1.922218	10	NCBP1 , LSM11 , SNRPD3, SNRPF , ZNF473, NCBP2 , SNRPG , SNRPE
Transport of Ribonucleoproteins into the Host Nucleus	0.00e+00	0.0000000	0.8140358	0.7776911	2.642174	31	KPNB1 , NDC1 , NUP153 , NUP155 , TPR , RANBP2 , NUP160 , POM121 , NUP50 , NUP205 , NUP54 , NUP93 , NUP188 , KPNA1 , NUP214 , NUP107 , NUP98 , NUP85 , SEH1L , NUP43 , NUP133 , POM121C, NUP62 , NUP88 , RAE1
Impaired BRCA2 binding to PALB2	0.00e+00	0.0000008	0.7195128	0.7767850	2.452612	24	XRCC2 , BLM , BRIP1 , BRCA1 , RAD51AP1, RAD51D , WRN , BRCA2 , EXO1 , PALB2 , DNA2 , BARD1 , RAD51 , RBBP8 , MRE11 , RAD51C
Export of Viral Ribonucleoproteins from Nucleus	0.00e+00	0.0000000	0.8140358	0.7760569	2.636622	31	NDC1 , NUP153 , NUP155 , TPR , RANBP2 , NUP160 , POM121 , NUP50 , XPO1 , NUP205 , NUP54 , NUP93 , NUP188 , NUP214 , NUP107 , NUP98 , RAN , NUP85 , SEH1L , NUP43 , NUP133 , POM121C, NUP62 , NUP88 , RAE1

head(fgsea_dn,n_pw) %>%
  kbl(row.names = FALSE, caption="FGSEA:Top downregulated pathways in Aza treatment") %>%
  kable_paper("hover", full_width = F)

FGSEA:Top downregulated pathways in Aza treatment
pathway	pval	padj	log2err	ES	NES	size	leadingEdge
Endosomal/Vacuolar pathway	0.0000010	0.0000131	0.6435518	-0.8981315	-2.123361	12	HLA-B, HLA-C, B2M , HLA-A, HLA-E, CTSS , HLA-H, HLA-F, HLA-G
Interferon alpha/beta signaling	0.0000000	0.0000000	1.1330899	-0.8564838	-2.792262	55	IFI27 , MX1 , OAS3 , HLA-B , HLA-C , OAS2 , OAS1 , IFIT1 , IFITM3, IFIT3 , STAT1 , HLA-A , ISG15 , IRF7 , STAT2 , IFI6 , IFIT5 , IFIT2 , IFI35 , BST2 , IFITM1, IRF5 , RSAD2 , GBP2 , XAF1 , USP18 , HLA-E , ADAR , MX2 , HLA-H , OASL , PSMB8 , SAMHD1, HLA-F , HLA-G , TYK2 , PTPN6 , ISG20
TNFs bind their physiological receptors	0.0010794	0.0061733	0.4550599	-0.8166001	-1.840523	10	TNFSF13B, CD70 , TNFRSF1B, EDA2R , TNFSF13 , TNFSF9
Interleukin-20 family signaling	0.0003349	0.0022963	0.4984931	-0.8141392	-1.924787	12	STAT1 , STAT2 , STAT5A, JAK3 , JAK2 , TYK2 , IL10RB, STAT3 , STAT5B, STAT4
Regulation of IFNA/IFNB signaling	0.0018159	0.0096850	0.4550599	-0.8037272	-1.811509	10	STAT1, STAT2, USP18, TYK2 , PTPN6
Maturation of nucleoprotein_9683610	0.0009719	0.0056305	0.4772708	-0.7952852	-1.835324	11	PARP9 , PARP14, PARP10
Interferon gamma signaling	0.0000000	0.0000000	1.0073180	-0.7934868	-2.647956	63	OAS3 , HLA-B , HLA-C , OAS2 , OAS1 , B2M , STAT1 , HLA-A , IRF7 , TRIM22, IRF5 , GBP2 , SMAD7 , TRIM25, HLA-E , TRIM14, HLA-H , OASL , TRIM5 , GBP3 , PTAFR , TRIM8 , SP100 , HLA-F , TRIM21, PML , JAK2 , HLA-G , IFNGR2, GBP5 , PTPN6 , IFNGR1, MT2A , NCAM1 , GBP4 , PRKCD , TRIM26
MyD88 deficiency (TLR2/4)	0.0004203	0.0027967	0.4984931	-0.7914249	-1.904508	13	MYD88 , S100A9, S100A8, CD36 , TLR6 , CD14 , TLR2
RUNX3 regulates p14-ARF	0.0068087	0.0272115	0.2320568	-0.7641348	-1.722272	10	CDKN2A, CCND1 , RUNX3 , RUNX1
Transcriptional regulation of pluripotent stem cells	0.0019331	0.0102287	0.4550599	-0.7636281	-1.805368	12	EPAS1, EOMES, PBX1 , NR5A1, STAT3, KLF4 , HHEX

fgsea_up <- head(fgsea_up,n_pw)
fgsea_up_vec <- fgsea_up$ES
names(fgsea_up_vec) <- fgsea_up$pathway
fgsea_up_vec <- rev(fgsea_up_vec)

fgsea_dn <- head(fgsea_dn,n_pw)
fgsea_dn_vec <- fgsea_dn$ES * -1
names(fgsea_dn_vec) <- fgsea_dn$pathway

Make a barplot

par(mar=c(5,20,5,1))

cols <- c(rep("blue",n_pw),rep("red",n_pw))

barplot(c(fgsea_dn_vec,fgsea_up_vec),
  horiz=TRUE,las=1,cex.names=0.65,col=cols,
  main="top DE Reactomes",
  xlab="ES",
  cex.main=0.9)

mtext("FCS test")

Compare pathway results

v1 <- list("ORA up"=ora_ups, "ORA fix up"=orafix_ups, "FCS up"=fgsea_ups)
v2 <- list("ORA dn"=ora_dns, "ORA fix dn"=orafix_dns, "FCS dn"=fgsea_dns)
v3 <- list("ORA"=union(ora_ups,ora_dns),
  "ORA fix"=union(orafix_ups,orafix_dns),
  "FCS"=union(fgsea_ups,fgsea_dns))

par(mar=c(10,10,10,10))
par(mfrow=c(2,1))
plot(euler(v1),quantities = list(cex = 2), labels = list(cex = 2),main="upregulated genes")

plot(euler(v2),quantities = list(cex = 2), labels = list(cex = 2),main="downregulated genes")

plot(euler(v3),quantities = list(cex = 2), labels = list(cex = 2),main="up- and down-regulated genes")

Jaccard index comparing ORA and FCS.

jaccard <- function(a, b) {
    intersection = length(intersect(a, b))
    union = length(a) + length(b) - intersection
    return (intersection/union)
}

ora <- c(ora_ups,ora_dns)
fcs <- c(fgsea_ups,fgsea_dns)

jaccard(ora,fcs)

## [1] 0.441247

Jaccard index comparing ORA and ORA fix.

ora <- c(ora_ups,ora_dns)
orafix <- c(orafix_ups,orafix_dns)

jaccard(ora,orafix)

## [1] 0.6566667

Jaccard index comparing ORA fix and FCS.

orafix <- c(orafix_ups,orafix_dns)
fcs <- c(fgsea_ups,fgsea_dns)

jaccard(orafix,fcs)

## [1] 0.5438596

Drill in to a specific pathway

In this case the top upregulated and downregulated pathways

upname <- head(fgsea_up,1)$pathway
plotEnrichment(gsets[[upname]], stat, gseaParam = 1, ticksSize = 0.2)

dnname <- head(fgsea_dn,1)$pathway
plotEnrichment(gsets[[dnname]], stat, gseaParam = 1, ticksSize = 0.2)

Now a heatmap of these gene sets.

# reads per million normalisation
rpm <- apply(mxf,2, function(x) { x/sum(x) * 1000000} )
colnames(rpm) <- c("UNTR1","UNTR2","UNTR3","TRT1","TRT2","TRT3")
gnames_up <- gsets[[which(names(gsets) == upname)]]
gnames_dn <- gsets[[which(names(gsets) == dnname)]]
gene_ids <- rownames(mxf)
gene_names <- sapply(strsplit(gene_ids," "),"[[",2)
rpm_up <- rpm[which(gene_names %in% gnames_up),]
rownames(rpm_up) <- sapply(strsplit(rownames(rpm_up)," "),"[[",2)
rpm_dn <- rpm[which(gene_names %in% gnames_dn),]
rownames(rpm_dn) <- sapply(strsplit(rownames(rpm_dn)," "),"[[",2)
colsidecols <- c("blue","blue","blue","red","red","red")
heatmap.2(rpm_up,scale="row",trace="none",margins=c(10,15),main=upname,ColSideColors=colsidecols)

heatmap.2(rpm_dn,scale="row",trace="none",margins=c(10,15),main=dnname,ColSideColors=colsidecols)

Session information

For reproducibility

Click HERE to show session info

sessionInfo()

## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.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: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] gplots_3.1.3                eulerr_7.0.0               
##  [3] fgsea_1.22.0                clusterProfiler_4.4.4      
##  [5] kableExtra_1.3.4            DESeq2_1.36.0              
##  [7] SummarizedExperiment_1.26.1 Biobase_2.56.0             
##  [9] MatrixGenerics_1.8.1        matrixStats_0.63.0         
## [11] GenomicRanges_1.48.0        GenomeInfoDb_1.32.2        
## [13] IRanges_2.30.0              S4Vectors_0.34.0           
## [15] BiocGenerics_0.42.0         getDEE2_1.6.0              
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3     rstudioapi_0.14        jsonlite_1.8.4        
##   [4] magrittr_2.0.3         farver_2.1.1           rmarkdown_2.21        
##   [7] zlibbioc_1.42.0        vctrs_0.6.2            memoise_2.0.1         
##  [10] RCurl_1.98-1.12        ggtree_3.4.1           webshot_0.5.4         
##  [13] htmltools_0.5.5        gridGraphics_0.5-1     sass_0.4.5            
##  [16] KernSmooth_2.23-20     bslib_0.4.2            plyr_1.8.8            
##  [19] cachem_1.0.7           igraph_1.4.2           lifecycle_1.0.3       
##  [22] pkgconfig_2.0.3        Matrix_1.5-4           R6_2.5.1              
##  [25] fastmap_1.1.1          GenomeInfoDbData_1.2.8 digest_0.6.31         
##  [28] aplot_0.1.10           enrichplot_1.16.1      colorspace_2.1-0      
##  [31] patchwork_1.1.2        AnnotationDbi_1.58.0   geneplotter_1.74.0    
##  [34] RSQLite_2.3.1          labeling_0.4.2         fansi_1.0.4           
##  [37] httr_1.4.5             polyclip_1.10-4        compiler_4.3.0        
##  [40] bit64_4.0.5            withr_2.5.0            downloader_0.4        
##  [43] BiocParallel_1.30.3    viridis_0.6.2          DBI_1.1.3             
##  [46] highr_0.10             ggforce_0.4.1          MASS_7.3-59           
##  [49] DelayedArray_0.22.0    caTools_1.18.2         gtools_3.9.4          
##  [52] tools_4.3.0            DO.db_2.9              scatterpie_0.1.9      
##  [55] ape_5.7-1              glue_1.6.2             nlme_3.1-162          
##  [58] GOSemSim_2.22.0        polylabelr_0.2.0       shadowtext_0.1.2      
##  [61] grid_4.3.0             reshape2_1.4.4         generics_0.1.3        
##  [64] gtable_0.3.3           tidyr_1.3.0            data.table_1.14.8     
##  [67] tidygraph_1.2.3        xml2_1.3.4             utf8_1.2.3            
##  [70] XVector_0.36.0         ggrepel_0.9.3          pillar_1.9.0          
##  [73] stringr_1.5.0          yulab.utils_0.0.6      genefilter_1.78.0     
##  [76] splines_4.3.0          dplyr_1.1.2            tweenr_2.0.2          
##  [79] treeio_1.20.1          lattice_0.21-8         survival_3.5-5        
##  [82] bit_4.0.5              annotate_1.74.0        tidyselect_1.2.0      
##  [85] GO.db_3.15.0           locfit_1.5-9.7         Biostrings_2.64.0     
##  [88] knitr_1.42             gridExtra_2.3          svglite_2.1.1         
##  [91] xfun_0.39              graphlayouts_0.8.4     stringi_1.7.12        
##  [94] lazyeval_0.2.2         ggfun_0.0.9            yaml_2.3.7            
##  [97] evaluate_0.20          codetools_0.2-19       ggraph_2.1.0          
## [100] tibble_3.2.1           qvalue_2.28.0          ggplotify_0.1.0       
## [103] cli_3.6.1              xtable_1.8-4           systemfonts_1.0.4     
## [106] munsell_0.5.0          jquerylib_0.1.4        Rcpp_1.0.10           
## [109] png_0.1-8              XML_3.99-0.14          parallel_4.3.0        
## [112] ggplot2_3.4.2          blob_1.2.4             DOSE_3.22.0           
## [115] htm2txt_2.2.2          bitops_1.0-7           tidytree_0.4.2        
## [118] viridisLite_0.4.1      scales_1.2.1           purrr_1.0.1           
## [121] crayon_1.5.2           rlang_1.1.1            fastmatch_1.1-3       
## [124] KEGGREST_1.36.3        rvest_1.0.3

Examining background gene lists in clusterProfiler

Mark Ziemann & Anusuiya Bora

2023-06-21

Introduction

Data quality control

Differential expression analysis

ORA with Clusterprofiler custom background with otherwise default analysis

ORA with Clusterprofiler with custom background and special gene set list to fix the potential bug

Enrichment analysis with functional class scoring

Compare pathway results

Drill in to a specific pathway

Session information