What hcocena does

hcocena builds co-expression networks and reads biology off their structure. The idea is that genes which vary together across samples tend to belong to the same programme, so a gene-gene correlation network, cut into communities, gives you a handful of modules instead of thousands of individual genes. Each module can then be characterised – which conditions drive it, which functions it is enriched for, which genes sit at its centre.

What makes it hCoCena (“horizontal CoCena”) is that several datasets can be analysed as one network. Different tissues, platforms or cohorts each contribute their own correlation network; those are merged, and the modules are found on the merged graph. A module is therefore a programme that exists across layers, and its behaviour can be compared between them.

This vignette runs the whole analysis end to end on a small simulated dataset, so every number and figure below is produced when the vignette is built. The steps that need an external service – ChEA3, Enrichr, DoRothEA, Cytoscape, an LLM endpoint – cannot run here; see Beyond this vignette for the full templates that cover them.

Installation

hcocena is distributed through Bioconductor and is installed with BiocManager:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}
BiocManager::install("hcocena")

Several optional steps rely on packages listed under Suggests (for example limma, decoupleR, RCy3 or ellmer); the functions that need them say so and name the package to install.

library(hcocena)
library(MultiAssayExperiment)
library(SummarizedExperiment)

The example data

Two layers of 40 genes and 32 samples each. The genes are real marker symbols for four blood cell types, and the simulation makes each set co-vary as a block, so the analysis should recover four modules that correspond to the four cell types. That is the point of using a simulation here: we know the answer in advance and can check that the pipeline finds it.

The samples are 16 donors measured at two timepoints, split across three conditions, with a batch label and an age covariate – enough structure for the statistics later on. inst/scripts/make-fixtures.R in the installed package contains the generating code.

extdir <- normalizePath(system.file("extdata", package = "hcocena"), winslash = "/")
counts_preview <- read.delim(file.path(extdir, "toy_layer1_counts.tsv"), nrows = 3)
counts_preview[, 1:5]
#>   SYMBOL L1_D01_T1 L1_D01_T2 L1_D02_T1 L1_D02_T2
#> 1   CD3D    22.303    25.599    26.800    22.430
#> 2   CD3E    23.631    25.702    27.998    24.200
#> 3   CD3G    22.243    26.985    25.278    21.879

anno_preview <- read.delim(file.path(extdir, "toy_layer1_anno.tsv"))
head(anno_preview, 3)
#>    SampleID   group batch donor timepoint age
#> 1 L1_D01_T1 control    B1   D01        T1  22
#> 2 L1_D01_T2    mild    B1   D01        T2  22
#> 3 L1_D02_T1  severe    B2   D02        T1  41

Setting up the analysis

Every analysis starts from an empty HCoCenaExperiment and is configured by passing it through a series of functions. The object carries the data, the settings and every intermediate result, so there is no hidden state to keep track of.

outdir <- file.path(tempdir(), "hcocena-vignette-output")
dir.create(outdir, recursive = TRUE, showWarnings = FALSE)

hc <- hc_init()
hc <- hc_set_paths(
  hc,
  dir_count_data = paste0(extdir, "/"),
  dir_annotation = paste0(extdir, "/"),
  dir_reference_files = paste0(extdir, "/"),
  dir_output = outdir
)
hc <- hc_init_save_folder(hc, name = "")
#> Using output directory directly (no additional save subfolder): /tmp/Rtmp1kBIuz/hcocena-vignette-output

hc_define_layers() names the layers and says which files belong to them. The names are used in file names and plot labels from here on.

hc <- hc_define_layers(
  hc,
  data_sets = list(
    Layer1 = c("toy_layer1_counts.tsv", "toy_layer1_anno.tsv"),
    Layer2 = c("toy_layer2_counts.tsv", "toy_layer2_anno.tsv")
  )
)

