CD8+ T cell Dual scRNA/TCR-seq analysis

scRNA/TCR-seq workflow for Maxwell et al Figures

This vignette reproduces figure panels from from Maxwell et al, 2024 Figure 2/S2 on CD8+ T cell scRNA and TCR sequencing using T cells sorted from wildtype (sgScramble) or ARID1A deficient (sgArid1a) B16F10 tumors.

Load packages.

#load packages with package manager package called pacman
if (!require("pacman")) install.packages("pacman")
pacman::p_load(here,  
               Seurat, #Primary package for scRNA analysis
               scDataviz, #UMAP density plot
               cowplot, #arranging density plots on same plane
               ggplot2, #plotting
               dplyr, #data wrangling
               ggpubfigs, #colorblind friendly color palettes 
               scRepertoire, #scTCR analysis
               SingleCellExperiment, #Make sce objects
               scater, #single cell toolkit
               Nebulosa, #Marker gene expression density
               stringr,#Reformatting strings
               forcats, #factors
               WebGestaltR) #GSEA

knitr::opts_chunk$set(echo = TRUE) #default setting is that code chunks output are visible 
set.seed(123)  # Setting seed for reproducibility

CD8+ T cell UMAP (Figure 2D)

Import data and generate UMAP using Seurat DimPlot function with color blind friendly palette from the ggpubfigs package.

#import data
dat.integrated <- readRDS("scRNA_CD8_B16F10.rds")

#Seurat Dimplot function for UMAP
DimPlot(dat.integrated, reduction = "umap", cols = friendly_pal("nickel_five"), label = TRUE, label.box = TRUE) + 
  ggtitle("CD8+ T cell clusters") + theme(plot.title = element_text(hjust = 0.5))

Cluster Marker Gene Heatmap (Figure S2B)

Generating marker gene heatmap of top 10 marker genes from FindMarkers Seurat function for each cluster. The heatmap will look blurry when viewed in Rstudio, but has much higher resolution when viewed in a PDF viewer such as Adobe or Drawboard PDF.

marker_genes <- read.delim("scRNA_marker_genes.txt", header = F)
marker_genes <- marker_genes[, 1]
  
# Generating the heatmap
heatmap_plot <- DoHeatmap(object = dat.integrated, features = marker_genes, group.colors = friendly_pal("nickel_five")) +
  scale_fill_viridis_c(limits = c(-0.5, 1.5), oob = scales::squish) +
  scale_y_discrete(position = "right")

# Save the heatmap with specified dimensions
ggsave("Marker_genes_heatmap.pdf", plot = heatmap_plot, width = 10, height = 14, device = "pdf")
#View heatmap in PDF viewer

Marker Gene Expression Density (Figure S2C)

To visualize cluster marker gene expression, we use Neubulosa package for gene weighted density visualization. This visualization is more sensitive to low expression than Seurat’s FeaturePlot function and the plots are beuatiful. Below are Nebulosa density plots for C5 proliferation cluster genes Mki67 and Top2a.

p5 <- plot_density(dat.integrated, c("Mki67", "Top2a"))
print(p5)

Density Plot of CD8+ T cell clusters (Figure 2E)

Visualize CD8+ T cell states from sgScramble to sgArid1a tumors using scDataViz package.

#subset Seurat object by treatment for density plot
sgScr.int=subset(x=dat.integrated,subset=Treatment=="sgScramble")
sgArid1a.int=subset(x=dat.integrated,subset=Treatment=="sgArid1a")

#save sce object for sgScr
sgScr_sce <- as.SingleCellExperiment(sgScr.int) 
#colData function from Scater single cell package
metadata(sgScr_sce) <- data.frame(colData(sgScr_sce))

#save sce object for sgArid1a
sgArid1a_sce <- as.SingleCellExperiment(sgArid1a.int) 
#colData function from Scater single cell package
metadata(sgArid1a_sce) <- data.frame(colData(sgArid1a_sce))

