zenith performs gene set analysis on the result of differential expression using linear (mixed) modeling with dream by considering the correlation between gene expression traits. This package implements the camera method from the limma package proposed by Wu and Smyth (2012). zenith() is a simple extension of camera() to be compatible with linear (mixed) models implemented in dream().

Standard workflow

# Load packages
library(zenith)
library(edgeR)
library(variancePartition)
library(tweeDEseqCountData)
library(kableExtra)

# Load RNA-seq data from LCL's
data(pickrell)
geneCounts = exprs(pickrell.eset)
df_metadata = pData(pickrell.eset)

# Filter genes
# Note this is low coverage data, so just use as code example
dsgn = model.matrix(~ gender, df_metadata)
keep = filterByExpr(geneCounts, dsgn, min.count=5)

# Compute library size normalization
dge = DGEList(counts = geneCounts[keep,])
dge = calcNormFactors(dge)

# Estimate precision weights using voom
vobj = voomWithDreamWeights(dge, ~ gender, df_metadata)

# Apply dream analysis
fit = dream(vobj, ~ gender, df_metadata)
fit = eBayes(fit)

# Load get_MSigDB database, Hallmark genes
# use gene 'SYMBOL', or 'ENSEMBL' id
# use get_GeneOntology() to load Gene Ontology
msdb.gs = get_MSigDB("H", to="ENSEMBL")
   
# Run zenith analysis, and specific which coefficient to evaluate
res.gsa = zenith_gsa(fit, msdb.gs, 'gendermale', progressbar=FALSE )

# Show top gene sets: head(res.gsa)
kable_styling(kable(head(res.gsa), row.names=FALSE))
coef Geneset NGenes Correlation delta se p.less p.greater PValue Direction FDR
gendermale M5890_HALLMARK_TNFA_SIGNALING_VIA_NFKB 119 0.01 -1.0019645 0.1618263 0.0000000 1.0000000 0.0000000 Down 0.0000019
gendermale M5892_HALLMARK_CHOLESTEROL_HOMEOSTASIS 37 0.01 -1.0915363 0.2296170 0.0000055 0.9999945 0.0000110 Down 0.0002582
gendermale M5932_HALLMARK_INFLAMMATORY_RESPONSE 93 0.01 -0.7819360 0.1722474 0.0000120 0.9999880 0.0000241 Down 0.0003769
gendermale M5947_HALLMARK_IL2_STAT5_SIGNALING 105 0.01 -0.6493291 0.1672385 0.0001194 0.9998806 0.0002389 Down 0.0028066
gendermale M5930_HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 53 0.01 -0.6136470 0.2031832 0.0017867 0.9982133 0.0035735 Down 0.0335905
gendermale M5925_HALLMARK_E2F_TARGETS 109 0.01 0.4817527 0.1659106 0.9975051 0.0024949 0.0049898 Up 0.0355160 
# for each cell type select 3 genesets with largest t-statistic
# and 1 geneset with the lowest
# Grey boxes indicate the gene set could not be evaluted because
#    to few genes were represented
plotZenithResults(res.gsa)

Session Info

## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [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] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] org.Hs.eg.db_3.16.0       AnnotationDbi_1.60.1     
##  [3] IRanges_2.32.0            S4Vectors_0.36.2         
##  [5] msigdbr_7.5.1             kableExtra_1.3.4         
##  [7] tweeDEseqCountData_1.36.0 Biobase_2.58.0           
##  [9] BiocGenerics_0.44.0       variancePartition_1.28.7 
## [11] BiocParallel_1.32.5       ggplot2_3.4.1            
## [13] edgeR_3.40.2              limma_3.54.2             
## [15] zenith_1.0.2              knitr_1.42               
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.4.1             BiocFileCache_2.6.1        
##   [3] systemfonts_1.0.4           plyr_1.8.8                 
##   [5] GSEABase_1.60.0             remaCor_0.0.11             
##   [7] splines_4.2.2               GenomeInfoDb_1.34.9        
##   [9] digest_0.6.31               foreach_1.5.2              
##  [11] htmltools_0.5.4             fansi_1.0.4                
##  [13] magrittr_2.0.3              memoise_2.0.1              
##  [15] doParallel_1.0.17           aod_1.3.2                  
##  [17] Biostrings_2.66.0           annotate_1.76.0            
##  [19] matrixStats_0.63.0          svglite_2.1.1              
##  [21] prettyunits_1.1.1           colorspace_2.1-0           
##  [23] blob_1.2.3                  rvest_1.0.3                
##  [25] rappdirs_0.3.3              rbibutils_2.2.13           
##  [27] xfun_0.37                   dplyr_1.1.0                
##  [29] crayon_1.5.2                RCurl_1.98-1.10            
##  [31] jsonlite_1.8.4              graph_1.76.0               
##  [33] lme4_1.1-31                 iterators_1.0.14           
##  [35] glue_1.6.2                  gtable_0.3.1               
##  [37] zlibbioc_1.44.0             XVector_0.38.0             
##  [39] webshot_0.5.4               DelayedArray_0.24.0        
##  [41] EnrichmentBrowser_2.28.1    Rgraphviz_2.42.0           
##  [43] RcppZiggurat_0.1.6          scales_1.2.1               
##  [45] mvtnorm_1.1-3               DBI_1.1.3                  
##  [47] Rcpp_1.0.10                 viridisLite_0.4.1          
##  [49] xtable_1.8-4                progress_1.2.2             
##  [51] bit_4.0.5                   httr_1.4.5                 
##  [53] gplots_3.1.3                ellipsis_0.3.2             
##  [55] pkgconfig_2.0.3             XML_3.99-0.13              
##  [57] farver_2.1.1                sass_0.4.5                 
##  [59] dbplyr_2.3.1                locfit_1.5-9.7             
##  [61] utf8_1.2.3                  labeling_0.4.2             
##  [63] tidyselect_1.2.0            rlang_1.0.6                
##  [65] reshape2_1.4.4              munsell_0.5.0              
##  [67] tools_4.2.2                 cachem_1.0.7               
##  [69] cli_3.6.0                   generics_0.1.3             
##  [71] RSQLite_2.3.0               broom_1.0.3                
##  [73] evaluate_0.20               stringr_1.5.0              
##  [75] fastmap_1.1.1               yaml_2.3.7                 
##  [77] babelgene_22.9              RhpcBLASctl_0.23-42        
##  [79] bit64_4.0.5                 caTools_1.18.2             
##  [81] purrr_1.0.1                 KEGGREST_1.38.0            
##  [83] nlme_3.1-162                KEGGgraph_1.58.3           
##  [85] xml2_1.3.3                  compiler_4.2.2             
##  [87] pbkrtest_0.5.2              rstudioapi_0.14            
##  [89] filelock_1.0.2              curl_5.0.0                 
##  [91] png_0.1-8                   clusterGeneration_1.3.7    
##  [93] tibble_3.2.0                bslib_0.4.2                
##  [95] stringi_1.7.12              highr_0.10                 
##  [97] lattice_0.20-45             Matrix_1.5-3               
##  [99] nloptr_2.0.3                vctrs_0.5.2                
## [101] pillar_1.8.1                lifecycle_1.0.3            
## [103] RUnit_0.4.32                Rdpack_2.4                 
## [105] jquerylib_0.1.4             bitops_1.0-7               
## [107] GenomicRanges_1.50.2        R6_2.5.1                   
## [109] KernSmooth_2.23-20          codetools_0.2-19           
## [111] boot_1.3-28.1               MASS_7.3-58.3              
## [113] gtools_3.9.4                SummarizedExperiment_1.28.0
## [115] withr_2.5.0                 GenomeInfoDbData_1.2.9     
## [117] parallel_4.2.2              hms_1.1.2                  
## [119] grid_4.2.2                  tidyr_1.3.0                
## [121] minqa_1.2.5                 rmarkdown_2.20             
## [123] Rfast_2.0.7                 MatrixGenerics_1.10.0