hc <- hc_read_data(
  hc,
  sep_counts = "\t",
  sep_anno = "\t",
  gene_symbol_col = "SYMBOL",
  sample_col = "SampleID",
  count_has_rn = FALSE,
  anno_has_rn = FALSE,
  auto_setup_output = FALSE
)

Global settings

variable_of_interest is the annotation column that defines the groups everything is later compared across; control_keyword names the reference level. range_GFC caps the group-fold-change scale introduced below.

hc <- hc_set_global_settings(
  hc,
  organism = "human",
  control_keyword = "control",
  variable_of_interest = "group",
  min_nodes_number_for_network = 1,
  min_nodes_number_for_cluster = 1,
  range_GFC = 2,
  layout_algorithm = "layout_with_stress",
  data_in_log = FALSE
)

Layer settings

top_var is the number of most variable genes per layer that enter the correlation. Genes that barely move across samples carry no co-expression signal but do add noise and computation, so the network is built on the variable ones. In a real analysis this is typically a few thousand; here the dataset only has 40 genes, so we take all of them.

min_corr and range_cutoff_length define the grid of correlation cutoffs that the next step scans.

hc <- hc_set_layer_settings(
  hc,
  top_var = c(40, 40),
  min_corr = c(0.5, 0.5),
  range_cutoff_length = c(50, 50),
  print_distribution_plots = c(FALSE, FALSE)
)

Part I: correlation and choosing a cutoff

hc_run_expression_analysis_1() computes all pairwise gene-gene correlations per layer and then, for every candidate cutoff on the grid, records what the resulting network would look like.

hc <- hc_run_expression_analysis_1(hc, corr_method = "pearson")

Why a cutoff has to be chosen at all. A correlation matrix is a complete graph – every gene is connected to every other one, just with different weights. To get a network with structure you keep only edges above some correlation. Set it too low and everything is connected to everything, so there are no modules to find; set it too high and the network falls apart into fragments and you lose most of your genes.

hcocena picks the cutoff by two criteria at once:

Scale-free topology (R²). Biological networks tend to have few highly connected hubs and many sparsely connected genes, a degree distribution that follows a power law. R² measures how well the network at a given cutoff fits that expectation. Higher is more biologically plausible.
How much network is left. The number of remaining nodes and the mean degree. A cutoff with a beautiful R² over eleven surviving genes is useless.

hc <- hc_set_cutoff(hc, auto = TRUE)

The statistics behind the choice can be inspected per layer – R², node count and mean degree along the cutoff grid:

hc_plot_cutoffs(hc, interactive = FALSE)