#save plot for sgScr
#contourPlot() function from scDataviz package
sgScramble <- contourPlot(sgScr_sce, contour=NA, reducedDim = 'UMAP',dimColnames = c('UMAP_1','UMAP_2'), title = "sgScramble") + 
  scale_fill_viridis_c(option="magma") +
  theme_classic() #creates cell density plot split by treatment

#save plot for sgArid1a
#contourPlot() function from scDataviz package
sgArid1a <- contourPlot(sgArid1a_sce, contour=NA, reducedDim = 'UMAP',dimColnames = c('UMAP_1','UMAP_2'), title = "sgArid1a") + 
  scale_fill_viridis_c(option="magma") + 
   theme_classic() #creates cell density plot split by treatment

#Run plot_grid() function from cowplot package
#Plot density UMAPs on same plane
plot_grid(sgScramble,sgArid1a)

CD8+ T cell cluster proportion quantification (Figure 2F)

Visualize cluster populations among all CD8+ T cells within each group. Data wrangling with dplyr and plotting with ggplot2.

#Extract cluster, genotype, and cell number from Seurat object
clusts=table(dat.integrated@active.ident,dat.integrated@meta.data$Treatment)
clusts=as.data.frame(clusts)

#Reshape dataframe and create proportion column using dplyr package
clusts <- clusts %>%
  rename(Cluster = Var1, Genotype = Var2) %>%
  group_by(Genotype) %>%
  mutate(prop = Freq/sum(Freq))

#Make Genotype a factor and set desired order for plotting
clusts$Genotype <- as.factor(clusts$Genotype)
clusts$Genotype <- factor(clusts$Genotype, levels = c("sgScramble", "sgArid1a"))

#color blind friendly palette
mypalette <- friendly_pal("nickel_five")

#barplot
ggplot(clusts,aes(x=Genotype,y=prop,fill=Cluster))+
  geom_bar(stat="identity", color="black")+
  theme_classic()+ scale_fill_manual(values=mypalette) + 
  ggtitle(substitute(paste(bold("Cluster Density")))) + 
  theme(plot.title = element_text(size = 18, hjust = 0.5),
        axis.text.x = element_text(size = 14),
        axis.title = element_blank(),
        axis.title.y = element_text(size = 14)) +
  labs(y = "Proportion of UMAP")

Anti-PD1 Response Gene Signature (Figure 2G)

Visualize expression of a previously defined CD8+ T cell specific anti-PD1 response signature Kumar et al in T cells from both sgScramble and sgArid1a tumors.

#Read in gene signatures from Maxwell et al
gene_signatures <- read.delim("Maxwell_et_al_Gene_signatures.txt", header = TRUE)

#Select gene signature from Kumar et al
anti_PD1_response <- gene_signatures %>%
  select("CD8_anti_PD1_response_sig")

#Convert from dataframe to vector
anti_PD1_response <- as.vector(anti_PD1_response)

#AddModuleScore function from Seurat to add a custom gene signature to object
dat.integrated <- AddModuleScore(object = dat.integrated, features = anti_PD1_response, ctrl = 5, name = "Anti_PD1_Response", seed = 12)

#Visualize expression of custom gene signature 
FeaturePlot(dat.integrated,features = "Anti_PD1_Response1",label=FALSE,split.by = "Treatment",order=TRUE, pt.size = .75) & theme(legend.position="right",text=element_text(size=10)) & scale_color_viridis_c(option="magma", limits=c(0.5, 1.1),oob=scales::squish)

Cluster 2 GSEA Enrichments (Figure 2H)

Perform differential gene expression (DGE) followed by GSEA comparing sgScramble to sgArid1a for all five clusters. DGE is performed using Seurat FindMarkers function and GSEA is performed using WebGestaltR package.

