Introduction to single-cell RNA-seq analysis
Differential Abundance
April 2024
Acknowledgments: much of the material in this section has been
derived from the chapters on differential expression and abundance in
the OSCA
book and the Hemberg
Group course materials. Additional material concerning
miloR has been based on the demonstration
from the Marioni Lab.
0.1 Differential Abundance Analysis
In the previous section we discussed how we can perform differential expression analysis using pseudo-bulk samples, when multiple biological replicates are available in a multi-condition experiment. Differential expression analysis addresses the question of whether a gene has a different average expression between the two groups being compared. However, another useful question to address is whether the cell composition also differs between conditions. We could imagine, for example, that natural killer cells are more abundant in leukemia samples than in healthy samples.
This type of analysis is referred to as differential abundance analysis.
0.2 Setup
We will continue working with our PBMMC and ETV6-RUNX1 samples, as we did in the differential expression section. Here is a reminder of the code to load the packages and read the data in.
# load packages
library(BiocParallel)
library(scran)
library(scater)
library(miloR)
library(tidyverse)
library(patchwork)
# load the SCE object
sce <- readRDS("R_objects/Caron_clustered.PBMMCandETV6RUNX1.rds")
# check the contents of the object
sce## class: SingleCellExperiment
## dim: 33102 30461
## metadata(1): Samples
## assays(3): counts logcounts reconstructed
## rownames(33102): MIR1302-2HG ENSG00000238009 ... ENSG00000276017
## ENSG00000278817
## rowData names(4): ID Symbol Type Chromosome
## colnames(30461): 1_AAACCTGAGACTTTCG-1 1_AAACCTGGTCTTCAAG-1 ...
## 11_TTTGTCATCAGTTGAC-1 11_TTTGTCATCTCGTTTA-1
## colData names(13): Sample Barcode ... k.60_cluster.fun.leiden label
## reducedDimNames(3): corrected TSNE_corrected UMAP_corrected
## mainExpName: NULL
## altExpNames(0):# plot UMAP done on the batch-corrected data
plotReducedDim(sce, dimred = "UMAP_corrected",
colour_by = "label",
text_by = "label")
0.3 Differential abundance between conditions
# Differential abundance with Milo ----One method to perform DA analysis is to simply count how many cells occur in each label + condition combination and then test for differences in the counts between the two conditions.
table(sce$label, sce$SampleName)##
## ETV6-RUNX1_1 ETV6-RUNX1_2 ETV6-RUNX1_3 ETV6-RUNX1_4
## B (c1) 1761 4939 1791 2331
## B (c2) 650 1023 251 1274
## B (c3) 244 97 17 27
## T (c4) 15 193 1897 283
## Erythrocytes (c5) 1 0 71 178
## CD20+ B (c6) 6 115 545 139
## B (c7) 83 389 91 95
## NK T (c8) 9 44 96 59
## Erythrocytes (c9) 3 13 95 159
## Mono (c10) 1 6 19 43
## B (c11) 0 17 9 5
## CD20+ B (c12) 0 2 0 0
## CD20+ B (c13) 0 1 0 0
## T (c14) 0 5 23 4
## Erythrocytes (c15) 0 1 89 1019
## Erythrocytes (c16) 0 0 26 79
## Mono (c17) 0 0 0 0
##
## PBMMC_1 PBMMC_2 PBMMC_3
## B (c1) 40 69 187
## B (c2) 16 39 113
## B (c3) 31 110 144
## T (c4) 129 1213 1143
## Erythrocytes (c5) 7 483 18
## CD20+ B (c6) 33 418 584
## B (c7) 6 21 5
## NK T (c8) 26 388 280
## Erythrocytes (c9) 26 541 119
## Mono (c10) 140 405 690
## B (c11) 4 32 28
## CD20+ B (c12) 124 80 253
## CD20+ B (c13) 251 211 780
## T (c14) 0 1 5
## Erythrocytes (c15) 0 788 39
## Erythrocytes (c16) 0 11 5
## Mono (c17) 13 42 37As these are count data, statistical methods used in flow-cytometry have been adapted to test for differences in cell abundance from such a matrix of counts. This is the approach used in the Differential Abundance Analysis chapter in the OSCA book.
However, this approach relies on a fixed set of clusters determined by us beforehand, which can be limiting in cases where the changes in abundance are more continuous (e.g. along a developmental trajectory, or where cell states are not completely discrete). To address this limitation, Dann et al. (2022) developed a method where cell abundance differences are tested along neighbourhoods of cells on a KNN graph. This means that the results aren’t dependent on our clustering results, and are instead treated in a more “continuous” way.
This method has been implemented by the authors in the R/Bioconductor
package MiloR,
which we cover in this section. Starting from a graph that faithfully
recapitulates the biology of the cell population, the Milo
analysis consists of 3 steps:
- Building a k-nearest neighbours (KNN) graph
- Sampling representative neighbourhoods in the graph (for computational efficiency)
- Testing for differential abundance of conditions in all neighbourhoods
- Accounting for multiple hypothesis testing using a weighted FDR procedure that accounts for the overlap of neighbourhoods
The first step in the analysis is to turn our single cell object into
a Milo object. This is very similar to the
SingleCellExperiment object we’ve been working with so far,
but also includes information about the neighbourhood graphs that are
built during the analysis.
# create the Milo object can be simply converted from a SCE
milo <- Milo(sce)
milo## class: Milo
## dim: 33102 30461
## metadata(1): Samples
## assays(3): counts logcounts reconstructed
## rownames(33102): MIR1302-2HG ENSG00000238009 ... ENSG00000276017
## ENSG00000278817
## rowData names(4): ID Symbol Type Chromosome
## colnames(30461): 1_AAACCTGAGACTTTCG-1 1_AAACCTGGTCTTCAAG-1 ...
## 11_TTTGTCATCAGTTGAC-1 11_TTTGTCATCTCGTTTA-1
## colData names(13): Sample Barcode ... k.60_cluster.fun.leiden label
## reducedDimNames(3): corrected TSNE_corrected UMAP_corrected
## mainExpName: NULL
## altExpNames(0):
## nhoods dimensions(2): 1 1
## nhoodCounts dimensions(2): 1 1
## nhoodDistances dimension(1): 0
## graph names(0):
## nhoodIndex names(1): 0
## nhoodExpression dimension(2): 1 1
## nhoodReducedDim names(0):
## nhoodGraph names(0):
## nhoodAdjacency dimension(2): 1 1## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 8299229 443.3 14391954 768.7 NA 14391954 768.7
## Vcells 172193273 1313.8 242807464 1852.5 16384 185283113 1413.6Notice how there are now several slots with prefix
nhoods, which we will explore as we go through the
analysis.
0.4 Construct KNN graph
# Build KNN graph ----The first step in our analysis is to build a KNN graph from our cells. This is very similar to what we did earlier in the clustering session, except now the graph will be stored inside the object:
# add KNN graph to Milo object
milo <- buildGraph(milo,
k = 60,
d = 50,
reduced.dim = "corrected",
BPPARAM = MulticoreParam(7))A couple of things to note about this:
kis the number of nearest neighbours to build the graph. This can be adjusted depending on how much of a fine resolution you want. If you use too high a number, you will end up loosing resolution as a higher diversity of cells will be connected in the graph. On the other hand, if you use too low a number, you may increase noise in the data and loose statistical power, as only very few cells will be connected to each other in a neighbourhood. The author’s recommendation is to use a value of k as you did for clustering and UMAP visualisation - this is what we’ve done in this case.dis the number of dimensions from our dimensionality reduction matrix to use. In this case we’re using the MNN-corrected matrix, which contains 50 dimensions, and we use all of them (50 is also the default, so we could have left this option out).
The object now has a new object inside the graph slot,
which a standard object type from the igraph package:
# the graph is stored as a standard igraph object
graph(milo)## IGRAPH 90b7933 U--- 30461 1412281 --
## + edges from 90b7933:
## [1] 1-- 138 1-- 351 1-- 462 1-- 536 1-- 569 1-- 597 1-- 600 1-- 606
## [9] 1-- 760 1-- 780 1-- 990 1-- 1542 1-- 1557 1-- 1674 1-- 1711 1-- 1936
## [17] 1-- 1964 1-- 1973 1-- 2013 1-- 2277 1-- 2538 1-- 2634 1-- 2817 1-- 3123
## [25] 1-- 3483 1-- 3854 1-- 5903 1-- 6022 1-- 6184 1-- 7528 1-- 7925 1-- 8095
## [33] 1-- 9222 1-- 9538 1--10135 1--11901 1--12550 1--12612 1--12783 1--14466
## [41] 1--15031 1--15089 1--15163 1--15209 1--15249 1--15611 1--15614 1--15698
## [49] 1--16125 1--16332 1--16552 1--16795 1--17338 1--17623 1--18057 1--18718
## [57] 1--18745 1--19261 1--19340 1--19630 2-- 27 2-- 85 2-- 216 2-- 239
## [65] 2-- 308 2-- 310 2-- 320 2-- 333 2-- 393 2-- 411 2-- 412 2-- 457
## + ... omitted several edges0.5 Define neighbourhoods
The next step in the analysis is to define cell neighbourhoods. This essentially consists of picking a focal cell and counting how many other cells it is connected to in the graph (and how many come from each group under comparison). However, if we did this for every single cell in the data, it could get computationally quite intractable. Instead, Milo implements a sampling step, where so-called “index cells” are sampled from the larger graph, and those will be the neighbourhoods used in the DA analysis.
This is done with the makeNhodds() function:
# sample index cells to define neighbourhoods
milo <- makeNhoods(milo,
prop = 0.1,
k = 60,
d = 50,
reduced_dims = "corrected")
# check our object again
milo## class: Milo
## dim: 33102 30461
## metadata(1): Samples
## assays(3): counts logcounts reconstructed
## rownames(33102): MIR1302-2HG ENSG00000238009 ... ENSG00000276017
## ENSG00000278817
## rowData names(4): ID Symbol Type Chromosome
## colnames(30461): 1_AAACCTGAGACTTTCG-1 1_AAACCTGGTCTTCAAG-1 ...
## 11_TTTGTCATCAGTTGAC-1 11_TTTGTCATCTCGTTTA-1
## colData names(13): Sample Barcode ... k.60_cluster.fun.leiden label
## reducedDimNames(3): corrected TSNE_corrected UMAP_corrected
## mainExpName: NULL
## altExpNames(0):
## nhoods dimensions(2): 30461 1574
## nhoodCounts dimensions(2): 1 1
## nhoodDistances dimension(1): 0
## graph names(1): graph
## nhoodIndex names(1): 1574
## nhoodExpression dimension(2): 1 1
## nhoodReducedDim names(0):
## nhoodGraph names(0):
## nhoodAdjacency dimension(2): 1 1Some things to note:
propis the proportion of cells to sample (the authors advise 0.1 - 0.2).kanddshould be set the same as the k value used to build the original KNN graph. These values will be used for the graph-sampling algorithm used behind the scenes.
We can also see that our Milo object now has the nhoods
slots populated. In this case, it indicates that the sampling algorithm
picked 1574 index cells to form neighbourhoods for DA analysis.
One good QC metric to check at this stage is to create a histogram of cell counts per neighbourhood. Like we said earlier, when we set the value of k to define our KNN graph, we want to make sure the value is not too low, such that most neighbourhoods have very few cells, nor for it to be so high that we have very large (and presumably heterogenous) neighbourhoods.
Conveniently, MiloR has a plotting function just for
this purpose:
# distribution of neighbourhood sizes
plotNhoodSizeHist(milo) +
geom_vline(xintercept = 100, col = "salmon")
The authors of Milo have stated several different parameters of
deciding if your histogram is correct. Either ‘peaking above 20’,
‘peaking between 50 and 100’ or ‘an average neighbourhood size over 5 x
N_samples’. Realistically, all of these statements give numbers in the
same ballpark and so we can make a judgement on our data. In our case,
we can see that our histogram peaks at ~100, suggesting a good
neighbourhood size. If this was not the case, we could re-run the
analysis from the start, making a KNN graph with a higher/lower
k value.
0.6 Counting cells
After defining our sample of neighbourhoods, the next step consists
of counting how many cells there are in each neighbourhood for each
sample replicate. In this step we need to define which column of our
metadata we want to use to do the counting. We should count cells at the
biological replicate leve, which in our case is stored in the
SampleName column.
# count cells in each neighbourhood
milo <- countCells(milo,
meta.data = colData(milo),
samples = "SampleName")
# Milo now has a counts matrix
head(nhoodCounts(milo))## 6 x 7 sparse Matrix of class "dgCMatrix"
## ETV6-RUNX1_1 ETV6-RUNX1_2 ETV6-RUNX1_3 ETV6-RUNX1_4 PBMMC_1 PBMMC_2 PBMMC_3
## 1 . . 3 1 1 52 10
## 2 7 313 5 56 . . 2
## 3 4 37 233 21 20 44 121
## 4 . 26 451 49 2 39 118
## 5 . 15 327 21 1 22 99
## 6 14 216 263 17 2 3 12The dimensions of the counts matrix corresponds to the number of neighbourhoods (rows) and samples (columns), in our case a 1574 by 7 matrix.
0.7 Neighbourhood connectivity
There is one more step that we need to do before we are ready to run our DA analysis. It consists of calculating the distances between each neighbourhood. As we said, the neighbourhoods on our graph may partially overlap, so when Milo corrects our p-values for multiple testing, it takes into account the spatial dependence of our tests. For example, neighbourhoods that are closer to each other may have similar p-values, and so we should avoid penalising them too much, otherwise we will sacrifice statistical power.
# calculate distances between neighbourhoods - for p-value correction
milo <- calcNhoodDistance(milo, d = 50, reduced.dim = "corrected")As before, the value of d (number of dimensions to
consider from our MNN-corrected matrix) should be the same that was used
for building our initial graph.
0.8 Running DA tests
Finally, we are ready to run the actual differential abundance step.
Similarly to the pseudo-bulk appproach for differential expression,
MiloR takes advantage of the edgeR package and
its negative binomial linear model implementation. This provides a
suitable statistical model to account for over-dispersed count data, as
is common with these kind of datasets.
First, we need to define a simple table with information about our samples - we will use this table to define our model formula (similarly to what we did in the differential expression step). In this case, we want to detect DA between our PBMMC and ETV6-RUNX1 sample groups, so we define a table based on those two groups. We could also include a batch column in this table, but to keep the demonstration simple we will skip this. We would normally do this if we know there is a batch effect that we want to account for it in DA testing.
# define a table for our model design
sample_info <- unique(colData(milo)[,c("SampleName", "SampleGroup")])
rownames(sample_info) <- sample_info$SampleName
sample_info## DataFrame with 7 rows and 2 columns
## SampleName SampleGroup
## <character> <character>
## ETV6-RUNX1_1 ETV6-RUNX1_1 ETV6-RUNX1
## ETV6-RUNX1_2 ETV6-RUNX1_2 ETV6-RUNX1
## ETV6-RUNX1_3 ETV6-RUNX1_3 ETV6-RUNX1
## ETV6-RUNX1_4 ETV6-RUNX1_4 ETV6-RUNX1
## PBMMC_1 PBMMC_1 PBMMC
## PBMMC_2 PBMMC_2 PBMMC
## PBMMC_3 PBMMC_3 PBMMCNow we can do the DA test, explicitly defining our experimental
design. In this case as discussed we will test the difference between
sample groups. The testNhoods function calculates a
Fold-change and corrected P-value for each neighbourhood, which
indicates whether there is significant differential abundance between
sample groups.
# run DA test
da_results <- testNhoods(milo,
design = ~ SampleGroup,
design.df = sample_info,
reduced.dim = "corrected")
# results are returned as a data.frame
da_results %>%
arrange(SpatialFDR) %>%
head()## logFC logCPM F PValue FDR Nhood SpatialFDR
## 48 10.10456 10.02440 24.72729 0.0002901171 0.006444996 48 0.007941485
## 59 10.44691 10.36420 23.67404 0.0003490528 0.006444996 59 0.007941485
## 62 10.25417 10.17140 24.99681 0.0002769508 0.006444996 62 0.007941485
## 108 10.24200 10.15622 24.62901 0.0002950979 0.006444996 108 0.007941485
## 129 10.39473 10.31135 24.40224 0.0003069756 0.006444996 129 0.007941485
## 148 10.60814 10.52343 24.18342 0.0003189665 0.006444996 148 0.0079414850.9 Inspecting DA results
A good diagnostic plot to make after running our analysis is a histogram of p-values. We should expect a distribution with a peak close to 0 (corresponding to differentially abundant neighbourhoods) and tailing off towards 1. This blog article explains the different p-value histogram profiles you may see and what they can mean.
# p-value histogram
ggplot(da_results, aes(PValue)) +
geom_histogram(bins = 50)
Next, we can get an overview of our results as a volcano plot, marking a significance threshold of 10% FDR:
# volcano plot
# each point in this plot corresponds to a neighbourhood (not a cell)
ggplot(da_results, aes(logFC, -log10(SpatialFDR))) +
geom_point(aes(colour = FDR < 0.1)) +
geom_hline(yintercept = 1) 
As we can see, several neighbourhoods fall below our FDR threshold, indicating significant differential abundance of cells between PBMMC and ETV6-RUNX1 samples.
There is an unsual discontinuity in our logFC axis, suggesting a sudden change in abundance in some of our neighbourhoods. We can get a better idea of the fold changes by visualising them as a graph for our neighbourhoods (rather than the original single-cell graph, which would be too big to display). We can then super-impose this graph on our original UMAP (or t-SNE), to compare with our original analysis.
# build neighbourhood graph embedding
milo <- buildNhoodGraph(milo)# our original UMAP with our previously annotated cell labels
umap_plot <- plotReducedDim(milo,
dimred = "UMAP_corrected",
colour_by = "label",
text_by = "label")
# the neighbourhood map adjusted to match UMAP embedding
nh_graph_plot <- plotNhoodGraphDA(milo,
da_results,
layout = "UMAP_corrected",
alpha = 0.05)
# the two plots together side-by-side
umap_plot + nh_graph_plot +
plot_layout(guides="collect")
On the left we have our original UMAP, with cell/cluster annotations we did previously (by standard clustering and using known marker genes to manually annotate our cells). On the right we have the Milo neighbourhood graph, where each node represents a neighbourhood, coloured by the log fold-change in abundance between PBMMC and ETV6-RUNX1 samples (positive values represent higher abundance in ETV6-RUNX1, and vice-versa). We can see, for example, a cluster of cells with negative log fold changes around our large B cell cluster, which likely explains the discontinuity in values we saw earlier in our volcano plot.
Although we have annotations for our cells, these annotations at the moment are not present in the differential abundance table:
head(da_results)## logFC logCPM F PValue FDR Nhood SpatialFDR
## 1 6.9219691 9.878859 12.0517888 0.003289915 0.02538396 1 0.0291073
## 2 -3.6154326 8.867091 5.0096231 0.040323653 0.09989293 2 0.1059409
## 3 1.9703218 11.539023 1.1480863 0.300403783 0.38039868 3 0.3872336
## 4 0.8197139 11.672708 0.1529064 0.701113049 0.74014215 4 0.7420397
## 5 0.8856981 11.221669 0.1835949 0.674215900 0.71800800 5 0.7205472
## 6 -1.9537987 10.317513 1.3559960 0.261899793 0.34238395 6 0.3495426We can transfer our cell labels to the neighbourhood DA results, by simply counting how many cells within that neighbourhood share the same label, and assign the label with the highest counts.
# annotate our neighbourhood DA results with our cell labels
da_results <- annotateNhoods(milo, da_results, coldata_col = "label")
head(da_results)## logFC logCPM F PValue FDR Nhood SpatialFDR
## 1 6.9219691 9.878859 12.0517888 0.003289915 0.02538396 1 0.0291073
## 2 -3.6154326 8.867091 5.0096231 0.040323653 0.09989293 2 0.1059409
## 3 1.9703218 11.539023 1.1480863 0.300403783 0.38039868 3 0.3872336
## 4 0.8197139 11.672708 0.1529064 0.701113049 0.74014215 4 0.7420397
## 5 0.8856981 11.221669 0.1835949 0.674215900 0.71800800 5 0.7205472
## 6 -1.9537987 10.317513 1.3559960 0.261899793 0.34238395 6 0.3495426
## label label_fraction
## 1 Erythrocytes (c9) 1
## 2 B (c1) 1
## 3 T (c4) 1
## 4 T (c4) 1
## 5 T (c4) 1
## 6 B (c1) 1The result includes the fraction of cells in the neighbourhood that shared that label. We can look at the distribution of this fraction as a quality check:
# histogram of fraction of cells in the neighbourhood with the same label
ggplot(da_results, aes(label_fraction)) +
geom_histogram(bins = 50)
We should expect the plot to be very biased towards 1, as is seen. This makes sense, since Milo’s neighbourhoods should be very similar to our previously-defined clusters. After all, both our clustering and Milo used KNN graphs as the starting point, and we’ve used the same settings for the k and d parameters in both steps of our analysis.
Despite most neighbourhoods being homogenous, some seem to have mixed labels. We can highlight these in our results table:
# add "mixed" label to neighbourhoods with less 70% consistency
da_results$label <- ifelse(da_results$label_fraction < 0.7,
"Mixed",
da_results$label)
head(da_results)## logFC logCPM F PValue FDR Nhood SpatialFDR
## 1 6.9219691 9.878859 12.0517888 0.003289915 0.02538396 1 0.0291073
## 2 -3.6154326 8.867091 5.0096231 0.040323653 0.09989293 2 0.1059409
## 3 1.9703218 11.539023 1.1480863 0.300403783 0.38039868 3 0.3872336
## 4 0.8197139 11.672708 0.1529064 0.701113049 0.74014215 4 0.7420397
## 5 0.8856981 11.221669 0.1835949 0.674215900 0.71800800 5 0.7205472
## 6 -1.9537987 10.317513 1.3559960 0.261899793 0.34238395 6 0.3495426
## label label_fraction
## 1 Erythrocytes (c9) 1
## 2 B (c1) 1
## 3 T (c4) 1
## 4 T (c4) 1
## 5 T (c4) 1
## 6 B (c1) 1Finally, we can visualise the distribution of DA fold changes in different labels.
# distribution of logFC across neighbourhood labels
plotDAbeeswarm(da_results, group.by = "label")
We can see that several clusters are enriched and some are depleted between our sample groups. There are not many DA neighbourhoods with a mixed label either. And we can see that B cells cluster 1 has quite a discontinuity in log fold changes - this is what was causing the odd distribution in the volcano plot, and we saw earlier is likely due to a sub-population of B cells with different representation across our samples.
It is also worth noting that this analysis may also reveal uncorrected batch effects, which would lead to an artificial separation of our samples according to their group. As always, we should always proceed in our analysis with care and back our conclusions with further confirmatory experiments.
There are further downstream analysis possible to identify genes that differentiate between our different neighbourhoods. This is similar in spirit to the marker gene analysis we did earlier, but based on the neighbourhood graph constructed by Milo. You can consult the package’s main vignette to see examples of these downstream analysis.
0.10 Exercise
:::exercise
We want to achieve the following:
- Rerun the DA analysis this time changing the ‘k’ and ‘d’ parameters, does this alter the results?
- Which cell type cluster has the most changing neighbourhoods? Does this make sense from what we know of the biology?
0.11 Session information
sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.7.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/London
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] patchwork_1.2.0 lubridate_1.9.3
## [3] forcats_1.0.0 stringr_1.5.1
## [5] dplyr_1.1.4 purrr_1.0.2
## [7] readr_2.1.5 tidyr_1.3.1
## [9] tibble_3.2.1 tidyverse_2.0.0
## [11] miloR_1.10.0 edgeR_4.0.16
## [13] limma_3.58.1 scater_1.30.1
## [15] ggplot2_3.5.0 scran_1.30.2
## [17] scuttle_1.12.0 SingleCellExperiment_1.24.0
## [19] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [21] GenomicRanges_1.54.1 GenomeInfoDb_1.38.8
## [23] IRanges_2.36.0 S4Vectors_0.40.2
## [25] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [27] matrixStats_1.3.0 BiocParallel_1.36.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 gridExtra_2.3
## [3] rlang_1.1.3 magrittr_2.0.3
## [5] compiler_4.3.2 DelayedMatrixStats_1.24.0
## [7] vctrs_0.6.5 pkgconfig_2.0.3
## [9] crayon_1.5.2 fastmap_1.1.1
## [11] XVector_0.42.0 labeling_0.4.3
## [13] ggraph_2.2.1 utf8_1.2.4
## [15] rmarkdown_2.26 tzdb_0.4.0
## [17] ggbeeswarm_0.7.2 xfun_0.43
## [19] bluster_1.12.0 zlibbioc_1.48.2
## [21] cachem_1.0.8 beachmat_2.18.1
## [23] jsonlite_1.8.8 highr_0.10
## [25] DelayedArray_0.28.0 tweenr_2.0.3
## [27] irlba_2.3.5.1 parallel_4.3.2
## [29] cluster_2.1.6 R6_2.5.1
## [31] stringi_1.8.3 RColorBrewer_1.1-3
## [33] bslib_0.7.0 jquerylib_0.1.4
## [35] Rcpp_1.0.12 knitr_1.45
## [37] splines_4.3.2 timechange_0.3.0
## [39] Matrix_1.6-5 igraph_2.0.3
## [41] tidyselect_1.2.1 rstudioapi_0.16.0
## [43] abind_1.4-5 yaml_2.3.8
## [45] viridis_0.6.5 codetools_0.2-20
## [47] lattice_0.22-6 withr_3.0.0
## [49] evaluate_0.23 polyclip_1.10-6
## [51] pillar_1.9.0 generics_0.1.3
## [53] RCurl_1.98-1.14 hms_1.1.3
## [55] sparseMatrixStats_1.14.0 munsell_0.5.1
## [57] scales_1.3.0 gtools_3.9.5
## [59] glue_1.7.0 metapod_1.10.1
## [61] tools_4.3.2 BiocNeighbors_1.20.2
## [63] ScaledMatrix_1.10.0 locfit_1.5-9.9
## [65] graphlayouts_1.1.1 cowplot_1.1.3
## [67] tidygraph_1.3.1 grid_4.3.2
## [69] colorspace_2.1-0 GenomeInfoDbData_1.2.11
## [71] beeswarm_0.4.0 BiocSingular_1.18.0
## [73] ggforce_0.4.2 vipor_0.4.7
## [75] cli_3.6.2 rsvd_1.0.5
## [77] fansi_1.0.6 S4Arrays_1.2.1
## [79] viridisLite_0.4.2 gtable_0.3.4
## [81] sass_0.4.9 digest_0.6.35
## [83] SparseArray_1.2.4 ggrepel_0.9.5
## [85] dqrng_0.3.2 farver_2.1.1
## [87] memoise_2.0.1 htmltools_0.5.8.1
## [89] lifecycle_1.0.4 statmod_1.5.0
## [91] MASS_7.3-60.0.1