# Advanced CSOA

# Introduction

This tutorial showcases some advanced functions of CSOA.

# Prerequisites

In addition to CSOA, you need to install patchwork, scRNAseq, scuttle and stringr for this tutorial.

# Refining gene signatures: I

First, we will repeat the setup from the CSOA tutorial. We need the `BaronPancreasData` dataset from `scRNAseq`, and we load the PanglaoDB-generated list of acinar markers.

    library(CSOA)
    library(ggplot2)
    library(henna)
    library(patchwork)
    library(scRNAseq)
    library(scuttle)
    library(stringr)
    library(Seurat)

    sceObj <- BaronPancreasData('human')
    sceObj <- logNormCounts(sceObj)

    seuratObj <- as.Seurat(sceObj)

    acinarMarkers <- c('PRSS1', 'KLK1', 'CTRC', 'PNLIP', 'AKR1C3', 'CTRB1',
                       'DUOXA2', 'ALDOB', 'REG3A', 'SERPINA3', 'PRSS3', 'REG1B',
                       'CFB',   'GDF15',    'MUC1','ANPEP', 'ANGPTL4', 'OLFM4',
                       'GSTA1', 'LGALS2', 'PDZK1IP1', 'RARRES2', 'CXCL17',
                       'UBD', 'GSTA2', 'LYZ', 'RBPJL', 'PTF1A', 'CELA3A',
                       'SPINK1', 'ZG16', 'CEL', 'CELA2A', 'CPB1', 'CELA1',
                       'PNLIPRP1', 'RNASE1', 'AMY2B', 'CPA2','CPA1', 'CELA3B',
                       'CTRB2', 'PLA2G1B', 'PRSS2', 'CLPS', 'REG1A', 'SYCN')

A distinctive feature of CSOA is that the calculation of per-cell scores does not necessarily use all the genes in the input gene signature. If a gene did not register any top overlap, it will not feature in the calculation of CSOA scores.

We now run CSOA with the option of calculating the contribution of individual gene pairs to the CSOA score toggled on using `pairFileName`:

    seuratObj <- runCSOA(seuratObj, list(CSOA_acinar = acinarMarkers),
                         pairFileTemplate='pairs')
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 194 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Saving pair scores file: pairsCSOA_acinar.qs2...
    #> Computing per-cell gene signature scores...

Next, we inspect the resulting data frame using `qGrab`, a wrapper around `qs_read` from the `qs2` package which deletes the file after reading it:

    acinarPairs <- qGrab('pairsCSOA_acinar.qs2')
    head(acinarPairs)
    #>      gene1 gene2 overlapScore overlapRank pairScore pairRank revCumsum
    #> 507   CPA1 PRSS2    0.8279889          22  1.643546        1 100.00000
    #> 662   CTRC PRSS2    0.8279889          22  1.403776        2  98.35645
    #> 815   KLK1 PRSS2    0.8886735          12  1.379073        3  96.95268
    #> 327 CELA2A PRSS2    0.7633825          34  1.332878        4  95.57361
    #> 977  PRSS1 PRSS2    0.5808436          82  1.309812        5  94.24073
    #> 572   CPB1 PRSS2    0.7633825          34  1.287813        6  92.93092

Using the `overlapGenes` function, we find that only `27` out of the `47` input genes contributed directly to the CSOA score:

    genes <- overlapGenes(acinarPairs)
    nGenes <- length(genes)
    nGenes
    #> [1] 27
    length(acinarMarkers)
    #> [1] 47
    setdiff(acinarMarkers, genes)
    #>  [1] "AKR1C3"   "DUOXA2"   "SERPINA3" "GDF15"    "MUC1"     "ANPEP"
    #>  [7] "ANGPTL4"  "OLFM4"    "LGALS2"   "PDZK1IP1" "CXCL17"   "UBD"
    #> [13] "LYZ"      "RBPJL"    "PTF1A"    "ZG16"     "CELA1"    "RNASE1"
    #> [19] "AMY2B"    "CPA2"

Next, we run CSOA with only the `27` genes as input:

    seuratObj <- runCSOA(seuratObj, list(CSOA_acinar2 = genes))
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 86 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Computing per-cell gene signature scores...

