Last updated: 2024-06-12

Checks: 7 0

Knit directory: spatialsnippets/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20231017) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 74dbb4e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    renv/library/
    Ignored:    renv/staging/

Untracked files:
    Untracked:  analysis/e_spatiallyRestrictedGenes.Rmd

Unstaged changes:
    Modified:   analysis/d_cosmxAlzPlaque.Rmd
    Modified:   analysis/index.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/e_DEPseudobulk_insitu.Rmd) and HTML (docs/e_DEPseudobulk_insitu.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 74dbb4e swbioinf 2024-06-12 wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd")
html 9cd7910 swbioinf 2024-06-11 Build site.
Rmd 561981e swbioinf 2024-06-11 wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd")
html 61d4499 swbioinf 2024-05-31 Build site.
Rmd cf7833a swbioinf 2024-05-31 wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd")
html 6978154 swbioinf 2024-05-30 Build site.
html 97393ec swbioinf 2024-05-28 Build site.
Rmd dd0f089 swbioinf 2024-05-28 wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd")
html d6d841a swbioinf 2024-05-22 Build site.
Rmd 9c53965 swbioinf 2024-05-22 wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd")
html 58ee2ec swbioinf 2024-05-22 Build site.
Rmd 36dd228 swbioinf 2024-05-22 wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd")
html 5a9d7e9 swbioinf 2024-05-16 Build site.
Rmd c67ef86 swbioinf 2024-05-16 wflow_publish("analysis/e_DEPseudobulk_insitu.Rmd")
html c77f76c swbioinf 2024-05-07 Build site.
Rmd 7b0be93 swbioinf 2024-05-07 wflow_publish(c("analysis/index.Rmd", "analysis/e_DEPseudobulk_insitu.Rmd",

Overview

Once we have identified cell types present in the samples, its typical to test how gene expression changes between experimental conditions, within each different cell type. Some cell types may be dramatically affected by the experimental conditions, while others are not. Likewise some genes may change only in a specific cell type, whereas others show a more general difference.

This document describes how to apply a pseudobulk approach to test for differences between groups, accounting for biological replicates. This is very similar to how a non-spatial single cell experiment may be analysed.

In a pseudobulk approach counts are obtaing by pooling together groups of cells; in this case cells from the of the same type from the same fov. These pooled counts can then be analysed more like a bulk RNAseq experiment.

Note that there are other approaches to calculate differential expression in this kind of data - including those that make use of individual cells (Soneson and Robinson 2018).

This requires:

  • Biological replicates for each group
  • Assigned cell types
  • [Optionally] Multiple fovs measured per sample

For example:

  • What genes are differentially expressed in epithelial cells in Crohn’s disease vs healthy individuals?
  • How do genes change with treatment in each different cell type in my sample?

Steps:

  1. Calculate pseudobulk
  2. Filter to testable pseudobulk groups (enough cells to pool)
  3. Filter to testable genes (enough expression to see changes)
  4. Test for changes in gene expression
  5. Plot DE results and individual genes.

Worked example

How does gene expression change within each cell type between Ulcerative colitis or Crohn’s disease, and Healthy controls?

Using data from (Garrido-Trigo et al. 2023)

Load libraries and data

library(Seurat)
library(speckle)
library(tidyverse)
library(limma)
library(DT)
library(edgeR)
data_dir              <- file.path("~/projects/spatialsnippets/datasets/GSE234713_IBDcosmx_GarridoTrigo2023/processed_data") 
seurat_file_01_loaded <- file.path(data_dir, "GSE234713_CosMx_IBD_seurat_01_loaded.RDS")
so <- readRDS(seurat_file_01_loaded)

Experimental design

There are three individuals per condition (one tissue sample from each individual). With multiple fovs on each physical tissue sample.

sample_table <- select(as_tibble(so@meta.data), condition, individual_code, fov_name) %>%
  unique() %>% 
  group_by(condition, individual_code) %>% 
  summarise(n_fovs= n(), item = str_c(fov_name, collapse = ", "))

DT::datatable(sample_table)

Count how many cells of each type in your data

Using a pseudobulk approach.

  • Need at least x reads in a cell to include it
  • Need at least x cells of a celltype within an fov to include that
  • Can only test where we have at least 2 samples on each side of a contrast.
min_reads_per_cell <- 200

ggplot(so@meta.data, aes(x=nCount_RNA)) +
  geom_density() +
  geom_vline(xintercept = min_reads_per_cell, lty=3) +
  scale_x_log10() +
  theme_bw()+
  ggtitle("How many reads per cell?")

Version Author Date
c77f76c swbioinf 2024-05-07
so<- so[,so$nCount_RNA >= min_reads_per_cell]

We will pool each celltype within each fov (cluster_group). But there needs to be a certain number of cells for that to work - less than a certain number of cells and a pseudobulk pool will be excluded.

Note there are much fewer t-cells overall, but given that we have a high number of samples, there should still be enough to include. Its typical that some of the less common cell types are difficult or impossible to reliably test.

min_cells_per_fovcluster <- 20

so$fov_cluster <- paste0(so$fov_name,"_", so$celltype_subset)

celltype_summary_table <- so@meta.data %>% 
  group_by(condition, group, individual_code, fov_name, celltype_subset, fov_cluster) %>%
  summarise(cells=n(), .groups = 'drop')
DT::datatable(celltype_summary_table)
ggplot(celltype_summary_table, aes(x=cells, col=celltype_subset)) +
  geom_density() +
  geom_vline(xintercept=min_cells_per_fovcluster, lty=3) +
  geom_rug() +
  scale_x_log10() +
  theme_bw() +
  ggtitle("How many cells per fov-cluster?")

Version Author Date
6978154 swbioinf 2024-05-30
c77f76c swbioinf 2024-05-07

Record the names of those fov_clusters that contain enough cells to be used. Will use this later to filter.

passed_fov_clusters <- celltype_summary_table$fov_cluster[celltype_summary_table$cells >= min_cells_per_fovcluster]

Calculate pseudobulk

Now use the PseudobulkExpression() function to sum up each gene’s gene expression across each fov_cluster.

pseudobulk_counts <- PseudobulkExpression(so, assays = "RNA", layer="counts",  method = 'aggregate', group.by = 'fov_cluster')
pseudobulk_counts_matrix <- pseudobulk_counts[["RNA"]]

# Change - back to _. Ideally we'd have neither, but - will cause problems later
colnames(pseudobulk_counts_matrix)<-gsub("-","_",colnames(pseudobulk_counts_matrix))

Now instead of counts for each individual cell, we have our pseodubulk matrix with the pooled sum of counts for each celltype within each fov region. So the numbers are much higher, fewer zeros. But each pool is a different size, which we will adress later with normalisation.

pseudobulk_counts_matrix[1:10,1:4]
10 x 4 sparse Matrix of class "dgCMatrix"
      CD_a_001_epi CD_a_001_myeloids CD_a_001_plasmas CD_a_001_stroma
AATK            38                 4               16              13
ABL1            65                 1                6              18
ABL2            33                 2               11              10
ACE             28                 .                4               6
ACE2            59                 .               17              18
ACKR1           16                 1               13              21
ACKR3           37                 4               12              15
ACKR4           36                 2                9               7
ACTA2           23                 1                7              45
ACTG2           32                 2               15              27

Next, filter to only those passed fov_clusters we saw earlier - this is what we’ll use for calculating differential expression.

We also need an annotation table to tell us what is in each fov_cluster (sample, celltype, condition…). Build this by filtering the celltype summary table made earlier. If we filter it to the same samples, in the same order as the pseudubulk matrix, we can then use those columns of data to build our differential expression model in the next step.

# Filter the pseudobulk table
pseudobulk_counts_matrix <- pseudobulk_counts_matrix[,passed_fov_clusters]

# And filter the celltype summary into an annotation table for only those fov_clusters
# pull in relevant annotation in a matched order
pseudobulk_anno_table <- celltype_summary_table
match_order <- match(passed_fov_clusters, pseudobulk_anno_table$fov_cluster)
pseudobulk_anno_table <- pseudobulk_anno_table[match_order,]

# Double check that the order of samples in the annotation table matches the pseudobulk table. 
# If this is wrong then the results will be nonsense!
stopifnot(all(colnames(pseudobulk_counts_matrix) == pseudobulk_anno_table$fov_cluster  ))

Calculate Differential Expression

We now have a counts matrix, and an annotation table that describes its samples. This looks very much like a bulk RNAseq experiment.

pseudobulk_counts_matrix[1:10,1:4]
10 x 4 sparse Matrix of class "dgCMatrix"
      CD_a_001_epi CD_a_001_plasmas CD_a_001_stroma CD_a_002_myeloids
AATK            38               16              13                 8
ABL1            65                6              18                 3
ABL2            33               11              10                11
ACE             28                4               6                 3
ACE2            59               17              18                 4
ACKR1           16               13              21                 4
ACKR3           37               12              15                 2
ACKR4           36                9               7                 3
ACTA2           23                7              45                 6
ACTG2           32               15              27                 3
head(pseudobulk_anno_table)
# A tibble: 6 × 7
  condition     group individual_code fov_name celltype_subset fov_cluster cells
  <fct>         <fct> <chr>           <fct>    <fct>           <chr>       <int>
1 Crohn's dise… CD    CD_a            CD_a_001 epi             CD_a_001_e…   279
2 Crohn's dise… CD    CD_a            CD_a_001 plasmas         CD_a_001_p…    99
3 Crohn's dise… CD    CD_a            CD_a_001 stroma          CD_a_001_s…    78
4 Crohn's dise… CD    CD_a            CD_a_002 myeloids        CD_a_002_m…    28
5 Crohn's dise… CD    CD_a            CD_a_002 stroma          CD_a_002_s…   230
6 Crohn's dise… CD    CD_a            CD_a_002 tcells          CD_a_002_t…    21

We will use a fairly standard limma differential expression analysis.

We will process each celltype one at a time as follows:

  1. Subset the counts matrix and annotation table to the celltype
  2. Check that we have enough replicates for that celltype to run the test(s). We need to have fovs from at least 2 biological samples on each side of our test (e.g. 2treatment, 2 control). Even though we have enough tissue samples, if its a less common cell type we might have filtered them away!
  3. Build a model. In our case, a simple one that looks for variouat by disease group, blocked on individual.
  4. Run the test.

How to build the model and run the test will need to be customised foryour experimental design. Online resources that describe approaches for bulk RNAseq analyses can be applied to pseubulk analyses - see the ‘More Information’ section of this document for suggestions.

In this case the line model.matrix( ~0 + group) builds a model where gene expression is expected to vary by group (UC/CD/HC). The ~0 part indicates a ‘intercept’, which is convenient becase it means we can include all 3 groups in our contrasts, rather than one of them being treated as the baseline. ( NB: Mathematically, you could use ~group alone without intercept, and define contrasts accordingly to get the same result. )

But, we also need to tae into acoount our individuals. We do this by blocking on individual, and calculageing a intra duplicate correlation with duplicateCorrelation().

If instead we didn’t have multiple fovs per replicate, we might use a model like ~0 + group + individual; no blocking or use of duplicateCorrelation(). That is more typical of a typical of a bulk rnaseq or single cell RNAseq pseuobulk analysis.

min_samples_to_calc <- 2 # require 2 samples on on either side of contrast

de_result_list <- list()

# celltype_subset is a matrix
for (the_celltype in levels(so$celltype_subset)) {
  
  anno_table.this   <- pseudobulk_anno_table[pseudobulk_anno_table$celltype_subset == the_celltype,]
  count_matrix.this <- pseudobulk_counts_matrix[,anno_table.this$fov_cluster]
 
  print(the_celltype)
  
  # skip clusters with nothing 
  if( nrow(anno_table.this) < 1 ) {next}
  
  # Setup objects for limma
  dge <- DGEList(count_matrix.this)
  dge <- calcNormFactors(dge)
  
  # Build model
  group           <- anno_table.this$group
  individual_code <- anno_table.this$individual_code

  # To do any calculations, we need at least 2 pseudobulk groups per contrast.
  # there are plenty in this experiemnt, but with less replicates and rare cell types
  # it may be neccesary to check and skip certain contrasts

  # Model design 
  design    <- model.matrix( ~0 + group)

  vm  <- voom(dge, design = design, plot = FALSE)
  
  # Adding dupliate correlation to use individual fovs, rather than pooled per biosample
  corrfit <- duplicateCorrelation(vm, design, block=individual_code)
  fit     <- lmFit(vm, design, correlation = corrfit$consensus, block=individual_code) 

  # Then fit contrasts and run ebayes
  contrasts <- makeContrasts(UCvHC  = groupUC - groupHC, 
                             CDvHC  = groupCD - groupHC,
                           levels=coef(fit))

  fit <- contrasts.fit(fit, contrasts)
  fit <- eBayes(fit)

  
  for ( the_coef in colnames(contrasts) ) {
    de_result.this <- topTable(fit, n = Inf, adjust.method = "BH", coef = the_coef) %>%
      rownames_to_column("target") %>%
      mutate(contrast=the_coef,
             celltype=the_celltype) %>%
      select(celltype,contrast,target,everything()) %>%
      arrange(P.Value)
    
    
      de_result_list[[paste(the_celltype, the_coef, sep="_")]] <- de_result.this
    
  }

  

 
}
[1] "epi"
[1] "myeloids"
[1] "plasmas"
[1] "stroma"
[1] "tcells"
# Join together results for all celltypes, and pull out those with a singificant adjusted p-value
de_results_all <- bind_rows(de_result_list)
de_results_sig <- filter(de_results_all, adj.P.Val < 0.01)

Table of significant results.

DT::datatable(de_results_sig)

DE plots

The below plots show the logFC calculate for each gene versus its average expression across all samples. This is a useful diagnostic plot to evaluate your differential expression results.

  • Only big changes are significant at lower expressions, simply becuase higher expression means more statistical confidence.
  • A lopsided plot might indicate some oddness around the normalisation of your data (e.g. extremely differnet cell counts).

[[ADD pvalue plot?]]

library(ggrepel) # gg_repel, For non-overlapping gene labels


make_ma_style_plot <- function(res_table, pval_threshold = 0.01, n_genes_to_label = 10) {
  p <- ggplot(res_table, aes(x=AveExpr, y=logFC, col=adj.P.Val < pval_threshold) ) +
    geom_hline(yintercept = c(0), col='grey80') +
    geom_point(pch=3) +
    geom_text_repel(data    = head(arrange(filter(res_table , adj.P.Val < pval_threshold ), P.Value), n=5),
                    mapping = aes(label=target), col="red" ) +
    theme_bw() +
    geom_hline(yintercept = c(-1,1), lty=3) +
    scale_colour_manual(values = c('FALSE'="black", 'TRUE'="red")) +
    theme(legend.position = 'none')
  return(p)
}
#res_table.UCvHC.epi <- filter(de_results_all, contrast == "UCvHC", celltype=="epi")

make_ma_style_plot(res_table = filter(de_results_all, contrast == "UCvHC", celltype=="epi"))

Version Author Date
6978154 swbioinf 2024-05-30
5a9d7e9 swbioinf 2024-05-16
make_ma_style_plot(res_table = filter(de_results_all, contrast == "UCvHC", celltype=="tcells"))

make_ma_style_plot(res_table = filter(de_results_all, contrast == "UCvHC", celltype=="stroma"))

Check some examples

Its always worth visualising how the expression of your differentially expressed genes really looks, with respect to your experimental design. How best to do this depends on your experiment.

The results suggests that TNFRSF18 was significantly DE between individuals with Ulcerative Colitis and Healthy Controls in plasma cells. There’s some very convenient seurat plots below;

VlnPlot(subset(so, celltype_subset == "plasmas"), features = "TNFRSF18", group.by = 'group', alpha = 0.1)

Version Author Date
6978154 swbioinf 2024-05-30
5a9d7e9 swbioinf 2024-05-16
FeaturePlot(so, "TNFRSF18", split.by = "group")

But it gets difficult to summarise data at the single cell level. We can also use the the normalised pseudobulk expression to see how gene expression varies within each fov,individual,celltype and condition - The plot below shows an overview of TNFRSF18 across the entire experiment.

NB: We have to plot normalised expression instead of the raw counts as there are vastly different numbers of cells in each fov+celltype grouping.

# Get tmm normalised coutns for all pseudobulk
# WHen we did the DE we calculated this a celltype at a time, so values might differ slightly!
dge <- DGEList(pseudobulk_counts_matrix)
dge <- calcNormFactors(dge)
norm_pseudobulk <- cpm(dge , log=TRUE) # uses tmm normalisation

# Plot expression for TNFRSF18
plottable <- cbind(pseudobulk_anno_table, expression = norm_pseudobulk["TNFRSF18",])
ggplot(plottable, aes(x=individual_code, y=expression, col=condition )) +
  geom_boxplot(outlier.shape = NA) +
  geom_point() + 
  theme_bw() + 
  theme(axis.text.x=element_text(angle = -90, hjust = 0)) + 
  facet_wrap(~celltype_subset) 

And you can compare that with IGHG1, which was flagged at differentially expressed across multiple cell types.

plottable <- cbind(pseudobulk_anno_table, expression = norm_pseudobulk["IGHG1",])
ggplot(plottable, aes(x=individual_code, y=expression, col=condition )) +
  geom_boxplot(outlier.shape = NA) +
  geom_point() + 
  theme_bw() + 
  theme(axis.text.x=element_text(angle = -90, hjust = 0)) + 
  facet_wrap(~celltype_subset) 

Code Snippet

library(Seurat)
library(edgeR)
library(limma)


min_reads_per_cell       <- 200
min_cells_per_fovcluster <- 20
min_samples_to_calc      <- 2 # require 2 samples on on either side of contrast


# Remove cells with too few counts
so <- so[,so$nCount_RNA >= min_reads_per_cell]


# Define fov+cluster groups, with all relevant sample annotation
# remove those with too few cells.
so$fov_cluster <- paste0(so$fov_name,"_", so$celltype_subset)
celltype_summary_table <- so@meta.data %>% 
  group_by(condition, group, individual_code, fov_name, celltype_subset, fov_cluster) %>%
  summarise(cells=n(), .groups = 'drop')


# Calculate pseudobulk for each fov+cluster group
pseudobulk_counts <- PseudobulkExpression(so, assays = "RNA", layer="counts",  method = 'aggregate', group.by = 'fov_cluster')
pseudobulk_counts_matrix <- pseudobulk_counts[["RNA"]]
# Change - back to _. Ideally have neither and skip this step
colnames(pseudobulk_counts_matrix)<-gsub("-","_",colnames(pseudobulk_counts_matrix))


# Determin fov_clusters with eough cells
# Filter both pseudobulk matrix and pseudobulk annotation 
passed_fov_clusters <- celltype_summary_table$fov_cluster[celltype_summary_table$cells >= min_cells_per_fovcluster]
pseudobulk_counts_matrix <- pseudobulk_counts_matrix[,passed_fov_clusters]
pseudobulk_anno_table    <- celltype_summary_table[passed_fov_clusters,]


# Calculate DE across every celltype
de_result_list <- list()

for (the_celltype in levels(so$celltype_subset)) {
  
  anno_table.this   <- pseudobulk_anno_table[pseudobulk_anno_table$celltype_subset == the_celltype,]
  count_matrix.this <- pseudobulk_counts_matrix[,anno_table.this$fov_cluster]
 
  print(the_celltype)
  
  # skip clusters with nothing 
  if( nrow(anno_table.this) < 1 ) {next}
  
  # Setup objects for limma
  dge <- DGEList(count_matrix.this)
  dge <- calcNormFactors(dge)
  
  # Build model
  group           <- anno_table.this$group
  individual_code <- anno_table.this$individual_code

  # To do any calculations, we need at least 2 pseudobulk groups per contrast.
  # there are plenty in this experiemnt, but with less replicates and rare cell types
  # it may be neccesary to check and skip certain contrasts

  # Model design 
  design    <- model.matrix( ~0 + group)

  vm  <- voom(dge, design = design, plot = FALSE)
  
  # Adding dupliate correlation to use individual fovs, rather than pooled per biosample
  corrfit <- duplicateCorrelation(vm, design, block=individual_code)
  fit     <- lmFit(vm, design, correlation = corrfit$consensus, block=individual_code) 

  # Then fit contrasts and run ebayes
  contrasts <- makeContrasts(UCvHC  = groupUC - groupHC
                           levels=coef(fit))

  fit <- contrasts.fit(fit, contrasts)
  fit <- eBayes(fit)

  
  for ( the_coef in colnames(contrasts) ) {
    de_result.this <- topTable(fit, n = Inf, adjust.method = "BH", coef = the_coef) %>%
      rownames_to_column("target") %>%
      mutate(contrast=the_coef,
             celltype=the_celltype) %>%
      select(celltype,contrast,target,everything()) %>%
      arrange(P.Value)
    
    
      de_result_list[[paste(the_celltype, the_coef, sep="_")]] <- de_result.this
    
  }

}

de_results_all <- bind_rows(de_result_list)
de_results_sig <- filter(de_results_all, adj.P.Val < 0.01)

Results

DT::datatable(head(de_results_all))

This table is the typical output of limma tests; With a couple of extra columns added by our code.

  • celltype: The celltype being tested (Added by example code)
  • contrast: The contrast being tested (Added by example code)
  • target : The gene name (Added by example code, is the rowname in limma output)
  • rownames : The tested cell types
  • logFC : Log 2 fold change between tested groups. For a test of Test-Con;
    • At logFC +1, A is doubled B.
    • At logFC -1, A is half of B.
    • A logFC 0 indicates no change.
  • AveExpr : Average expression of a gene across all replicates.
  • t : Moderated T-statistic. See Limma documentation.
  • P.Value : P.value
  • adj.P.Val : A multiple-hypothesis corrected p-value
  • B : B statistic (rarely used). See Limma documentation.

More Information

References

Garrido-Trigo, Alba, Ana M. Corraliza, Marisol Veny, Isabella Dotti, Elisa Melón-Ardanaz, Aina Rill, Helena L. Crowell, et al. 2023. “Macrophage and Neutrophil Heterogeneity at Single-Cell Spatial Resolution in Human Inflammatory Bowel Disease.” Nature Communications 14 (1): 4506. https://doi.org/10.1038/s41467-023-40156-6.
Ritchie, Matthew E., Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, and Gordon K. Smyth. 2015. “Limma Powers Differential Expression Analyses for RNA-sequencing and Microarray Studies.” Nucleic Acids Research 43 (7): e47. https://doi.org/10.1093/nar/gkv007.
Soneson, Charlotte, and Mark D Robinson. 2018. “Bias, Robustness and Scalability in Single-Cell Differential Expression Analysis.” Nature Methods 15 (4): 255–61. https://doi.org/10.1038/nmeth.4612.

sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.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.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
 [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8   
 [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] ggrepel_0.9.5      edgeR_4.0.16       DT_0.33            limma_3.58.1      
 [5] lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
 [9] purrr_1.0.2        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[13] ggplot2_3.5.0      tidyverse_2.0.0    speckle_1.2.0      Seurat_5.0.3      
[17] SeuratObject_5.0.1 sp_2.1-3           workflowr_1.7.1   

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22            splines_4.3.2              
  [3] later_1.3.2                 bitops_1.0-7               
  [5] polyclip_1.10-6             fastDummies_1.7.3          
  [7] lifecycle_1.0.4             rprojroot_2.0.4            
  [9] globals_0.16.3              processx_3.8.4             
 [11] lattice_0.22-6              MASS_7.3-60.0.1            
 [13] crosstalk_1.2.1             magrittr_2.0.3             
 [15] plotly_4.10.4               sass_0.4.9                 
 [17] rmarkdown_2.26              jquerylib_0.1.4            
 [19] yaml_2.3.8                  httpuv_1.6.15              
 [21] sctransform_0.4.1           spam_2.10-0                
 [23] spatstat.sparse_3.0-3       reticulate_1.35.0          
 [25] cowplot_1.1.3               pbapply_1.7-2              
 [27] RColorBrewer_1.1-3          abind_1.4-5                
 [29] zlibbioc_1.48.2             Rtsne_0.17                 
 [31] GenomicRanges_1.54.1        BiocGenerics_0.48.1        
 [33] RCurl_1.98-1.14             git2r_0.33.0               
 [35] GenomeInfoDbData_1.2.11     IRanges_2.36.0             
 [37] S4Vectors_0.40.2            irlba_2.3.5.1              
 [39] listenv_0.9.1               spatstat.utils_3.0-4       
 [41] goftest_1.2-3               RSpectra_0.16-1            
 [43] spatstat.random_3.2-3       fitdistrplus_1.1-11        
 [45] parallelly_1.37.1           leiden_0.4.3.1             
 [47] codetools_0.2-20            DelayedArray_0.28.0        
 [49] tidyselect_1.2.1            farver_2.1.1               
 [51] matrixStats_1.2.0           stats4_4.3.2               
 [53] spatstat.explore_3.2-7      jsonlite_1.8.8             
 [55] progressr_0.14.0            ggridges_0.5.6             
 [57] survival_3.5-8              tools_4.3.2                
 [59] ica_1.0-3                   Rcpp_1.0.12                
 [61] glue_1.7.0                  gridExtra_2.3              
 [63] SparseArray_1.2.4           xfun_0.43                  
 [65] MatrixGenerics_1.14.0       GenomeInfoDb_1.38.8        
 [67] withr_3.0.0                 BiocManager_1.30.22        
 [69] fastmap_1.1.1               fansi_1.0.6                
 [71] callr_3.7.6                 digest_0.6.35              
 [73] timechange_0.3.0            R6_2.5.1                   
 [75] mime_0.12                   colorspace_2.1-0           
 [77] scattermore_1.2             tensor_1.5                 
 [79] spatstat.data_3.0-4         utf8_1.2.4                 
 [81] generics_0.1.3              renv_1.0.5                 
 [83] data.table_1.15.4           httr_1.4.7                 
 [85] htmlwidgets_1.6.4           S4Arrays_1.2.1             
 [87] whisker_0.4.1               uwot_0.1.16                
 [89] pkgconfig_2.0.3             gtable_0.3.4               
 [91] lmtest_0.9-40               SingleCellExperiment_1.24.0
 [93] XVector_0.42.0              htmltools_0.5.8.1          
 [95] dotCall64_1.1-1             scales_1.3.0               
 [97] Biobase_2.62.0              png_0.1-8                  
 [99] knitr_1.45                  rstudioapi_0.16.0          
[101] tzdb_0.4.0                  reshape2_1.4.4             
[103] nlme_3.1-164                cachem_1.0.8               
[105] zoo_1.8-12                  KernSmooth_2.23-22         
[107] parallel_4.3.2              miniUI_0.1.1.1             
[109] pillar_1.9.0                grid_4.3.2                 
[111] vctrs_0.6.5                 RANN_2.6.1                 
[113] promises_1.2.1              xtable_1.8-4               
[115] cluster_2.1.6               evaluate_0.23              
[117] cli_3.6.2                   locfit_1.5-9.9             
[119] compiler_4.3.2              rlang_1.1.3                
[121] crayon_1.5.2                future.apply_1.11.2        
[123] labeling_0.4.3              ps_1.7.6                   
[125] getPass_0.2-4               plyr_1.8.9                 
[127] fs_1.6.3                    stringi_1.8.3              
[129] viridisLite_0.4.2           deldir_2.0-4               
[131] munsell_0.5.1               lazyeval_0.2.2             
[133] spatstat.geom_3.2-9         Matrix_1.6-5               
[135] RcppHNSW_0.6.0              hms_1.1.3                  
[137] patchwork_1.2.0             future_1.33.2              
[139] statmod_1.5.0               shiny_1.8.1.1              
[141] highr_0.10                  SummarizedExperiment_1.32.0
[143] ROCR_1.0-11                 igraph_2.0.3               
[145] bslib_0.7.0