1 Introduction
Since the introduction of microarrays (Schena et al. 1995; Nguyen et al. 2002), high-throughput omics technologies have advanced substantially, leading to the development of more advanced technologies such as RNA-seq (Chu and Corey 2012; Van den Berge et al. 2019), single-cell sequencing (Jovic et al. 2022), and other recent ones, which have transformed modern biology and medicine (Negi et al. 2022; Lee 2023). These high-throughput experiments allow us to study the behaviour of hundreds or thousands of genes (or proteins, metabolites, etc., often referred to as “features” in a general way) under different experimental conditions, resulting in lists of these features (e.g., a list of differentially expressed genes in a health condition) with critical information to solve relevant problems, such as the detection of bio-markers, drug development, molecular sub-classification of cancer, clinical applications, personalised medicine, among many others (Jain 2004; Pedersen and Kanigan 2016; Climente-González et al. 2017; Cieślik and Chinnaiyan 2018).
Although obtaining feature lists is an important first step, a comprehensive study should also include a functional analysis to elucidate the biological significance of these data. External biological knowledge databases, such as Gene Ontology GO (Ashburner et al. 2000; Consortium et al. 2023), Kyoto Encyclopedia of Genes and Genomes KEGG (Kanehisa and Goto 2000; Kanehisa et al. 2025), Reactome (Jassal et al. 2020), etc., serve as foundational resources for this purpose, enabling researchers to associate features in the lists with biological concepts or pathways and thereby gain meaningful insights into underlying mechanisms.
In particular, the goSorensen R package focuses on the functional comparison of gene lists using the Gene Ontology (GO) as a framework for representing biological knowledge. The GO is a structured, standardised representation of biological knowledge in which biological concepts, known as GO terms, are connected through formally defined relationships. It is designed to be species-agnostic, allowing gene products, such as the proteins or functional RNAs encoded by genes, from different organisms to be annotated using a common vocabulary and facilitating the comparison and integration of functional information across biological databases.
The GO is organised into three aspects, or sub-ontologies: Molecular Function (MF), Cellular Component (CC), and Biological Process (BP). Molecular Function terms represent molecular-level activities performed by gene products, such as catalytic or binding activities. Cellular Component terms capture the cellular location where a molecular function takes place, including cellular structures, membrane-enclosed compartments, and protein-containing complexes. Biological Process terms represent larger biological processes or programs accomplished by the concerted action of multiple molecular activities. Thus, the three GO aspects provide complementary descriptions of the activities performed by gene products, the cellular locations in which these activities occur, and the broader biological processes in which they participate.
Within each aspect, GO terms are organised hierarchically. The GO is structured as a graph in which terms are represented as nodes and the relationships between them as edges. Child terms are generally more specialised than their parent terms, although a term may have more than one parent. Consequently, the GO hierarchy is not a simple tree but a directed acyclic graph (DAG). Terms representing more general concepts connected to a given term are referred to as its ancestors, whereas more specialised terms below it in the hierarchy are its descendants.
Each of the three GO aspects is represented by a separate root term. GO defines several types of relationships between terms; in this work, particular attention is given to the is a relation, which represents subclass relationships and defines the basic hierarchical structure used by the method. Since a GO term may have more than one parent, several paths can connect a term to the corresponding root.
Throughout this work, we define the GO level of a term as the length of the shortest path composed exclusively of is a relationships connecting that term to the root of its corresponding GO aspect. Under this operational definition, larger GO levels generally correspond to more specific biological concepts.
Figure 1: An example of GO terms in the BP ontology and their relationship to their ancestor terms.
For example, Figure 1 illustrates the multiple inheritance of the GO term GO:0042446 in the Biological Process ontology through is a relationships with its ancestors. The shortest is a path connecting this term to the root has length 4; accordingly, the term is assigned to GO level 4.
Gene products are associated with GO terms through annotations, which link a gene product to a biological concept represented in the ontology. Given a gene list obtained from an experiment, enrichment analysis can be used to identify GO terms that are represented more frequently than expected under an appropriate reference model (Hosack et al. 2003). These enriched terms provide a functional summary of the biological characteristics represented in the gene list and constitute the basic information used by goSorensen to compare gene lists from a functional perspective.
When the goal is to identify similarities across two or more gene lists, it is typically more informative to compare the biological concepts associated with the genes (annotations) rather than the genes themselves as simple elements. Some methods developed for this purpose include the Semantic Similarity (Pesquita 2017) or the PANTHER web tool (Mi et al. 2017). These approaches often use comparison metrics, such as the Resnik, Lin, or Jaccard indices, to quantify similarity based on the overlap of annotations between lists. However, these methods generally overlook that the gene lists represent a sample derived from an experiment (e.g., from RNA-Seq), which may vary across repetitions or samplings. So, studying the sampling distribution associated with these metrics can provide a valuable complement, allowing one to assess whether an observed biological similarity is statistically significant or simply the result of chance.
One of the few methods introduced in the scientific literature to address these statistical aspects to provide significant functional comparison results is goProfiles. For each gene list, this method establishes a vector of annotation frequencies across all GO terms at a given GO level, referred to as an annotation profile. A distance metric is then defined between annotation profiles, and its sampling distribution is evaluated to develop tests for significant functional similarity. While goProfiles represents a robust statistical framework, it also has certain limitations. In particular, including GO terms with low or zero frequency, or assigning equal weight to both enriched and non-enriched terms, may introduce unnecessary noise into the distance metric.
To address these issues, we have developed a method for comparing feature lists, available in Flores (2025) and Flores et al. (2022), which we hereafter refer to as the goSorensen method. The fundamental guiding principle underlying this method is that two lists can be considered biologically similar if there is a significant prevalence of shared enriched GO terms in both lists compared to non-shared enriched GO terms. We assess the prevalence of this proportion using a dissimilarity measure derived from the Sorensen-Dice index (Sorensen 1948), whose sample distribution was studied and implemented, allowing for more general conclusions regarding the biological similarity between lists. The goSorensen R package (Flores and Ocana 2026), implements this novel method within the R/Bioconductor ecosystem. Its design integrates enrichment analysis, statistical testing, and visualisation in a unified workflow, ensuring full compatibility with widely used packages such as clusterProfiler (Wu et al. 2021), goProfiles (Sanchez et al. 2026), and different databases of functional annotation. By formalising the entire analytical process into a coherent, object-oriented structure and offering extensive documentation and vignettes, goSorensen enables researchers to perform statistically grounded and reproducible functional comparisons of gene lists efficiently.
This article provides a comprehensive overview of how the goSorensen package manages the data flow (genes, GO annotations and enrichment analysis, statistical analysis, etc.) involved in the goSorensen method to compare gene lists, providing the user with an easy-to-understand and use environment. The following section briefly explains the theoretical framework of the goSorensen method. The third section describes the general design of the package, elucidating its structure and leading advantages. In section four, we present an application of the package using real gene lists, and finally, section five is devoted to the relevant conclusions.
2 The goSorensen method
Figure 2: Sequential flow of the goSorensen method to detect biological similarity between feature lists
Figure 2 illustrates an overview of the sequential steps involved in the goSorensen method, which are detailed in Flores et al. (2022) and Flores (2025). As discussed above, the first step in comparing the \(s\) gene lists derived from a high-throughput experiment involves annotating each element of these lists in the \(n\) GO terms of a specific ontology (BP, CC, or MF) and level for performing an analysis to detect the enriched and non-enriched GO terms by each of the lists. The enrichment matrix compiles the results of this enrichment analysis. In Figure 2 the enriched terms of the matrix are shown with a number \(1\) and the non-enriched terms with a \(0\) but, in fact, in package goSorensen this is a logical matrix, of TRUE/FALSE values.
The \(2 \times 2\) enrichment contingency table, also depicted in Figure 2, summarises the joint and the marginal enrichment patterns of any two specific lists from the enrichment matrix, where:
\(n_{11}\) denotes the number of GO terms enriched in both lists, representing the joint enrichment.
\(n_{10}\) denotes the number of GO terms enriched exclusively in the first list \(L_i\), and non-enriched in the second list \(L_j\).
\(n_{01}\) represents the number of GO terms enriched only in the second list \(L_j\), and non-enriched in the first list \(L_i\).
\(n_{00}\) represents the number of GO terms not enriched in either of the two gene lists.
The next step is to assess the degree of prevalence of the joint enrichment (\(n_{11}\)) over the marginal enrichment (\(n_{01} + n_{10}\)). For this purpose, we use the following measure, based on the Sorensen-Dice index, which we refer to as the Sorensen dissimilarity \(\widehat{d}_S\): \[\begin{equation*} \widehat{d}_S = 1 - \dfrac{2n_{11}}{2n_{11} + n_{10} + n_{01}}, \\ \end{equation*}\] or, in terms of relatives frequencies:
\[\begin{equation*} \begin{aligned} \widehat{d}_S = & 1 - \dfrac{2\widehat{p}_{11}}{2\hat{p}_{11} + \widehat{p}_{10} + \widehat{p}_{01}}, \hspace{0.75cm} \text{with} \hspace{0.2cm} \widehat{p}_{ij} = n_{ij} / n. \end{aligned} \end{equation*}\]
Considering that the frequencies (\(n_{11}, n_{01}, n_{10}, n_{00}\)) come from a multinomial distribution with parameters \((n, (p_{11}, p_{01}, p_{10}, p_{00}))\), where \(p_{11}\) denotes the joint enrichment probability and (\(p_{01}, p_{10}\)) denotes the marginal enrichment probabilities, it is possible to prove that the sample distribution of \(\widehat{d}_S\) can be approximated by a normal distribution with mean given by the theoretical dissimilarity \(d_S\), and variance given by \(\widehat \sigma _S^2\):
\[\begin{equation*} \begin{aligned} \widehat{d}_S \thickapprox \text{Norm} \left(d_S, \frac{\widehat{\sigma}_s}{\sqrt{n}}\right), \hspace{0.5cm} \text{when } n \to \infty. \end{aligned} \end{equation*}\]
being
\(d_S = 1 - \cfrac{2p_{11}}{2p_{11} + p_{10} + p_{01}} \quad \text{and} \quad \widehat{\sigma}_S^2 = \cfrac{4\widehat{p}_{11}\left(\widehat{p}_{10} + \widehat{p}_{01}\right)\left(\widehat{p}_{11} + \widehat{p}_{10} + \widehat{p}_{01}\right)}{\left(2\widehat{p}_{11} + \widehat{p}_{10} + \widehat{p}_{01}\right)^4}\)
So, using a significance level \(\alpha\), one can obtain the following upper bound for a confidence interval of \(d_S\):
\[\widehat{d}_S - z_\alpha \cfrac{\widehat \sigma_S}{\sqrt{n}},\] where \(z_\alpha\) is the \(\alpha\) quantile associated with the studentised statistic \(Z_S=\cfrac{\widehat{d}_S - d_S}{\widehat \sigma_S/\sqrt{n}}\), which asymptotically follows a standard normal distribution \(\text{Norm}(0, 1)\):
\[\begin{equation*} \begin{aligned} Z_S = \cfrac{\widehat{d}_S - d_S}{\widehat \sigma_S/\sqrt{n}} \thickapprox \text{Norm}(0, 1), \hspace{0.5cm} \text{when } n \to \infty. \end{aligned} \end{equation*}\]
While the standard normal distribution adequately fits the true distribution of the studentized statistic \(Z_S\), its accuracy decreases at low enrichment levels, which results in an increased type I error probability (TIEP) during inference processes. Figure 3 interactively illustrates the results of a simulation process, which demonstrates that for low enrichment levels, the distribution derived from the bootstrap methodology fits the true distribution of \(Z_S\) better than the theoretical normal distribution \(\text{Norm}(0, 1)\).
Figure 3: The true distribution (approximated) of the studentized statistic \(Z_S\) and its fit degree with its bootstrap estimation \(Z_S^*\) and the normal distribution \(\text{Norm}(0, 1)\), for three different enrichment levels: Low \((p_{11} = 0.0125, p_{01} = 0.005, p_{10} = 0.005)\), Moderate \((p_{11} = 0.125, p_{01} = 0.05, p_{10} = 0.05)\) and High \((p_{11} = 0.4, p_{01} = 0.2, p_{10} = 0.2)\).
Both approximations (normal and bootstrap) to the true sample distribution play a crucial role in the sequential flow illustrated in Figure 2. They provide the basis for developing the statistical analysis to assess biological similarities between gene lists. These statistical analysis include:
A test for the equivalence hypotheses: \[\begin{equation*} \begin{split} H_0: d_S \ge d_0 \\ H_1: d_S < d_0 \end{split} \end{equation*}\] to establish the irrelevance of biological dissimilarities, up to a predefined equivalence limit \(d_0\). This test can be extended to the equivalence comparison of \(s > 2\) lists.
A new index of dissimilarity between lists, which is based on the irrelevance-threshold determining whether the lists are significantly equivalent. This dissimilarity index can be employed in subsequent analyses, such as clustering or multidimensional scaling (MDS), where the grouping of biologically similar gene lists is grounded in the theory of equivalence testing.
See (Flores et al. 2022) and (Flores 2025) for more details on these statistical methods, as well as for other issues such as criteria for setting equivalence limits or the desirability of performing equivalence tests (Wellek 2002) versus the typical (and inadvisable, in our opinion) difference testing approach \(H_0: d_S=0\) versus \(H_1:d_S>0\).
3 goSorensen: A Bioconductor R-package
goSorensen is the informatics tool designed to implement the goSorensen method. This software package is hosted on Bioconductor, an open-source project based on the R programming language, specifically designed for bioinformatics research (Gentleman et al. 2004). This platform provides goSorensen with multiple advantages, including ease of visualisation, access, community support, documentation, reproducibility, and, most importantly, the ability to integrate with dependency packages essential for performing calculations that are part of the goSorensen method. Key dependency packages include software packages such as clusterProfiler (Yu et al. 2012; Wu et al. 2021) and goProfiles (Sanchez et al. 2026) for developing the enrichment analysis at different GO levels, and the GO.db annotation package (Carlson 2026a), which provides extensive information on ontologies, terms, and annotations of the Gene Ontology, which is constantly updated to reflect changes and expansions to this biological knowledge database. Furthermore, Bioconductor hosts numerous genomic databases annotation packages that enable goSorensen to expand its analysis to several species or organisms, such as org.Hs.eg.db for the Homo sapiens (Carlson 2026c), org.Mm.eg.db for the Mus musculus (Carlson 2026d), org.Dm.eg.db for the Drosophila melanogaster (Carlson 2026b), and others. As new annotation resources become available in Bioconductor, goSorensen can be easily extended to additional organisms, ensuring scalability and broad applicability.
3.1 Design and structure of the goSorensen package.
Managing the extensive data flow associated with the goSorensen method is a considerable challenge. This data flow mainly includes the gene lists, which often contain hundreds or even thousands of elements; the annotations of each of these elements in each of the thousands of GO terms located in a specific ontology and level; the enrichment analysis for these GO terms; and the results of the statistical analyses, such as hypothesis tests, irrelevance-threshold matrices of dissimilarities, relevant statistics, and plots. Figure 2 summarises this data flow, from obtaining gene lists to the enrichment and statistical analyses that ultimately assess biological similarity. Each step in this process is represented in goSorensen by specific data structures: Objects of class “list” correspond to the gene lists used as inputs; objects of class “matrix” and “tableList” represent the enrichment analyses; and objects of class “equivSDhtestList” and “distList” capture the statistical analysis and final results.
Figure 4: Organisational structure of the main classes, functions, and documentation elements in the goSorensen package, designed to manage the extensive data flow associated with the goSorensen method.
The components of the goSorensen package are structured to handle this large data flow efficiently, ensuring consistency between intermediate objects and facilitating integration with the accompanying documentation and vignettes. Figure 4 illustrates this organisation, showing the most relevant classes, functions, and documentation elements that define the design of goSorensen. To support their use, the package includes several example objects representing instances of these classes. These examples, accessible through the standard R function data (e.g., data("allOncoGeneLists")), are used throughout the vignettes and help files to demonstrate the package workflow. Below, we describe the most relevant classes, functions and documentation elements depicted in Figure 4:
The standard R class “list” is used to represent a collection of several gene lists (e.g., the gene lists obtained by different laboratories that have carried out a similar experiment). Each one of the components of this “list” is an object of class “character”, containing a vector of gene identifiers1. This structure provides a flexible and intuitive way to manage multiple gene lists as input to the goSorensen method. The object allOncoGeneLists exemplifies this data organisation and serves as a reference dataset for demonstration purposes:
[1] "list"length(allOncoGeneLists)[1] 7names(allOncoGeneLists)[1] "atlas" "cangenes" "cis" "miscellaneous"
[5] "sanger" "Vogelstein" "waldman" class(allOncoGeneLists$atlas)[1] "character"# First ENTREZ gene identifiers in gene list atlas:
head(allOncoGeneLists$atlas)[1] "11186" "27086" "239" "7764" "5934" "246" In goSorensen, the R class “matrix” is used to store the results obtained from the enrichment analysis. Each row corresponds to a GO term, and each column represents one of the analysed gene lists, where the logical values indicate whether the GO term is enriched (TRUE) or not (FALSE). This representation offers a simple and efficient way to summarise enrichment results across several gene lists, making it easier to apply the subsequent statistical steps of the goSorensen method. As an example, the object enrichedInBP4 illustrates this structure for the ontology BP at the GO level 4.
data("enrichedInBP4")class(enrichedInBP4)[1] "matrix" "array" dim(enrichedInBP4)[1] 473 7# Only those GO terms enriched in almost one gene list are included as rows of
# enrichedInBP4. The true number of terms in this ontology and level is:
attr(enrichedInBP4, "nTerms")[1] 3381colnames(enrichedInBP4)[1] "atlas" "cangenes" "cis" "miscellaneous"
[5] "sanger" "Vogelstein" "waldman" [1] "GO:0001649" "GO:0030278" "GO:0030279" "GO:0030282" "GO:0045778"
[6] "GO:0048755"The class “tableList” is specific to goSorensen and corresponds to the typical output of the function buildEnrichTable. It is an inheritance specialisation of “list”, which stores the enrichment contingency tables obtained from all pairwise comparisons between gene lists. Each element of an object “tableList” represents one of the analysed lists and contains the contingency tables associated with its comparisons to the other lists. Although this nested structure may appear complex at first, it offers a practical way to organise the large number of pairwise enrichment results produced by the method in a clear and consistent format. The object cont_all_BP4 illustrates this structure for the BP ontology at GO level 4:
[1] "tableList" "list" class(cont_all_BP4$sanger)[1] "list"class(cont_all_BP4$sanger$atlas)[1] "table"The class “equivSDhtestList” is also specific to goSorensen and corresponds to the typical output of the function equivTestSorensen. It extends the S3 class structure to store all hypothesis tests performed between pairs of gene lists for a given ontology and GO level. An object of this class is a list whose elements correspond to the analysed gene lists, and within each element are the results of the equivalence tests comparing that list with all others. Each of these inner elements is an object of class “equivSDhtest”, which contains the detailed results of a single pairwise test, including the test statistic, Sorensen dissimilarity, standard error, confidence interval, and p-value. The object eqTest_all_BP4 illustrates this structure for the BP ontology at GO level 4:
The class “distList” is also specific to goSorensen and corresponds to the output of the function allSorenThreshold. It extends the R class “list” to store the irrelevance-threshold matrices of dissimilarities obtained for different GO levels and ontologies. Each element of a “distList” object represents a particular GO ontology and level, and contains an object of class “dist” with the corresponding irrelevance-threshold matrix of dissimilarities. This structure provides an organised and consistent way to compare functional similarity between gene lists across multiple ontologies and levels. The object allDissMatrx illustrates this class, encompassing the three ontologies (BP, CC, and MF) and GO levels 3 through 10:
The goSorensen functions are the core components of the package. Figure 4 depicts these functions highlighted in turquoise, where the functions enclosed in the box outlined in turquoise with solid borders comprise the main functions of the package since they are directly related to the sequential flow of the goSorensen method depicted in Figure 2. Among these main functions are those devoted to the enrichment analysis, which exclusively accept objects of class “list” as input. Like most functions of goSorensen, these are implemented as generic functions under the S3 object-oriented programming paradigm, providing a flexible yet consistent framework for defining outputs, ensuring that results from the different stages of the method can be stored, accessed, and extended in a coherent and standardised way. For example, Figure 5 illustrates how the function buildEnrichTable generates different outputs based on the type of input provided by the user. If the input is a pair of vectors of class “character”, the function returns an object of class “table” containing the enrichment contingency table between these two compared gene lists. Whereas, if the input is an object of class “list” with two or more gene lists, the function returns an object of class “tableList” containing all the possible enrichment contingency tables between the compared gene lists.
Figure 5: buildEnrichTable: a generic function to compute the enrichment contingency table from two or more gene lists.
The remaining main functions of goSorensen complete the last step (statistical analysis) in the sequential flow depicted in Figure 2 to detect biological similarity through the goSorensen method. These functions can generate their outputs based on either the normal or the bootstrap sample distribution of the Sorensen dissimilarity, depending on the user’s choice.
As illustrated in Figure 4, these functions (highlighted in the turquoise dashed box) encompass three main tasks: estimating relevant statistics, performing the results for the equivalence hypothesis test to prove functional similarity, and computing the irrelevance-threshold matrix of dissimilarities.
These generic functions accept as input either objects of class “list” (raw data, i.e., gene lists) or objects of class “tableList” (enrichment contingency tables). While both inputs yield the same outputs, the computational efficiency is significantly enhanced when the user provides an object of class “tableList” directly, as the enrichment analysis step has been previously performed. Figure 6 illustrates this behaviour using the function equivTestSorensen, as an example. This figure shows that an object of class “equivSDhtest”, which contains the results of an equivalence test between two gene lists, can be generated using either two “character” vectors representing the gene lists or a “matrix” object containing the corresponding enrichment contingency table as input.
Moreover, when comparing more than two gene lists, the same function, equivTestSorensen, accepts either an object of class “list” containing all the gene lists to be compared or an object of class “tableList” comprising all the possible contingency tables between the lists as input. Regardless of the input, the output is an object of class equivSDhtestList, which includes the results of all pairwise hypothesis tests performed between the compared gene lists.
Figure 6: equivTestSorensen: a generic function to compute the outputs of the equivalence hypothesis tests to detect biological similarity between two or more gene lists.
In addition to the main functions already described, goSorensen includes two extended functions designed to operate across multiple GO ontologies and levels. The function allBuildEnrichTable generalises buildEnrichTable to generate enrichment contingency tables simultaneously for several GO levels and ontologies, producing objects of class “allTableList.” Similarly, the function allEquivTestSorensen extends equivTestSorensen to perform equivalence tests over multiple ontologies and levels, generating objects of class “AllEquivSDhtest”. These functions ensure that the package can scale the goSorensen method to more complex analyses, maintaining consistency and reproducibility across the full GO structure.
Finally, goSorensen includes a set of complementary functions designed to facilitate the management and interpretation of results. Accessor functions allow users to extract specific components from the outputs of equivalence tests, such as sample dissimilarities, p-values, standard errors, contingency tables, or confidence interval limits. These functions are particularly useful when handling large-scale analyses involving numerous pairwise tests, as they enable efficient retrieval and post-processing of relevant statistics, including adjustments for multiple testing. In addition, an upgrade function that enables users to update existing test results when modifying key parameters, such as the significance level, irrelevance threshold, or sampling distribution. This approach improves computational efficiency by recalculating only the affected components rather than repeating the full analysis. The package also provides plotting functions that generate objects to plot dendrograms to visualise hierarchical clustering of the irrelevance-threshold matrices of similarities.
According to Figure 4, the last component of the goSorensen package is its documentation files. In addition to the standard documentation that every R package provides, such as help pages, descriptions, illustrative examples, and a Reference Manual, goSorensen includes two detailed vignettes that serve as comprehensive user guides. The first vignette, ‘An Introduction to goSorensen R-Package’, offers a comprehensive overview of the package’s structure, its main classes and functions. It explains, through reproducible examples, how to perform the enrichment and equivalence analyses and how to interpret their results. The second vignette, ‘Working with the Irrelevance-threshold Matrix of Dissimilarities’, focuses on the statistical analyses that can be carried out using the dissimilarity matrices produced by the package, showing how these can be explored to study biological similarity across gene lists and GO levels. Together, these documents provide a clear and practical introduction to the package, bridging methodological concepts with applied examples and supporting reproducibility and transparency in its use
4 An example application
All examples in this section were generated using Bioconductor 3.23, the current Bioconductor release at the time of revision. Running the examples with a different Bioconductor release may yield different numerical results because the underlying annotation resources, including GO.db and org.Hs.eg.db, evolve between releases. A reproducibility analysis comparing Bioconductor releases 3.21, 3.22 and 3.23, together with the corresponding code and computational environments, is available at https://github.com/ASPresearch/goSorensen-Bioc-Reproducibility.
The objective of this section is to illustrate the application of the main functions of goSorensen and, consequently, how to obtain the different classes of objects described in the previous section. For this purpose, the package includes the object allOncoGeneLists, a “list” containing seven gene lists derived from a compilation of data carried out by the Bushman Lab on the human species (https://www.med.upenn.edu/bushmanlab/resources.html). These lists represent sets of transcribed genes with functional evidence of expression at different stages of cancer progression.
Before applying the goSorensen functions to the data contained in allOncoGeneLists, it is necessary to define a vector representing the gene universe available in the annotation databases for the species under study. In the case of Homo sapiens, this vector can be obtained from the annotation package org.Hs.eg.db as follows:
The object humanEntrezIDs is a vector of class “character” containing the ENTREZ identifiers that define the gene universe for the human species. The ENTREZ ID is a system of gene identifiers recognised by most biological databases, including GO, and will be used as the reference identifier throughout this section. However, the gene universe can also be represented using other identifier systems, such as ENSEMBL IDs, GENE SYMBOL IDs, or NCBI IDs.
4.1 Enrichment analysis
Enrichment matrix
The first step of the analysis is to compute the enrichment matrix. For example, for the ontology BP, at GO level 4, this can be done as follows::
enrichedInBP4 <- enrichedIn(allOncoGeneLists,
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
onto = "BP", GOLevel = 4) By default, this matrix includes only GO terms enriched in at least one list. Setting onlyEnriched = FALSE (the default is TRUE) returns the full matrix, and the attribute nTerms indicates the total number of GO terms evaluated. Since the full matrix is mostly composed of non-enriched terms, using the reduced version is generally more practical.
Analysing a large number of gene lists, which would involve repeated enrichment analyses, could potentially create a computational bottleneck. In this case, it might be advisable to use the enrichedIn argument parallel = TRUE. This is certainly not the case in the present example – parallelisation, particularly on Windows, can result in a certain amount of computational overhead, especially during initialisation. All functions whose execution ultimately depends on enrichedIn have this argument, which is set to FALSE by default.
Enrichment contingency table
Based on the total number of GO terms annotated in the selected ontology and level, the enrichment contingency tables are obtained as the cross-frequency between pairs of lists from the enrichment matrix. For two specific lists (e.g., atlas and sanger) and a given ontology-level combination (e.g., BP, level 4), the function buildEnrichTable internally computes the corresponding enrichment matrix to generate the contingency table between both lists. If the argument storeEnrichedIn is set to TRUE (value by default), the enrichment matrix is accessible in the attribute enriched. The code routine to obtain this example contingency table is as follows:
cont_atlas.sanger_BP4 <- buildEnrichTable(allOncoGeneLists$atlas,
allOncoGeneLists$sanger,
listNames = c("atlas", "sanger"),
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
onto = "BP", GOLevel = 4)
cont_atlas.sanger_BP4 Enriched in sanger
Enriched in atlas TRUE FALSE
TRUE 204 186
FALSE 27 2964To obtain all pairwise tables at once, one can run:
cont_all_BP4 <- buildEnrichTable(allOncoGeneLists,
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
onto = "BP", GOLevel = 4)and for several ontologies and GO levels:
allContTabs <- allBuildEnrichTable(allOncoGeneLists,
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
ontos = c("BP", "CC", "MF"),
GOLevels = 3:10)4.2 Statistical analysis
Hypothesis tests
The next step is to evaluate functional similarity through equivalence testing on the Sorensen dissimilarity (\(H_0: d_S \geq d_0\) Vs. \(H_1: d_S < d_0\)). The irrelevance threshold \(d_0\) defines the maximum dissimilarity at which two lists are still considered functionally equivalent. Following the approach of Flores et al. (2022), \(d_0\) can be derived from a ratio \(\rho = 2p_{11}/(p_{01}+p_{10})\) of joint versus marginal enrichment probability, which can be expressed as:
\[ d_0 = \cfrac{1}{1+\rho}. \]
Following conventional criteria of robustness in statistical model building, which accept variations of up to \(20\%\) as practically irrelevant (Box 1979), and in line with the standard bioequivalence range of \([80\%\) \(–\) \(125\%]\) adopted by international regulatory guidelines (FDA 1992; EMA 2010), Flores et al. (2022) propose a ratio \(\rho = 1/0.8=1.25\). This value represents the upper bound of acceptable deviation and corresponds to an equivalence limit of \(d_0=0.4444\) in the context of the goSorensen method.
For illustration:
eqTest_atlas.sanger_BP4 <- equivTestSorensen(allOncoGeneLists$atlas,
allOncoGeneLists$sanger,
listNames = c("atlas", "sanger"),
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
onto = "BP", GOLevel = 4,
d0 = 0.4444, conf.level = 0.95)eqTest_atlas.sanger_BP4
Normal asymptotic test for 2x2 contingency tables based on the
Sorensen-Dice dissimilarity
data: cont_atlas.sanger_BP4
(d - d0) / se = -4.5899, p-value = 2.217e-06
alternative hypothesis: true equivalence limit d0 is less than 0.4444444
95 percent confidence interval:
0.0000000 0.3793505
sample estimates:
Sorensen dissimilarity
0.3429952
attr(,"se")
standard error
0.02210248 In this example, the test does not reject the null hypothesis of non-equivalence at \(\alpha = 0.05\). Therefore, there is insufficient statistical evidence to conclude that the atlas and sanger lists are functionally equivalent under the threshold \(d0 = 0.4444\), for the BP ontology at GO level 4
Similarly, for a specific ontology and GO level, when all pairwise tests are needed, one can use:
eqTest_all_BP4 <- equivTestSorensen(allOncoGeneLists,
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
onto = "BP", GOLevel = 4,
d0 = 0.4444, conf.level = 0.95)If contingency tables are pre-computed (e.g., cont_all_BP4), the same output can be generated more efficiently:
eqTest_all_BP4 <- equivTestSorensen(cont_all_BP4,
d0 = 0.4444, conf.level = 0.95)Across multiple ontologies and GO levels, the function allEquivTestSorensen generalises this procedure. Users can also apply bootstrap inference by setting boot = TRUE, or update existing results with the upgrade function without rerunning the full analysis.
See the vignette An Introduction to goSorensen R-Package for more details about the use of the functions described above.
Irrelevance-threshold matrix of dissimilarities
Each element of this matrix represents the irrelevance-threshold dissimilarity for the corresponding pair of gene lists; equivalence is statistically supported for irrelevance thresholds larger than this value. For the ontology BP at GO level 4:
dissMatrx_BP4 <- sorenThreshold(allOncoGeneLists,
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
onto = "BP", GOLevel = 4,
trace = FALSE)
dissMatrx_BP4 atlas cis miscellaneous sanger Vogelstein
cis 0.7135865
miscellaneous 0.5667299 0.6286086
sanger 0.4606465 0.6589836 0.3384844
Vogelstein 0.3461926 0.6896145 0.5109424 0.2489031
waldman 0.3051649 0.6742598 0.2756020 0.4153017 0.3744204Using pre-computed tables (e.g., cont_all_BP4) makes this step faster. Multiple ontologies and levels can be processed simultaneously using:
allDissMatrx <- allSorenThreshold(allOncoGeneLists,
geneUniverse = humanEntrezIDs,
orgPackg = "org.Hs.eg.db",
ontos = c("BP", "CC", "MF"),
GOLevels = 3:10)These matrices are obtained using the normal approach; if one wants to use the bootstrap approach instead, the argument boot = TRUE, which is set to FALSE by default, can be set.
Hierarchical clustering of these matrices can be performed directly with the hclustThreshold function, which extends the standard hclust framework to objects generated by goSorensen. For “distList” objects containing matrices across several ontologies and GO levels, the function allHclustThreshold generalises this procedure. The resulting clustering objects can be visualised using the standard plot method. In contrast, MDS can be performed directly using standard R functions, since the methodological contribution of goSorensen lies in the construction of the irrelevance-threshold matrix of dissimilarities rather than in introducing a new multidimensional scaling procedure. Worked examples for both visualisations, including all the R code required to reproduce the dendrogram and MDS biplot, are provided in the vignette Working with the Irrelevance-threshold Matrix of Dissimilarities.
Although the irrelevance-threshold matrix of dissimilarities, which can be conveniently visualised through dendrograms or MDS biplots, provides an alternative way of assessing statistical equivalence between gene lists from that provided by p-values, both results are available to users and can be examined together. For any irrelevance-threshold matrix computed by goSorensen, such as dissMatrx_BP4, the enrichment contingency tables used in its computation are retained in the all2x2Tables attribute. These tables can therefore be recovered and used to perform the corresponding equivalence tests without repeating the enrichment analysis. The accessor functions provided by goSorensen can then be used to retrieve the p-values, as well as other outputs of the equivalence tests. For example, the pairwise p-values associated with dissMatrx_BP4 can be obtained as follows:
# Extract the contingency tables
contTable <- attr(dissMatrx_BP4, "all2x2Tables")
# Get the p-values
pvals <- getPvalue(equivTestSorensen(contTable))
Figure 7: MDS – biplot and dendrogram representing the clustered lists resulting from the irrelevance-threshold matrix of dissimilarities generated for the lists included in allOncoGeneLists, for the ontology BP and the GO level 4.
Figure 7 presents the two-dimensional MDS biplot and the dendrogram for the matrix contained in the object dissMatrx_BP4. Since these dissimilarities have an inferential basis, the MDS plot and dendrogram provide useful exploratory representations of patterns of functional similarity. However, proximity or clustering in these visualisations should not by itself be interpreted as evidence of pairwise statistical equivalence; the corresponding equivalence tests should be considered for inferential conclusions. For example, we can see that the lists waldman, vogelstein, sanger, and miscellaneous are minimally distanced from each other, which suggests a group of functionally similar lists; the corresponding pairwise equivalence tests can then be examined to assess the statistical evidence supporting these relationships.
The abscissa and ordinate axes of the MDS - biplot represent the two most representative dimensions that capture the variability from the original space. The characterisation of these dimensions constitutes an interesting component of the statistical analysis, as it allows us to detect which GO terms are associated with the formation of these dimensions. In other words, a characterisation of the biplot dimensions enables us to detect which biological concepts contribute to the patterns of functional similarity represented in the biplot.
The main idea behind the characterisation of the MDS - Biplot dimensions is that a GO term with as different and stable enrichment as possible in the lists located at the extremes of dimensions helps explain the formation of clusters.
The following interactive table presents the results obtained by applying this characterisation to the object dissMatrx_BP4:
For specific details of the irrelevance-threshold matrix of dissimilarities, its clustering and visualisation, and the identification of the most relevant GO terms explaining biological similarity through the characterisation of MDS biplot dimensions, readers are referred to Chapter 5 of Flores (2025) and to the vignette Working with the Irrelevance-threshold Matrix of Dissimilarities.
4.3 Interactive visualisation dashboard
Since it is not feasible to observe all the results depicted in this section, we have created an interactive Shiny application to help users explore all of them. This dashboard, included as supplementary material, is designed exclusively for visualisation purposes and is not intended to serve as an additional analysis tool for goSorensen. It is accessible directly at https://pablof1988.shinyapps.io/app_allOnco/, and the source code with a brief guide to the application is available from the GitHub repository located at https://github.com/pablof1988/allOncoGeneLists_viewer
Figure 8: Screenshots of the tabs in the Shiny app. An interactive dashboard to explore the results of applying the goSorensen functions to the gene lists contained in the object allOncoGeneLists.
Figure 8 illustrates interactively the screenshots of the tabs in the Shiny application. The application is organised into three main tabs following the sequential steps of the goSorensen method:
GENE LISTS displays the ENTREZ identifiers of the selected lists.
ENRICHMENT ANALYSIS shows the enrichment matrix for all lists and the contingency table for the selected pair of lists, both adjusted to the chosen ontology and GO level.
STATISTICAL ANALYSIS presents in two subtabs the results of the equivalence hypothesis tests and the irrelevance-threshold matrix of dissimilarities, including visualisations such as dendrograms and MDS biplots. Users can dynamically modify key parameters (e.g., irrelevance limit \(d_0\), confidence level, or sampling distribution) to update the visual outputs.
Additionally, the app includes a tool to identify the GO terms contributing to the biological similarity observed in each MDS dimension.
5 Conclusions and discussions
This article introduces goSorensen, an R package developed to implement the goSorensen method for assessing biological similarity between gene lists. The package integrates enrichment analysis, equivalence testing, and visualisation, allowing the entire process to be carried out in a coherent and reproducible way, promoting flexible workflows and facilitating integration with other Bioconductor packages. Built under the object-oriented programming paradigm, goSorensen allows the user to manipulate and combine objects of different classes within its functions. The statistical calculations used by goSorensen to determine significant equivalence between lists are mainly based on a normal approximation to the sampling distribution of the Sorensen dissimilarity, which quantifies their functional similarity. Although this approximation provides an efficient and straightforward approach for inference, its accuracy may decrease at low enrichment levels, where a bootstrap alternative—also implemented in the package—offers more reliable results but at the expense of a higher computational cost.
Since the goSorensen package is a direct implementation of the goSorensen method, future work is closely linked to improving, optimising, and scaling this approach. A first potential improvement is to develop an analysis procedure that does not depend on selecting a specific GO level. Because GO terms within the same level differ in their degree of specificity, difficulties may arise when functional similarity is inferred from overly general or uninformative concepts. In some cases, a term may fail to appear enriched at a general level even though enrichment is evident at more specific levels—or conversely, enrichment may only be observed for broader ancestral terms. Developing a procedure that integrates information across levels would therefore provide a more consistent representation of biological similarity, reducing the influence of arbitrary level selection.
A limitation of the current implementation arises from the components of the Sorensen dissimilarity. Enriched GO terms—either joint or marginal—carry an inherent level of uncertainty associated with their detection. As is well known, the assessment of GO term enrichment involves a certain degree of arbitrariness since it depends on criteria established by the researcher’s perception, such as cut-off points. In addition, the identification of enrichment is itself the outcome of a hypothesis test, which depends on earlier decisions, such as the choice of significance levels or the method adopted to address multiple testing. Developing adaptive approaches that explicitly model the uncertainty derived from the enrichment stage could control a possible overall inflation of type I error.
Another important challenge is related to the management of large volumes of information. When the analysis involves a great number of genes or lists, the resulting amount of data can substantially increase computational demand. In some situations, the required resources may exceed the capacity of standard computing equipment, limiting the practical applicability of the method. Future developments should therefore consider strategies to reduce the computational burden.
Since goSorensen relies on external databases to provide biological meaning, another direction for future work is to explore the behaviour of the method using additional resources beyond the Gene Ontology, for example, the Disease Ontology, the Drug Ontology, or even databases based on biological pathways, such as REACTOME or KEGG. Furthermore, the conceptual foundation underlying goSorensen can be generalised to other similarity measures beyond the Sorensen–Dice index, such as the Jaccard coefficient or other metrics to assess joint enrichment overlap. In addition, the framework can be extended to other omics domains, such as proteomics, metabolomics, or transcriptomics, where assessing the equivalence between feature lists is essential for interpreting the biological knowledge associated with more complex studies.
In summary, goSorensen offers a statistically grounded and reproducible framework for evaluating functional equivalence between gene lists within R and Bioconductor. Its continued development will aim to refine the underlying method, enhance computational efficiency, and expand its applicability to other biological domains. In this way, goSorensen will continue to contribute to the integration of rigorous statistical methodology with practical tools for interpreting complex biological data.
Acknowledgements
This work was supported by the Agencia Estatal de Investigación (AEI), grant PID2023-148013OB-C22, and by the Instituto de Salud Carlos III through CIBERFES (CB16/10/00269).
5.1 Supplementary materials
Supplementary materials are available in addition to this article. It can be downloaded at RJ-2026-051.zip
5.2 Bioconductor packages used
goSorensen, goProfiles, clusterProfiler
At this point, it may be important to reiterate the distinction we make between the term (commonly used in bioinformatics) gene list, represented in goSorensen by an object of class “character”, and the class “list” , which is used here to designate several gene lists.↩︎