¶ EBAII n1 2024
Differential expression analysis ¶

thibault.dayris@gustaveroussy.fr bastien.job@gustaveroussy.fr

1 Forewords

I need three teams for this session: team Wilcoxon, team Student, and team ROC.

I will also need your help because I can’t make DESeq2 work correctly. But I’m sure, that we will solve my issue: you’re in the best session here at EBAII.

1.1 TLDR: R command lines

In this presentation, there will be screen captures for you to follow the lesson. There will also be every single R command lines. Do not take care of the command lines if you find them too challenging. Our goal here, is to understand the main mechanism of Differential Expression Analysis. R is just a tool.

Below are the libraries we need to perform this whole session:

base::library(package = "BiocParallel")    # Optionally multithread some steps
base::library(package = "DT")              # Display nice table in HTML
base::library(package = "ggplot2")         # Draw graphs and plots
base::library(package = "ggpubr")          # Draw nicer graphs
base::library(package = "rstatix")         # Base R statistics
base::library(package = "knitr")           # Build this presentation
base::library(package = "dplyr")           # Handle big tables
base::library(package = "Seurat")          # Handle SingleCell analyses
base::library(package = "SeuratObject")    # Handle SingleCell objects
base::library(package = "SingleCellExperiment") # Handle SingleCell file formats
base::library(package = "UpSetR")          # Nice venn-like graphs
base::library(package = "EnhancedVolcano") # Draw Volcano plot

First, we load Seurat object:

sobj <- base::readRDS(
  # Path to the RDS file
  file = "/shared/projects/ebaii_sc_teachers/SC_TD/06_Integration/RESULTS/12_TD3A.TDCT_S5_Integrated_12926.3886.RDS"
)

Then we join layers:

Seurat::Idents(sobj) <- sobj$orig.ident
sobj <- SeuratObject::JoinLayers(sobj)

Then we perform differential expression analysis:

sobj_de <- Seurat::FindMarkers(
  # Object on which to perform DEA
  object = sobj,
  # Name of factor in condition 1
  ident.1 = "TD3A",
  # Name of factor in condition 2
  ident.2 = "TDCT"
)

And that’s all. Our goal is to understand these lines, being able to write them is a bonus.

1.2 Purpose of this session

Up to now, we have:

1. Identified to which cell each sequenced reads come from
2. Identified to which gene each read come from
3. Identified possible bias in gene expression for each cell
4. Filtered and corrected these bias as well as we can

We would like to identify the list of genes that characterize differences between cell cycle phases G1 and S groups.

At the end of this session you will know:

1. how to select a differential analysis method
2. how to select the correct normalization (if any?) that must be provided to your differential analysis method
3. How to read differential expression results

1.3 Load RDS dataset

You already have a dataset loaded ? Good. Keep on going with it ! You don’t have one ? Use mine:

# The function `readRDS` from package `base`.
sobj <- base::readRDS(
  # Path to the RDS file
  file = "/shared/projects/ebaii_sc_teachers/SC_TD/06_Integration/RESULTS/12_TD3A.TDCT_S5_Integrated_12926.3886.RDS"
)

The command above, uses the function readRDS from base R package. and git it the path to the RDS file we built in previous session.

1.4 Insight

We are wondering what’s in our dataset. Let’s have a look, with the function print form the package base.

base::print(
  # Object to display
  x = sobj
)

An object of class Seurat 
12926 features across 3886 samples within 1 assay 
Active assay: RNA (12926 features, 2000 variable features)
 5 layers present: counts.TD3A, counts.TDCT, data.TD3A, data.TDCT, scale.data
 6 dimensional reductions calculated: pca, CCAIntegration, umap, RPCAIntegration, HarmonyIntegration, HarmonyStandalone

We have 12508 features (aka genes), across 2018 samples (aka cells) in the TD3A condition. We have 12254 features (aka genes), across 1868 samples (aka cells) in the TD3A condition.

Let us have a look at the RNA counts for 10 cells and their annotation, with the function head from the package utils.

utils::head(
  # Object to visualize
  x = sobj,
  # Number of lines to display
  n = 10
)

There is no counts, normalized or not. Where are they ?