We aim to compare the results with those obtained using the original `47` acinar markers. To this end, we build a function that computes and prints several statistics of interest here—Pearson correlation, and mean, median, maximum and minimum of the absolute-value differences between the two sets of CSOA scores:

    printStats <- function(seuratObj, col1, col2){
        paste0('Pearson correlation: ', cor(seuratObj[[]][[col1]],
                                            seuratObj[[]][[col2]]))
        absDiffs <- abs(seuratObj[[]][[col1]] - seuratObj[[]][[col2]])
        absDiffsStats <- sapply(list(mean, median, max, min),
                                function(fun) fun(absDiffs))
        names(absDiffsStats) <- paste0(c('mean', 'median', 'max', 'min'),
                                       'AbsDiff')
        absDiffsStats
    }

The results are quite close to those obtained using all the `47` genes:

    printStats(seuratObj, 'CSOA_acinar', 'CSOA_acinar2')
    #>   meanAbsDiff medianAbsDiff    maxAbsDiff    minAbsDiff
    #>  0.0047680035  0.0001103177  0.1434652424  0.0000000000

It is instructive to think why the results between the two runs were not identical. This is because the additional overlaps computed for the `47` acinar markers impacted each individual component of the overlap rank (p-value rank, ratio rank), inducing changes in the overlap ranking and subsequent selection and scoring.

# Refining gene signatures: II

We can further refine gene signatures by eliminating overlaps which account only for a small percent of the total score, and subsequently setting the corresponding genes as new input genes.

We return to the `acinarPairs` data frame. The `revCumsum` column shows that `90%` of the total CSOA score is accounted for by the top-scoring `128` gene pairs:

    acinarPairsSub <- subset(acinarPairs, revCumsum > 10)
    nrow(acinarPairsSub)
    #> [1] 128

The `overlapGenes` function finds the `24` genes that correspond to these pairs:

    genes90 <- overlapGenes(acinarPairsSub)
    genes90
    #>  [1] "CPA1"     "CTRC"     "KLK1"     "CELA2A"   "PRSS1"    "CPB1"
    #>  [7] "PNLIP"    "CELA3B"   "PRSS2"    "CELA3A"   "PNLIPRP1" "CTRB1"
    #> [13] "CTRB2"    "REG1B"    "PLA2G1B"  "PRSS3"    "GSTA2"    "GSTA1"
    #> [19] "CLPS"     "CEL"      "SYCN"     "REG1A"    "REG3A"    "RARRES2"
    nGenes90 <- length(genes90)
    nGenes90
    #> [1] 24

Next, we run CSOA using solely these `24` genes as input genes. The results remain close to those obtained using all acinar markers:

    seuratObj <- runCSOA(seuratObj, list(CSOA_acinar3 = genes90))
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 97 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Computing per-cell gene signature scores...

    printStats(seuratObj, 'CSOA_acinar', 'CSOA_acinar3')
    #>   meanAbsDiff medianAbsDiff    maxAbsDiff    minAbsDiff
    #>  0.0068344918  0.0001244179  0.1725938747  0.0000000000

# Addressing noisy gene sets

This section will demonstrate CSOA’s ability to extract the most biologically meaningful genes from an input gene set. First, as random numbers are going to be involved, we set a seed to ensure reproducibility:

    set.seed(123)

We will proceed by replacing randomly chosen genes among the top `24` acinar markers with random genes from the Seurat object, and running CSOA using the resulting gene set.

First, we build a function for this task:

    runCSOAReplace <- function(seuratObj, markers, nReplacedGenes,
                               pairFileTemplate = 'pairs'){
        genesComplement <- setdiff(rownames(seuratObj), markers)
        geneSet <- list(c(sample(markers, length(markers) - nReplacedGenes),
                          sample(genesComplement, nReplacedGenes)))
        names(geneSet) <- paste0('CSOA_replace', nReplacedGenes)
        seuratObj <- runCSOA(seuratObj, geneSet, pairFileTemplate=pairFileTemplate)
        return(seuratObj)
    }

Now we run CSOA on a set in which half of top `24` acinar markers (randomly chosen) have been replaced by random genes. The results are very close to those obtained using the original `24` markers, though with an increase in the maximum absolute difference—CSOA no longer recognizes a very few acinar cells as such, because of the lower number of input acinar markers.

    nReplacedGenes <- ceiling(nGenes90 / 2)
    seuratObj <- runCSOAReplace(seuratObj, genes90, nReplacedGenes)
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 33 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Saving pair scores file: pairsCSOA_replace12.qs2...
    #> Computing per-cell gene signature scores...

    newCol <- paste0('CSOA_replace', nReplacedGenes)
    printStats(seuratObj, 'CSOA_acinar3', newCol)
    #>   meanAbsDiff medianAbsDiff    maxAbsDiff    minAbsDiff
    #>  0.0113040121  0.0002535866  0.2785675871  0.0000000000
    VlnPlot(seuratObj, newCol, group.by='label') + NoLegend()

