Example gene set analysis

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

Intro

Here we are performing an analysis of some gene expression data to demonstrate the difference between ORA and FCS methods and to highlight the differences caused by improper background gene set use.

The dataset being used is SRP128998 and we are comparing the cells grown in normal glucose condition (control) to the high glucose condition (case).

Data are obtained from http://dee2.io/

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

dir.create("images")

## Warning in dir.create("images"): 'images' already exists

Get expression data and make an MDS plot

I’m using some RNA-seq data looking at the effect of hyperglycemia on hepatocytes.

name="SRP128998"
mdat<-getDEE2Metadata("hsapiens")
samplesheet <- mdat[grep("SRP128998",mdat$SRP_accession),]
samplesheet<-samplesheet[order(samplesheet$SRR_accession),]
samplesheet$trt<-as.factor(c(1,1,1,1,1,1,0,0,0,0,0,0))
samplesheet$VPA<-as.factor(c(0,0,0,1,1,1,0,0,0,1,1,1))
s1 <- subset(samplesheet,VPA==0)

s1 %>% kbl(caption = "sample sheet") %>% kable_paper("hover", full_width = F)

sample sheet
	SRR_accession	QC_summary	SRX_accession	SRS_accession	SRP_accession	Sample_name	GEO_series	trt
475895	SRR6467479	PASS	SRX3557428	SRS2830728	SRP128998	GSM2932791	GSE109140	1
475896	SRR6467480	PASS	SRX3557429	SRS2830730	SRP128998	GSM2932792	GSE109140	1
475897	SRR6467481	PASS	SRX3557430	SRS2830729	SRP128998	GSM2932793	GSE109140	1
475901	SRR6467485	PASS	SRX3557434	SRS2830733	SRP128998	GSM2932797	GSE109140	0
475902	SRR6467486	PASS	SRX3557435	SRS2830734	SRP128998	GSM2932798	GSE109140	0
475903	SRR6467487	PASS	SRX3557436	SRS2830735	SRP128998	GSM2932799	GSE109140	0

w<-getDEE2("hsapiens",samplesheet$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

x<-Tx2Gene(w)
x<-x$Tx2Gene

# save the genetable for later
gt<-w$GeneInfo[,1,drop=FALSE]
gt$accession<-rownames(gt)

# counts 
x1<-x[,which(colnames(x) %in% s1$SRR_accession)]

Here show the number of genes in the annotation set, and those detected above the detection threshold.

# filter out lowly expressed genes
x1<-x1[which(rowSums(x1)/ncol(x1)>=(10)),]
nrow(x)

## [1] 39297

nrow(x1)

## [1] 15635

Now multidimensional scaling (MDS) plot to show the correlation between the datasets. If the control and case datasets are clustered separately, then it is likely that there will be many differentially expressed genes with FDR<0.05.

plot(cmdscale(dist(t(x1))), xlab="Coordinate 1", ylab="Coordinate 2", pch=19, col=s1$trt, main="MDS")

Differential expression

Now run DESeq2 for control vs case.

y <- DESeqDataSetFromMatrix(countData = round(x1), colData = s1, design = ~ trt)

## converting counts to integer mode

y <- DESeq(y)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

de <- results(y)
de<-as.data.frame(de[order(de$pvalue),])
rownames(de)<-sapply(strsplit(rownames(de),"\\."),"[[",1)
head(de) %>% kbl() %>% kable_paper("hover", full_width = F)

	baseMean	log2FoldChange	lfcSE	stat
ENSG00000145050	5839.731	-2.753692	0.1518338	-18.13623
ENSG00000149131	1346.633	2.161115	0.1427218	15.14215
ENSG00000044574	124889.027	-2.033391	0.1343040	-15.14021
ENSG00000128228	1676.368	-2.836358	0.1895729	-14.96183
ENSG00000179218	78785.663	-2.227516	0.1586383	-14.04148
ENSG00000090520	6751.044	-2.138112	0.1538129	-13.90073

Now let’s have a look at some of the charts showing differential expression. In particular, an MA plot and volcano plot.

maplot <- function(de,contrast_name) {
  sig <-subset(de, padj < 0.05 )
  up <-rownames(subset(de, padj < 0.05 & log2FoldChange > 0))
  dn <-rownames(subset(de, padj < 0.05 & log2FoldChange < 0))
  GENESUP <- length(up)
  GENESDN <- length(dn)
  DET=nrow(de)
  SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down", DET, "detected")
  ns <-subset(de, padj > 0.05 )
  plot(log2(de$baseMean),de$log2FoldChange, 
       xlab="log2 basemean", ylab="log2 foldchange",
       pch=19, cex=0.5, col="dark gray",
       main=contrast_name, cex.main=0.7)
  points(log2(sig$baseMean),sig$log2FoldChange,
         pch=19, cex=0.5, col="red")
  mtext(SUBHEADER,cex = 0.7)
}

make_volcano <- function(de,name) {
    sig <- subset(de,padj<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,log2FoldChange>0))
    N_DN=nrow(subset(sig,log2FoldChange<0))
    DET=nrow(de)
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
    plot(de$log2FoldChange,-log10(de$padj),cex=0.5,pch=19,col="darkgray",
        main=name, xlab="log2 FC", ylab="-log10 pval", xlim=c(-6,6))
    mtext(HEADER)
    grid()
    points(sig$log2FoldChange,-log10(sig$padj),cex=0.5,pch=19,col="red")
}