In order to explore the content of sobj, use the function str from the package utils:

utils::str(
  # Object to explore
  object = sobj@assays
)

Alternatively, in RStudio, you can click on the object pane and explore manually the content of the object. If we explore the slot assays, then we find the counts.

You can access them with:

utils::head(
  # Object to visualize
  x = SeuratObject::GetAssayData(object = sobj, assay = "RNA", layer = "data.TD3A"),
  # Number of rows to display
  n = 10
)

We have one gene per line, one cell per column, and RNA counts in each row.

For the sake of this session, we won’t compare the whole dataset. It would take up to 15 minutes to complete. During the rest of this session, we will compare clusters 8 and 10, as annotated with Harmony.

We need to re-annotate layers to do so. This is done in two steps:

1. redefine Idents in order to be able to call cells by their cluster names.
2. JoinLayers to have a single count table with both TD3A and TDCT cells from clusters 8 and 10 together.

Seurat::Idents(sobj) <- sobj$HarmonyStandalone_clusters
sobj <- SeuratObject::JoinLayers(sobj)

Let’s check if our cells are joint:

base::print(sobj)

An object of class Seurat 
12926 features across 3886 samples within 1 assay 
Active assay: RNA (12926 features, 2000 variable features)
 3 layers present: data, counts, scale.data
 6 dimensional reductions calculated: pca, CCAIntegration, umap, RPCAIntegration, HarmonyIntegration, HarmonyStandalone

We have one gene per line, one cell per column, and RNA counts in each row. Are these counts normalized ? Are they scaled ? Are they filtered ? Are they corrected ?

Answer

These counts are normalized, scaled, filtered. This information is available in the seurat object itself, within the slot commands. See an example below:

names(sobj@commands)

 [1] "NormalizeData.RNA"                   
 [2] "FindVariableFeatures.RNA"            
 [3] "ScaleData.RNA"                       
 [4] "RunPCA.RNA"                          
 [5] "RunUMAP.RNA.CCAIntegration"          
 [6] "FindNeighbors.RNA.CCAIntegration"    
 [7] "RunUMAP.RNA.RPCAIntegration"         
 [8] "FindNeighbors.RNA.RPCAIntegration"   
 [9] "RunUMAP.RNA.HarmonyIntegration"      
[10] "FindNeighbors.RNA.HarmonyIntegration"
[11] "RunUMAP.RNA.HarmonyStandalone"       
[12] "FindNeighbors.RNA.HarmonyStandalone" 
[13] "FindClusters"                        

However, please be aware that counts in the slot count are raw counts. Normalized counts are in the slot data and scaled data are in the slot scaled.data. And it you do not find that clear, I totally agree with you.

Is it normal that we have so many zeros ? And what about surch low counts, is it because we downsampled the sequenced reads for this session ?

Answer

The large number of null counts is completely normal. In maths/stats we talk about matrix sparcity, i.e a table with lots of zeros. If the data were to be downsampled, we had had done this by cropping over a small chromosome, and not reducing the general amount of reads.

2 Select a DE method

2.1 Available methods

Seurat let us use multiple differential analysis methods with its function FindMarkers.

1. wilcox: The wilcoxon test tests the mean of expression and looks for a difference in these means.
2. MAST: This tool has been built for Single Cell. It is based on a statistical model called “Hurdle Model”, which excells with data that contains lots of zeros (which is our case in Single Cell RNA-Seq: most of the genes are not expressed.)
3. DESeq2: This tool has originally been built for bulk RNA-Seq but now includes specific functions for Single Cell. It performs well when counts are highly variable or when you wand to compare a handful of cells.
4. t-test: The t-test uses a comparison of means to assess differential expression probability.
5. roc: An AUC value of 1 means that expression values for this gene alone can perfectly classify the two groupings (i.e. Each of the cells in cells.1 exhibit a higher level than each of the cells in cells.2). An AUC value of 0 also means there is perfect classification, but in the other direction. A value of 0.5 implies that the gene has no predictive power to classify the two groups.

The main question now is how to choose the right test: spoilers, there are no option better than another in all ways.