The noise introduced by random genes did not directly affect CSOA results at all. Using the `overlapGenes` function, we can verify that none of the random genes participated in the `23` top overlaps:

    length(setdiff(overlapGenes(qGrab('pairsCSOA_replace12.qs2')), genes90))
    #> [1] 0

We now replace all but `2` of the top `24` acinar markers with random genes and run CSOA again. The mean, median and minimum absolute differences remain very small, and the correlation with the original results remains very high. However, the maximum absolute difference has become substantial, indicating that an additional small number of cells could no longer be detected as acinar cells. Noise now makes an effect, though a small one: a single random gene (`PTGR1`) introduces two overlaps into CSOA scoring, which together account for `27.8%` of the total CSOA score.

    nReplacedGenes <- nGenes90 - 2
    seuratObj <- runCSOAReplace(seuratObj, genes90, nReplacedGenes)
    #> Computing percentile sets...
    #> Warning in percentileSets(geneSetExp, percentile): 4 gene(s) had no top cells
    #> at the indicated percentile. These are now excluded from the gene signature.
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 3 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Saving pair scores file: pairsCSOA_replace22.qs2...
    #> Computing per-cell gene signature scores...

    newCol <- paste0('CSOA_replace', nReplacedGenes)
    printStats(seuratObj, 'CSOA_acinar3', newCol)
    #>   meanAbsDiff medianAbsDiff    maxAbsDiff    minAbsDiff
    #>    0.01975495    0.00000000    0.54736806    0.00000000
    dfReplace <- qGrab('pairsCSOA_replace22.qs2')
    dfReplace
    #>      gene1 gene2 overlapScore overlapRank pairScore pairRank revCumsum
    #> 70  CELA2A CTRB1    1.0000000           1  72.20135        1 100.00000
    #> 152  CTRB1 PTGR1    0.6201145           2  16.47713        2  27.79865
    #> 82  CELA2A PTGR1    0.6201145           2  11.32151        3  11.32151
    setdiff(overlapGenes(dfReplace), genes90)
    #> [1] "PTGR1"
    VlnPlot(seuratObj, newCol, group.by='label') + NoLegend()

These results have implications upon the construction of input gene sets for CSOA. Due to CSOA’s efficient noise-filtering ability, it is generally better to err on the side of including more genes in the input, even though this will incur higher computational costs.

# Operating with gene sets that characterize multiple subpopulations

Cell signatures do not necessarily characterize a single cell type or state. Therefore, it is desirable that a gene set analysis method is able to recognize the different cell types or states characterized by different genes in the input set.

To illustrate this functionality of CSOA, we will first construct a mixed signature using acinar and stellate cell markers selected from PanglaoDB:

    stellateMarkers <- c('COL6A1', 'RGS5', 'COL1A1', 'MMP11', 'THY1', 'COL6A3',
                         'SFRP2', 'COL1A2', 'TNFAIP6', 'TIMP3', 'SPARC', 'COL3A1',
                         'MGP', 'COL6A2', 'COL4A1', 'FN1', 'SPON2', 'TIMP1',
                         'TGFB1', 'INHBA', 'PDGFRA', 'NDUFA4L2', 'MMP14', 'CTGF',
                         'CYGB', 'KRT10', 'PDGFRB', 'DYNLT1', 'GEM')

    mixedMarkers <- union(acinarMarkers, stellateMarkers)

Now we run CSOA with the mixed markers as input. As expected, we see enrichment of CSOA scores among both acinar and stellate cells:

    seuratObj <- runCSOA(seuratObj, list(CSOA_mixed = mixedMarkers),
                         pairFileTemplate='pairs')
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 249 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Saving pair scores file: pairsCSOA_mixed.qs2...
    #> Computing per-cell gene signature scores...
    VlnPlot(seuratObj, 'CSOA_mixed', group.by='label') + NoLegend()