#Add a "cluster_treatment" vector to object metadata (e.g., Naive_sgArid1a or Naive_sgScramble)
dat.integrated@meta.data$cluster_treatment <- paste(as.vector(dat.integrated$cell.type), as.vector(dat.integrated$Treatment), sep = "_")
table(dat.integrated$cluster_treatment)

## 
##        Effector_sgArid1a      Effector_sgScramble      Exhaustion_sgArid1a 
##                     8276                     3669                     3568 
##    Exhaustion_sgScramble          Innate_sgArid1a        Innate_sgScramble 
##                     1052                     2632                      679 
##           Naive_sgArid1a         Naive_sgScramble   Proliferating_sgArid1a 
##                     7220                     8943                     1740 
## Proliferating_sgScramble 
##                      386

#Make cluster_treatment object ident 
Idents(dat.integrated) <- dat.integrated$cluster_treatment
table(Idents(dat.integrated))

## 
##         Naive_sgScramble      Effector_sgScramble    Exhaustion_sgScramble 
##                     8943                     3669                     1052 
##        Innate_sgScramble Proliferating_sgScramble           Naive_sgArid1a 
##                      679                      386                     7220 
##        Effector_sgArid1a          Innate_sgArid1a      Exhaustion_sgArid1a 
##                     8276                     2632                     3568 
##   Proliferating_sgArid1a 
##                     1740

##########
#part1: DE analysis with Seurat FindMarkers function
# outputs: .csv files for DE analysis

#set seed for reproducible analysis
set.seed(12345)
#cluster names vector
name <- as.vector(unique(dat.integrated$cell.type))
#cluster_genotype vectors
sgScramble <- paste0(name, "_sgScramble")
sgArid1a <- paste0(name, "_sgArid1a")

#DGE analysis
for (i in 1:length(name)){
  outname <- paste0("response_", name[i], ".csv")
  temp <- FindMarkers(dat.integrated, ident.1 = sgArid1a[i], ident.2 = sgScramble[i], logfc.threshold = 0, min.pct = 0.05)
  write.csv(temp, outname, quote=F)
}

##########
#part2: generate GSEA input file - ranked gene list with scores
# outputs: /GSEA/*.rnk files for GSEA
##########
#gene rank

#DGE .csv files for each cluster
csv <- list.files(pattern="*.csv")

#Generate rank list for GSEA
for(i in 1:length(csv)){
  # Extract the base file name without the extension and add the new extension
  base_name <- gsub(pattern="\\.csv$", replacement="", x=csv[i])
  rank <- paste0("./GSEA/", base_name, ".rnk") # Construct the output path correctly
  
  # Read the CSV file
  deg <- read.csv(csv[i], header=TRUE)
  
  # Calculate the rank score
  rank.score <- -log10(deg$p_val_adj)*sign(deg$avg_log2FC)
  
  # Handle infinite and NA values
  rank.score[is.na(rank.score)] <- 0
  rank.score[rank.score == Inf] <- max(rank.score[which(rank.score < Inf)])+50
  rank.score[rank.score == -Inf] <- min(rank.score[which(rank.score > -Inf)])-50
  
  # Check for NA and infinite values before proceeding
  if(!any(is.na(rank.score)) && !any(is.infinite(rank.score))){
    deg <- cbind(deg$X, rank.score)
    write.table(deg, rank, quote=FALSE, row.names = FALSE, col.names = FALSE, sep="\t")
  }
}


fnms <- list.files(pattern="*.csv")

for(i in 1:length(fnms)){
  # Extract the base file name without the extension and add the new extension
  base_name <- gsub(pattern="\\.csv$", replacement="", x=fnms[i])
  onms <- paste0("./GSEA/", base_name, ".rnk") # Construct the output path correctly
  
  # Read the CSV file
  deg <- read.csv(fnms[i], header=TRUE)
  
  # Calculate the rank score
  rank.score <- -log10(deg$p_val_adj)*sign(deg$avg_log2FC)
  
  # Handle infinite and NA values
  rank.score[is.na(rank.score)] <- 0
  rank.score[rank.score == Inf] <- max(rank.score[which(rank.score < Inf)])+50
  rank.score[rank.score == -Inf] <- min(rank.score[which(rank.score > -Inf)])-50
  
  # Check for NA and infinite values before proceeding
  if(!any(is.na(rank.score)) && !any(is.infinite(rank.score))){
    deg <- cbind(deg$X, rank.score)
    write.table(deg, onms, quote=FALSE, row.names = FALSE, col.names = FALSE, sep="\t")
  }
}