From Soneson & Robinson (2018) Nature Methods:

de_tools

Here, researchers have designed an artificial dataset where they knew in advance the list of differentially expressed genes. They have used all these algorithms and consigned the results.

1. DESeq2, Limma seems to have a higher number of false positives (genes called differentially expressed while they were not.)
2. Wilcoxon seems to be better in general
3. Mast performs well in absolute

ANOVA was not present in this study.

2.2 Case: Plac8

The question is now to guess whether this gene is differnetially expressed or not.

2.2.1 Cell observations

Let’s have a look at the gene named ‘Plac8’, involved in positive regulation of cold-induced thermogenesis and positive regulation of transcription by RNA polymerase II. In order to plot its expression across all cells, we are going to use the function VlnPlot from Seurat package. The input object is obviously contained in the sobj variable, since it is our only Seurat object. In addition, we are going to select the feature Plac8, and split the graph according to the clusters we annotated earlier.

Seurat::VlnPlot(
  # A subset of the Seurat object
  # limited to clusters 8 and 10, 
  # or else we will plot all the clusters
  object = subset(sobj, HarmonyStandalone_clusters %in% c("8", "10")),
  # The name of the gene of interest (feature = gene)
  features = "Plac8",
  # The name of the Seurat cell annotation
  split.by = "HarmonyStandalone_clusters",
  # Change color for presentation
  cols = c("darkslategray3", "olivedrab")
)

Using your ‘intuition’, is this gene differentially expressed between cluster 8 and cluster 10 ?

Answer

In cluster 10, the violin plot highlights almost no cells with low or zero Plac8 expression. The highest density of cells has a Plac8 normalized expression aroung 1.5.

In Cluster 8, cells seems to have no expression of Plac8 at all.

IMHO, and you can disagree, the expression of the gene Plac8 differs between cluster 8 and cluster 10. This is purely intuitive.

Using your ‘intuition’, is this gene differentially expressed between G1 and S phases ?

Answer

Seurat::VlnPlot(
  # The Seurat object
  object = sobj,
  # The name of the gene of interest (feature = gene)
  features = "Plac8",
  # The name of the Seurat cell-cycle annotation
  group.by = "CC_Seurat_Phase",
  # Change color for presentation
  cols = c("darkslategray3", "olivedrab", "orangered3")
)

Most of the expression is null, some cells express the gene.

Okay, let’s have some informations about these distributions.

# Store counts in a variable
counts <- base::as.data.frame(
  # The matrix to reformat into a dataframe
  x = SeuratObject::GetAssayData(object = sobj, assay = "RNA", layer = "data")
)
# Rename cells with cell harmony cluster
base::colnames(counts) <- paste(
  # The names of the cell cluster for each cell
  sobj$CC_Seurat_Phase,
  # The names of the cells themselves
  colnames(sobj),
  sep = "_"
)

We have 2838 cells within the G1 group:

countsG1 <- select(counts, matches("^G1."))

Plac8G1 <- base::as.data.frame(base::t(countsG1["Plac8", ]))

base::summary(Plac8G1)

     Plac8        
 Min.   :0.00000  
 1st Qu.:0.00000  
 Median :0.00000  
 Mean   :0.04823  
 3rd Qu.:0.00000  
 Max.   :3.48145  

We have 711 cells withing the S group:

countsS <- select(counts, matches("^S."))

Plac8S <- base::as.data.frame(base::t(countsS["Plac8", ]))

base::summary(Plac8S)

     Plac8       
 Min.   :0.0000  
 1st Qu.:0.0000  
 Median :0.0000  
 Mean   :0.2097  
 3rd Qu.:0.0000  
 Max.   :3.5292  

2.2.2 From biology to statistics

Okay, let us resort on statistics to evaluate our chances to be guess correctly, or our risks to guess wrong.

We have lots of observations: 2838 cells within the G1 phase, and 711 cells withing the S phase.Statisticians really like to have a lot of observations! Ideally, statisticians always want to have more observation than tests. We have a total of 3549 observations and we are testing 1 gene. For them, this is a dream come true!

Are our cells supposed to be interacting with each others ? Are they independent from each others ? This is very important, and usually, it requires a discussion.