We can interpret the overlap data frame as a graph, with genes as vertices and overlaps as edges. We can build and visualize the graph as follows:

    dfNetwork <- qGrab('pairsCSOA_mixed.qs2')
    overlapNetworkPlot(dfNetwork, rankCol='overlapRank')
    #> Warning in element_text(hjust = 0.5, ...): `...` must be empty.
    #> ✖ Problematic argument:
    #> • rankCol = "overlapRank"

We remark that the graph is not connected: it has no vertex from which we can reach all other vertices by travelling on the graph edges. Instead, some subgraphs are connected. We refer to these as connected components. Their significance is that they can be interpreted as gene modules—groups of genes linked by the highly significant overlaps of some of their corresponding cell sets. Using the `connectedComponents` function from `henna`, we find that the graph contains two connected components, `1` and `2`:

    dfNetwork <- connectedComponents(dfNetwork)
    dplyr::count(dfNetwork, component)
    #>   component   n
    #> 1         1 133
    #> 2         2 116

Next, we are interested in scoring the major gene modules in order to ascertain the identity of the cells they characterize. We can score the major components separately as follows:

    components <- c(1, 2)
    seuratObj <- scoreModules(seuratObj, dfNetwork, components)
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 101 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Computing per-cell gene signature scores...
    #> 87 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Computing per-cell gene signature scores...
    p1 <- VlnPlot(seuratObj, 'CSOA_component1', group.by='label') + NoLegend() +
        theme(axis.text.y = element_blank())
    p2 <- VlnPlot(seuratObj, 'CSOA_component2', group.by='label') + NoLegend() +
        theme(axis.text.y = element_blank())
    p1 / p2

Thus, component `1` is an acinar one and component `2` is a stellate one. The genes defining each component can be extracted from the overlap matrix after `connectedComponents` has been run:

    genes1 <- overlapGenes(subset(dfNetwork, component == 1))
    genes1
    #>  [1] "CPA1"     "PNLIP"    "CTRC"     "CELA2A"   "CPB1"     "KLK1"
    #>  [7] "PRSS2"    "CELA3B"   "CELA3A"   "PRSS1"    "CTRB1"    "CTRB2"
    #> [13] "GSTA2"    "PLA2G1B"  "CLPS"     "GSTA1"    "PRSS3"    "REG1B"
    #> [19] "CEL"      "RARRES2"  "SYCN"     "REG1A"    "PNLIPRP1" "REG3A"
    #> [25] "SPINK1"
    genes2 <- overlapGenes(subset(dfNetwork, component == 2))
    genes2
    #>  [1] "COL6A2" "COL6A1" "COL1A1" "COL3A1" "COL1A2" "MMP14"  "COL6A3" "FN1"
    #>  [9] "SPARC"  "INHBA"  "CYGB"   "PDGFRB" "SFRP2"  "MMP11"  "SPON2"  "MGP"
    #> [17] "TIMP3"  "THY1"   "TIMP1"  "PDGFRA"

It may happen, though, that the network of overlaps generated by CSOA for a mixed signature is connected. Let us construct a mixed signature of acinar and ductal markers, run CSOA, and identity the connected components of the graph determined by the top overlaps:

    ductalMarkers <- c('CFTR', 'SERPINA5', 'SLPI', 'TFF1', 'CFB', 'LGALS4',
                       'CTSH',  'PERP', 'PDLIM3',   'WFDC2', 'SLC3A1', 'AQP1',
                       'ALDH1A3', 'VTCN1',  'KRT19', 'TFF2', 'KRT7', 'CLDN4',
                       'LAMB3', 'TACSTD2', 'CCL2', 'DCDC2','CXCL2', 'CLDN10',
                       'HNF1B', 'KRT20', 'MUC1', 'ONECUT1', 'AMBP', 'HHEX',
                       'ANXA4', 'SPP1', 'PDX1', 'SERPINA3', 'GDF15', 'AKR1C3',
                       'MMP7', 'DEFB1', 'SERPING1', 'TSPAN8', 'CLDN1', 'S100A10',
                       'PIGR')

    mixedMarkers <- union(acinarMarkers, ductalMarkers)
    seuratObj <- runCSOA(seuratObj, list(CSOA_mixed = mixedMarkers),
                         pairFileTemplate='pairs')
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 760 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Saving pair scores file: pairsCSOA_mixed.qs2...
    #> Computing per-cell gene signature scores...
    VlnPlot(seuratObj, 'CSOA_mixed', group.by='label') + NoLegend()

    dfNetwork <- qGrab('pairsCSOA_mixed.qs2')
    dfNetwork <- connectedComponents(dfNetwork)
    dplyr::count(dfNetwork, component)
    #>   component   n
    #> 1         1 760

