library(goSorensen)

1 Introduction

The goal of goSorensen is to implement the equivalence test introduced in Flores, P., Salicrú, M., Sánchez-Pla, A. and Ocaña, J.(2022) “An equivalence test between features lists, based on the Sorensen - Dice index and the joint frequencies of GO node enrichment”, BMC Bioinformatics, 2022 23:207.

Given two gene lists, \(L_1\) and \(L_2\), (the data) and a given set of n Gene Ontology (GO) terms (the frame of reference for biological information in these lists), the test is devoted to answer the following question (quite informally stated for the moment): The dissimilarity between the biological information in both lists, is it negligible? To measure the dissimilarity we use the Sorensen-Dice index:

\[ \hat d_{12} = d(L_1,L_2) = \frac{2n_{11}}{2n_{11} + n_{10} + n_{01}} \]

where \(n_{11}\) corresponds to the number of GO terms (among the n GO terms under consideration) which are enriched in both gene lists, \(n_{10}\) corresponds to the GO terms enriched in \(L_1\) but not in \(L_2\) and \(n_{01}\) the reverse, those enriched in \(L_2\) but not in \(L_1\). For notation completeness, \(n_{00}\) would correspond to those GO terms not enriched in both lists; it is not considered by the Sorensen-Dice index but would be necessary in some computations. Obviously, \(n = n_{11} + n_{10} + n_{01} + n_{00}\).

More precisely, the above problem can be restated as follows: Given a negligibility threshold \(d_0\) for the Sorensen-Dice values, to decide negligibility corresponds to rejecting the null hypothesis \(H_0: d \ge d_0\) in favor of the alternative \(H_1: d < d_0\), where \(d\) stands for the “true” value of the Sorensen-Dice dissimilarity (\(L_1\) and \(L_2\) are samples, and the own process of declaring enrichment of a GO term is random, so \(\hat d = d(L_1,L_2)\) is an estimate of \(d\)). Then, a bit more precise statement of the problem is “The dissimilarity between the biological information in two gene lists, is it negligible up to a degree \(d_0\)?” Where this information is expressed by means of the Sorensen-Dice dissimilarity measured on the degree of coincidence and non-coincidence in GO terms enrichment among a given set of GO terms.

For the moment, the reference set of GO terms can be only all those GO terms in a given level of one GO ontology, either BP, CC or MF.

2 Installation

goSorensen package has to be installed with a working R version (>=4.2.0). Installation could take a few minutes on a regular desktop or laptop. Package can be installed from Bioconductor or devtools package, then it needs to be loaded using library(goSorensen)

To install from Bioconductor (recommended):

## Only if BiocManager is not previosly installed:
install.packages("BiocManager")

## otherwise, directly:
BiocManager::install("goSorensen")

To install from Github

devtools::install_github("pablof1988/goSorensen", build_vignettes = TRUE)

3 Data.

The dataset used in this vignette, allOncoGeneLists, is based on the gene lists compiled at http://www.bushmanlab.org/links/genelists, a comprehensive set of gene lists related to cancer. The package goSorensen loads this dataset by means of data(allOncoGeneLists):

data("allOncoGeneLists")

It is a “list” object of length 7. Each one of its elements is a “character” with the gene identifiers of a gene list related to cancer.

4 Performing the equivalence test

4.1 From a previously built contingency table of join enrichment

Function equivTestSorensen performs the equivalence test. One possibility is to build first the mutual enrichment contingency table by means of the function buildEnrichTable and then to perform the equivalence test:

data("humanEntrezIDs")
length(allOncoGeneLists)
## [1] 7
sapply(allOncoGeneLists, length)
##         atlas      cangenes           cis miscellaneous        sanger 
##           991           189           613           187           450 
##    Vogelstein       waldman 
##           419           426
# First 20 gene identifiers of gene lists Vogelstein and sanger:
allOncoGeneLists[["Vogelstein"]][1:20]
##  [1] "10006"  "25"     "27"     "23305"  "91"     "4299"   "3899"   "27125" 
##  [9] "207"    "238"    "139285" "324"    "367"    "23092"  "23365"  "8289"  
## [17] "57492"  "196528" "405"    "79058"
allOncoGeneLists[["sanger"]][1:20]
##  [1] "25"     "27"     "2181"   "57082"  "10962"  "51517"  "27125"  "10142" 
##  [9] "207"    "208"    "217"    "238"    "57714"  "324"    "23365"  "399"   
## [17] "8289"   "405"    "79058"  "171023"
# Build the enrichment contingency table between gene lists Vogelstein and 
# sanger for the MF ontology at GO level 5:
enrichTab <- buildEnrichTable(allOncoGeneLists[["Vogelstein"]],
                              allOncoGeneLists[["sanger"]],
                              geneUniverse = humanEntrezIDs, orgPackg = "org.Hs.eg.db",
                              onto = "MF", GOLevel = 5, listNames = c("Vogelstein", "sanger"))
enrichTab
##                       Enriched in sanger
## Enriched in Vogelstein TRUE FALSE
##                  TRUE    29    11
##                  FALSE    1  2170
# Equivalence test for an equivalence (or negligibility) limit 0.2857
testResult <- equivTestSorensen(enrichTab, d0 = 0.2857)
testResult
## 
##  Normal asymptotic test for 2x2 contingency tables based on the
##  Sorensen-Dice dissimilarity
## 
## data:  enrichTab
## (d - d0) / se = -2.315, p-value = 0.01031
## alternative hypothesis: true equivalence limit d0 is less than 0.2857
## 95 percent confidence interval:
##  0.000000 0.252619
## sample estimates:
## Sorensen dissimilarity 
##              0.1714286 
## attr(,"se")
## standard error 
##     0.04936026

4.2 Directly from the gene lists.

To perform the test directly from the gene lists (internally building the contingency table) is also possible:

equivTestSorensen(allOncoGeneLists[["Vogelstein"]], allOncoGeneLists[["sanger"]], d0 = 0.2857,
                              geneUniverse = humanEntrezIDs, orgPackg = "org.Hs.eg.db",
                              onto = "MF", GOLevel = 5, listNames = c("Vogelstein", "sanger"))
## 
##  Normal asymptotic test for 2x2 contingency tables based on the
##  Sorensen-Dice dissimilarity
## 
## data:  tab
## (d - d0) / se = -2.315, p-value = 0.01031
## alternative hypothesis: true equivalence limit d0 is less than 0.2857
## 95 percent confidence interval:
##  0.000000 0.252619
## sample estimates:
## Sorensen dissimilarity 
##              0.1714286 
## attr(,"se")
## standard error 
##     0.04936026

To save computing time, the first option (building the contingency table separately, first) may be preferable: buildEnrichTable may take some time (many enrichment tests) and it would be advantageous to have the contingency table ready for further computations.

The above tests use a standard normal approximation to the sample distribution of the \((\hat d - d) / \widehat {se}\) statistic, where \(\widehat {se}\) stands for the standard error of the sample dissimilarity, \(\hat d\).

4.3 Using a bootstrap aproximation.

Alternatively, it is possible to estimate this distribution by means of bootstrap:

boot.testResult <- equivTestSorensen(enrichTab, d0 = 0.2857, boot = TRUE)
boot.testResult
## 
##  Bootstrap test for 2x2 contingency tables based on the Sorensen-Dice
##  dissimilarity (10000 bootstrap replicates)
## 
## data:  enrichTab
## (d - d0) / se = -2.315, p-value = 0.0349
## alternative hypothesis: true equivalence limit d0 is less than 0.2857
## 95 percent confidence interval:
##  0.0000000 0.2718094
## sample estimates:
## Sorensen dissimilarity 
##              0.1714286 
## attr(,"se")
## standard error 
##     0.04936026

For low frequencies in the contingency table, bootstrap is a more conservative but preferable approach, with better type I error control.

5 Accessing to specific fields

To access specific fields of the test result:

getDissimilarity(testResult)
## Sorensen dissimilarity 
##              0.1714286 
## attr(,"se")
## standard error 
##     0.04936026
getSE(testResult)
## standard error 
##     0.04936026
getPvalue(testResult)
##    p-value 
## 0.01030512
getTable(testResult)
##                       Enriched in sanger
## Enriched in Vogelstein TRUE FALSE
##                  TRUE    29    11
##                  FALSE    1  2170
getUpper(testResult)
##   dUpper 
## 0.252619
# In the bootstrap approach, only these differ:
getPvalue(boot.testResult)
##    p-value 
## 0.03489651
getUpper(boot.testResult)
##    dUpper 
## 0.2718094
# (Only available for bootstrap tests) efective number of bootstrap resamples:
getNboot(boot.testResult)
## [1] 10000

7 All pairwise tests (or other computations)

For objects of class list, all these functions (equivTestSorensen, dSorensen, seSorensen, duppSorensen) assume a list of character objects containing gene identifiers and all pairwise computations are performed. For example, to obtain the matrix of all pairwise Sorensen-Dice dissimilarities:

dSorensen(allOncoGeneLists, onto = "MF", GOLevel = 5, 
          geneUniverse = humanEntrezIDs, orgPackg = "org.Hs.eg.db")
##                   atlas cangenes       cis miscellaneous    sanger Vogelstein
## atlas         0.0000000        1 0.7358491     0.5000000 0.3866667  0.3647059
## cangenes      1.0000000        0 1.0000000     1.0000000 1.0000000  1.0000000
## cis           0.7358491        1 0.0000000     0.6774194 0.6842105  0.7083333
## miscellaneous 0.5000000        1 0.6774194     0.0000000 0.4339623  0.5238095
## sanger        0.3866667        1 0.6842105     0.4339623 0.0000000  0.1714286
## Vogelstein    0.3647059        1 0.7083333     0.5238095 0.1714286  0.0000000
## waldman       0.3647059        1 0.7083333     0.3333333 0.4571429  0.4750000
##                 waldman
## atlas         0.3647059
## cangenes      1.0000000
## cis           0.7083333
## miscellaneous 0.3333333
## sanger        0.4571429
## Vogelstein    0.4750000
## waldman       0.0000000

Similarly, the following code performs all pairwise tests:

allTests <- equivTestSorensen(allOncoGeneLists, d0 = 0.2857, 
                              onto = "MF", GOLevel = 5, 
                              geneUniverse = humanEntrezIDs, 
                              orgPackg = "org.Hs.eg.db")
getPvalue(allTests)
##           cangenes.atlas.p-value                cis.atlas.p-value 
##                              NaN                       0.99999999 
##             cis.cangenes.p-value      miscellaneous.atlas.p-value 
##                              NaN                       0.99786392 
##   miscellaneous.cangenes.p-value        miscellaneous.cis.p-value 
##                              NaN                       0.99979410 
##             sanger.atlas.p-value          sanger.cangenes.p-value 
##                       0.93448380                              NaN 
##               sanger.cis.p-value     sanger.miscellaneous.p-value 
##                       0.99996940                       0.96363878 
##         Vogelstein.atlas.p-value      Vogelstein.cangenes.p-value 
##                       0.90045435                              NaN 
##           Vogelstein.cis.p-value Vogelstein.miscellaneous.p-value 
##                       0.99999944                       0.99879102 
##        Vogelstein.sanger.p-value            waldman.atlas.p-value 
##                       0.01030512                       0.90045435 
##         waldman.cangenes.p-value              waldman.cis.p-value 
##                              NaN                       0.99999944 
##    waldman.miscellaneous.p-value           waldman.sanger.p-value 
##                       0.75372379                       0.99090749 
##       waldman.Vogelstein.p-value 
##                       0.99717965
getDissimilarity(allTests, simplify = FALSE)
##                   atlas cangenes       cis miscellaneous    sanger Vogelstein
## atlas         0.0000000        1 0.7358491     0.5000000 0.3866667  0.3647059
## cangenes      1.0000000        0 1.0000000     1.0000000 1.0000000  1.0000000
## cis           0.7358491        1 0.0000000     0.6774194 0.6842105  0.7083333
## miscellaneous 0.5000000        1 0.6774194     0.0000000 0.4339623  0.5238095
## sanger        0.3866667        1 0.6842105     0.4339623 0.0000000  0.1714286
## Vogelstein    0.3647059        1 0.7083333     0.5238095 0.1714286  0.0000000
## waldman       0.3647059        1 0.7083333     0.3333333 0.4571429  0.4750000
##                 waldman
## atlas         0.3647059
## cangenes      1.0000000
## cis           0.7083333
## miscellaneous 0.3333333
## sanger        0.4571429
## Vogelstein    0.4750000
## waldman       0.0000000