maplot(de,name)

make_volcano(de,name)

Gene sets from Reactome

In order to perform gene set analysis, we need some gene sets.

if (! file.exists("ReactomePathways.gmt")) {
  download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", 
    destfile="ReactomePathways.gmt.zip")
  unzip("ReactomePathways.gmt.zip")
}
genesets<-gmt_import("ReactomePathways.gmt")

FCS with Mitch

Mitch uses rank-ANOVA statistics for enrichment detection.

m <- mitch_import(de,DEtype = "DEseq2", geneTable = gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 15635

## Note: no. genes in output = 14533

## Note: estimated proportion of input genes in output = 0.93

mres <- mitch_calc(m,genesets = genesets)

## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.

m_up <- subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist > 0)[,1]
m_dn <- subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist < 0)[,1]
message(paste("Number of up-regulated pathways:",length(m_up) ))

## Number of up-regulated pathways: 95

message(paste("Number of down-regulated pathways:",length(m_dn) ))

## Number of down-regulated pathways: 316

head(mres$enrichment_result,10)  %>% kbl() %>% kable_paper("hover", full_width = F)

	set	setSize	s.dist
154	Cellular responses to stimuli	642	-0.2205691
155	Cellular responses to stress	634	-0.2150514
506	Infectious disease	632	-0.2030281
1035	SRP-dependent cotranslational protein targeting to membrane	109	-0.4728714
74	Asparagine N-linked glycosylation	252	-0.3060943
86	Axon guidance	424	-0.2324048
627	Metabolism of proteins	1550	-0.1176239
705	Nervous system development	444	-0.2086060
153	Cellular response to starvation	150	-0.3551396
570	L13a-mediated translational silencing of Ceruloplasmin expression	108	-0.4039887

ORA with clusterprofiler

Clusterprofiler uses a hypergeometric test. Firstly I will conduct the analysis separately for up and down regulated genes and with the correct backgound (as intended by the developers).

genesets2 <- read.gmt("ReactomePathways.gmt")

de_up <- rownames(subset(de,log2FoldChange>0,padj<0.05))
de_up <- unique(gt[which(rownames(gt) %in% de_up),1])

de_dn <- rownames(subset(de,log2FoldChange<0,padj<0.05))
de_dn <- unique(gt[which(rownames(gt) %in% de_dn),1])

de_bg <- rownames(de)
de_bg <- unique(gt[which(rownames(gt) %in% de_bg),1])

c_up <- as.data.frame(enricher(gene = de_up, universe = de_bg,  maxGSSize = 5000, TERM2GENE = genesets2))
c_up <- rownames(subset(c_up, p.adjust < 0.05))
       
c_dn <- as.data.frame(enricher(gene = de_dn, universe = de_bg,  maxGSSize = 5000, TERM2GENE = genesets2))
c_dn <- rownames(subset(c_dn, p.adjust < 0.05))

Now performing ORA with clusterprofiler combining up and down.

de_de <- rownames(subset(de,padj<0.05))
de_de <- unique(gt[which(rownames(gt) %in% de_de),1])

d_de <- as.data.frame(enricher(gene = de_de, universe = de_bg,  maxGSSize = 5000, TERM2GENE = genesets2))
d_de <- rownames(subset(d_de, p.adjust < 0.05))

Now performing ORA with clusterprofiler with whole genome background list

de_bg <- w$GeneInfo$GeneSymbol

f_up <- as.data.frame(enricher(gene = de_up, universe = de_bg,  maxGSSize = 5000, TERM2GENE = genesets2))
f_up <- rownames(subset(f_up, p.adjust < 0.05))
       
f_dn <- as.data.frame(enricher(gene = de_dn, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2))
f_dn <- rownames(subset(f_dn, p.adjust < 0.05))

Now performing ORA (combining up and down gene lists) with clusterprofiler with whole genome background list

