library(goSorensen)

Introduction.

The majority of functions of goSorensen are devoted to implementing a statistical hypothesis test to detect equivalence between two or more gene lists. This method was 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.

Following this method, we develop the irrelevance-threshold matrix of dissimilarities \(\mathfrak{D}\):

\[\begin{equation} \begin{aligned} & \begin{matrix} \hspace{0.6cm} L_1 & L_2 & \hspace{0.1cm} L_3 & \hspace{0.1cm} \ldots \hspace{0.1cm} & L_{s-1} \end{matrix} \\ \mathfrak{D} = \begin{matrix} L_2 \\ L_3 \\ L_4 \\ \vdots \\ L_s \end{matrix} & \begin{pmatrix} \mathfrak{d}_{21} & & & & \\ \mathfrak{d}_{31} & \mathfrak{d}_{32} & & & \\ \mathfrak{d}_{41} & \mathfrak{d}_{42} & \mathfrak{d}_{43} & & \\ \vdots & \vdots & \vdots & & \\ \mathfrak{d}_{s1} & \mathfrak{d}_{s2} & \mathfrak{d}_{s3} & \ldots & \mathfrak{d}_{s(s-1)} \\ \end{pmatrix} \end{aligned} \label{it-diss} \end{equation}\]

The goal is to determine a dissimilarity measure \(\mathfrak{d}_{ij}\) for comparing two lists (\(L_i, L_j\)). This dissimilarity is not only a descriptive measure, but it is based on the irrelevance threshold determining whether two lists (\(L_i, L_j\)) are equivalent, which ensures that the dissimilarity \(\mathfrak{d}_{ij}\) is directly associated with the statistically significant equivalence declaration for a specific ontology (BP, MF, or CC) and GO level.

The dissimilarity measures between all the compared \(s\) feature lists can be collectively obtained in the symmetric matrix \(\mathfrak{D}\) using the following algorithm:

  • Step 1. Establish \(h = s(s-1)\) equivalence hypothesis tests. Each test \(ET_I\) (with \(I = 1, \ldots, h\)) compares some specific pair of lists \(L_{i}, L_{j}\) (with \(i, j = 1, 2, \ldots, s\)).

  • Step 2. Let \(\mathfrak{d}_h\) be the smallest irrelevance limit that makes all the null hypotheses corresponding to the \(h\) equivalence tests be rejected (All the \(s\) lists are equivalent): \[\mathfrak{d}_h = \text{min} \left(\mathfrak{d}: P_v(\mathfrak{d})_{(I)} \leq \frac{\alpha}{h+1-I}; \hspace{0.5cm} I=1, 2, \ldots,h \right).\]

  • Step 3. Obtain \(\mathfrak{d}_h\) and use it as the irrelevance limit between the lists \(L_i, L_j\) corresponding to the last position in a vector of ordered p-values \(P_v(\mathfrak{d}_{h})_{(I)}\), i.e. \(\mathfrak{d}_{ij} = \mathfrak{d}_{h}.\)

  • Step 4. Set \(h = h-1\), excluding the comparison between the lists \(L_i, L_j\) from Step \(3\), and go back to Step \(2\) until \(h=0\).

Data used in this vignette.

We use a specific database to provide a clear and precise explanation of the computation, visualization, and interpretation of the \(\mathfrak{D}\) matrix. Specifically, we will employ the set of PBTs lists, which is available as a resource in the goSorensen package. Access to this database can be obtained as follows:

data("pbtGeneLists")

PBTs is an object of class list, containing 14 character vectors with the ENTREZ gene identifiers of a gene list related to cell graft rejection:

sapply(pbtGeneLists, length)
ABMR_RATs       BAT       CT1     ENDAT     GRIT2     GRIT3    IRITD1    IRITD3 
      108       105       948       118       131       205       181       321 
   IRITD5       KT1       LT1       LT3  Rej_RATs TCMR_RATs 
      215       562       913       693       107       121