oes the expression in our cells follow a normal distribution ? It’s easy to check. Let’s draw the expression like we did above.

First, we use the function rbind from base package. Be careful, the function rbind also exists in DelayedArray, data.table, and BiocGenerics packages. We want to use the basic one. Here is an example of odd behaviors that occurs when not using the package name before a function call.

# Add column idientifiers in each count dataframes
Plac8G1$phase <- "G1"
Plac8S$phase <- "S"
# Paste the rows beneith each other
Plac8 <- base::rbind(
  #  variable pointing to G1 counts
  Plac8G1,
  #  variable pointing to S counts
  Plac8S,
  # A Boolean, requesting that strings/characters
  # should not be casted as `logical`. It breaks graphs.
  stringsAsFactors = FALSE
)

Secondly, we use the function gghistogram from package ggpubr in order to display relative abundance of gene expression:

ggpubr::gghistogram(
  Plac8,
  x = "Plac8",
  y = '..density..',
  fill = "steelblue",
  bins = 15,
  add_density = TRUE
)

Our distribution doesn’t seems to be normal, nor binomial we will have to rely on non-parametric tests.

Let’s run a non-parametric test base on the mean of distributions, since it’s the clothest to our ‘intuitive’ approach. Let there be Wilcoxon test.

In R, it’s quite straightforward: we have the function wilcoxon_test to perform the test, then we can plot the result.

# On the expression table stored in the varialbe `Plac8`,
# first apply the function `wilcox_test` from package `rstatix`,
# then we apply the function `add_significance` from package `rstatix`
stat.test <- Plac8 %>% rstatix::wilcox_test(Plac8 ~ phase) %>% rstatix::add_significance()
# While doing so, we usually also compute effect size
# eff.size <- Plac8 %>% rstatix::wilcox_effsize(Plac8 ~ phase)

Wilcoxon test says: the distributions are different, with a 3^{-13} % of risks of being wrong. The gene Plac8 can safely be said differentially expressed. We can compute a fold change and conclude.

Just out of curiosity, be aware that t_test from rstatix package, gives the same answer. However, DESeq gives a p-value of 0.05 and an adjusted p-value equal to 1.

Depending on the test used, this gene is, or is not differentially expressed.

2.3 Conclusion

In conclusion, to select your method, use the following:

1. If you have already done a study with one of these methods, keep using the same. This is crutial if you ever want to compare your new result with the old ones.
2. If you want to compare your study with a published one, use the same method.
3. If you have no idea, use Wilcoxon.
4. If you have bulk data analyzed with DESeq2/Limma, use DESeq2/Limma. It will be easier to take both works in consideration.

Please, never use a simple Wilcoxon on bulk RNA-Seq data.

3 Select a dataset

3.1 Dataset depends on selected method

There it is quite easier:

3.2 FindMarkers

With the function FindMarkers from package Seurat, we want to make three groups:

1. One using wilcoxon to perform DEA between clusters “8” and “10”.
2. One using t-test to perform DEA between clusters “8” and “10”.
3. One using ROC to perform DEA between “8” and “10”.

We will observe the results and compare our gene lists.

Hey, why are you looking at me? It’s your turn to work! Use the all the notions seen above to select the right counts (slot), the right input object, and the right arguments.

10 minutes practice !

Answers

Here are the code for each team:

sobj_wilcox <- Seurat::FindMarkers(
  # The variable that contains Seurat Object
  object = sobj,
  # Name of condition 1
  ident.1 = "8",
  # Name of condition 2
  ident.2 = "10",
  # Factor name in the Seurat Object
  group.by = "HarmonyStandalone_clusters",
  # Differential analysis method
  test.use = "wilcox"
)

sobj_t <- Seurat::FindMarkers(
  # The variable that contains Seurat Object
  object = sobj,
  # Name of condition 1
  ident.1 = "8",
  # Name of condition 2
  ident.2 = "10",
  # Factor name in the Seurat Object
  group.by = "HarmonyStandalone_clusters",
  # Differential analysis method
  test.use = "t"
)