e_de <- as.data.frame(enricher(gene = de_de, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2))
e_de <- rownames(subset(e_de, p.adjust < 0.05))

Venn diagram comparison

The Venn (or Euler to be more correct) diagram is useful to visualise the overlaps between sets.

par(cex.main=0.5)

par(mar=c(2,2,2,2))

v0 <- list("ORA up"=c_up,"ORA dn"=c_dn,
           "ORA comb" = d_de)

plot(euler(v0),quantities = TRUE, edges = "gray", main="effect of combining up and down regulated genes")

v1 <- list("FCS up"=m_up, "FCS dn"=m_dn,
           "ORA up"=c_up,"ORA dn"=c_dn)
  
plot(euler(v1),quantities = TRUE, edges = "gray", main="FCS compared to ORA")

v2 <- list("ORA up"=c_up,"ORA dn"=c_dn, 
           "ORA* up"=f_up,"ORA* dn"=f_dn )

plot(euler(v2),quantities = TRUE, edges = "gray", main="Effect of inappropriate background* (whole genome)")

vx <- list("ORA up"=c_up,"ORA dn"=c_dn,
           "ORA comb" = d_de, "ORA* comb" = e_de)

plot(euler(vx),quantities = TRUE, edges = "gray", main="combining up and down genes and whole genome bg*")

v3 <- list("ORA up"=c_up,"ORA dn"=c_dn, 
           "ORA* up"=f_up,"ORA* dn"=f_dn ,
           "FCS up"=m_up, "FCS dn"=m_dn)

png("images/fcs_ora1.png")
plot(euler(v1),quantities = TRUE, edges = "gray", main="FCS vs ORA")
dev.off()

## png 
##   2

png("images/orabg1.png")
plot(euler(v2),quantities = TRUE, edges = "gray", main="Effect of inappropriate background* (whole genome)")
dev.off()

## png 
##   2

png("images/oracomb1.png")
plot(euler(vx),quantities = TRUE, main="combining up and down genes and whole genome bg*")
dev.off()

## png 
##   2

pdf("images/fcs_ora1.pdf",width=4,height=4)
plot(euler(v1),quantities = TRUE, edges = "gray", main="FCS vs ORA")
dev.off()

## png 
##   2

pdf("images/orabg1.pdf",width=4,height=4)
plot(euler(v2),quantities = TRUE, edges = "gray", main="Effect of inappropriate background* (whole genome)")
dev.off()

## png 
##   2

pdf("images/oracomb1.pdf",width=4,height=4)
plot(euler(vx),quantities = TRUE, edges = "gray", main="combining up and down genes and whole genome bg*")
dev.off()

## png 
##   2

Jaccard calculation

# ORA vs ORA combined
dc <- length(intersect(d_de, c(c_up,c_dn))) / length(union(d_de, c(c_up,c_dn)))

# ORA vs ORA* combined
ec <- length(intersect(e_de, c(c_up,c_dn))) / length(union(e_de, c(c_up,c_dn)))

# FCS vs ORA
cm <- length(intersect(c(c_up,c_dn), c(m_up,m_dn))) / length(union(c(c_up,c_dn), c(m_up,m_dn)))

m_up <- gsub("^","up ",m_up)
m_dn <- gsub("^","dn ",m_dn)
m_de <- union(m_up,m_dn)

c_up <- gsub("^","up ",c_up)
c_dn <- gsub("^","dn ",c_dn)
c_de <- union(c_up,c_dn)

f_up <- gsub("^","up ",f_up)
f_dn <- gsub("^","dn ",f_dn)
f_de <- union(f_up,f_dn)

# ORA vs ORA*
cf <- length(intersect(c_de, f_de )) / length(union(c_de, f_de))

# FCS vs ORA*
mf <- length(intersect(m_de, f_de )) / length(union(m_de, f_de))

dat <- c("FCS vs ORA"=cm,"ORA vs ORA*"=cf,"FCS vs ORA*"=mf, "ORA vs ORA comb"=dc, "ORA vs ORA* comb"=ec)

dat

##       FCS vs ORA      ORA vs ORA*      FCS vs ORA*  ORA vs ORA comb 
##       0.65647059       0.41193595       0.52594670       0.03960396 
## ORA vs ORA* comb 
##       0.07854985

barplot(dat,ylab="jaccard metric",horiz=TRUE,las=1)

saveRDS(dat,file = "ex1dat.rds")

