Case Study: Gene Set Enrichment Analysis (Human Data Set)
Alejandro Reyes and Keegan Korthauer
October 30, 2018
1 Summary
The objective of this document is to compare different methods for multiple testing correction in the context of Gene Set Enrichment Analysis (GSEA). We use RNA-seq data consisting of cortex and cerebellum samples of a subset of the GTEx individuals. Shortly, samples were downloaded using the Short Read Archive Toolkit and mapped to the human reference genome version GRCh38 using STAR v2.4.2a. htseq-count was used to tabulate the number of uniquely mapping reads for each gene. We used the DESeq2 package to format the data into a DESeqDataSet object.
We use the fgsea Bioconductor package to implement the GSEA method. This is a Functional Class Scoring approach, which does not require setting an arbitrary threshold for Differential Expression, but instead relies on the gene’s rank (here we rank by DESeq2 test statistic).
2 Workspace Setup
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(scales)
library(DESeq2)## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind,
## colMeans, colnames, colSums, dirname, do.call, duplicated,
## eval, evalq, Filter, Find, get, grep, grepl, intersect,
## is.unsorted, lapply, lengths, Map, mapply, match, mget, order,
## paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind,
## Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which, which.max,
## which.min##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges##
## Attaching package: 'IRanges'## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:dplyr':
##
## exprs## Loading required package: DelayedArray## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians## The following object is masked from 'package:dplyr':
##
## count## Loading required package: BiocParallel##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following objects are masked from 'package:base':
##
## aperm, applylibrary(EnsDb.Hsapiens.v86)## Loading required package: ensembldb## Loading required package: GenomicFeatures## Loading required package: AnnotationDbi##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: AnnotationFilter##
## Attaching package: 'ensembldb'## The following object is masked from 'package:dplyr':
##
## filter## The following object is masked from 'package:stats':
##
## filterlibrary(biomaRt)
library(fgsea)## Loading required package: Rcpp## load helper functions
for (f in list.files("../R", "\\.(r|R)$", full.names = TRUE)) {
source(f)
}
## project data/results folders
datdir <- "data"
resdir <- "results"
sbdir <- "../../results/GSEA"
dir.create(datdir, showWarnings = FALSE)
dir.create(resdir, showWarnings = FALSE)
dir.create(sbdir, showWarnings = FALSE)
## intermediary files we create below
count_file <- file.path(datdir, "human-counts.rds")
deseq_file <- file.path(datdir, "human-deseq.rds")
goset_file <- file.path(datdir, "human-gosets.rds")
result_file <- file.path(resdir, "human-results.rds")
bench_file <- file.path(sbdir, "human-benchmark.rds")
bench_file_uninf <- file.path(sbdir, "human-uninf-benchmark.rds")3 Data Preparation
We download the DESeqDataSet from zenodo, if not present locally, and modify the design to test for differences in gene expression between cerebellums and cortex.
if (!file.exists(count_file)) {
download.file("https://zenodo.org/record/1475409/files/gsea-human-counts.rds?download=1",
destfile = count_file)
}
dsdObject <- readRDS(count_file)
design( dsdObject ) <- ~ tissueTo keep running times short, we only perform differential tesing on protein coding genes, as specified in Ensembl release 86.
gnType <- genes(EnsDb.Hsapiens.v86, column = "gene_biotype")
protGenes <- gnType$gene_id[gnType$gene_biotype == "protein_coding"]
dsdObject <- dsdObject[rownames(dsdObject) %in% protGenes, ]4 Data Analysis
4.1 Enrichment Analysis
Next, we run DESeq2 to retrieve a list of differentially expressed genes at a FDR of 10%.
if (!file.exists(deseq_file)) {
dsdObject <- DESeq(dsdObject)
res <- results(dsdObject, independentFiltering = FALSE)
saveRDS(res, file = deseq_file)
} else {
res <- readRDS(deseq_file)
}
genes <- as.numeric(res$padj < 0.1)
names(genes) <- rownames(res)
sum(genes, na.rm=TRUE)## [1] 9853Next we’ll use the biomaRt package to download and associate GO annotations for each gene.
if (!file.exists(goset_file)) {
mart <- useMart("ensembl", "hsapiens_gene_ensembl")
goSets <- getBM(c("ensembl_gene_id", "go_id"), mart = mart,
filters = "ensembl_gene_id", values = rownames(res))
goSets <- goSets[!nchar( goSets$go_id ) == 0, ]
goSets <- with(goSets, split(go_id, ensembl_gene_id))
saveRDS(goSets, file = goset_file)
} else {
goSets <- readRDS(goset_file)
}Now we use the fgsea package to perform the gene set enrichment analysis and obtain a enrichment p-value for each pathway.
# invert the list so each item is a pathway instead of a gene
goSets <- split(rep(names(goSets), lengths(goSets)), unlist(goSets))
stats <- res$stat
names(stats) <- rownames(res)
stats <- stats[!is.na(stats)]
if (!file.exists(result_file)) {
goRes <- fgsea(goSets,
stats,
nperm=10000,
maxSize=500,
minSize=5)
saveRDS(goRes, file = result_file)
} else {
goRes <- readRDS(result_file)
}Add a random (uninformative covariate) to the dataset.
## Add random (uninformative) covariate
set.seed(66778)
goRes$rand_covar <- rnorm(nrow(goRes))4.2 Covariate Diagnostics
Here, we want to check whether the size of the gene set is actually informative and independent under the null.
4.2.1 Gene Set Size
We will explore whether the size of the gene set can be used as a covariate for modern multiple-testing correction methods in the context of GSEA. In the plot below, the log10 of the p-values is plotted as a function of the size of the gene set. There is a pattern in which gene sets with a higher number of genes tend to have smaller p-values, which is indicative that gene set size is an informative covariate.
rank_scatter(dat = goRes, pval = "pval",
covariate = "size", bins = 50, funx = log2,
funfill = log10_trans()) +
ggtitle("Enriched gene sets") +
xlab(expression(log[10]~"(# of genes)")) +
ylab(expression(-log[10]~"(p-value)") ) 
We can also explore if the covariate seems to be independent under the null.
strat_hist(goRes, pval="pval", covariate="size", maxy=11)
4.2.2 Random
We will explore whether the random covariate can be used as a covariate for modern multiple-testing correction methods in the context of GSEA. In the plot below, the log10 of the p-values is plotted as a function of the random covariate. This covariate looks independent of the p-values.
rank_scatter(dat = goRes, pval = "pval",
covariate = "rand_covar", bins = 50,
funfill = log10_trans()) +
ggtitle("Enriched gene sets") +
ylab(expression(-log[10]~"(p-value)") ) 
We can also explore if the covariate seems to be independent under the null.
strat_hist(goRes, pval="pval", covariate="rand_covar", maxy=11)
4.3 Multiple-Testing Correction
We then execute the benchDesign and generate a SummarizedBenchmark object containing multiple-testing corrections using several methods.
## rename columns and prepare for benchmarking
res <- dplyr:::select(goRes, c("pval", "size", "rand_covar")) %>%
dplyr:::rename(ind_covariate = size)
## generate default BenchDesign
bd <- initializeBenchDesign()We don’t include ashq, fdrreg-e and fdrreg-t from the analysis because the necessary assumptions are not met in the current case study. Namely, effect sizes and standard errors are not available for ASH, and the test statistics are not normally distributed under the null and alternative, as required by Scott’s FDR regression methods.
if (!file.exists(bench_file)) {
sGSEA <- buildBench(bd, res, ftCols = "ind_covariate")
saveRDS(sGSEA, file = bench_file)
} else {
sGSEA <- readRDS(bench_file)
}We’ll also compare the results to an uninformative (random) covariate.
if (!file.exists(bench_file_uninf)) {
res$ind_covariate <- res$rand_covar
sGSEA_rand <- buildBench(bd, res, ftCols = "ind_covariate")
saveRDS(sGSEA_rand, file = bench_file_uninf)
} else {
sGSEA_rand <- readRDS(bench_file_uninf)
}4.4 Benchmark Metrics
4.4.1 Gene Set Size
assayNames(sGSEA) <- "qvalue"
sGSEA <- addDefaultMetrics(sGSEA)
sGSEA <- estimatePerformanceMetrics(sGSEA, addColData=TRUE)rejections_scatter(sGSEA, as_fraction=FALSE, supplementary=FALSE)## Found already estimated performance metrics, replacing these
plotFDRMethodsOverlap(sGSEA, alpha=0.1, supplementary=FALSE, order.by="freq", nsets=100)
covariateLinePlot(sGSEA, alpha = 0.1, covname = "ind_covariate", trans = "log1p")## Warning: Removed 25 rows containing missing values (geom_path).
4.4.2 Random
assayNames(sGSEA_rand) <- "qvalue"
sGSEA_rand <- addDefaultMetrics(sGSEA_rand)
sGSEA_rand <- estimatePerformanceMetrics(sGSEA_rand, addColData=TRUE)rejections_scatter(sGSEA_rand, as_fraction=FALSE, supplementary=FALSE)## Found already estimated performance metrics, replacing these
plotFDRMethodsOverlap(sGSEA_rand, alpha=0.1, supplementary=FALSE, order.by="freq", nsets=100)
covariateLinePlot(sGSEA_rand, alpha = 0.1, covname = "ind_covariate", trans = "log1p")## Warning: Removed 25 rows containing missing values (geom_path).
4.5 Covariate comparison
Here we compare the method ranks for the two covariates at alpha = 0.10.
plotMethodRanks(c(bench_file, bench_file_uninf),
colLabels = c("Set Size", "Random"),
alpha = 0.10, xlab = "Covariate",
excludeMethods = NULL)
5 Session Info
sessionInfo()## R version 3.5.0 (2018-04-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS: /usr/lib64/libblas.so.3.4.2
## LAPACK: /usr/lib64/liblapack.so.3.4.2
##
## 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
##
## attached base packages:
## [1] splines parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] adaptMT_1.0.0 FDRreg_0.2-1
## [3] mvtnorm_1.0-8 BayesLogit_0.6
## [5] fda_2.4.7 Matrix_1.2-14
## [7] fdrtool_1.2.15 swfdr_1.6.0
## [9] qvalue_2.12.0 ashr_2.2-7
## [11] IHW_1.8.0 SummarizedBenchmark_0.99.1
## [13] mclust_5.4 stringr_1.3.1
## [15] rlang_0.2.1 UpSetR_1.3.3
## [17] tidyr_0.8.1 bindrcpp_0.2.2
## [19] magrittr_1.5 cowplot_0.9.2
## [21] hexbin_1.27.2 fgsea_1.6.0
## [23] Rcpp_0.12.17 biomaRt_2.36.1
## [25] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.4.1
## [27] AnnotationFilter_1.4.0 GenomicFeatures_1.32.0
## [29] AnnotationDbi_1.42.1 DESeq2_1.20.0
## [31] SummarizedExperiment_1.10.1 DelayedArray_0.6.1
## [33] BiocParallel_1.14.2 matrixStats_0.53.1
## [35] Biobase_2.40.0 GenomicRanges_1.32.3
## [37] GenomeInfoDb_1.16.0 IRanges_2.14.10
## [39] S4Vectors_0.18.3 BiocGenerics_0.26.0
## [41] scales_0.5.0 ggplot2_3.0.0
## [43] dplyr_0.7.5
##
## loaded via a namespace (and not attached):
## [1] colorspace_1.3-2 rprojroot_1.3-2
## [3] htmlTable_1.12 XVector_0.20.0
## [5] ggdendro_0.1-20 base64enc_0.1-3
## [7] rstudioapi_0.7 bit64_0.9-7
## [9] mosaic_1.2.0 codetools_0.2-15
## [11] mnormt_1.5-5 pscl_1.5.2
## [13] doParallel_1.0.11 geneplotter_1.58.0
## [15] knitr_1.20 Formula_1.2-3
## [17] Rsamtools_1.32.0 broom_0.4.4
## [19] annotate_1.58.0 cluster_2.0.7-1
## [21] compiler_3.5.0 httr_1.3.1
## [23] backports_1.1.2 assertthat_0.2.0
## [25] lazyeval_0.2.1 acepack_1.4.1
## [27] htmltools_0.3.6 prettyunits_1.0.2
## [29] tools_3.5.0 gtable_0.2.0
## [31] glue_1.2.0 GenomeInfoDbData_1.1.0
## [33] reshape2_1.4.3 fastmatch_1.1-0
## [35] slam_0.1-43 Biostrings_2.48.0
## [37] nlme_3.1-137 rtracklayer_1.40.3
## [39] iterators_1.0.9 psych_1.8.4
## [41] mosaicCore_0.5.0 XML_3.98-1.11
## [43] zlibbioc_1.26.0 MASS_7.3-50
## [45] ProtGenerics_1.12.0 RColorBrewer_1.1-2
## [47] yaml_2.1.19 curl_3.2
## [49] mosaicData_0.16.0 memoise_1.1.0
## [51] gridExtra_2.3 rpart_4.1-13
## [53] latticeExtra_0.6-28 stringi_1.2.2
## [55] RSQLite_2.1.1 SQUAREM_2017.10-1
## [57] genefilter_1.62.0 foreach_1.4.4
## [59] checkmate_1.8.5 truncnorm_1.0-8
## [61] pkgconfig_2.0.1 bitops_1.0-6
## [63] evaluate_0.10.1 lattice_0.20-35
## [65] lpsymphony_1.8.0 purrr_0.2.5
## [67] bindr_0.1.1 GenomicAlignments_1.16.0
## [69] htmlwidgets_1.2 labeling_0.3
## [71] bit_1.1-14 tidyselect_0.2.4
## [73] ggformula_0.7.0 plyr_1.8.4
## [75] R6_2.2.2 Hmisc_4.1-1
## [77] DBI_1.0.0 pillar_1.2.3
## [79] foreign_0.8-70 withr_2.1.2
## [81] survival_2.41-3 RCurl_1.95-4.10
## [83] nnet_7.3-12 tibble_1.4.2
## [85] rmarkdown_1.10 progress_1.1.2
## [87] locfit_1.5-9.1 grid_3.5.0
## [89] data.table_1.11.4 blob_1.1.1
## [91] digest_0.6.15 xtable_1.8-2
## [93] munsell_0.4.3