1 Performing the Irrelevance - Threshold Matrix of Dissimilarities.

1.1 For an Specific Ontology and GO level.

One method to compute the dissimilarity matrix \(\mathfrak{D}\) for a given ontology and level (e.g., ontology BP, level 5) is to directly input the gene lists into the sorenThreshold function:

# Previously load the genomic annotation package for the studied specie:
library(org.Hs.eg.db)
humanEntrezIDs <- keys(org.Hs.eg.db, keytype = "ENTREZID")

# Computing the irrelevance-threshold matrix of dissimilarities
dismatBP5 <- sorenThreshold(pbtGeneLists, onto = "BP", GOLevel = 5,
                            geneUniverse = humanEntrezIDs, 
                            orgPackg = "org.Hs.eg.db")

 Wed Jul 31 17:39:37 2024   Building all enrichment contingency tables...

 Wed Jul 31 17:43:04 2024   Computing all threshold equivalence distances...
round(dismatBP5, 2)
          ABMR_RATs  BAT  CT1 ENDAT GRIT2 GRIT3 IRITD1 IRITD3 IRITD5  KT1  LT1
BAT            0.98                                                           
CT1            1.00 1.00                                                      
ENDAT          0.89 1.00 1.00                                                 
GRIT2          0.71 0.97 1.00  0.97                                           
GRIT3          0.47 0.98 1.00  0.97  0.35                                     
IRITD1         1.00 1.00 1.00  0.99  1.00  1.00                               
IRITD3         0.92 1.00 1.00  0.79  0.96  0.96   1.00                        
IRITD5         0.87 1.00 1.00  0.88  0.90  0.87   1.00   0.90                 
KT1            1.00 1.00 0.81  0.95  1.00  1.00   1.00   1.00   1.00          
LT1            1.00 1.00 1.00  1.00  1.00  1.00   1.00   1.00   1.00 1.00     
LT3            1.00 1.00 1.00  1.00  1.00  1.00   1.00   1.00   1.00 1.00 0.29
Rej_RATs       0.32 0.98 1.00  0.95  0.64  0.38   1.00   0.94   0.85 1.00 1.00
TCMR_RATs      0.52 0.96 1.00  0.98  0.58  0.42   1.00   0.93   0.86 1.00 1.00
           LT3 Rej_RATs
BAT                    
CT1                    
ENDAT                  
GRIT2                  
GRIT3                  
IRITD1                 
IRITD3                 
IRITD5                 
KT1                    
LT1                    
LT3                    
Rej_RATs  1.00         
TCMR_RATs 1.00     0.39

Each value in this matrix represents dissimilarities determined by the equivalence threshold which determines whether the compared lists are equivalent. In addition, it is essential to consider that when the Sorensen dissimilarity value is irrelevantly equal to zero (equivalent to zero), it suggests biological similarity. Based on this, we can deduce that the smallest values (near zero) in the matrix \(\mathfrak{D}\) indicate equivalence between lists, whereas the largest values (close to 1) indicate that the lists cannot be considered equivalent.

Another way is to compute all the pairwise \(2 \times 2\) enrichment contingency tables previously. From this table, we can obtain the dissimilarity matrix as follows:

# Previously compute the 2x2 contingency tables:
ctableBP5 <- buildEnrichTable(pbtGeneLists, onto = "BP", GOLevel = 5,
                            geneUniverse = humanEntrezIDs, 
                            orgPackg = "org.Hs.eg.db")

# Computing the irrelevance-threshold matrix of dissimilarities
dismatBP5 <- sorenThreshold(ctableBP5)

The result is exactly the same matrix od dissimilarities above showed.

1.2 For More than One Ontology and GO Level Defined by the User.

Similarly, we can directly perform the computation using the gene lists dataset:

allDismat <- allSorenThreshold (pbtGeneLists, geneUniverse = humanEntrezIDs, 
                            orgPackg = "org.Hs.eg.db", ontos = c("BP", "CC"), 
                            GOLevels = 4:7)