#> An object of class "HCoCenaExperiment"
#> Slot "mae":
#> A MultiAssayExperiment object of 2 listed
#>  experiments with user-defined names and respective classes.
#>  Containing an ExperimentList class object of length 2:
#>  [1] set1: SummarizedExperiment with 40 rows and 32 columns
#>  [2] set2: SummarizedExperiment with 40 rows and 32 columns
#> Functionality:
#>  experiments() - obtain the ExperimentList instance
#>  colData() - the primary/phenotype DataFrame
#>  sampleMap() - the sample coordination DataFrame
#>  `$`, `[`, `[[` - extract colData columns, subset, or experiment
#>  *Format() - convert into a long or wide DataFrame
#>  assays() - convert ExperimentList to a SimpleList of matrices
#>  exportClass() - save data to flat files
#> 
#> Slot "config":
#> An object of class "HCoCenaConfig"
#> Slot "global":
#> DataFrame with 1 row and 9 columns
#>   save_folder    organism     control         voi min_nodes_number_for_network
#>   <character> <character> <character> <character>                    <numeric>
#> 1                   human     control       group                            1
#>   min_nodes_number_for_cluster range_GFC   layout_algorithm data_in_log
#>                      <numeric> <numeric>        <character>   <logical>
#> 1                            1         2 layout_with_stress       FALSE
#> 
#> Slot "layer":
#> DataFrame with 2 rows and 9 columns
#>      layer_id  layer_name          count_source   annotation_source   top_var
#>   <character> <character>           <character>         <character> <numeric>
#> 1        set1      Layer1 toy_layer1_counts.tsv toy_layer1_anno.tsv        40
#> 2        set2      Layer2 toy_layer2_counts.tsv toy_layer2_anno.tsv        40
#>    min_corr range_cutoff_length print_distribution_plots    cutoff
#>   <numeric>           <numeric>                <logical> <numeric>
#> 1       0.5                  50                    FALSE     0.604
#> 2       0.5                  50                    FALSE     0.809
#> 
#> Slot "paths":
#> DataFrame with 1 row and 4 columns
#>           dir_count_data         dir_annotation    dir_reference_files
#>              <character>            <character>            <character>
#> 1 /tmp/RtmpHAEeNr/Rins.. /tmp/RtmpHAEeNr/Rins.. /tmp/RtmpHAEeNr/Rins..
#>               dir_output
#>              <character>
#> 1 /tmp/Rtmp1kBIuz/hcoc..
#> 
#> 
#> Slot "references":
#> An object of class "HCoCenaReferences"
#> Slot "registry":
#> DataFrame with 0 rows and 0 columns
#> 
#> Slot "data":
#> List of length 0
#> 
#> 
#> Slot "layer_results":
#> List of length 2
#> names(2): set1 set2
#> 
#> Slot "integration":
#> An object of class "HCoCenaIntegration"
#> Slot "combined_edgelist":
#> DataFrame with 0 rows and 0 columns
#> 
#> Slot "graph":
#> NULL
#> 
#> Slot "gfc":
#> DataFrame with 0 rows and 0 columns
#> 
#> Slot "cluster":
#> List of length 0
#> 
#> 
#> Slot "satellite":
#> List of length 1
#> names(1): cutoff_selection
#> 
#> Slot "provenance":
#> DataFrame with 1 row and 2 columns
#>               created_at          source
#>              <character>     <character>
#> 1 2026-10-04 00:55:58... legacy_hcobject

Part II: group fold changes

hc_run_expression_analysis_2() filters each layer’s network to the chosen cutoff and computes group fold changes (GFCs).

A GFC summarises, for one gene and one group of samples, how far that group sits from the overall mean of the gene, on a scale capped at range_GFC. It is what makes modules comparable across layers: absolute expression differs between platforms and tissues, but “up in severe relative to this gene’s own average” is the same quantity everywhere.

hc <- hc_run_expression_analysis_2(hc, plot_HM = FALSE)

Integration

Each layer now has its own filtered network. hc_build_integrated_network() merges them into one graph.

mode = "u" (union) keeps an edge present in any layer – the inclusive choice, and the usual one.
mode = "i" (intersection) keeps only edges present in every layer, which is stricter and needs the layers to measure largely the same genes.
multi_edges says what to do when an edge exists in several layers with different weights: "min" is conservative, "mean" and "max" are the alternatives.

hc <- hc_build_integrated_network(hc, mode = "u", multi_edges = "min")

g <- hc_graph(hc)
c(nodes = igraph::vcount(g), edges = igraph::ecount(g))
#> nodes edges 
#>    40   180

Module detection

Communities are found with the Leiden algorithm. resolution controls granularity: higher values give more, smaller modules.

Leiden is stochastic, so hcocena runs it no_of_iterations times and keeps the consensus. A gene ends up in the module it was assigned to most often, and genes that never settle are dropped – which is why a module assignment from this package is more stable than a single Leiden run.

hc <- hc_cluster_calculation(
  hc,
  cluster_algo = "cluster_leiden",
  no_of_iterations = 5,
  resolution = 1
)

gtc <- hc_gene_to_cluster(hc)
table(gtc$color)
#> 
#>      coral       gold lightgreen  steelblue 
#>         10         10         10         10

Four modules of ten genes, as the simulation intended. At this point the modules are identified by colour; the readable M1-style labels are assigned when the module heatmap is drawn, which also shows their GFC profile across the groups of both layers:

hc <- hc_plot_cluster_heatmap(hc, file_name = FALSE)

Module membership per gene is the table most downstream steps work from:

module_genes <- as.data.frame(hc_satellite(hc, "module_gene_list"))
split(module_genes$genes, module_genes$module)[["M1"]]
#>  [1] "SPON2"  "FGFBP2" "KLRD1"  "KLRF1"  "NKG7"   "GZMB"   "GNLY"   "KLRC1" 
#>  [9] "NCR1"   "PRF1"

Functional enrichment

Modules are tested for over-representation of gene sets. In a normal analysis these come from GO, KEGG, Hallmark or Reactome via gene_sets =; here we write a small GMT of the four marker panels instead, so that the vignette stays offline and self-contained.

The background matters. universe = "all_genes" (the default) tests each module against every measured gene, answering “is this module enriched relative to the transcriptome?”. universe = "network" tests against the genes that entered the network – the genes that could have landed in a module at all – answering “is this module enriched relative to the other modules?”. The second removes the bias that top-variance selection and the correlation cutoff introduce, and is the stricter question.

gmt <- file.path(tempdir(), "celltype_markers.gmt")
writeLines(c(
  paste(c("T cells", "-", "CD3D", "CD3E", "CD3G", "CD2", "CD28",
          "LCK", "ZAP70", "IL7R", "CD7", "TRAC"), collapse = "\t"),
  paste(c("B cells", "-", "CD19", "MS4A1", "CD79A", "CD79B", "BLNK",
          "PAX5", "CR2", "FCRL1", "TNFRSF13B", "VPREB3"), collapse = "\t"),
  paste(c("Monocytes", "-", "CD14", "LYZ", "FCN1", "VCAN", "S100A8",
          "S100A9", "CSF1R", "ITGAM", "CD68", "FCGR3A"), collapse = "\t"),
  paste(c("NK cells", "-", "NKG7", "GNLY", "KLRD1", "KLRF1", "PRF1",
          "GZMB", "NCR1", "KLRC1", "FGFBP2", "SPON2"), collapse = "\t")
), gmt)

hc <- hc_functional_enrichment(
  hc,
  gene_sets = character(0),
  custom_gmt_files = c(CellTypes = gmt),
  universe = "network"
)
#> Warning in calculate_qvalue(ora_res$pvalue): qvalue::qvalue() failed, returning
#> NA for qvalue. Error: missing values and NaN's not allowed if 'na.rm' is FALSE
#> Warning in calculate_qvalue(ora_res$pvalue): qvalue::qvalue() failed, returning
#> NA for qvalue. Error: missing values and NaN's not allowed if 'na.rm' is FALSE
#> Warning in calculate_qvalue(ora_res$pvalue): qvalue::qvalue() failed, returning
#> NA for qvalue. Error: missing values and NaN's not allowed if 'na.rm' is FALSE
#> Warning in calculate_qvalue(ora_res$pvalue): qvalue::qvalue() failed, returning
#> NA for qvalue. Error: missing values and NaN's not allowed if 'na.rm' is FALSE

enr <- hc_satellite(hc, "enrichments")$significant_enrichments_all_dbs
as.data.frame(enr)[, c("module_label", "term", "GeneRatio", "qvalue")]
#>   module_label      term GeneRatio       qvalue
#> 1           M1  NK cells     10/10 1.179718e-09
#> 2           M2   B cells     10/10 1.179718e-09
#> 3           M3   T cells     10/10 1.179718e-09
#> 4           M4 Monocytes     10/10 1.179718e-09

Each module recovers exactly one cell type with all ten of its markers, which is the answer the simulation was built to have.

Characterising the modules

Which conditions drive a module

hc_module_condition_significance() reduces each module to one value per sample – the mean expression of its genes – and tests that against the conditions. Wilcoxon (or Kruskal-Wallis for more than two groups) is rank-based, limma fits a linear model with pairwise contrasts, and an optional mixed model handles repeated measures via a donor random intercept.