sobj_roc <- Seurat::FindMarkers(
  # The variable that contains Seurat Object
  object = sobj,
  # Name of condition 1
  ident.1 = "8",
  # Name of condition 2
  ident.2 = "10",
  # Factor name in the Seurat Object
  group.by = "HarmonyStandalone_clusters",
  # Differential analysis method
  test.use = "roc"
)

In the function argument, there is a FoldChange threshold. Should we filter gene expression based on FoldChange? In case of positive answer, how much should that threshold be?

Answer

About thresholds on FDR (False Discovery Rate) and Log2(FC) (Log of the Fold Change), there are many discussions.

Remember, here the threshold on Fold Change is Logged. A log2(1) =0. And keep in mind the following:

1. If one selects a fold change threshold above 1.7, then their study concludes that smoking is not related to lung cancer.
2. If one selects a fold change threshold above 1, then their study concludes that a fast-food based diet does not lead to weight gain.
3. If one selects a fold change threshold above 1/25 000 000, then their study concludes: it is a complete hazard that mice have fetal malformation when in presence of Bisphanol.

In conclusion, there one, and only one reason to filter on fold change: in my experience, a fold change below 0.7 will be hard to see/verify on wet-lab (qRT).

If you need to reduce a too large number of differentially expressed genes, then reduce the FDR to 0.01, or even better, to 0.001. With that, you reduce your number of false claims.

Can you help me with DEseq2?

When I run the following command line, I have an error :

sobj_deseq2 <- Seurat::FindMarkers(
  # The variable that contains Seurat Object
  object = sobj,
  # Name of condition 1
  ident.1 = 8,
  # Name of condition 2
  ident.2 = 10,
  # Factor name in the Seurat Object
  group.by = "HarmonyStandalone_clusters",
  # Differential analysis method
  test.use = "deseq2",
  # Use non-normalized data with DESeq2
  slot = "counts"
)