The violin plot shows enrichment in both the acinar and ductal groups, but the overlap graph is connected. To divide the mixed signature into acinar and ductal components, we will run CSOA with a non-null value set for the `jaccardCutoff` parameter:

    seuratObj <- runCSOA(seuratObj, list(CSOA_mixed = mixedMarkers),
                         jaccardCutoff=1/3, pairFileTemplate='pairs')
    #> Computing percentile sets...
    #> Generating cell set overlaps...
    #> Processing overlaps...
    #> 760 overlaps have been selected for Jaccard-based filtering.
    #> 221 edges with low neighbor Jaccard  scores have been removed.
    #> 42 edges with low neighbor Jaccard  scores have been removed.
    #> 3 edges with low neighbor Jaccard  scores have been removed.
    #> 3 edges with low neighbor Jaccard  scores have been removed.
    #> 5 edges with low neighbor Jaccard  scores have been removed.
    #> 486 overlaps will be used in the calculation of CSOA scores.
    #> Normalizing expression matrix by rows...
    #> Computing per-cell scores for gene pairs...
    #> Saving pair scores file: pairsCSOA_mixed.qs2...
    #> Computing per-cell gene signature scores...

This function has reduced the number of overlaps from 760 to 486 by selectively removing overlaps for which the neighbors of the corresponding genes overlapped little, as calculated using the Jaccard score. The resulting graph has two connected components:

    dfNetwork <- qGrab('pairsCSOA_mixed.qs2')
    overlapNetworkPlot(dfNetwork, rankCol='overlapRank')
    #> Warning in element_text(hjust = 0.5, ...): `...` must be empty.
    #> ✖ Problematic argument:
    #> • rankCol = "overlapRank"

These correspond to acinar and ductal cells, respectively:

    dfNetwork <- connectedComponents(dfNetwork)
    components <- unique(dfNetwork$component)
    seuratObj <- scoreModules(seuratObj, dfNetwork, components)
    p1 <- VlnPlot(seuratObj, 'CSOA_component1', group.by='label') + NoLegend() +
        theme(axis.text.y=element_blank())
    p2 <- VlnPlot(seuratObj, 'CSOA_component2', group.by='label') + NoLegend() +
        theme(axis.text.y=element_blank())
    p1 / p2

# Visualizing gene participation in top overlaps

To display gene participation in top overlaps, CSOA provides the `geneRadialPlot` function:

    geneRadialPlot(acinarPairs,
                   title='Top 100 overlap genes - Acinar markers')

Additionally, we can visualize and compare the genes involved in the top overlaps across multiple gene sets. We will reload the alpha and gamma markers from the Getting started with CSOA tutorial, and run again CSOA with four gene sets. To facilitate visualization, we will select the top `15` overlaps for each gene set:

    alphaMarkers <- c('GCG', 'TTR', 'PCSK2', 'FXYD5', 'LDB2', 'MAFB',
                      'CHGA', 'SCGB2A1', 'GLS', 'FAP', 'DPP4', 'GPR119',
                      'PAX6', 'NEUROD1', 'LOXL4', 'PLCE1', 'GC', 'KLHL41',
                      'FEV', 'PTGER3', 'RFX6', 'SMARCA1', 'PGR', 'IRX1',
                      'UCP2', 'RGS4', 'KCNK16', 'GLP1R', 'ARX', 'POU3F4',
                      'RESP18', 'PYY', 'SLC38A5', 'TM4SF4', 'CRYBA2', 'SH3GL2',
                      'PCSK1', 'PRRG2', 'IRX2', 'ALDH1A1','PEMT', 'SMIM24',
                      'F10', 'SCGN', 'SLC30A8')

    gammaMarkers <- c('PPY', 'ABCC9', 'FGB', 'ZNF503', 'MEIS1', 'LMO3', 'EGR3',
                      'CHN2', 'PTGFR', 'ENTPD2', 'AQP3', 'THSD7A', 'CARTPT',
                      'ISL1', 'PAX6', 'NEUROD1', 'APOBEC2', 'SEMA3E', 'SLITRK6',
                      'SERTM1', 'PXK', 'PPY2P', 'ETV1', 'ARX', 'CMTM8', 'SCGB2A1',
                      'FXYD2', 'SCGN')

    geneSets <- list(acinarMarkers, alphaMarkers, ductalMarkers, gammaMarkers)
    names(geneSets) <- c('CSOA_acinar', 'CSOA_alpha', 'CSOA_ductal', 'CSOA_gamma')
    seuratObj <- runCSOA(seuratObj, geneSets, pairFileTemplate='pairs',
                                 keepOverlapOrder=TRUE)
    #> Warning in percentileSets(geneSetExp, percentile): 1 gene(s) had no top cells
    #> at the indicated percentile. These are now excluded from the gene signature.

    geneRadialPlot(lapply(paste0('pairs', names(geneSets), '.qs2'), qGrab),
                   groupLegendTitle='Gene set',
                   groupNames=str_remove(names(geneSets), 'CSOA_'),
                   cutoff=15,
                   title='Top 15 overlap genes',
                   extraCircles=2)