#########
# part3: perform GSEA with WebGestaltR
#########

outdir <- "./GSEA/hallmark/"

runEnrich <- function(rankGene, outfnm, outdir) {
  fdr <- 1
  enrichTestGSEA <- WebGestaltR(enrichMethod = "GSEA",
                                organism="mmusculus", 
                                #enrichDatabaseFile="./GSEA/hallmark_mouse.gmt",
                                enrichDatabaseFile="./GSEA/hallmark_mouse.gmt", #change gmt file
                                enrichDatabaseType="genesymbol",
                                interestGene = rankGene,
                                interestGeneType = "genesymbol",
                                referenceSet="genome",
                                minNum=3,maxNum=2000,
                                perNum=10000,
                                fdrMethod="BH", sigMethod="top", 
                                topThr = 51, reportNum = 51,
                                isOutput=T,
                                outputDirectory=outdir,
                                projectName=outfnm,
                                saveRawGseaResult=T,
                                nThreads=20)
  
}

csv <- list.files("./GSEA/", pattern="*.rnk")
for (i in 1:length(csv)){
  tmn <- paste0("./GSEA/", csv[i])
  rankGene <- read.delim(tmn, header=F)
  outfnm <- strsplit(tmn, "response_|.rnk")[[1]][2]
  runEnrich(rankGene, outfnm, outdir)
} 

## Loading the functional categories...
## Loading the ID list...
## Summarizing the uploaded ID list by GO Slim data...

## Performing the enrichment analysis...
## 10000 permutations of score complete...

## Begin affinity propagation...
## End affinity propagation...
## Begin weighted set cover...
## End weighted set cover...
## Generate the final report...
## Results can be found in the ./GSEA/hallmark//Project_Effector!
## Loading the functional categories...
## Loading the ID list...
## Summarizing the uploaded ID list by GO Slim data...

## Performing the enrichment analysis...
## 10000 permutations of score complete...

## Begin affinity propagation...
## End affinity propagation...
## Begin weighted set cover...
## End weighted set cover...
## Generate the final report...
## Results can be found in the ./GSEA/hallmark//Project_Exhaustion!
## Loading the functional categories...
## Loading the ID list...
## Summarizing the uploaded ID list by GO Slim data...

## Performing the enrichment analysis...
## 10000 permutations of score complete...

## Begin affinity propagation...
## End affinity propagation...
## Begin weighted set cover...
## End weighted set cover...
## Generate the final report...
## Results can be found in the ./GSEA/hallmark//Project_Innate!
## Loading the functional categories...
## Loading the ID list...
## Summarizing the uploaded ID list by GO Slim data...

## Performing the enrichment analysis...
## 10000 permutations of score complete...

## Begin affinity propagation...
## End affinity propagation...
## Begin weighted set cover...
## End weighted set cover...
## Generate the final report...
## Results can be found in the ./GSEA/hallmark//Project_Naive!
## Loading the functional categories...
## Loading the ID list...
## Summarizing the uploaded ID list by GO Slim data...

## Performing the enrichment analysis...
## 10000 permutations of score complete...

## Begin affinity propagation...
## End affinity propagation...
## Begin weighted set cover...
## End weighted set cover...
## Generate the final report...
## Results can be found in the ./GSEA/hallmark//Project_Proliferating!

#end.