Or, previously computing all the pairwise contingency tables for the ontologies (BP, CC, MF) and GO levels defined by the user:

allTabs <- allBuildEnrichTable(pbtGeneLists, geneUniverse = humanEntrezIDs, 
                            orgPackg = "org.Hs.eg.db", ontos = c("BP", "CC"), 
                            GOLevels = 4:7)
allDismat <- allSorenThreshold(allTabs)

Given the substantial quantity of outcomes obtained, exhibiting all of these results is not feasible. The result is a collection of dissimilarity matrices for the different ontologies and GO levels defined in the ontos and GOLevels arguments.

Suppose the ontos and GOLevels arguments are not provided. In that case, the allSorenThreshold function automatically computes the dissimilarity matrices for all three ontologies (BP, CC and MF) and GO levels ranging from 3 to 10.

2 Representing the Dissimilarity Matrix in Statistic Graphs.

2.1 Dendrogram.

Below, we graph a dendrogram based on the irrelevance-threshold matrix of dissimilarities for the particular case of the MF ontology, GO level 5, in our PBTs dataset. The dissimilarity matrix for this particular case was computed above and it is saved in the object dismatBP5.

clust.threshold <- hclustThreshold(dismatBP5)
plot(clust.threshold)

The dendrogram displays the lists of genes that are most similar and most distant. For instance, when considering a dissimilarity of 0.6, chosen arbitrarily, we can observe the formation of three groups. Within each group, the lists can be regarded as similar but different from those in other groups. There are three groups of lists that can be considered equivalent. The first group includes LT3 and LT1. The second group contains GRIT2 and GRIT3. The third group comprises the lists TCMR_RATs, ABMR_RATs, and Rej_RATs.

2.2 MDS-Biplot.

We initiate the process by performing a multidimensional scaling MDS analysis. MDS transforms the initial dissimilarity matrix into Euclidean distances while maintaining the maximum possible variation of the original dissimilarities. In this instance, we compute the first two dimensions, which account for the most relevant proportion of the original variability.

mds <- as.data.frame(cmdscale(dismatBP5, k = 2))
mds
                   V1          V2
ABMR_RATs  0.38854334  0.04392676
BAT       -0.16327259 -0.09079187
CT1       -0.26862731 -0.24698285
ENDAT     -0.16612448 -0.28420320
GRIT2      0.35204607  0.07836793
GRIT3      0.45673252  0.10300851
IRITD1    -0.22078619 -0.12783291
IRITD3    -0.13587210 -0.24607399
IRITD5     0.01069298 -0.13659312
KT1       -0.27602226 -0.27344008
LT1       -0.41625933  0.50991290
LT3       -0.41625933  0.50991290
Rej_RATs   0.43787856  0.08027021
TCMR_RATs  0.41733012  0.08051881

Below we plot the biplot:

library(ggplot2)
library(ggrepel)
graph <- ggplot() +
  geom_point(aes(mds[,1], mds[,2]), color = "blue", size = 3) +
  geom_text_repel(aes(mds[,1], mds[,2], label = attr(dismatBP5, "Labels")), 
                  color = "black", size = 3) +
    xlab("Dim 1") +
    ylab("Dim 2") +
    theme_light()
graph

The MDS-Biplot suggests that the formation of groups is consistent with what was observed in the dendrogram. It is important to remember that the distances observed between lists in the biplot are defined based on the dissimilarity threshold that establishes the equivalence of the compared lists. For instance, the proximity between the Rej_RATs and GRIT3 lists implies that both are equivalent.

However, it is more interesting to identify the GO terms most strongly associated with the dimensions generated in the biplot rather than just determining if the lists are equivalent. These GO terms elucidate the biological functions associated with the observed distances between lists, explaining the formation of groups and the biological reasons for detecting equivalences between lists. In the next section, we develop a 5-step procedure for identifying these GO terms:

3 Characterizing MDS-Biplot Dimensions.

The following 5-step procedure helps identify the GO terms that biologically elucidate the creation of groups in each dimension of the Biplot. In this specific instance, we will execute the technique on Dimension 1. However, the approach remains identical if we were to explain Dimension 2.

  • STEP 1: Split the axis (Dimension 1) in three parts to identify the lists located at the extremes.
# Split the axis 20% to the left, 60% to the middle and 20% to the right:
prop <- c(0.2, 0.6, 0.2) 
# Sort according  dimension 1:
sorted <- mds[order(mds[, 1]), ] 
# Determine the range for dimension 1.
range <- sorted[, 1][c(1, nrow(mds))] 
# Find the cutpoints to split the axis:
cutpoints <- cumsum(prop) * diff(range) + range[1]
cutpoints <- cutpoints[-length(cutpoints)]

# Identify lists to the left:  
lleft <- rownames(sorted[sorted[, 1] < cutpoints[1], ])
# Identify lists to the right
lright <- rownames(sorted[sorted[, 1] > cutpoints[2], ]) 

lleft
[1] "LT1" "LT3" "KT1" "CT1"
lright
[1] "GRIT2"     "ABMR_RATs" "TCMR_RATs" "Rej_RATs"  "GRIT3"    
graph +
  geom_vline(xintercept = cutpoints, color = "red", linetype = "dashed", 
             linewidth = 0.75)

The basic idea is that if the enrichment of a GO term in these two groups of lists is different, then we can assume that this GO term discriminates the groups of lists, consequently explaining the formation of groups in that particular dimension.

  • STEP 2: Identify the enriched and non-enriched GO terms for each GO term and each list selected at the extremes, differentiating the lists by the groups to which they belong.
tableleft <- enrichedIn(pbtGeneLists[lleft], onto = "BP", GOLevel = 5,
                            geneUniverse = humanEntrezIDs, 
                            orgPackg = "org.Hs.eg.db")

tableright <- enrichedIn(pbtGeneLists[lright], onto = "BP", GOLevel = 5,
                            geneUniverse = humanEntrezIDs, 
                            orgPackg = "org.Hs.eg.db")

Due to the large size of the results, we only show the first 15 rows of each matrix:

tableleft[1:15, ]
             LT1   LT3   KT1   CT1
GO:0000055 FALSE FALSE FALSE FALSE
GO:0000056 FALSE FALSE FALSE FALSE
GO:2000200 FALSE FALSE FALSE FALSE
GO:2000201 FALSE FALSE FALSE FALSE
GO:2000202 FALSE FALSE FALSE FALSE
GO:1903137 FALSE FALSE FALSE FALSE
GO:1903138 FALSE FALSE FALSE FALSE
GO:1903139 FALSE FALSE FALSE FALSE
GO:0051229 FALSE FALSE FALSE FALSE
GO:0051232 FALSE FALSE FALSE FALSE
GO:0090306 FALSE FALSE FALSE FALSE
GO:1990395 FALSE FALSE FALSE FALSE
GO:0001578  TRUE  TRUE FALSE FALSE
GO:0007020 FALSE FALSE FALSE FALSE
GO:0007051 FALSE FALSE FALSE FALSE
tableright[1:15, ]
           GRIT2 ABMR_RATs TCMR_RATs Rej_RATs GRIT3