# Session information

    sessionInfo()
    #> R version 4.5.2 (2025-10-31)
    #> Platform: x86_64-pc-linux-gnu
    #> Running under: Ubuntu 24.04.4 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.26.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=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: Etc/UTC
    #> tzcode source: system (glibc)
    #>
    #> attached base packages:
    #> [1] stats4    stats     graphics  grDevices utils     datasets  methods
    #> [8] base
    #>
    #> other attached packages:
    #>  [1] Seurat_5.4.0                SeuratObject_5.3.0
    #>  [3] sp_2.2-1                    stringr_1.6.0
    #>  [5] scuttle_1.21.0              scRNAseq_2.25.0
    #>  [7] SingleCellExperiment_1.33.1 SummarizedExperiment_1.41.1
    #>  [9] Biobase_2.71.0              GenomicRanges_1.63.1
    #> [11] Seqinfo_1.1.0               IRanges_2.45.0
    #> [13] S4Vectors_0.49.0            BiocGenerics_0.57.0
    #> [15] generics_0.1.4              MatrixGenerics_1.23.0
    #> [17] matrixStats_1.5.0           patchwork_1.3.2
    #> [19] henna_0.7.5                 ggplot2_4.0.2
    #> [21] CSOA_1.1.3                  BiocStyle_2.39.0
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] ProtGenerics_1.43.0      spatstat.sparse_3.1-0    bitops_1.0-9
    #>   [4] httr_1.4.8               RColorBrewer_1.1-3       tools_4.5.2
    #>   [7] sctransform_0.4.3        alabaster.base_1.11.2    R6_2.6.1
    #>  [10] HDF5Array_1.39.0         lazyeval_0.2.2           uwot_0.2.4
    #>  [13] rhdf5filters_1.23.3      withr_3.0.2              gridExtra_2.3
    #>  [16] progressr_0.18.0         cli_3.6.5                spatstat.explore_3.7-0
    #>  [19] fastDummies_1.7.5        prismatic_1.1.2          labeling_0.4.3
    #>  [22] alabaster.se_1.11.0      sass_0.4.10              S7_0.2.1
    #>  [25] spatstat.data_3.1-9      ggridges_0.5.7           pbapply_1.7-4
    #>  [28] Rsamtools_2.27.1         parallelly_1.46.1        RSQLite_2.4.6
    #>  [31] BiocIO_1.21.0            ica_1.0-3                spatstat.random_3.4-4
    #>  [34] dplyr_1.2.0              wesanderson_0.3.7        Matrix_1.7-4
    #>  [37] abind_1.4-8              lifecycle_1.0.5          yaml_2.3.12
    #>  [40] rhdf5_2.55.15            SparseArray_1.11.10      BiocFileCache_3.1.0
    #>  [43] Rtsne_0.17               paletteer_1.7.0          grid_4.5.2
    #>  [46] blob_1.3.0               promises_1.5.0           ExperimentHub_3.1.0
    #>  [49] crayon_1.5.3             miniUI_0.1.2             lattice_0.22-9
    #>  [52] beachmat_2.27.3          cowplot_1.2.0            GenomicFeatures_1.63.1
    #>  [55] cigarillo_1.1.0          KEGGREST_1.51.1          sys_3.4.3
    #>  [58] maketools_1.3.2          pillar_1.11.1            knitr_1.51
    #>  [61] abdiv_0.2.0              rjson_0.2.23             bayesbio_1.0.0
    #>  [64] future.apply_1.20.2      codetools_0.2-20         glue_1.8.0
    #>  [67] spatstat.univar_3.1-6    data.table_1.18.2.1      gypsum_1.7.0
    #>  [70] vctrs_0.7.1              png_0.1-8                spam_2.11-3
    #>  [73] gtable_0.3.6             kernlab_0.9-33           rematch2_2.1.2
    #>  [76] cachem_1.1.0             xfun_0.56                S4Arrays_1.11.1
    #>  [79] mime_0.13                tidygraph_1.3.1          survival_3.8-6
    #>  [82] fitdistrplus_1.2-6       ROCR_1.0-12              liver_1.27
    #>  [85] nlme_3.1-168             bit64_4.6.0-1            alabaster.ranges_1.11.0
    #>  [88] filelock_1.0.3           RcppAnnoy_0.0.23         GenomeInfoDb_1.47.2
    #>  [91] bslib_0.10.0             irlba_2.3.7              KernSmooth_2.23-26
    #>  [94] otel_0.2.0               DBI_1.3.0                tidyselect_1.2.1
    #>  [97] bit_4.6.0                compiler_4.5.2           curl_7.0.0
    #> [100] httr2_1.2.2              h5mread_1.3.1            textshape_1.7.5
    #> [103] DelayedArray_0.37.0      plotly_4.12.0            stringfish_0.18.0
    #> [106] rtracklayer_1.71.3       scales_1.4.0             lmtest_0.9-40
    #> [109] ggeasy_0.1.6             rappdirs_0.3.4           digest_0.6.39
    #> [112] goftest_1.2-3            spatstat.utils_3.2-2     alabaster.matrix_1.11.0
    #> [115] rmarkdown_2.30           XVector_0.51.0           htmltools_0.5.9
    #> [118] pkgconfig_2.0.3          ensembldb_2.35.0         dbplyr_2.5.2
    #> [121] fastmap_1.2.0            UCSC.utils_1.7.1         rlang_1.1.7
    #> [124] htmlwidgets_1.6.4        shiny_1.13.0             farver_2.1.2
    #> [127] jquerylib_0.1.4          zoo_1.8-15               jsonlite_2.0.0
    #> [130] BiocParallel_1.45.0      RCurl_1.98-1.17          magrittr_2.0.4
    #> [133] dotCall64_1.2            Rhdf5lib_1.33.2          Rcpp_1.1.1
    #> [136] kerntools_1.2.1          ggnewscale_0.5.2         viridis_0.6.5
    #> [139] reticulate_1.45.0        alabaster.schemas_1.11.0 stringi_1.8.7
    #> [142] ggalluvial_0.12.6        ggraph_2.2.2             MASS_7.3-65
    #> [145] AnnotationHub_4.1.0      plyr_1.8.9               parallel_4.5.2
    #> [148] listenv_0.10.1           ggrepel_0.9.7            deldir_2.0-4
    #> [151] Biostrings_2.79.5        graphlayouts_1.2.3       splines_4.5.2
    #> [154] tensor_1.5.1             igraph_2.2.2             spatstat.geom_3.7-0
    #> [157] RcppHNSW_0.6.0           buildtools_1.0.0         reshape2_1.4.5
    #> [160] qs2_0.1.7                BiocVersion_3.23.1       XML_3.99-0.22
    #> [163] evaluate_1.0.5           RcppParallel_5.1.11-2    BiocManager_1.30.27
    #> [166] tweenr_2.0.3             httpuv_1.6.16            RANN_2.6.2
    #> [169] tidyr_1.3.2              purrr_1.2.1              polyclip_1.10-7
    #> [172] alabaster.sce_1.11.0     future_1.69.0            scattermore_1.2
    #> [175] ggforce_0.5.0            xtable_1.8-8             AnnotationFilter_1.35.0
    #> [178] restfulr_0.0.16          RSpectra_0.16-2          later_1.4.8
    #> [181] viridisLite_0.4.3        tibble_3.3.1             memoise_2.0.1
    #> [184] AnnotationDbi_1.73.0     GenomicAlignments_1.47.0 cluster_2.1.8.2
    #> [187] globals_0.19.0