The CD8+ T cell immunotherapy response signature indicated the highest expression is in C2 effector like CD8+ T cells and this is increased in sgArid1a tumor group. To investigate what gene signatures are enriched in sgArid1a C2 T cells vs sgScrmable C2 T cells, we pulled data from GSEA results for C2 Effector like cluster and plot enriched gene sets with ggplot2.

rawGseaResult <- readRDS("./GSEA/hallmark/Project_Effector/Project_Effector_GSEA/rawGseaResult.rds")


Gene_set <- names(rawGseaResult$Items_in_Set)
#replace underscores with spaces
Gene_set <- gsub("_", " ", Gene_set) 
#Capitalize only first letter of each word, stringr package
Gene_set <- str_to_title(Gene_set)
#delete hallmark text, too reduant for later plotting
Gene_set <- gsub("Hallmark ", "", Gene_set)

NES <- rawGseaResult$Enrichment_Results$NES
FDR <- rawGseaResult$Enrichment_Results$fdr

#make gsea_df for plotting
C2_gsea_df <- data.frame(Gene_set, NES, FDR)

#Take top and bottom 4 gene sets, dplyr package
C2_top4 <- C2_gsea_df %>% 
  arrange(NES) %>%
  slice_head(n=4)
C2_bottom4 <- C2_gsea_df %>%
  arrange(NES) %>%
  slice_tail(n=4)

#join the top enrichments for sgArid1a and sgScr C2 T cells
C2_top_df <- rbind(C2_top4, C2_bottom4)

#GSEA vertical bars (log10pvalue scale, y=gene set, x=NES) #ggplot2 and forcats packages
ggplot(C2_top_df, aes(NES, fct_reorder(Gene_set, NES), fill=FDR)) + geom_col() + 
         scale_fill_gradientn(colors=c("#8b0000","#FEFEBE","#9FE2BF","slateblue"))+theme_classic()+theme(panel.border = element_rect(color="black",fill=NA,size=1.5))+ggtitle("GSEA")+labs(fill="FDR")+geom_vline(xintercept = c(0),size=1) +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_text(size = 14),
        plot.title = element_text(size = 18, face = "bold")) +
  labs(title = "C2 Effector-like \nGSEA", 
       x = "NES \n(sgArid1a/sgScramble)") 

CD8+ T cell ISG Metagene (Figure 2I)

Since IFN responses are prominently enriched in C2 effector like T cells, we created an interferon stimulated gene (ISG) metagene from leading edge genes in IFN Alpha/Gamma response signatures enriched in any of the five clusters from sgArid1a tumors to visualize this ISG signature in CD8+ T cells.

#Select CD8+ T cell metagene signature
CD8_ISG_score <- gene_signatures %>%
  select("CD8_ISG_Metagene")

#Convert from dataframe to vector
CD8_ISG_score <- as.vector(CD8_ISG_score)

#AddModuleScore function from Seurat to add a custom gene signature to object
dat.integrated <- AddModuleScore(object = dat.integrated, features = CD8_ISG_score, ctrl = 5, name = "ISG_Metagene", seed = 12)

#Visualize expression of custom gene signature 
FeaturePlot(dat.integrated,features = "ISG_Metagene1",label=FALSE,split.by = "Treatment",order=TRUE, pt.size = .75) & theme(legend.position="right",text=element_text(size=10)) & scale_color_viridis_c(option="magma", limits=c(0, 1), oob=scales::squish)

TCR clonotype UMAP (Figure 2J)

To further understand the T cell response in sgArid1a tumors, we investigated the clonal expansion of T cells via 10X VDJ kit for TCR-sequencing. This was primarily via the scRepertoire package.

#Read in VDJ filtered annotation files from 10X's CellRanger
sgArid1a_1=read.csv("./TCR_files/sgArid1a_1_filtered_contig_annotations.csv") 
sgArid1a_2=read.csv("./TCR_files/sgArid1a_2_filtered_contig_annotations.csv")
sgScr_1=read.csv("./TCR_files/sgScr_1_filtered_contig_annotations.csv")
sgScr_2=read.csv("./TCR_files/sgScr_2_filtered_contig_annotations.csv")