Error in PerformDE(object = object, cells.1 = cells.1, cells.2 = cells.2, : Unknown test: deseq2

Answer

Oh! My fault, it was a typo in my command! Thank you all for your help!

sobj_deseq2 <- Seurat::FindMarkers(
  # The variable that contains Seurat Object
  object = sobj,
  # Name of condition 1
  ident.1 = "8",
  # Name of condition 2
  ident.2 = "10",
  # Factor name in the Seurat Object
  group.by = "HarmonyStandalone_clusters",
  # Differential analysis method
  test.use = "DESeq2",
  # Use non-normalized data with DESeq2
  slot = "counts"
)

Remark: by doing surch modification, some fold changes have been modified: remember the gene Atad2 with a mean expression of 0.08 in G1, and 0.2 in S phases? Mean expressions are now around 1.08 for G1, and 1.2 for S phases. This may be the reason why it was not differentially expressed in DESeq2, while Wilcoxon and T-test returned adjusted pvalues far below 0.05.

4 Explore results

4.1 Big tables

Let us have a look at the results:

We have the following columns:

1. p_val: Ignore this column. Always ignore raw p-values. Look at corrected ones, and if they are missing, then compute them.
2. avg_log2FC: Average Log2(FoldChange). Illustrates how much a gene is differentially expressed between samples in each condition.
3. pct.1: Percent of cells with gene expression in condition one, here in cluster 8.
4. pct.2: Percent of cells with gene expression in condition two, here in cluster 10.
5. p_val_adj: The confidence we have in the result. The closer to 0, the lesser is the risk of error.

We have the following columns:

1. p_val: Ignore this column. Always ignore raw p-values. Look at corrected ones, and if they are missing, then compute them.
2. avg_log2FC: Average Log2(FoldChange). Illustrates how much a gene is differentially expressed between samples in each condition.
3. pct.1: Percent of cells with gene expression in condition one, here in cluster 8.
4. pct.2: Percent of cells with gene expression in condition two, here in cluster 10.
5. p_val_adj: The confidence we have in the result. The closer to 0, the lesser is the risk of error.

We have the following columns:

1. myAUC: The area under the curve
2. avg_diff: Average difference in AUC
3. power: abs(AUC-0.5) * 2, useful to sort genes based on AUC
4. pct.1: Percent of cells with gene expression in condition one, here in cluster 8.
5. pct.2: Percent of cells with gene expression in condition two, here in cluster 10.

We have the following columns:

1. p_val: Ignore this column. Always ignore raw p-values. Look at corrected ones, and if they are missing, then compute them.
2. avg_log2FC: Average Log2(FoldChange). Illustrates how much a gene is differentially expessed between samples in each condition.
3. pct.1: Percent of cells with gene expression in condition one, here in cluster 8.
4. pct.2: Percent of cells with gene expression in condition two, here in cluster 10.
5. p_val_adj: The confidence we have in the result. The closer to 0, the lesser is the risk of error.

4.2 Filter DEA results

What kind of threshold should be used to filter each results?

If we must label certain scores as good or bad, we can reference the following rule of thumb:

0.5 = No discrimination 0.5-0.7 = Poor discrimination 0.7-0.8 = Acceptable discrimination 0.8-0.9= Excellent discrimination 0.9 = Outstanding discrimination

Hosmer and Lemeshow in Applied Logistic Regression (p. 177)

4.3 Add results to Seurat objects

We’d like to store the results of differential expression analysis in the Seurat object.

sobj@misc$wilcox <- sobj_wilcox

4.4 Common results

Now we can plot intersections in an up-set graph. It is like a Venn diagram:

UpSetR::upset(
  data = UpSetR::fromList(data),
  order.by = "freq"
)

4.5 Heatmap

We’d like to display the expression of genes identified by FindMarkers. Then we use the function DoHeatmap from the package Seurat.

In order to limit the graph to differentially expressed reads, we use the function rownames from R base package on the DEA result table. In this example, I use the results of wilcoxon, but you shall use any of the results you previously obtained.

Seurat::DoHeatmap(
  # variable pointing to seurat object
  object = sobj,
  # name of DE genes
  features = base::rownames(sobj_wilcox),
  # Cluster annotation
  group.by = "HarmonyStandalone_clusters",
)

4.6 Volcano plot

A Volcano plot is usefull to identify differnetial expression analysis bias.

The package EnhancedVolcano has an eponym function for that:

EnhancedVolcano::EnhancedVolcano(
  #  variable pointing to the DEA results
  toptable = sobj_wilcox,
  # Gene names
  lab = rownames(sobj_wilcox),
  # Column in which to find Fold Change
  x = "avg_log2FC",
  # Column in which to find confidence interval
  y = "p_val_adj",
  # Lower fold-change cut-off
  FCcutoff = 0.2
)

4.7 Session Info

This list of all packages used while you work should be included in each en every R presentation:

utils::sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: x86_64-conda-linux-gnu
Running under: Ubuntu 20.04.6 LTS

Matrix products: default
BLAS/LAPACK: /shared/ifbstor1/software/miniconda/envs/r-4.4.1/lib/libopenblasp-r0.3.27.so;  LAPACK version 3.12.0

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

time zone: Europe/Paris
tzcode source: system (glibc)

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

other attached packages:
 [1] EnhancedVolcano_1.22.0      ggrepel_0.9.6              
 [3] UpSetR_1.4.0                SingleCellExperiment_1.26.0
 [5] SummarizedExperiment_1.34.0 Biobase_2.64.0             
 [7] GenomicRanges_1.56.2        GenomeInfoDb_1.40.1        
 [9] IRanges_2.38.1              S4Vectors_0.42.1           
[11] BiocGenerics_0.50.0         MatrixGenerics_1.16.0      
[13] matrixStats_1.4.1           Seurat_5.1.0               
[15] SeuratObject_5.0.2          sp_2.1-4                   
[17] dplyr_1.1.4                 knitr_1.48                 
[19] rstatix_0.7.2               ggpubr_0.6.0               
[21] ggplot2_3.5.1               DT_0.33                    
[23] BiocParallel_1.38.0        

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22        splines_4.4.1           later_1.3.2            
  [4] tibble_3.2.1            polyclip_1.10-7         rpart_4.1.23           
  [7] fastDummies_1.7.4       lifecycle_1.0.4         globals_0.16.3         
 [10] lattice_0.22-6          MASS_7.3-61             crosstalk_1.2.1        
 [13] backports_1.5.0         magrittr_2.0.3          limma_3.60.6           
 [16] Hmisc_5.1-3             plotly_4.10.4           sass_0.4.9             
 [19] rmarkdown_2.28          jquerylib_0.1.4         yaml_2.3.10            
 [22] httpuv_1.6.15           sctransform_0.4.1       spam_2.11-0            
 [25] spatstat.sparse_3.1-0   reticulate_1.39.0       cowplot_1.1.3          
 [28] pbapply_1.7-2           RColorBrewer_1.1-3      abind_1.4-8            
 [31] zlibbioc_1.50.0         Rtsne_0.17              purrr_1.0.2            
 [34] nnet_7.3-19             GenomeInfoDbData_1.2.12 irlba_2.3.5.1          
 [37] listenv_0.9.1           spatstat.utils_3.1-0    goftest_1.2-3          
 [40] RSpectra_0.16-2         spatstat.random_3.3-2   fitdistrplus_1.2-1     
 [43] parallelly_1.38.0       DelayedArray_0.30.1     leiden_0.4.3.1         
 [46] codetools_0.2-20        tidyselect_1.2.1        UCSC.utils_1.0.0       
 [49] farver_2.1.2            base64enc_0.1-3         spatstat.explore_3.3-2 
 [52] jsonlite_1.8.9          progressr_0.14.0        Formula_1.2-5          
 [55] ggridges_0.5.6          survival_3.7-0          tools_4.4.1            
 [58] ica_1.0-3               Rcpp_1.0.13             glue_1.8.0             
 [61] SparseArray_1.4.8       gridExtra_2.3           DESeq2_1.44.0          
 [64] xfun_0.48               withr_3.0.1             fastmap_1.2.0          
 [67] fansi_1.0.6             digest_0.6.37           R6_2.5.1               
 [70] mime_0.12               colorspace_2.1-1        scattermore_1.2        
 [73] tensor_1.5              spatstat.data_3.1-2     utf8_1.2.4             
 [76] tidyr_1.3.1             generics_0.1.3          data.table_1.16.2      
 [79] S4Arrays_1.4.1          httr_1.4.7              htmlwidgets_1.6.4      
 [82] uwot_0.2.2              pkgconfig_2.0.3         gtable_0.3.5           
 [85] lmtest_0.9-40           XVector_0.44.0          htmltools_0.5.8.1      
 [88] carData_3.0-5           dotCall64_1.2           scales_1.3.0           
 [91] png_0.1-8               spatstat.univar_3.0-1   rstudioapi_0.17.0      
 [94] reshape2_1.4.4          checkmate_2.3.1         nlme_3.1-165           
 [97] cachem_1.1.0            zoo_1.8-12              stringr_1.5.1          
[100] KernSmooth_2.23-24      vipor_0.4.7             parallel_4.4.1         
[103] miniUI_0.1.1.1          foreign_0.8-86          ggrastr_1.0.2          
[106] pillar_1.9.0            grid_4.4.1              vctrs_0.6.5            
[109] RANN_2.6.2              promises_1.3.0          car_3.1-3              
[112] xtable_1.8-4            cluster_2.1.6           beeswarm_0.4.0         
[115] htmlTable_2.4.2         evaluate_1.0.1          locfit_1.5-9.9         
[118] cli_3.6.3               compiler_4.4.1          crayon_1.5.3           
[121] rlang_1.1.4             future.apply_1.11.2     ggsignif_0.6.4         
[124] labeling_0.4.3          ggbeeswarm_0.7.2        plyr_1.8.9             
[127] stringi_1.8.4           viridisLite_0.4.2       deldir_2.0-4           
[130] munsell_0.5.1           lazyeval_0.2.2          spatstat.geom_3.3-3    
[133] Matrix_1.7-1            RcppHNSW_0.6.0          patchwork_1.3.0        
[136] future_1.34.0           statmod_1.5.0           shiny_1.9.1            
[139] highr_0.11              ROCR_1.0-11             igraph_2.1.1           
[142] broom_1.0.7             bslib_0.8.0