8 Alternative representations of contingency tables of joint enrichment

Besides admitting objects of class table, character and list, functions equivTestSorensen, dSorensen, seSorensen and duppSorensen are also adequate for contingency tables represented as a plain matrix or a numeric:

enrichMat <- matrix(c(20, 1, 9, 2149), nrow = 2)
enrichMat
##      [,1] [,2]
## [1,]   20    9
## [2,]    1 2149
dSorensen(enrichMat)
## [1] 0.2
enrichVec <- c(20, 1, 9, 2149)
equivTestSorensen(enrichVec)
## 
##  Normal asymptotic test for 2x2 contingency tables based on the
##  Sorensen-Dice dissimilarity
## 
## data:  enrichVec
## (d - d0) / se = -3.8784, p-value = 5.257e-05
## alternative hypothesis: true equivalence limit d0 is less than 0.4444444
## 95 percent confidence interval:
##  0.0000000 0.3036703
## sample estimates:
## Sorensen dissimilarity 
##                    0.2 
## attr(,"se")
## standard error 
##     0.06302709
equivTestSorensen(enrichVec, boot = TRUE)
## 
##  Bootstrap test for 2x2 contingency tables based on the Sorensen-Dice
##  dissimilarity (10000 bootstrap replicates)
## 
## data:  enrichVec
## (d - d0) / se = -3.8784, p-value = 0.0043
## alternative hypothesis: true equivalence limit d0 is less than 0.4444444
## 95 percent confidence interval:
##  0.0000000 0.3323029
## sample estimates:
## Sorensen dissimilarity 
##                    0.2 
## attr(,"se")
## standard error 
##     0.06302709
len3Vec <- c(20, 1, 9)
dSorensen(len3Vec)
## [1] 0.2
seSorensen(len3Vec)
## [1] 0.06302709
duppSorensen(len3Vec)
## [1] 0.3036703
# Error, bootstrapping requires the full (4 values) contingency table:
try(duppSorensen(len3Vec, boot = TRUE), TRUE)

Session information