standardize_modules = TRUE z-scores each module within a layer first. Raw means carry the module’s absolute expression level, which distorts limma’s variance moderation and makes effect sizes incomparable between modules. Applied per layer it also removes the offset between layers, which matters because the model formula itself carries no layer term.

hc <- hc_module_condition_significance(
  hc,
  condition_col = "group",
  run_limma = FALSE,
  standardize_modules = TRUE,
  export_excel = FALSE
)

sig <- hc_satellite(hc, "module_condition_significance")$wilcox
as.data.frame(sig)[, c("cluster", "test", "p", "p_adj")]
#>      cluster             test         p     p_adj
#> 1       gold kruskal_fallback 0.8188524 0.8188524
#> 2      coral kruskal_fallback 0.6230004 0.8188524
#> 3 lightgreen kruskal_fallback 0.2910916 0.8188524
#> 4  steelblue kruskal_fallback 0.5062805 0.8188524

Nothing is significant here, and that is the correct answer: the simulation assigned the condition labels at random, so there is no group effect to find.

Correlating modules with metadata

Numeric covariates are correlated with module expression directly; categorical ones are correlated across the groups of the variable of interest, which needs at least three of them.

hc_meta_correlation_num(hc, set = 1, meta = "age")

Module GFC tables

The numbers behind the module heatmap, per group and per sample:

hc <- hc_module_gfc_tables(hc, export_xlsx = FALSE)

head(as.data.frame(hc_satellite(hc, "module_gfc_tables")$by_group), 4)
#>   module       mild     severe
#> 1     M2  0.1953234  0.4271183
#> 2     M4  0.2121004  1.0714531
#> 3     M1  1.0583028 -0.8238035
#> 4     M3 -0.4097938 -0.2087299

Hub genes

Hubs are the genes most central to their own module, by weighted degree, closeness and betweenness. They are the natural candidates when a module has to be reduced to a few representatives.

hc <- hc_find_hubs(hc, top = 5)

vapply(hc_satellite(hc, "hub_out"), function(x) paste(x[[1]], collapse = ", "),
       character(1))
#>                                  M2                                  M4 
#> "CD79A, VPREB3, FCRL1, CD79B, BLNK"  "VCAN, CD68, S100A8, S100A9, CD14" 
#>                                  M1                                  M3 
#> "SPON2, NKG7, KLRC1, FGFBP2, KLRD1"       "CD3D, LCK, CD28, ZAP70, CD7"

The hubs of each module are exactly the markers of the corresponding cell type, which is another way of seeing that the modules are clean.

Module scores

A module score is the fraction of a gene’s neighbours that lie inside its own module – a measure of how cleanly the community was separated. Scores near 1 mean the module is well isolated in the graph.

hc <- hc_get_module_scores(hc)

The object that carries all of this

Everything above is stored in the same HCoCenaExperiment. The assays live in a MultiAssayExperiment:

names(experiments(hc_mae(hc)))
#> [1] "set1" "set2"
vapply(experiments(hc_mae(hc)), function(x) dim(x)[[1]], integer(1))
#> set1 set2 
#>   40   40

Settings stay attached, so an analysis can be reconstructed from the object alone:

as.data.frame(methods::slot(hc_config(hc), "global"))[, c("organism", "voi", "range_GFC")]
#>   organism   voi range_GFC
#> 1    human group         2

Downstream results are collected in the satellite slot, which hc_satellite() returns:

names(hc_satellite(hc))
#> [1] "cutoff_selection"              "module_gene_list"             
#> [3] "module_gfc_means"              "enrichments"                  
#> [5] "module_condition_significance" "module_gfc_tables"            
#> [7] "hub_out"                       "module_scores"

Beyond this vignette

Several parts of hcocena need services that cannot be reached while a vignette is being built, and are therefore not shown above:

transcription-factor enrichment through ChEA3 (hc_tf_overrep_module(), hc_check_tf()),
cell-type annotation from Enrichr libraries (hc_celltype_annotation(), hc_celltype_activity_decoupler()),
upstream regulator and pathway inference via CollecTRI/DoRothEA/PROGENy (hc_upstream_inference()),
network layout round-trips with Cytoscape (hc_export_to_cytoscape()),
LLM-based module interpretation (hc_llm_enrichment()) through the API of Claude, OpenAI, Gemini or a local vLLM server – or, without an API key, by pasting the request into any provider’s web chat and reading the answer back with hc_llm_import().

The complete walkthroughs are installed with the package as R Markdown templates:

dir(system.file("scripts", "workflows", package = "hcocena"), pattern = "[.]Rmd$")
#> [1] "hcocena_main.Rmd"      "hcocena_satellite.Rmd"

hcocena_main.Rmd is the full analysis; hcocena_satellite.Rmd collects the optional downstream steps. Copy one next to your data and edit the paths at the top:

file.copy(
  system.file("scripts", "workflows", "hcocena_main.Rmd", package = "hcocena"),
  "hcocena_main.Rmd"
)

They are templates rather than reproducible documents: they point at your own count and annotation files and contact the services listed above.

Citing hcocena

The method – per-layer co-expression networks, group fold changes, cutoff selection by scale-free topology, and integration of the layer networks – is described in Oestreich et al. (2022), Bioinformatics 38(20):4727-4734, doi:10.1093/bioinformatics/btac589. A step-by-step protocol for a full analysis is given in Holsten et al. (2024), STAR Protocols 5(1):102922, doi:10.1016/j.xpro.2024.102922. Both papers describe the earlier function names without the hc_ prefix; each of those functions has an hc_ counterpart here (for example functional_enrichment() is now hc_functional_enrichment()).