Session information

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
##  [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8   
##  [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] eulerr_6.1.1                kableExtra_1.3.4           
##  [3] mitch_1.4.1                 clusterProfiler_4.0.5      
##  [5] DESeq2_1.32.0               SummarizedExperiment_1.22.0
##  [7] Biobase_2.52.0              MatrixGenerics_1.4.3       
##  [9] matrixStats_0.61.0          GenomicRanges_1.44.0       
## [11] GenomeInfoDb_1.28.4         IRanges_2.26.0             
## [13] S4Vectors_0.30.2            BiocGenerics_0.38.0        
## [15] getDEE2_1.2.0              
## 
## loaded via a namespace (and not attached):
##   [1] shadowtext_0.0.9       fastmatch_1.1-3        systemfonts_1.0.3     
##   [4] plyr_1.8.6             igraph_1.2.8           lazyeval_0.2.2        
##   [7] polylabelr_0.2.0       splines_4.1.2          BiocParallel_1.26.2   
##  [10] ggplot2_3.3.5          digest_0.6.28          yulab.utils_0.0.4     
##  [13] htmltools_0.5.2        GOSemSim_2.18.1        viridis_0.6.2         
##  [16] GO.db_3.13.0           fansi_0.5.0            magrittr_2.0.1        
##  [19] memoise_2.0.0          Biostrings_2.60.2      annotate_1.70.0       
##  [22] graphlayouts_0.7.1     svglite_2.0.0          enrichplot_1.12.3     
##  [25] colorspace_2.0-2       rvest_1.0.2            blob_1.2.2            
##  [28] ggrepel_0.9.1          xfun_0.28              dplyr_1.0.7           
##  [31] crayon_1.4.2           RCurl_1.98-1.5         jsonlite_1.7.2        
##  [34] scatterpie_0.1.7       genefilter_1.74.1      survival_3.2-13       
##  [37] ape_5.5                glue_1.5.0             polyclip_1.10-0       
##  [40] gtable_0.3.0           zlibbioc_1.38.0        XVector_0.32.0        
##  [43] webshot_0.5.2          htm2txt_2.1.1          DelayedArray_0.18.0   
##  [46] scales_1.1.1           DOSE_3.18.3            DBI_1.1.1             
##  [49] GGally_2.1.2           Rcpp_1.0.7             viridisLite_0.4.0     
##  [52] xtable_1.8-4           gridGraphics_0.5-1     tidytree_0.3.6        
##  [55] bit_4.0.4              htmlwidgets_1.5.4      httr_1.4.2            
##  [58] fgsea_1.18.0           gplots_3.1.1           RColorBrewer_1.1-2    
##  [61] ellipsis_0.3.2         pkgconfig_2.0.3        reshape_0.8.8         
##  [64] XML_3.99-0.8           farver_2.1.0           sass_0.4.0            
##  [67] locfit_1.5-9.4         utf8_1.2.2             ggplotify_0.1.0       
##  [70] tidyselect_1.1.1       rlang_0.4.12           reshape2_1.4.4        
##  [73] later_1.3.0            AnnotationDbi_1.54.1   munsell_0.5.0         
##  [76] tools_4.1.2            cachem_1.0.6           downloader_0.4        
##  [79] generics_0.1.1         RSQLite_2.2.8          evaluate_0.14         
##  [82] stringr_1.4.0          fastmap_1.1.0          yaml_2.2.1            
##  [85] ggtree_3.0.4           knitr_1.36             bit64_4.0.5           
##  [88] tidygraph_1.2.0        caTools_1.18.2         purrr_0.3.4           
##  [91] KEGGREST_1.32.0        ggraph_2.0.5           nlme_3.1-153          
##  [94] mime_0.12              aplot_0.1.1            xml2_1.3.2            
##  [97] DO.db_2.9              rstudioapi_0.13        compiler_4.1.2        
## [100] beeswarm_0.4.0         png_0.1-7              treeio_1.16.2         
## [103] tibble_3.1.6           tweenr_1.0.2           geneplotter_1.70.0    
## [106] bslib_0.3.1            stringi_1.7.5          highr_0.9             
## [109] lattice_0.20-45        Matrix_1.3-4           vctrs_0.3.8           
## [112] pillar_1.6.4           lifecycle_1.0.1        jquerylib_0.1.4       
## [115] data.table_1.14.2      cowplot_1.1.1          bitops_1.0-7          
## [118] httpuv_1.6.3           patchwork_1.1.1        qvalue_2.24.0         
## [121] R6_2.5.1               promises_1.2.0.1       KernSmooth_2.23-20    
## [124] echarts4r_0.4.2        gridExtra_2.3          gtools_3.9.2          
## [127] MASS_7.3-54            assertthat_0.2.1       GenomeInfoDbData_1.2.6
## [130] grid_4.1.2             ggfun_0.0.4            tidyr_1.1.4           
## [133] rmarkdown_2.11         ggforce_0.3.3          shiny_1.7.1