GO:0000055 FALSE     FALSE     FALSE    FALSE FALSE
GO:0000056 FALSE     FALSE     FALSE    FALSE FALSE
GO:2000200 FALSE     FALSE     FALSE    FALSE FALSE
GO:2000201 FALSE     FALSE     FALSE    FALSE FALSE
GO:2000202 FALSE     FALSE     FALSE    FALSE FALSE
GO:1903137 FALSE     FALSE     FALSE    FALSE FALSE
GO:1903138 FALSE     FALSE     FALSE    FALSE FALSE
GO:1903139 FALSE     FALSE     FALSE    FALSE FALSE
GO:0051229 FALSE     FALSE     FALSE    FALSE FALSE
GO:0051232 FALSE     FALSE     FALSE    FALSE FALSE
GO:0090306 FALSE     FALSE     FALSE    FALSE FALSE
GO:1990395 FALSE     FALSE     FALSE    FALSE FALSE
GO:0001578 FALSE     FALSE     FALSE    FALSE FALSE
GO:0007020 FALSE     FALSE     FALSE    FALSE FALSE
GO:0007051 FALSE     FALSE     FALSE    FALSE FALSE
  • STEP 3: Compute the means \(\overline{X}_{GO_j}\) and variances \(S^2_{GO_j}\) for each of the GO terms analyzed at each extreme group:
lmnsd <- apply(tableleft, 1, 
                 function(x){c("meanLeft" = mean(x), "sdLeft" = sd(x))})
rmnsd <- apply(tableright, 1, 
                 function(x){c("meanRight" = mean(x), "sdRight" = sd(x))})
  • STEP 4: Establish a statistic to assess the degree of enrichment disparity between extreme groups for each GO term. This measure should consider that a GO term discriminates a dimension if: i) the lists show opposing patterns of enrichment in the GO term, indicated by a large difference between \(| _L\overline{X}_{GO_j} - _R \overline{X}_{GO_j}|\) (\(L\) represents Left and \(R\) Rigth), and ii) the behavior of the lists remains stable at each extreme, indicated by a small amount of \(_LS^2_{GO_j} + _RS^2_{GO_j}\).

Considering these specific details, the following statistic, is suggested:

\[st_j = \dfrac{| _L\overline{X}_{GO_j} - _R \overline{X}_{GO_j}|}{\sqrt{\dfrac{_LS^2_{GO_j}}{l} + \dfrac{_RS^2_{GO_j}}{r} + \epsilon }}\] with \(l\) the number of lists to the extreme left and \(r\) the number of lists to the extreme right. Using \(\epsilon= 0.00000001\) the statistic is computed as follows:

nl <- ncol(tableleft) 
nr <- ncol(tableright) 
stat <- abs(lmnsd[1, ] - rmnsd[1, ]) / 
  sqrt((((lmnsd[2, ] / nl) + (rmnsd[2, ] / nr))) + 0.00000001)
  • STEP 5: Select the GO terms with the highest statistic (those with a value as close as possible to \(1/\sqrt{\epsilon}\) ). These selected GO terms mainly contribute to forming groups in the analyzed dimension.
result <- stat[stat == max(stat)]
result
GO:0001819 GO:0002718 GO:0032649 GO:0002720 GO:0032729 GO:0001910 GO:0042267 
     10000      10000      10000      10000      10000      10000      10000 
GO:0002758 GO:0002460 GO:0002819 GO:0002821 GO:0002366 GO:0002285 GO:0001909 
     10000      10000      10000      10000      10000      10000      10000 
GO:0002449 GO:0002703 GO:0002704 GO:0002705 GO:0002685 GO:0002697 GO:1903706 
     10000      10000      10000      10000      10000      10000      10000 
GO:0050777 GO:0002696 GO:0002699 GO:0051250 GO:0051251 GO:0002700 GO:0002702 
     10000      10000      10000      10000      10000      10000      10000 
GO:0002757 GO:0002768 GO:0002832 GO:0002833 GO:0045088 GO:0045824 GO:0045089 
     10000      10000      10000      10000      10000      10000      10000 
GO:0060326 GO:0031348 GO:0070665 GO:0051607 GO:0002237 GO:0050900 GO:0002478 
     10000      10000      10000      10000      10000      10000      10000 
GO:0022407 GO:0045785 GO:0032102 GO:0050727 GO:0070663 GO:0002228 GO:0002274 
     10000      10000      10000      10000      10000      10000      10000 
