Dreamlet analysis of single cell RNA-seq
Linear (mixed) model analysis of pseudobulk data
Developed by Gabriel Hoffman
Run on 2024-02-05 16:08:14
Source:vignettes/dreamlet.Rmd
dreamlet.RmdIntroduction
As the scale of single cell/nucleus RNA-seq has increased, so has the complexity of study designs. Analysis of datasets with simple study designs can be performed using linear model as in the muscat package. Yet analysis of datsets with complex study designs such as repeated measures or many technical batches can benefit from linear mixed model analysis to model to correlation structure between samples. We previously developed dream to apply linear mixed models to bulk RNA-seq data using a limma-style workflow. Dreamlet extends the previous work of dream and muscat to apply linear mixed models to pseudobulk data. Dreamlet also supports linear models and facilitates application of 1) variancePartition to quantify the contribution of multiple variables to expression variation, and 2) zenith to perform gene set analysis on the differential expression signatures.
Installation
To install this package, start R (version “4.3”) and enter:
if (!require("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
# Select release #1 or #2
# 1) Bioconductor release
BiocManager::install("dreamlet")
# 2) Latest stable release
devtools::install_github("DiseaseNeurogenomics/dreamlet")Process single cell count data
Here we perform analysis of PBMCs from 8 individuals stimulated with interferon-β (Kang, et al, 2018, Nature Biotech). This is a small dataset that does not have repeated measures or high dimensional batch effects, so the sophisticated features of dreamlet are not strictly necessary. But this gives us an opportunity to walk through a standard dreamlet workflow.
Preprocess data
Here, single cell RNA-seq data is downloaded from ExperimentHub.
library(dreamlet)
library(muscat)
library(ExperimentHub)
library(zenith)
library(scater)
# Download data, specifying EH2259 for the Kang, et al study
eh <- ExperimentHub()
sce <- eh[["EH2259"]]
# only keep singlet cells with sufficient reads
sce <- sce[rowSums(counts(sce) > 0) > 0, ]
sce <- sce[, colData(sce)$multiplets == "singlet"]
# compute QC metrics
qc <- perCellQCMetrics(sce)
# remove cells with few or many detected genes
ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
sce <- sce[, !ol]
# set variable indicating stimulated (stim) or control (ctrl)
sce$StimStatus <- sce$stimAggregate to pseudobulk
Dreamlet, like muscat, performs analysis at the pseudobulk-level by
summing raw counts across cells for a given sample and cell type.
aggregateToPseudoBulk is substantially faster for large
on-disk datasets than muscat::aggregateData.
# Since 'ind' is the individual and 'StimStatus' is the stimulus status,
# create unique identifier for each sample
sce$id <- paste0(sce$StimStatus, sce$ind)
# Create pseudobulk data by specifying cluster_id and sample_id
# Count data for each cell type is then stored in the `assay` field
# assay: entry in assayNames(sce) storing raw counts
# cluster_id: variable in colData(sce) indicating cell clusters
# sample_id: variable in colData(sce) indicating sample id for aggregating cells
pb <- aggregateToPseudoBulk(sce,
assay = "counts",
cluster_id = "cell",
sample_id = "id",
verbose = FALSE
)
# one 'assay' per cell type
assayNames(pb)## [1] "B cells" "CD14+ Monocytes" "CD4 T cells"
## [4] "CD8 T cells" "Dendritic cells" "FCGR3A+ Monocytes"
## [7] "Megakaryocytes" "NK cells"Voom for pseudobulk
Apply voom-style normalization for pseudobulk counts within each cell cluster using voomWithDreamWeights to handle random effects (if specified).
# Normalize and apply voom/voomWithDreamWeights
res.proc <- processAssays(pb, ~StimStatus, min.count = 5)
# the resulting object of class dreamletProcessedData stores
# normalized data and other information
res.proc## class: dreamletProcessedData
## assays(8): B cells CD14+ Monocytes ... Megakaryocytes NK cells
## colData(4): ind stim multiplets StimStatus
## metadata(3): cell id cluster
## Samples:
## min: 13
## max: 16
## Genes:
## min: 164
## max: 5262
## details(7): assay n_retain ... n_errors error_initialprocessAssays() retains samples with at least
min.cells in a given cell type. While dropping a few
samples usually is not a problem, in some cases dropping sames can mean
that a variable included in the regression formula no longer has any
variation. For example, dropping all stimulated samples from analysis of
a given cell type would be mean the variable StimStatus has
no variation and is perfectly colinear with the intercept term. This
colinearity issue is detected internally and variables with these
problem are dropped from the regression formula for that particular cell
type. The number of samples retained and the resulting formula used in
each cell type can be accessed as follows. In this analysis, samples are
dropped from 3 cell types but the original formula remains valid in each
case.
# view details of dropping samples
details(res.proc)## assay n_retain formula formDropsTerms n_genes n_errors
## 1 B cells 16 ~StimStatus FALSE 1961 0
## 2 CD14+ Monocytes 16 ~StimStatus FALSE 3087 0
## 3 CD4 T cells 16 ~StimStatus FALSE 5262 0
## 4 CD8 T cells 16 ~StimStatus FALSE 1030 0
## 5 Dendritic cells 13 ~StimStatus FALSE 164 0
## 6 FCGR3A+ Monocytes 16 ~StimStatus FALSE 1160 0
## 7 Megakaryocytes 13 ~StimStatus FALSE 172 0
## 8 NK cells 16 ~StimStatus FALSE 1656 0
## error_initial
## 1 FALSE
## 2 FALSE
## 3 FALSE
## 4 FALSE
## 5 FALSE
## 6 FALSE
## 7 FALSE
## 8 FALSEHere the mean-variance trend from voom is shown for each cell type. Cell types with sufficient number of cells and reads show a clear mean-variance trend. While in rare cell types like megakaryocytes, fewer genes have sufficient reads and the trend is less apparent.
# show voom plot for each cell clusters
plotVoom(res.proc)
# Show plots for subset of cell clusters
# plotVoom( res.proc[1:3] )
# Show plots for one cell cluster
# plotVoom( res.proc[["B cells"]])Variance partitioning
The variancePartition package uses linear and linear mixed models to quanify the contribution of multiple sources of expression variation at the gene-level. For each gene it fits a linear (mixed) model and evalutes the fraction of expression variation explained by each variable.
Variance fractions can be visualized at the gene-level for each cell type using a bar plot, or genome-wide using a violin plot.
# run variance partitioning analysis
vp.lst <- fitVarPart(res.proc, ~StimStatus)
# Show variance fractions at the gene-level for each cell type
genes <- vp.lst$gene[2:4]
plotPercentBars(vp.lst[vp.lst$gene %in% genes, ])
# Summarize variance fractions genome-wide for each cell type
plotVarPart(vp.lst, label.angle = 60)
Differential expression
Since the normalized expression data and metadata are stored within
res.proc, only the regression formula remains to be
specified. Here we only included the stimulus status, but analyses of
larger datasets can include covariates and random effects. With formula
~ StimStatus, an intercept is fit and coefficient
StimStatusstim log fold change between simulated and
controls.
# Differential expression analysis within each assay,
# evaluated on the voom normalized data
res.dl <- dreamlet(res.proc, ~StimStatus)
# names of estimated coefficients
coefNames(res.dl)## [1] "(Intercept)" "StimStatusstim"
# the resulting object of class dreamletResult
# stores results and other information
res.dl## class: dreamletResult
## assays(8): B cells CD14+ Monocytes ... Megakaryocytes NK cells
## Genes:
## min: 164
## max: 5262
## details(7): assay n_retain ... n_errors error_initial
## coefNames(2): (Intercept) StimStatusstimVolcano plots
The volcano plot can indicate the strength of the differential expression signal with each cell type. Red points indicate FDR < 0.05.
plotVolcano(res.dl, coef = "StimStatusstim")
Gene-level heatmap
For each cell type and specified gene, show z-statistic from
dreamlet analysis. Grey indicates that insufficient reads
were observed to include the gene in the analysis.
genes <- c("ISG20", "ISG15")
plotGeneHeatmap(res.dl, coef = "StimStatusstim", genes = genes)
Extract results
Each entry in res.dl stores a model fit by dream(),
and results can be extracted using topTable() as in
limma by specifying the coefficient of interest. The
results shows the gene name, log fold change, average expression,
t-statistic, p-value, FDR (i.e. adj.P.Val).
# results from full analysis
topTable(res.dl, coef = "StimStatusstim")## DataFrame with 10 rows and 9 columns
## assay ID logFC AveExpr t P.Value
## <character> <character> <numeric> <numeric> <numeric> <numeric>
## 1 B cells ISG15 5.46383 10.19452 29.8807 3.82697e-22
## 2 CD4 T cells ISG20 3.08327 10.40321 37.7709 6.82462e-22
## 3 B cells ISG20 3.23736 11.38145 28.4318 1.39774e-21
## 4 CD4 T cells MT2A 4.32138 7.80016 36.4780 1.47906e-21
## 5 CD14+ Monocytes IL1RN 7.04042 7.78556 39.9678 1.48612e-21
## 6 CD4 T cells IFI6 5.78194 8.52278 33.9813 7.12211e-21
## 7 CD4 T cells HERC5 4.28329 6.72836 33.1920 1.19816e-20
## 8 CD4 T cells TNFSF10 4.48757 7.48063 32.2412 2.27814e-20
## 9 NK cells ISG15 4.60250 11.00679 26.9896 2.57986e-20
## 10 NK cells ISG20 3.69323 10.96970 26.7509 3.21371e-20
## adj.P.Val B z.std
## <numeric> <numeric> <numeric>
## 1 4.30737e-18 40.6512 29.8807
## 2 4.30737e-18 40.2095 37.7709
## 3 4.30737e-18 39.4170 28.4318
## 4 4.30737e-18 39.3739 36.4780
## 5 4.30737e-18 38.5477 39.9678
## 6 1.72023e-17 37.8585 33.9813
## 7 2.48053e-17 37.2367 33.1920
## 8 4.12685e-17 36.7118 32.2412
## 9 4.15415e-17 36.4711 26.9896
## 10 4.65731e-17 36.2740 26.7509
# only B cells
topTable(res.dl[["B cells"]], coef = "StimStatusstim")## logFC AveExpr t P.Value adj.P.Val B
## ISG15 5.463825 10.194524 29.88069 3.826968e-22 7.504683e-19 40.65119
## ISG20 3.237359 11.381452 28.43176 1.397743e-21 1.370487e-18 39.41704
## LY6E 4.199663 8.690971 23.65134 1.631076e-19 1.066180e-16 34.64028
## PLSCR1 4.070973 8.256238 22.90123 3.725818e-19 1.718635e-16 33.79523
## EPSTI1 3.652905 8.089555 22.75648 4.382038e-19 1.718635e-16 33.63648
## SAT1 2.140671 9.651797 21.92936 1.127706e-18 3.685718e-16 32.73085
## IRF7 3.655292 8.662936 21.37136 2.173109e-18 6.087810e-16 32.09392
## UBE2L6 2.688268 9.233683 21.17978 2.731666e-18 6.370928e-16 31.85761
## TRIM22 2.906534 7.160498 21.12311 2.923934e-18 6.370928e-16 31.70909
## SOCS1 2.427295 8.581121 21.00831 3.357573e-18 6.584201e-16 31.66236Forest plot
A forest plot shows the log fold change and standard error of a given gene across all cell types. The color indicates the FDR.
plotForest(res.dl, coef = "StimStatusstim", gene = "ISG20")
Box plot
Examine the expression of ISG20 between stimulation conditions within
CD14+ Monocytes. Use extractData() to create a
tibble with gene expression data and metadata from
colData() from one cell type.
# get data
df <- extractData(res.proc, "CD14+ Monocytes", genes = "ISG20")
# set theme
thm <- theme_classic() +
theme(aspect.ratio = 1, plot.title = element_text(hjust = 0.5))
# make plot
ggplot(df, aes(StimStatus, ISG20)) +
geom_boxplot() +
thm +
ylab(bquote(Expression ~ (log[2] ~ CPM))) +
ggtitle("ISG20")
Advanced used of contrasts
A hypothesis test of the difference between two or more
coefficients can be performed using contrasts. The contrast matrix is
evaluated for each cell type in the backend, so only the contrast string
must be supplied to dreamlet().
# create a contrasts called 'Diff' that is the difference between expression
# in the stimulated and controls.
# More than one can be specified
contrasts <- c(Diff = "StimStatusstim - StimStatusctrl")
# Evalaute the regression model without an intercept term.
# Instead estimate the mean expression in stimulated, controls and then
# set Diff to the difference between the two
res.dl2 <- dreamlet(res.proc, ~ 0 + StimStatus, contrasts = contrasts)
# see estimated coefficients
coefNames(res.dl2)## [1] "Diff" "StimStatusctrl" "StimStatusstim"
# Volcano plot of Diff
plotVolcano(res.dl2[1:2], coef = "Diff")
This new Diff variable can then be used downstream for
analysis asking for a coefficient. But note that since there is no
intercept term in this model, the meaning of StimStatusstim
changes here. When the formula is 0 + StimStatus then
StimStatusstim is the mean expression in stimulated
samples.
For further information about using contrasts see makeContrastsDream() and vignette.
Gene set analysis
While standard enrichment methods like Fishers exact test, requires
specifying a FDR cutoff to identify differentially expressed genes.
However, dichotomizing differential expression results is often too
simple and ignores the quantitative variation captured by the
differential expression test statistics. Here we use
zenith, a wrapper for limma::camera,
to perform gene set analysis using the full spectrum of differential
expression test statistics. zenith/camera is a competetive test
that compares the mean test statistic for genes in a given gene set, to
genes not in that set while accounting for correlation between
genes.
Here, zenith_gsa takes a dreamletResult
object, the coefficient of interest, and gene sets as a
GeneSetCollection object from GSEABase.
# Load Gene Ontology database
# use gene 'SYMBOL', or 'ENSEMBL' id
# use get_MSigDB() to load MSigDB
# Use Cellular Component (i.e. CC) to reduce run time here
go.gs <- get_GeneOntology("CC", to = "SYMBOL")
# Run zenith gene set analysis on result of dreamlet
res_zenith <- zenith_gsa(res.dl, coef = "StimStatusstim", go.gs)
# examine results for each ell type and gene set
head(res_zenith)## assay coef
## 1 B cells StimStatusstim
## 2 B cells StimStatusstim
## 3 B cells StimStatusstim
## 4 B cells StimStatusstim
## 5 B cells StimStatusstim
## 6 B cells StimStatusstim
## Geneset NGenes
## 1 GO0022626: cytosolic ribosome 74
## 2 GO0022625: cytosolic large ribosomal subunit 45
## 3 GO0031256: leading edge membrane 20
## 4 GO0090575: RNA polymerase II transcription regulator complex 37
## 5 GO0000151: ubiquitin ligase complex 34
## 6 GO0005839: proteasome core complex 13
## Correlation delta se p.less p.greater PValue Direction
## 1 0.01 -1.359918 0.4149538 0.002752843 0.997247157 0.005505685 Down
## 2 0.01 -1.374381 0.4847978 0.006618402 0.993381598 0.013236805 Down
## 3 0.01 -1.858272 0.6589135 0.006813596 0.993186404 0.013627192 Down
## 4 0.01 1.424376 0.5191995 0.992075627 0.007924373 0.015848745 Up
## 5 0.01 1.392016 0.5355516 0.989495380 0.010504620 0.021009240 Up
## 6 0.01 1.993645 0.7922932 0.987661570 0.012338430 0.024676860 Up
## FDR
## 1 0.2741364
## 2 0.4283927
## 3 0.4283927
## 4 0.4289412
## 5 0.4407344
## 6 0.4693180Heatmap of top genesets
# for each cell type select 5 genesets with largest t-statistic
# and 1 geneset with the lowest
# Grey boxes indicate the gene set could not be evaluted because
# to few genes were represented
plotZenithResults(res_zenith, 5, 1)
All gene sets with FDR < 30%
Here, show all genes with FDR < 5% in any cell type
# get genesets with FDR < 30%
# Few significant genesets because uses Cellular Component (i.e. CC)
gs <- unique(res_zenith$Geneset[res_zenith$FDR < 0.3])
# keep only results of these genesets
df <- res_zenith[res_zenith$Geneset %in% gs, ]
# plot results, but with no limit based on the highest/lowest t-statistic
plotZenithResults(df, Inf, Inf)
Comparing expression in clusters
Identifying genes that are differentially expressed between cell clusters incorporates a paired analysis design, since each individual is observed for each cell cluster.
# test differential expression between B cells and the rest of the cell clusters
ct.pairs <- c("CD4 T cells", "rest")
fit <- dreamletCompareClusters(pb, ct.pairs, method = "fixed")
# The coefficient 'compare' is the value logFC between test and baseline:
# compare = cellClustertest - cellClusterbaseline
df_Bcell <- topTable(fit, coef = "compare")
head(df_Bcell)## logFC AveExpr t P.Value adj.P.Val
## TYROBP -7.422120 7.903788 -90.89526 9.293731e-32 4.368054e-28
## FTL -4.940554 13.169295 -66.65308 1.450833e-28 3.409458e-25
## CD63 -5.083331 8.757645 -59.82930 1.867983e-27 2.926507e-24
## C15orf48 -7.254208 8.381008 -57.06268 5.719817e-27 6.720785e-24
## HLA-DRA_ENSG00000204287 -6.420675 8.859898 -52.69037 3.758019e-26 3.532538e-23
## CCL2 -7.278426 8.712962 -52.08993 4.925158e-26 3.858040e-23
## B
## TYROBP 61.17806
## FTL 55.57623
## CD63 52.88828
## C15orf48 51.66970
## HLA-DRA_ENSG00000204287 49.96082
## CCL2 49.63850Gene-cluster specificity
Evaluate the specificity of each gene for each cluster, retaining only highly expressed genes:
df_cts <- cellTypeSpecificity(pb)
# retain only genes with total CPM summed across cell type > 100
df_cts <- df_cts[df_cts$totalCPM > 100, ]
# Violin plot of specificity score for each cell type
plotViolin(df_cts)
Highlight expression fraction for most specific gene from each cell type:
genes <- rownames(df_cts)[apply(df_cts, 2, which.max)]
plotPercentBars(df_cts, genes = genes)
dreamlet::plotHeatmap(df_cts, genes = genes)
Session Info
## R version 4.3.0 (2023-04-21)
## Platform: x86_64-apple-darwin22.4.0 (64-bit)
## Running under: macOS Ventura 13.5
##
## Matrix products: default
## BLAS: /Users/gabrielhoffman/prog/R-4.3.0/lib/libRblas.dylib
## LAPACK: /usr/local/Cellar/r/4.3.0_1/lib/R/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: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] GO.db_3.17.0 org.Hs.eg.db_3.17.0
## [3] AnnotationDbi_1.62.1 muscData_1.14.0
## [5] scater_1.28.0 scuttle_1.10.1
## [7] zenith_1.4.1 ExperimentHub_2.8.1
## [9] AnnotationHub_3.8.0 BiocFileCache_2.8.0
## [11] dbplyr_2.3.2 muscat_1.14.0
## [13] dreamlet_1.1.5 SingleCellExperiment_1.22.0
## [15] SummarizedExperiment_1.30.1 Biobase_2.60.0
## [17] GenomicRanges_1.52.0 GenomeInfoDb_1.36.1
## [19] IRanges_2.34.1 S4Vectors_0.38.1
## [21] BiocGenerics_0.46.0 MatrixGenerics_1.12.0
## [23] matrixStats_1.0.0 variancePartition_1.33.2
## [25] BiocParallel_1.34.2 limma_3.56.2
## [27] ggplot2_3.4.4 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.2 bitops_1.0-7
## [3] httr_1.4.6 RColorBrewer_1.1-3
## [5] doParallel_1.0.17 Rgraphviz_2.44.0
## [7] numDeriv_2016.8-1.1 tools_4.3.0
## [9] sctransform_0.3.5 backports_1.4.1
## [11] utf8_1.2.3 R6_2.5.1
## [13] GetoptLong_1.0.5 withr_2.5.0
## [15] prettyunits_1.1.1 gridExtra_2.3
## [17] cli_3.6.1 textshaping_0.3.6
## [19] sandwich_3.0-2 labeling_0.4.2
## [21] sass_0.4.6 KEGGgraph_1.60.0
## [23] SQUAREM_2021.1 mvtnorm_1.2-2
## [25] blme_1.0-5 pkgdown_2.0.7
## [27] mixsqp_0.3-48 systemfonts_1.0.4
## [29] parallelly_1.36.0 invgamma_1.1
## [31] RSQLite_2.3.1 generics_0.1.3
## [33] shape_1.4.6 gtools_3.9.4
## [35] dplyr_1.1.2 Matrix_1.5-4.1
## [37] ggbeeswarm_0.7.2 fansi_1.0.4
## [39] abind_1.4-5 lifecycle_1.0.3
## [41] multcomp_1.4-23 yaml_2.3.7
## [43] edgeR_3.42.4 gplots_3.1.3
## [45] grid_4.3.0 blob_1.2.4
## [47] promises_1.2.0.1 crayon_1.5.2
## [49] lattice_0.21-8 beachmat_2.16.0
## [51] msigdbr_7.5.1 annotate_1.78.0
## [53] KEGGREST_1.40.0 pillar_1.9.0
## [55] knitr_1.43 ComplexHeatmap_2.16.0
## [57] rjson_0.2.21 boot_1.3-28.1
## [59] estimability_1.4.1 corpcor_1.6.10
## [61] future.apply_1.11.0 codetools_0.2-19
## [63] glue_1.6.2 data.table_1.14.8
## [65] vctrs_0.6.3 png_0.1-8
## [67] Rdpack_2.4 gtable_0.3.3
## [69] assertthat_0.2.1 cachem_1.0.8
## [71] xfun_0.39 mime_0.12
## [73] rbibutils_2.2.13 S4Arrays_1.0.4
## [75] Rfast_2.0.7 coda_0.19-4
## [77] survival_3.5-5 iterators_1.0.14
## [79] ellipsis_0.3.2 interactiveDisplayBase_1.38.0
## [81] TH.data_1.1-2 nlme_3.1-162
## [83] pbkrtest_0.5.2 bit64_4.0.5
## [85] filelock_1.0.2 progress_1.2.2
## [87] EnvStats_2.7.0 rprojroot_2.0.3
## [89] bslib_0.4.2 TMB_1.9.4
## [91] irlba_2.3.5.1 vipor_0.4.5
## [93] KernSmooth_2.23-21 colorspace_2.1-0
## [95] rmeta_3.0 DBI_1.1.3
## [97] DESeq2_1.40.1 tidyselect_1.2.0
## [99] emmeans_1.8.7 curl_5.0.0
## [101] bit_4.0.5 compiler_4.3.0
## [103] graph_1.78.0 BiocNeighbors_1.18.0
## [105] desc_1.4.2 DelayedArray_0.26.3
## [107] bookdown_0.34 scales_1.2.1
## [109] caTools_1.18.2 remaCor_0.0.17
## [111] rappdirs_0.3.3 stringr_1.5.0
## [113] digest_0.6.33 minqa_1.2.5
## [115] rmarkdown_2.22 aod_1.3.2
## [117] XVector_0.40.0 RhpcBLASctl_0.23-42
## [119] htmltools_0.5.5 pkgconfig_2.0.3
## [121] lme4_1.1-34 sparseMatrixStats_1.12.0
## [123] highr_0.10 mashr_0.2.69
## [125] fastmap_1.1.1 rlang_1.1.1
## [127] GlobalOptions_0.1.2 shiny_1.7.4
## [129] DelayedMatrixStats_1.22.0 farver_2.1.1
## [131] jquerylib_0.1.4 zoo_1.8-12
## [133] jsonlite_1.8.5 BiocSingular_1.16.0
## [135] RCurl_1.98-1.12 magrittr_2.0.3
## [137] GenomeInfoDbData_1.2.10 munsell_0.5.0
## [139] Rcpp_1.0.11 babelgene_22.9
## [141] viridis_0.6.3 EnrichmentBrowser_2.30.1
## [143] RcppZiggurat_0.1.6 stringi_1.7.12
## [145] zlibbioc_1.46.0 MASS_7.3-60
## [147] plyr_1.8.8 listenv_0.9.0
## [149] parallel_4.3.0 ggrepel_0.9.3
## [151] Biostrings_2.68.1 splines_4.3.0
## [153] hms_1.1.3 circlize_0.4.15
## [155] locfit_1.5-9.7 reshape2_1.4.4
## [157] ScaledMatrix_1.8.1 BiocVersion_3.17.1
## [159] XML_3.99-0.14 evaluate_0.21
## [161] BiocManager_1.30.20 httpuv_1.6.11
## [163] nloptr_2.0.3 foreach_1.5.2
## [165] tidyr_1.3.0 purrr_1.0.2
## [167] future_1.32.0 clue_0.3-64
## [169] scattermore_1.1 ashr_2.2-54
## [171] rsvd_1.0.5 broom_1.0.5
## [173] xtable_1.8-4 fANCOVA_0.6-1
## [175] later_1.3.1 viridisLite_0.4.2
## [177] ragg_1.2.5 truncnorm_1.0-9
## [179] tibble_3.2.1 lmerTest_3.1-3
## [181] glmmTMB_1.1.7 memoise_2.0.1
## [183] beeswarm_0.4.0 cluster_2.1.4
## [185] globals_0.16.2 GSEABase_1.62.0