All software and respective versions used to produce this document are listed below.

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
## 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_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       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] goSorensen_1.4.0 BiocStyle_2.30.0
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3            jsonlite_1.8.7               
##   [3] magrittr_2.0.3                farver_2.1.1                 
##   [5] rmarkdown_2.25                fs_1.6.3                     
##   [7] zlibbioc_1.48.0               vctrs_0.6.4                  
##   [9] memoise_2.0.1                 RCurl_1.98-1.12              
##  [11] CompQuadForm_1.4.3            ggtree_3.10.0                
##  [13] htmltools_0.5.6.1             AnnotationHub_3.10.0         
##  [15] curl_5.1.0                    gridGraphics_0.5-1           
##  [17] sass_0.4.7                    bslib_0.5.1                  
##  [19] plyr_1.8.9                    cachem_1.0.8                 
##  [21] igraph_1.5.1                  mime_0.12                    
##  [23] lifecycle_1.0.3               pkgconfig_2.0.3              
##  [25] gson_0.1.0                    Matrix_1.6-1.1               
##  [27] R6_2.5.1                      fastmap_1.1.1                
##  [29] GenomeInfoDbData_1.2.11       shiny_1.7.5.1                
##  [31] digest_0.6.33                 aplot_0.2.2                  
##  [33] enrichplot_1.22.0             colorspace_2.1-0             
##  [35] patchwork_1.1.3               AnnotationDbi_1.64.0         
##  [37] S4Vectors_0.40.0              RSQLite_2.3.1                
##  [39] org.Hs.eg.db_3.18.0           MPO.db_0.99.7                
##  [41] filelock_1.0.2                fansi_1.0.5                  
##  [43] httr_1.4.7                    polyclip_1.10-6              
##  [45] HPO.db_0.99.2                 compiler_4.3.1               
##  [47] bit64_4.0.5                   withr_2.5.1                  
##  [49] BiocParallel_1.36.0           viridis_0.6.4                
##  [51] DBI_1.1.3                     ggforce_0.4.1                
##  [53] MASS_7.3-60                   rappdirs_0.3.3               
##  [55] HDO.db_0.99.1                 tools_4.3.1                  
##  [57] scatterpie_0.2.1              ape_5.7-1                    
##  [59] interactiveDisplayBase_1.40.0 httpuv_1.6.12                
##  [61] glue_1.6.2                    nlme_3.1-163                 
##  [63] GOSemSim_2.28.0               promises_1.2.1               
##  [65] shadowtext_0.1.2              grid_4.3.1                   
##  [67] reshape2_1.4.4                fgsea_1.28.0                 
##  [69] generics_0.1.3                gtable_0.3.4                 
##  [71] tidyr_1.3.0                   data.table_1.14.8            
##  [73] tidygraph_1.2.3               utf8_1.2.4                   
##  [75] XVector_0.42.0                BiocGenerics_0.48.0          
##  [77] ggrepel_0.9.4                 BiocVersion_3.18.0           
##  [79] pillar_1.9.0                  stringr_1.5.0                
##  [81] goProfiles_1.64.0             yulab.utils_0.1.0            
##  [83] later_1.3.1                   splines_4.3.1                
##  [85] dplyr_1.1.3                   tweenr_2.0.2                 
##  [87] treeio_1.26.0                 BiocFileCache_2.10.0         
##  [89] lattice_0.22-5                bit_4.0.5                    
##  [91] tidyselect_1.2.0              GO.db_3.18.0                 
##  [93] Biostrings_2.70.0             knitr_1.44                   
##  [95] gridExtra_2.3                 bookdown_0.36                
##  [97] IRanges_2.36.0                stats4_4.3.1                 
##  [99] xfun_0.40                     graphlayouts_1.0.1           
## [101] Biobase_2.62.0                stringi_1.7.12               
## [103] lazyeval_0.2.2                ggfun_0.1.3                  
## [105] yaml_2.3.7                    evaluate_0.22                
## [107] codetools_0.2-19              ggraph_2.1.0                 
## [109] tibble_3.2.1                  qvalue_2.34.0                
## [111] BiocManager_1.30.22           ggplotify_0.1.2              
## [113] cli_3.6.1                     xtable_1.8-4                 
## [115] munsell_0.5.0                 jquerylib_0.1.4              
## [117] Rcpp_1.0.11                   GenomeInfoDb_1.38.0          
## [119] dbplyr_2.3.4                  png_0.1-8                    
## [121] parallel_4.3.1                ellipsis_0.3.2               
## [123] ggplot2_3.4.4                 blob_1.2.4                   
## [125] clusterProfiler_4.10.0        DOSE_3.28.0                  
## [127] bitops_1.0-7                  viridisLite_0.4.2            
## [129] tidytree_0.4.5                scales_1.2.1                 
## [131] purrr_1.0.2                   crayon_1.5.2                 
## [133] rlang_1.1.1                   cowplot_1.1.1                
## [135] fastmatch_1.1-4               KEGGREST_1.42.0

References

Appendix

Flores, P., Salicrú, M., Sánchez-Pla, A. et al. An equivalence test between features lists, based on the Sorensen–Dice index and the joint frequencies of GO term enrichment. BMC Bioinformatics 23, 207 (2022). https://doi.org/10.1186/s12859-022-04739-2