GO:0022409 GO:0030098 GO:0046651 GO:0002474 GO:0050867 GO:0031341 GO:0032943 
     10000      10000      10000      10000      10000      10000      10000 
GO:0071219 GO:0007159 GO:0032496 GO:1902105 
     10000      10000      10000      10000

Finally, we can identify the specific biological functions associated with the detected GO terms.

library(GO.db)
library(knitr)
# Previous function to get the description of the identified GO terms:
get_go_description <- function(go_id) {
  go_term <- Term(GOTERM[[go_id]])
  return(go_term)
}

# GO terms description:
kable(data.frame(Description = sapply(names(result), get_go_description, 
                                             USE.NAMES = TRUE)))
Description
GO:0001819 positive regulation of cytokine production
GO:0002718 regulation of cytokine production involved in immune response
GO:0032649 regulation of type II interferon production
GO:0002720 positive regulation of cytokine production involved in immune response
GO:0032729 positive regulation of type II interferon production
GO:0001910 regulation of leukocyte mediated cytotoxicity
GO:0042267 natural killer cell mediated cytotoxicity
GO:0002758 innate immune response-activating signaling pathway
GO:0002460 adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains
GO:0002819 regulation of adaptive immune response
GO:0002821 positive regulation of adaptive immune response
GO:0002366 leukocyte activation involved in immune response
GO:0002285 lymphocyte activation involved in immune response
GO:0001909 leukocyte mediated cytotoxicity
GO:0002449 lymphocyte mediated immunity
GO:0002703 regulation of leukocyte mediated immunity
GO:0002704 negative regulation of leukocyte mediated immunity
GO:0002705 positive regulation of leukocyte mediated immunity
GO:0002685 regulation of leukocyte migration
GO:0002697 regulation of immune effector process
GO:1903706 regulation of hemopoiesis
GO:0050777 negative regulation of immune response
GO:0002696 positive regulation of leukocyte activation
GO:0002699 positive regulation of immune effector process
GO:0051250 negative regulation of lymphocyte activation
GO:0051251 positive regulation of lymphocyte activation
GO:0002700 regulation of production of molecular mediator of immune response
GO:0002702 positive regulation of production of molecular mediator of immune response
GO:0002757 immune response-activating signaling pathway
GO:0002768 immune response-regulating cell surface receptor signaling pathway
GO:0002832 negative regulation of response to biotic stimulus
GO:0002833 positive regulation of response to biotic stimulus
GO:0045088 regulation of innate immune response
GO:0045824 negative regulation of innate immune response
GO:0045089 positive regulation of innate immune response
GO:0060326 cell chemotaxis
GO:0031348 negative regulation of defense response
GO:0070665 positive regulation of leukocyte proliferation
GO:0051607 defense response to virus
GO:0002237 response to molecule of bacterial origin
GO:0050900 leukocyte migration
GO:0002478 antigen processing and presentation of exogenous peptide antigen
GO:0022407 regulation of cell-cell adhesion
GO:0045785 positive regulation of cell adhesion
GO:0032102 negative regulation of response to external stimulus
GO:0050727 regulation of inflammatory response
GO:0070663 regulation of leukocyte proliferation
GO:0002228 natural killer cell mediated immunity
GO:0002274 myeloid leukocyte activation
GO:0022409 positive regulation of cell-cell adhesion
GO:0030098 lymphocyte differentiation
GO:0046651 lymphocyte proliferation
GO:0002474 antigen processing and presentation of peptide antigen via MHC class I
GO:0050867 positive regulation of cell activation
GO:0031341 regulation of cell killing
GO:0032943 mononuclear cell proliferation
GO:0071219 cellular response to molecule of bacterial origin
GO:0007159 leukocyte cell-cell adhesion
GO:0032496 response to lipopolysaccharide
GO:1902105 regulation of leukocyte differentiation