citation("hcocena")
#> To cite hcocena, please use the method paper and, where the protocol
#> was followed, the STAR Protocols paper:
#> 
#>   Oestreich M, Holsten L, Agrawal S, Dahm K, Koch P, Jin H, Becker M,
#>   Ulas T (2022). hCoCena: horizontal integration and analysis of
#>   transcriptomics datasets. Bioinformatics 38(20):4727-4734.
#>   doi:10.1093/bioinformatics/btac589
#> 
#>   Holsten L, Dahm K, Oestreich M, Becker M, Ulas T (2024). hCoCena: A
#>   toolbox for network-based co-expression analysis and horizontal
#>   integration of transcriptomic datasets. STAR Protocols 5(1):102922.
#>   doi:10.1016/j.xpro.2024.102922
#> 
#> To see these entries in BibTeX format, use 'print(<citation>,
#> bibtex=TRUE)', 'toBibtex(.)', or set
#> 'options(citation.bibtex.max=999)'.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] MultiAssayExperiment_1.39.1 SummarizedExperiment_1.43.0
#>  [3] Biobase_2.73.2              GenomicRanges_1.65.4       
#>  [5] Seqinfo_1.3.2               IRanges_2.47.5             
#>  [7] S4Vectors_0.51.10           BiocGenerics_0.59.12       
#>  [9] generics_0.1.4              MatrixGenerics_1.25.0      
#> [11] matrixStats_1.5.0           hcocena_0.99.12            
#> [13] BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] RColorBrewer_1.1-3      sys_3.4.3               jsonlite_2.0.0         
#>   [4] shape_1.4.6.1           tidydr_0.0.6            magrittr_2.0.5         
#>   [7] ggtangle_0.1.3          magick_2.9.1            farver_2.1.2           
#>  [10] rmarkdown_2.32          GlobalOptions_0.1.4     fs_2.1.0               
#>  [13] vctrs_0.7.3             memoise_2.0.1           Cairo_1.7-0            
#>  [16] ggtree_4.3.0            htmltools_0.5.9         S4Arrays_1.13.2        
#>  [19] BiocBaseUtils_1.15.1    SparseArray_1.13.4      gridGraphics_0.5-1     
#>  [22] sass_0.4.10             bslib_0.12.0            htmlwidgets_1.6.4      
#>  [25] plyr_1.8.9              httr2_1.3.0             cachem_1.1.0           
#>  [28] buildtools_1.0.0        igraph_2.3.4            lifecycle_1.0.5        
#>  [31] iterators_1.0.14        pkgconfig_2.0.3         gson_0.2.1             
#>  [34] Matrix_1.7-6            R6_2.6.1                fastmap_1.2.0          
#>  [37] clue_0.3-68             digest_0.6.39           aplot_0.3.2            
#>  [40] enrichplot_1.99.7       colorspace_2.1-3        ggnewscale_0.5.2       
#>  [43] patchwork_1.3.2         AnnotationDbi_1.75.2    aisdk_1.4.12           
#>  [46] ps_1.9.3                RSQLite_3.53.3          labeling_0.4.3         
#>  [49] polyclip_1.10-7         httr_1.4.9              abind_1.4-8            
#>  [52] compiler_4.6.1          fontquiver_0.2.1        bit64_4.8.6            
#>  [55] withr_3.0.3             doParallel_1.0.17       S7_0.2.2               
#>  [58] DBI_1.3.0               ggforce_0.5.0           MASS_7.3-66            
#>  [61] rappdirs_0.3.4          DelayedArray_0.39.8     rjson_0.2.23           
#>  [64] ggsci_5.2.0             gtools_3.9.5            tools_4.6.1            
#>  [67] otel_0.2.0              scatterpie_0.2.6        ape_5.8-1              
#>  [70] zip_3.0.2               glue_1.8.1              callr_3.8.0            
#>  [73] nlme_3.1-171            GOSemSim_2.39.3         grid_4.6.1             
#>  [76] cluster_2.1.8.3         reshape2_1.4.5          gtable_0.3.6           
#>  [79] tidyr_1.3.2             data.table_1.18.6.1     XVector_0.53.0         
#>  [82] ggrepel_0.9.8           foreach_1.5.2           pillar_1.11.1          
#>  [85] stringr_1.6.0           yulab.utils_0.2.5       splines_4.6.1          
#>  [88] circlize_0.4.18         tweenr_2.0.3            dplyr_1.2.1            
#>  [91] treeio_1.37.1           lattice_0.23-1          bit_4.6.0              
#>  [94] tidyselect_1.2.1        fontLiberation_0.1.0    GO.db_3.23.1           
#>  [97] ComplexHeatmap_2.29.0   maketools_1.3.2         Biostrings_2.81.9      
#> [100] knitr_1.52              fontBitstreamVera_0.1.1 xfun_0.61              
#> [103] graphlayouts_1.2.5      leidenbase_0.1.37       stringi_1.8.9          
#> [106] lazyeval_0.2.3          ggfun_0.2.1             yaml_2.3.12            
#> [109] evaluate_1.0.5          codetools_0.2-20        qvalue_2.45.0          
#> [112] gdtools_0.5.1           tibble_3.3.1            BiocManager_1.30.27    
#> [115] ggplotify_0.1.3         cli_3.6.6               systemfonts_1.3.2      
#> [118] processx_3.9.0          jquerylib_0.1.4         Rcpp_1.1.2             
#> [121] png_0.1-9               parallel_4.6.1          ggplot2_4.0.3          
#> [124] blob_1.3.0              clusterProfiler_4.21.3  DOSE_4.7.3             
#> [127] viridisLite_0.4.3       tidytree_0.4.8          ggiraph_0.9.6          
#> [130] enrichit_0.2.5          rlist_0.4.6.2           scales_1.4.0           
#> [133] openxlsx_4.2.9          purrr_1.2.2             crayon_1.5.3           
#> [136] GetoptLong_1.1.1        rlang_1.3.0             cowplot_1.2.0          
#> [139] KEGGREST_1.53.6

hcocena: a complete worked analysis