#Make list of TCR samples
contig_list=list(sgArid1a_1,sgArid1a_2,sgScr_1,sgScr_2)

#scRepertoire command to make object
combinedTCR=combineTCR(contig_list,samples=c("sgArid1a","sgArid1a","sgScr","sgScr"),ID=c("1","2","1","2"))

#Integrate TCR data with Seurat object
dat.integrated = combineExpression(combinedTCR,dat.integrated, cloneCall="gene", proportion = FALSE, cloneSize = c(Single=1,Small=5,Medium=20,Large=100,Hyperexpanded=500))

#Make new metadata feature for CloneType TCR expansion categories in Seurat object
slot(dat.integrated,"meta.data")$cloneType=factor(slot(dat.integrated,"meta.data")$cloneType, levels=c("Hyperexpanded (100 < X <= 500)","Large (20 < X <= 100)","Medium (5 < X <= 20)","Small (1 < X <= 5)","Single (0 < X <= 1)", "No TCR mapped"))

# Define custom color palette
custom_colors <- c("Hyperexpanded (100 < X <= 500)" = "#FF0000",
                   "Large (20 < X <= 100)" = "#FF7F00",
                   "Medium (5 < X <= 20)" = "#80A440",
                   "Small (1 < X <= 5)" = "#4A9BFF",
                   "Single (0 < X <= 1)" = "#AAEEFF",
                   "No TCR mapped" = "grey50" )


#Make sure any NA value is populated with "No TCR mapped"
dat.integrated@meta.data$cloneType[is.na(dat.integrated@meta.data$cloneType)] <- "No TCR mapped"


#DimPlot for CloneType 
DimPlot(dat.integrated, group.by = "cloneType", split.by = "Treatment", pt.size = 0.2,
        cols = custom_colors, order = c("Hyperexpanded (100 < X <= 500)",
                                                        "Large (20 < X <= 100)",
                                                        "Medium (5 < X <= 20)",
                                                        "Small (1 < X <= 5)",
                                                        "Single (0 < X <= 1)",
                                                        "No TCR mapped"))

TCR clonotype classes quantification (Figure 2J)

Visualizing the classes of clonotype expansion in each tumor genotype. Note: The TCR clonotype barchart groups are defined by a proportion of total TCR in their respective groups which are equal to the integer values seen on the TCR UMAP above.

sgArid1a_combine <- rbind(sgArid1a_1, sgArid1a_2)
sgScr_combine <- rbind(sgScr_1, sgScr_2)

#sgArid1a_comb=rbind(sgArid1a_1,sgArid1a_2)
contig_list=list(sgScr_combine, sgArid1a_combine)

#scRepetoire
combinedTCR=combineTCR(contig_list,samples=c("sgScr", "sgArid1a"))

mypalette2 <- rev(c("#FF0000", "#FF7F00", "#80A440", "#4A9BFF", "#AAEEFF"))

#scRepertoire function
clonalHomeostasis(combinedTCR, cloneSize = c(Single = 1e-04, Small = 0.0003, Medium = 0.0012, Large = 0.006, Hyperexpanded= 0.3),cloneCall = "gene",chain="both", )+scale_fill_manual(values=mypalette2) +theme(axis.title.x = element_blank())

TCR clonal scatterplot (Figure 2K)

Visualizing expansion of individual TCRs in each tumor genotype and whether they are expanded in both genotypes (dual expanded) or unique between genotypes. Colors for this figure were edited in illustrator. Interestingly, the top TCR in sgArid1a tumors is shared in sgScramble tumors, but at low frequency in sgScramble tumor group. This suggests increased antigen presentation of a common antigen such as a B16F10 tumor antigen (most likely), but could conceivably be reactive to other types of antigens such as an auto-antigen.