Session information.

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

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.20-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] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] knitr_1.48           GO.db_3.19.1         ggrepel_0.9.5       
##  [4] ggplot2_3.5.1        org.Hs.eg.db_3.19.1  AnnotationDbi_1.67.0
##  [7] IRanges_2.39.2       S4Vectors_0.43.2     Biobase_2.65.0      
## [10] BiocGenerics_0.51.0  goSorensen_1.7.1     BiocStyle_2.33.1    
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.2.3               gson_0.1.0              shadowtext_0.1.4       
##   [4] gridExtra_2.3           rlang_1.1.4             magrittr_2.0.3         
##   [7] DOSE_3.31.2             compiler_4.4.1          RSQLite_2.3.7          
##  [10] png_0.1-8               vctrs_0.6.5             reshape2_1.4.4         
##  [13] stringr_1.5.1           pkgconfig_2.0.3         crayon_1.5.3           
##  [16] fastmap_1.2.0           magick_2.8.4            XVector_0.45.0         
##  [19] labeling_0.4.3          ggraph_2.2.1            utf8_1.2.4             
##  [22] HDO.db_0.99.1           rmarkdown_2.27          enrichplot_1.25.0      
##  [25] UCSC.utils_1.1.0        tinytex_0.52            purrr_1.0.2            
##  [28] bit_4.0.5               xfun_0.46               zlibbioc_1.51.1        
##  [31] cachem_1.1.0            aplot_0.2.3             GenomeInfoDb_1.41.1    
##  [34] jsonlite_1.8.8          blob_1.2.4              highr_0.11             
##  [37] BiocParallel_1.39.0     tweenr_2.0.3            parallel_4.4.1         
##  [40] R6_2.5.1                RColorBrewer_1.1-3      bslib_0.8.0            
##  [43] stringi_1.8.4           jquerylib_0.1.4         GOSemSim_2.31.0        
##  [46] Rcpp_1.0.13             bookdown_0.40           goProfiles_1.67.0      
##  [49] Matrix_1.7-0            splines_4.4.1           igraph_2.0.3           
##  [52] tidyselect_1.2.1        viridis_0.6.5           qvalue_2.37.0          
##  [55] yaml_2.3.10             codetools_0.2-20        lattice_0.22-6         
##  [58] tibble_3.2.1            plyr_1.8.9              treeio_1.29.0          
##  [61] withr_3.0.1             KEGGREST_1.45.1         evaluate_0.24.0        
##  [64] gridGraphics_0.5-1      CompQuadForm_1.4.3      scatterpie_0.2.3       
##  [67] polyclip_1.10-7         Biostrings_2.73.1       ggtree_3.13.0          
##  [70] pillar_1.9.0            BiocManager_1.30.23     clusterProfiler_4.13.0 
##  [73] ggfun_0.1.5             generics_0.1.3          tidytree_0.4.6         
##  [76] munsell_0.5.1           scales_1.3.0            glue_1.7.0             
##  [79] lazyeval_0.2.2          tools_4.4.1             data.table_1.15.4      
##  [82] fgsea_1.31.0            fs_1.6.4                graphlayouts_1.1.1     
##  [85] fastmatch_1.1-4         tidygraph_1.3.1         cowplot_1.1.3          
##  [88] grid_4.4.1              ape_5.8                 tidyr_1.3.1            
##  [91] colorspace_2.1-1        nlme_3.1-165            patchwork_1.2.0        
##  [94] GenomeInfoDbData_1.2.12 ggforce_0.4.2           cli_3.6.3              
##  [97] fansi_1.0.6             viridisLite_0.4.2       dplyr_1.1.4            
## [100] gtable_0.3.5            yulab.utils_0.1.5       sass_0.4.9             
## [103] digest_0.6.36           ggplotify_0.1.2         farver_2.1.2           
## [106] memoise_2.0.1           htmltools_0.5.8.1       lifecycle_1.0.4        
## [109] httr_1.4.7              bit64_4.0.5             MASS_7.3-61