#clonotype scatter plot
clonalScatter(combinedTCR,cloneCall="gene",y.axis="sgArid1a",x.axis="sgScr",dot.size="total",graph="proportion", palette = "Berlin") #exportTable=TRUE exports data to table

Session Info

For reproducibility, my session info is listed below outlining the version of R and package versions that I’ve using for this vignette.

print(sessionInfo())

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.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: America/Los_Angeles
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] WebGestaltR_0.4.6           forcats_1.0.0              
##  [3] stringr_1.5.1               Nebulosa_1.12.1            
##  [5] patchwork_1.2.0             scater_1.31.2              
##  [7] scuttle_1.10.3              scRepertoire_2.0.0         
##  [9] ggpubfigs_0.0.1             dplyr_1.1.4                
## [11] ggplot2_3.5.0               cowplot_1.1.3              
## [13] scDataviz_1.12.0            SingleCellExperiment_1.22.0
## [15] SummarizedExperiment_1.30.2 Biobase_2.60.0             
## [17] GenomicRanges_1.52.1        GenomeInfoDb_1.36.4        
## [19] IRanges_2.34.1              MatrixGenerics_1.12.3      
## [21] matrixStats_1.2.0           S4Vectors_0.38.2           
## [23] BiocGenerics_0.46.0         Seurat_5.0.1               
## [25] SeuratObject_5.0.1          sp_2.1-3                   
## [27] here_1.0.1                  pacman_0.5.1               
## 
## loaded via a namespace (and not attached):
##   [1] spatstat.sparse_3.0-3     bitops_1.0-7             
##   [3] doParallel_1.0.17         httr_1.4.7               
##   [5] RColorBrewer_1.1-3        doRNG_1.8.6              
##   [7] tools_4.3.1               sctransform_0.4.1        
##   [9] utf8_1.2.4                R6_2.5.1                 
##  [11] lazyeval_0.2.2            uwot_0.1.16              
##  [13] withr_3.0.0               gridExtra_2.3            
##  [15] progressr_0.14.0          textshaping_0.3.7        
##  [17] quantreg_5.97             cli_3.6.2                
##  [19] Cairo_1.6-2               spatstat.explore_3.2-6   
##  [21] fastDummies_1.7.3         iNEXT_3.0.1              
##  [23] isoband_0.2.7             labeling_0.4.3           
##  [25] sass_0.4.8                mvtnorm_1.2-4            
##  [27] spatstat.data_3.0-4       readr_2.1.5              
##  [29] ggridges_0.5.6            pbapply_1.7-2            
##  [31] askpass_1.2.0             systemfonts_1.0.5        
##  [33] svglite_2.1.3             stringdist_0.9.12        
##  [35] parallelly_1.36.0         limma_3.56.2             
##  [37] flowCore_2.14.2           VGAM_1.1-10              
##  [39] rstudioapi_0.15.0         generics_0.1.3           
##  [41] vroom_1.6.5               ica_1.0-3                
##  [43] spatstat.random_3.2-2     Matrix_1.6-5             
##  [45] RProtoBufLib_2.14.1       ggbeeswarm_0.7.2         
##  [47] fansi_1.0.6               abind_1.4-5              
##  [49] lifecycle_1.0.4           whisker_0.4.1            
##  [51] yaml_2.3.8                Rtsne_0.17               
##  [53] grid_4.3.1                promises_1.2.1           
##  [55] crayon_1.5.2              miniUI_0.1.1.1           
##  [57] lattice_0.22-5            beachmat_2.16.0          
##  [59] pillar_1.9.0              knitr_1.45               
##  [61] rjson_0.2.21              future.apply_1.11.1      
##  [63] codetools_0.2-19          leiden_0.4.3.1           
##  [65] glue_1.7.0                data.table_1.15.0        
##  [67] vctrs_0.6.5               png_0.1-8                
##  [69] spam_2.10-0               gtable_0.3.4             
##  [71] cachem_1.0.8              ks_1.14.2                
##  [73] xfun_0.42                 S4Arrays_1.2.1           
##  [75] mime_0.12                 tidygraph_1.3.0          
##  [77] pracma_2.4.4              survival_3.5-7           
##  [79] iterators_1.0.14          cytolib_2.14.1           
##  [81] ellipsis_0.3.2            fitdistrplus_1.1-11      
##  [83] ROCR_1.0-11               nlme_3.1-164             
##  [85] bit64_4.0.5               RcppAnnoy_0.0.22         
##  [87] evd_2.3-6.1               rprojroot_2.0.4          
##  [89] bslib_0.6.1               irlba_2.3.5.1            
##  [91] vipor_0.4.7               KernSmooth_2.23-22       
##  [93] colorspace_2.1-0          ggrastr_1.0.2            
##  [95] tidyselect_1.2.0          bit_4.0.5                
##  [97] curl_5.2.0                compiler_4.3.1           
##  [99] BiocNeighbors_1.18.0      SparseM_1.81             
## [101] ggdendro_0.1.23           DelayedArray_0.26.7      
## [103] plotly_4.10.4             scales_1.3.0             
## [105] lmtest_0.9-40             apcluster_1.4.11         
## [107] digest_0.6.34             goftest_1.2-3            
## [109] presto_1.0.0              spatstat.utils_3.0-4     
## [111] rmarkdown_2.25            XVector_0.40.0           
## [113] htmltools_0.5.7           pkgconfig_2.0.3          
## [115] umap_0.2.10.0             sparseMatrixStats_1.12.2 
## [117] highr_0.10                fastmap_1.1.1            
## [119] rlang_1.1.3               htmlwidgets_1.6.4        
## [121] shiny_1.8.0               DelayedMatrixStats_1.22.6
## [123] farver_2.1.1              jquerylib_0.1.4          
## [125] zoo_1.8-12                jsonlite_1.8.8           
## [127] BiocParallel_1.34.2       mclust_6.1               
## [129] BiocSingular_1.16.0       RCurl_1.98-1.14          
## [131] magrittr_2.0.3            GenomeInfoDbData_1.2.10  
## [133] dotCall64_1.1-1           munsell_0.5.0            
## [135] Rcpp_1.0.12               evmix_2.12               
## [137] viridis_0.6.5             reticulate_1.35.0        
## [139] truncdist_1.0-2           stringi_1.8.3            
## [141] ggalluvial_0.12.5         ggraph_2.1.0             
## [143] zlibbioc_1.46.0           MASS_7.3-60.0.1          
## [145] plyr_1.8.9                parallel_4.3.1           
## [147] listenv_0.9.1             ggrepel_0.9.5            
## [149] deldir_2.0-2              graphlayouts_1.1.0       
## [151] splines_4.3.1             hash_2.2.6.3             
## [153] tensor_1.5                hms_1.1.3                
## [155] igraph_1.6.0              spatstat.geom_3.2-8      
## [157] cubature_2.1.0            rngtools_1.5.2           
## [159] RcppHNSW_0.6.0            reshape2_1.4.4           
## [161] ScaledMatrix_1.8.1        evaluate_0.23            
## [163] tzdb_0.4.0                foreach_1.5.2            
## [165] tweenr_2.0.2              httpuv_1.6.14            
## [167] MatrixModels_0.5-3        RANN_2.6.1               
## [169] tidyr_1.3.1               openssl_2.1.1            
## [171] purrr_1.0.2               polyclip_1.10-6          
## [173] future_1.33.1             scattermore_1.2          
## [175] ggforce_0.4.1             rsvd_1.0.5               
## [177] xtable_1.8-4              RSpectra_0.16-1          
## [179] later_1.3.2               ragg_1.2.7               
## [181] viridisLite_0.4.2         gsl_2.1-8                
## [183] tibble_3.2.1              beeswarm_0.4.0           
## [185] cluster_2.1.6             corrplot_0.92            
## [187] globals_0.16.2