B-ALL Immune Cell Reintegration
2026-02-20
Overview
This script reintegrates the non-leukaemia immune cell populations from the B-ALL scRNAseq dataset. Three objects are used:
- no_leuk_obj_anno (base object,
05_srt_no_leukaemia_clustered_scType_annotated.qs) — all 19,040 non-leukaemia cells with scType automated annotations in$scType_classification. Note: the$annotationscolumn in this object contains logicalFALSEvalues from an earlier filtering step and is overwritten with scType labels at load time. - no_leuk_gdt (
10_gdTcells_clustered.qs) — γδ T cell subset with refined annotations (Vγ6+Vδ4+vsγδ_like), functional bias scores (AntiTumorScore1,ProTumorScore1), andcluster_class - no_leuk_cd8
(
10_Teff_Tex_Tpex_subset.qs) — CD8 T cell subset with refined annotations (CD8_Tex,CD8_Tpex,T_eff)
Important note on batch correction: Unlike the leukaemia cells, which were transplanted from three donor sources (AD4_EM3, AD2.2_EM4, AD5_EM4) and required Harmony correction for transplant source, the immune cells are the endogenous immune system of each recipient mouse. They are entirely independent of transplant source. The relevant grouping variable for any batch consideration is Mouse_ID (six individual mice, each an independent biological replicate). Whether Harmony correction is needed at all is assessed from the standard PCA/UMAP — if Mouse_ID effects are not driving structure in the UMAP, no correction is applied.
The annotation pipeline: (1) doublets removed via scDblFinder; (2) BM
macrophages and CNS microglia-like macrophages separated by marker-based
scoring; (3) refined CD8 and γδ annotations transferred by barcode
matching; (4) all γδ T cells reclassified by V gene expression from the
count matrix; (5) ambiguous populations assessed before exclusion
decisions; (6) unified cell_annotation column built; (7)
standard UMAP computed with optional Harmony correction by Mouse_ID if
individual mouse effects are visible.
Setup
library(Seurat)## Loading required package: SeuratObject## Loading required package: sp##
## Attaching package: 'SeuratObject'## The following objects are masked from 'package:base':
##
## intersect, tlibrary(ggplot2)
library(patchwork)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(knitr)
library(qs2)## qs2 0.1.7library(viridis)## Loading required package: viridisLitelibrary(tidyr)
library(tibble)
library(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:viridis':
##
## viridis_pallibrary(harmony)## Loading required package: Rcpplibrary(ggrepel)tissue_cols <- c("BM" = "steelblue", "CNS" = "firebrick")
# Six individual mice — colour palette for Mouse_ID
mouse_cols <- c(
"1838161" = "#E69F00",
"1838162" = "#56B4E9",
"1838163" = "#009E73",
"1838171" = "#F0E442",
"1838279" = "#CC79A7",
"1838280" = "#D55E00"
)base_path <- "/exports/eddie/scratch/aduguid3/uploaded_data"
no_leuk_obj_anno <- qs_read(file.path(base_path,
"05_srt_no_leukaemia_clustered_scType_annotated.qs"))
no_leuk_gdt <- qs_read(file.path(base_path,
"10_gdTcells_clustered.qs"))
no_leuk_cd8 <- qs_read(file.path(base_path,
"10_Teff_Tex_Tpex_subset.qs"))
cat("Base object cells:", ncol(no_leuk_obj_anno), "\n")## Base object cells: 19040cat("γδ T cell cells: ", ncol(no_leuk_gdt), "\n")## γδ T cell cells: 358cat("CD8 T cell cells: ", ncol(no_leuk_cd8), "\n")## CD8 T cell cells: 499cat("\nBase object tissues:\n")##
## Base object tissues:print(table(no_leuk_obj_anno$Tissue))##
## BM CNS
## 13799 5241cat("\nCells per mouse:\n")##
## Cells per mouse:print(table(no_leuk_obj_anno$Mouse_ID, no_leuk_obj_anno$Tissue))##
## BM CNS
## 1838161 373 923
## 1838162 593 732
## 1838163 403 535
## 1838171 4662 1144
## 1838279 5173 756
## 1838280 2595 1151# $annotations contains logical FALSE — overwrite with scType labels
cat("\nscType_classification distribution:\n")##
## scType_classification distribution:print(sort(table(no_leuk_obj_anno$scType_classification), decreasing = TRUE))##
## Macrophages
## 3625
## Pre-B cells
## 3159
## Naive CD8+ T cells
## 1798
## γδ-T cells
## 1594
## Platelets
## 1512
## Neutrophils
## 1200
## Natural killer cells
## 999
## Plasma B cells
## 901
## Effector CD4+ T cells
## 872
## Progenitor cells
## 753
## Non-classical monocytes
## 616
## Basophils
## 588
## Effector CD8+ T cells
## 548
## Erythroid-like and erythroid precursor cells
## 315
## CD8+ NKT-like cells
## 308
## Naive B cells
## 129
## Myeloid Dendritic cells
## 123no_leuk_obj_anno$annotations <- as.character(
no_leuk_obj_anno$scType_classification)
cat("\nannotations set from scType_classification:",
length(unique(no_leuk_obj_anno$annotations)), "unique values\n")##
## annotations set from scType_classification: 17 unique valuesStep 1 — Doublet Check and Removal
cat("scDblFinder.class distribution:\n")## scDblFinder.class distribution:print(table(no_leuk_obj_anno$scDblFinder.class, useNA = "ifany"))##
## singlet doublet
## 17990 1050n_before <- ncol(no_leuk_obj_anno)
if ("doublet" %in% no_leuk_obj_anno$scDblFinder.class) {
no_leuk_obj_anno <- subset(no_leuk_obj_anno,
subset = scDblFinder.class == "singlet")
cat("\nDoublets removed via scDblFinder.class\n")
} else {
cat("\nNo doublets detected — object already filtered or none present\n")
}##
## Doublets removed via scDblFinder.classcat("Cells before:", n_before, "\n")## Cells before: 19040cat("Cells after: ", ncol(no_leuk_obj_anno), "\n")## Cells after: 17990cat("Removed: ", n_before - ncol(no_leuk_obj_anno), "\n")## Removed: 1050Step 2 — Macrophage Refinement
Macrophages classified by scType are refined into two subpopulations using marker-based scoring. Microglia markers (Cx3cr1, P2ry12, Tmem119, Sall1, Fcrls, Olfml3, Hexb) and BM macrophage markers (Csf1r, Adgre1, Cd68, Mrc1, Ccl2) are scored onto all cells via AddModuleScore. Macrophages scoring higher on the microglia signature are labelled “Microglia-like macrophages”.
These cells are CNS border-associated macrophages (likely meningeal or perivascular BAMs) rather than true parenchymal microglia. The term “Microglia-like macrophages” reflects this ambiguity — BAMs express canonical microglia markers but are ontogenetically and transcriptionally distinct from homeostatic parenchymal microglia.
Validation: Microglia-like macrophages were 95.9% CNS (mean microglia sig 3.60 vs BM sig 1.19). Macrophages were 44% CNS (mean BM sig 1.97 vs microglia sig 0.31). Score separation is robust.
micro_markers <- c("Cx3cr1", "P2ry12", "Tmem119", "Sall1",
"Fcrls", "Olfml3", "Hexb")
bm_macro_markers <- c("Csf1r", "Adgre1", "Cd68", "Mrc1", "Ccl2")
no_leuk_obj_anno <- AddModuleScore(no_leuk_obj_anno,
features = list(micro_markers),
name = "Microglia_sig",
seed = 42)
no_leuk_obj_anno <- AddModuleScore(no_leuk_obj_anno,
features = list(bm_macro_markers),
name = "BMMacro_sig",
seed = 42)
no_leuk_obj_anno$macro_refined <- ifelse(
no_leuk_obj_anno$scType_classification == "Macrophages" &
no_leuk_obj_anno$Microglia_sig1 > no_leuk_obj_anno$BMMacro_sig1,
"Microglia-like macrophages",
ifelse(
no_leuk_obj_anno$scType_classification == "Macrophages",
"Macrophages",
NA_character_
)
)
cat("Refined macrophage labels × Tissue:\n")## Refined macrophage labels × Tissue:print(table(no_leuk_obj_anno$macro_refined,
no_leuk_obj_anno$Tissue, useNA = "ifany"))##
## BM CNS
## Macrophages 1127 904
## Microglia-like macrophages 64 1476
## <NA> 11737 2682no_leuk_obj_anno@meta.data %>%
filter(scType_classification == "Macrophages") %>%
group_by(macro_refined) %>%
summarise(
n = n(),
mean_micro_sig = round(mean(Microglia_sig1), 4),
mean_bm_sig = round(mean(BMMacro_sig1), 4),
pct_CNS = round(mean(Tissue == "CNS") * 100, 1),
.groups = "drop"
) %>%
kable(caption = "Macrophage refinement — score summary and tissue distribution")| macro_refined | n | mean_micro_sig | mean_bm_sig | pct_CNS |
|---|---|---|---|---|
| Macrophages | 2031 | 0.3342 | 1.9996 | 44.5 |
| Microglia-like macrophages | 1540 | 3.5698 | 1.1658 | 95.8 |
Step 3 — Transfer Refined Annotations From Subset Objects
Refined γδ T cell and CD8 T cell annotations are transferred from
subset objects by exact barcode matching. The γδ subset
(no_leuk_gdt, 358 cells) contains fewer cells than the
scType γδ-T cells population (728 cells) — this discrepancy is resolved
in Step 4 by V gene reclassification of all γδ T cells.
gdt_meta <- no_leuk_gdt@meta.data %>%
rownames_to_column("barcode") %>%
select(barcode,
gdt_annotation = annotations,
gdt_bias = gamma_delta_bias,
gdt_class = cluster_class,
AntiTumorScore1,
ProTumorScore1)
base_barcodes <- rownames(no_leuk_obj_anno@meta.data)
gdt_match <- match(base_barcodes, gdt_meta$barcode)
no_leuk_obj_anno$gdt_annotation <- gdt_meta$gdt_annotation[gdt_match]
no_leuk_obj_anno$gdt_bias <- gdt_meta$gdt_bias[gdt_match]
no_leuk_obj_anno$gdt_class <- gdt_meta$gdt_class[gdt_match]
no_leuk_obj_anno$AntiTumorScore1 <- gdt_meta$AntiTumorScore1[gdt_match]
no_leuk_obj_anno$ProTumorScore1 <- gdt_meta$ProTumorScore1[gdt_match]
cat("γδ barcodes matched:", sum(!is.na(no_leuk_obj_anno$gdt_annotation)),
"of", nrow(gdt_meta), "\n")## γδ barcodes matched: 357 of 358cd8_meta <- no_leuk_cd8@meta.data %>%
rownames_to_column("barcode") %>%
select(barcode, cd8_annotation = annotations)
cd8_match <- match(base_barcodes, cd8_meta$barcode)
no_leuk_obj_anno$cd8_annotation <- cd8_meta$cd8_annotation[cd8_match]
cat("CD8 barcodes matched:", sum(!is.na(no_leuk_obj_anno$cd8_annotation)),
"of", nrow(cd8_meta), "\n")## CD8 barcodes matched: 497 of 499Step 4 — Population QC and γδ Reclassification
Pre-B Cells
Pre-B cells (2,893 BM / 205 CNS) — assessed to distinguish normal B cell progenitors from leukaemia escapees. High Cd79a (99.7%) and Cd19 (83.2%) confirm B lineage. Rag1 (12.3%) confirms active VDJ recombination. Cd34 (0.1%) and Dntt (0.5%) essentially absent — rules out leukaemia identity. Decision: retain.
preb_cells <- subset(no_leuk_obj_anno,
subset = scType_classification == "Pre-B cells")
cat("Pre-B cells:", ncol(preb_cells), "\n")## Pre-B cells: 3098print(table(preb_cells$Tissue))##
## BM CNS
## 2893 205preb_markers <- c("Cd19","Vpreb1","Igll1","Cd79a",
"Cd34","Flt3","Il7r",
"Dntt","Rag1","Rag2",
"Cd24","Cd43","Cd25")
preb_found <- preb_markers[preb_markers %in% rownames(preb_cells)]
sapply(preb_found, function(g) {
round(mean(GetAssayData(preb_cells, layer = "counts")[g, ] > 0) * 100, 2)
}) %>%
sort(decreasing = TRUE) %>%
data.frame(pct_detected = .) %>%
rownames_to_column("gene") %>%
kable(caption = "Pre-B cell marker detection rates")| gene | pct_detected |
|---|---|
| Cd79a | 99.71 |
| Cd19 | 83.21 |
| Il7r | 18.62 |
| Rag1 | 12.30 |
| Rag2 | 3.58 |
| Flt3 | 2.26 |
| Igll1 | 0.58 |
| Dntt | 0.52 |
| Cd34 | 0.10 |
Platelets — Reclassified as Stroma
ScType classified 1,442 cells as Platelets (719 BM / 723 CNS). The equal tissue split is biologically implausible for true platelets. Marker assessment confirmed stromal/endothelial identity: Col1a1 (54.8%), Fn1 (48.0%), Col3a1 (47.3%), Pdgfra (37.0%). Canonical platelet markers absent: Ppbp (1.6%), Pf4 (8.5%). Decision: relabel as Stroma and retain — relevant for downstream cell-cell interaction analysis.
plt_cells <- subset(no_leuk_obj_anno,
subset = scType_classification == "Platelets")
cat("Platelet cells:", ncol(plt_cells), "\n")## Platelet cells: 1442print(table(plt_cells$Tissue))##
## BM CNS
## 719 723plt_markers <- c("Ppbp","Pf4","Itga2b","Vwf",
"Col1a1","Col3a1","Pdgfra","Fn1",
"Pecam1","Cdh5","Eng",
"Acta2","Tagln")
plt_found <- plt_markers[plt_markers %in% rownames(plt_cells)]
sapply(plt_found, function(g) {
round(mean(GetAssayData(plt_cells, layer = "counts")[g, ] > 0) * 100, 2)
}) %>%
sort(decreasing = TRUE) %>%
data.frame(pct_detected = .) %>%
rownames_to_column("gene") %>%
kable(caption = "Platelet/stromal marker detection — scType 'Platelets' cells")| gene | pct_detected |
|---|---|
| Col1a1 | 54.79 |
| Fn1 | 47.99 |
| Col3a1 | 47.30 |
| Pdgfra | 36.96 |
| Eng | 27.67 |
| Pecam1 | 24.62 |
| Cdh5 | 23.93 |
| Itga2b | 20.60 |
| Tagln | 10.06 |
| Pf4 | 8.53 |
| Vwf | 8.18 |
| Acta2 | 2.98 |
| Ppbp | 1.60 |
Erythroid-like Cells
BM-restricted (289 BM / 2 CNS), consistent with sort impurity. Non-immune cells not relevant to immune microenvironment analysis. Decision: exclude.
ery_cells <- subset(no_leuk_obj_anno,
subset = scType_classification ==
"Erythroid-like and erythroid precursor cells")
cat("Erythroid cells:", ncol(ery_cells), "\n")## Erythroid cells: 291print(table(ery_cells$Tissue))##
## BM CNS
## 289 2ery_markers <- c("Hba-a1","Hbb-bt","Gypa","Klf1","Gata1","Epor","Tfrc")
ery_found <- ery_markers[ery_markers %in% rownames(ery_cells)]
sapply(ery_found, function(g) {
round(mean(GetAssayData(ery_cells, layer = "counts")[g, ] > 0) * 100, 2)
}) %>%
sort(decreasing = TRUE) %>%
data.frame(pct_detected = .) %>%
rownames_to_column("gene") %>%
kable(caption = "Erythroid marker detection rates")| gene | pct_detected |
|---|---|
| Klf1 | 87.63 |
| Tfrc | 86.25 |
| Gypa | 84.54 |
| Gata1 | 83.51 |
| Hba-a1 | 80.76 |
| Hbb-bt | 73.20 |
| Epor | 73.20 |
γδ T Cell Reclassification by V Gene Expression
Rather than relying on barcode matching between
no_leuk_gdt (358 cells) and the base object (728 scType
γδ-T cells), all γδ T cells are reclassified using V gene expression
directly from the count matrix. This is more robust as it is based on
transcriptional evidence present in the data and provides biologically
meaningful Vγ/Vδ labels for every γδ T cell regardless of object
provenance.
All seven Trgv genes and four Trdv genes were detected. The dominant CNS-enriched population was Vγ6Vδ4 (2 BM / 119 CNS), confirming CNS-restricted biology consistent with the role of Vγ6Vδ4 T cells as tissue-resident IL-17 producers at CNS barriers.
Labels consolidated into five groups: Vγ6Vδ4, Vδ4 (other Vγ) (Vδ4+ without confirmed Vγ6), Vδ5, γδ-T cells (Vδ unknown) (γ chain detected, δ chain not captured), and γδ-T cells (unclassified) (no V gene detected).
trgv_genes <- c("Trgv1","Trgv2","Trgv3","Trgv4","Trgv5","Trgv6","Trgv7")
trdv_genes <- c("Trdv1","Trdv2","Trdv3","Trdv4","Trdv5")
trgv_found <- trgv_genes[trgv_genes %in% rownames(no_leuk_obj_anno)]
trdv_found <- trdv_genes[trdv_genes %in% rownames(no_leuk_obj_anno)]
cat("Trgv genes found:", paste(trgv_found, collapse = ", "), "\n")## Trgv genes found: Trgv1, Trgv2, Trgv3, Trgv4, Trgv5, Trgv6, Trgv7cat("Trdv genes found:", paste(trdv_found, collapse = ", "), "\n")## Trdv genes found: Trdv1, Trdv3, Trdv4, Trdv5gdt_all <- subset(no_leuk_obj_anno,
subset = scType_classification == "γδ-T cells")
count_mat <- GetAssayData(gdt_all, layer = "counts")
trgv_max <- if (length(trgv_found) > 0) {
apply(count_mat[trgv_found, , drop = FALSE], 2, function(x) {
if (max(x) > 0) trgv_found[which.max(x)] else NA_character_
})
} else rep(NA_character_, ncol(gdt_all))
trdv_max <- if (length(trdv_found) > 0) {
apply(count_mat[trdv_found, , drop = FALSE], 2, function(x) {
if (max(x) > 0) trdv_found[which.max(x)] else NA_character_
})
} else rep(NA_character_, ncol(gdt_all))
vgene_label <- case_when(
!is.na(trgv_max) & !is.na(trdv_max) ~
paste0("Vγ", gsub("Trgv","",trgv_max), "Vδ", gsub("Trdv","",trdv_max)),
!is.na(trgv_max) & is.na(trdv_max) ~
paste0("Vγ", gsub("Trgv","",trgv_max), "Vδ?"),
is.na(trgv_max) & !is.na(trdv_max) ~
paste0("Vγ?Vδ", gsub("Trdv","",trdv_max)),
TRUE ~ "γδ-T cells (no V gene detected)"
)
cat("\nV gene distribution:\n")##
## V gene distribution:print(sort(table(vgene_label), decreasing = TRUE))## vgene_label
## γδ-T cells (no V gene detected) Vγ?Vδ4
## 388 153
## Vγ2Vδ? Vγ6Vδ4
## 139 121
## Vγ1Vδ? Vγ4Vδ?
## 98 66
## Vγ5Vδ? Vγ5Vδ4
## 39 24
## Vγ2Vδ4 Vγ1Vδ4
## 22 17
## Vγ4Vδ4 Vγ2Vδ5
## 5 4
## Vγ3Vδ? Vγ7Vδ?
## 3 3
## Vγ6Vδ? Vγ?Vδ5
## 2 1
## Vγ1Vδ1 Vγ4Vδ5
## 1 1
## Vγ7Vδ4
## 1cat("\nV gene labels × Tissue:\n")##
## V gene labels × Tissue:print(table(vgene_label, gdt_all$Tissue))##
## vgene_label BM CNS
## Vγ?Vδ4 61 92
## Vγ?Vδ5 1 0
## Vγ1Vδ? 61 37
## Vγ1Vδ1 0 1
## Vγ1Vδ4 6 11
## Vγ2Vδ? 87 52
## Vγ2Vδ4 1 21
## Vγ2Vδ5 4 0
## Vγ3Vδ? 1 2
## Vγ4Vδ? 39 27
## Vγ4Vδ4 1 4
## Vγ4Vδ5 0 1
## Vγ5Vδ? 20 19
## Vγ5Vδ4 9 15
## Vγ6Vδ? 0 2
## Vγ6Vδ4 2 119
## Vγ7Vδ? 2 1
## Vγ7Vδ4 0 1
## γδ-T cells (no V gene detected) 219 169# Store in staging column — applied to cell_annotation in Step 5
vgene_vec <- setNames(vgene_label, colnames(gdt_all))
gdt_main_idx <- rownames(no_leuk_obj_anno@meta.data) %in% names(vgene_vec)
no_leuk_obj_anno$vgene_label <- NA_character_
no_leuk_obj_anno$vgene_label[gdt_main_idx] <-
vgene_vec[rownames(no_leuk_obj_anno@meta.data)[gdt_main_idx]]
cat("\nvgene_label assigned to",
sum(!is.na(no_leuk_obj_anno$vgene_label)), "cells\n")##
## vgene_label assigned to 1088 cellsExclusion and Relabelling
# Relabel Platelets → Stroma
stroma_idx <- no_leuk_obj_anno$scType_classification == "Platelets"
no_leuk_obj_anno$annotations[stroma_idx] <- "Stroma"
cat("Platelets relabelled as Stroma:", sum(stroma_idx), "cells\n")## Platelets relabelled as Stroma: 1442 cells# Exclude Erythroid-like cells only
n_before <- ncol(no_leuk_obj_anno)
no_leuk_obj_anno <- subset(
no_leuk_obj_anno,
subset = scType_classification !=
"Erythroid-like and erythroid precursor cells"
)
cat("Cells before exclusion:", n_before, "\n")## Cells before exclusion: 17990cat("Cells after exclusion: ", ncol(no_leuk_obj_anno), "\n")## Cells after exclusion: 17699cat("Excluded (Erythroid): ", n_before - ncol(no_leuk_obj_anno), "\n")## Excluded (Erythroid): 291Step 5 — Build Unified Annotation
# Baseline — scType with Stroma relabelling applied
no_leuk_obj_anno$cell_annotation <- no_leuk_obj_anno$annotations
# Overwrite with macrophage refinement
mac_idx <- !is.na(no_leuk_obj_anno$macro_refined)
no_leuk_obj_anno$cell_annotation[mac_idx] <-
no_leuk_obj_anno$macro_refined[mac_idx]
# Overwrite with refined CD8 annotations
cd8_idx <- !is.na(no_leuk_obj_anno$cd8_annotation)
no_leuk_obj_anno$cell_annotation[cd8_idx] <-
no_leuk_obj_anno$cd8_annotation[cd8_idx]
# Apply V gene labels to γδ T cells
vgene_idx <- !is.na(no_leuk_obj_anno$vgene_label)
no_leuk_obj_anno$cell_annotation[vgene_idx] <-
no_leuk_obj_anno$vgene_label[vgene_idx]
# Consolidate γδ V gene labels
no_leuk_obj_anno$cell_annotation[
no_leuk_obj_anno$cell_annotation %in% c(
"Vγ?Vδ4","Vγ1Vδ4","Vγ2Vδ4",
"Vγ4Vδ4","Vγ5Vδ4","Vγ7Vδ4")] <- "Vδ4 (other Vγ)"
no_leuk_obj_anno$cell_annotation[
no_leuk_obj_anno$cell_annotation %in% c(
"Vγ2Vδ5","Vγ4Vδ5","Vγ?Vδ5")] <- "Vδ5"
no_leuk_obj_anno$cell_annotation[
no_leuk_obj_anno$cell_annotation %in% c(
"Vγ1Vδ?","Vγ2Vδ?","Vγ3Vδ?","Vγ4Vδ?",
"Vγ5Vδ?","Vγ6Vδ?","Vγ7Vδ?","Vγ1Vδ1")] <- "γδ-T cells (Vδ unknown)"
no_leuk_obj_anno$cell_annotation[
no_leuk_obj_anno$cell_annotation %in% c(
"γδ-T cells (no V gene detected)",
"γδ_like")] <- "γδ-T cells (unclassified)"
# Standardise remaining scType CD8 labels
no_leuk_obj_anno$cell_annotation[
no_leuk_obj_anno$cell_annotation == "Effector CD8+ T cells"] <-
"CD8_Teff (unclassified)"
no_leuk_obj_anno$cell_annotation[
no_leuk_obj_anno$cell_annotation == "Naive CD8+ T cells"] <- "CD8_Naive"
cat("Cells with cell_annotation:", sum(!is.na(no_leuk_obj_anno$cell_annotation)),
"of", ncol(no_leuk_obj_anno), "\n")## Cells with cell_annotation: 17699 of 17699cat("\nFull distribution:\n")##
## Full distribution:print(sort(table(no_leuk_obj_anno$cell_annotation), decreasing = TRUE))##
## Pre-B cells Macrophages
## 3098 2031
## CD8_Naive Microglia-like macrophages
## 1768 1540
## Stroma Neutrophils
## 1442 1145
## Natural killer cells Plasma B cells
## 979 844
## Effector CD4+ T cells Progenitor cells
## 814 714
## Non-classical monocytes Basophils
## 606 555
## γδ-T cells (unclassified) γδ-T cells (Vδ unknown)
## 388 351
## CD8+ NKT-like cells CD8_Tex
## 288 284
## Vδ4 (other Vγ) CD8_Tpex
## 222 143
## Vγ6Vδ4 Myeloid Dendritic cells
## 121 117
## Naive B cells CD8_Teff (unclassified)
## 97 80
## T_eff Vδ5
## 66 6cat("\ncell_annotation × Tissue:\n")##
## cell_annotation × Tissue:print(table(no_leuk_obj_anno$cell_annotation, no_leuk_obj_anno$Tissue))##
## BM CNS
## Basophils 523 32
## CD8_Naive 1639 129
## CD8_Teff (unclassified) 49 31
## CD8_Tex 183 101
## CD8_Tpex 113 30
## CD8+ NKT-like cells 173 115
## Effector CD4+ T cells 631 183
## Macrophages 1127 904
## Microglia-like macrophages 64 1476
## Myeloid Dendritic cells 76 41
## Naive B cells 94 3
## Natural killer cells 887 92
## Neutrophils 1103 42
## Non-classical monocytes 356 250
## Plasma B cells 754 90
## Pre-B cells 2893 205
## Progenitor cells 709 5
## Stroma 719 723
## T_eff 32 34
## Vγ6Vδ4 2 119
## Vδ4 (other Vγ) 78 144
## Vδ5 5 1
## γδ-T cells (unclassified) 219 169
## γδ-T cells (Vδ unknown) 210 141Step 6 — Standard PCA, UMAP and Clustering
A standard PCA and UMAP is computed first without any batch correction. Whether Mouse_ID effects are visible in the UMAP is assessed before deciding if Harmony correction is needed. Given these are endogenous immune cells from individual mice — not transplanted cells — Mouse_ID represents biological variation as much as technical batch, and overcorrection could remove genuine inter-individual immune variation.
standard_path <- "/exports/eddie/scratch/aduguid3/harmony_clustering/11_immune_standard.qs"
if (!file.exists(standard_path)) {
# Wipe any existing reductions from the 05_ object
no_leuk_obj_anno@reductions <- list()
no_leuk_obj_anno <- NormalizeData(no_leuk_obj_anno, verbose = FALSE)
no_leuk_obj_anno <- FindVariableFeatures(no_leuk_obj_anno,
selection.method = "vst",
nfeatures = 3000,
verbose = FALSE)
no_leuk_obj_anno <- ScaleData(no_leuk_obj_anno, verbose = FALSE)
no_leuk_obj_anno <- RunPCA(no_leuk_obj_anno, verbose = FALSE)
ElbowPlot(no_leuk_obj_anno, ndims = 50) +
ggtitle("PCA elbow plot — immune cells (standard)")
no_leuk_obj_anno <- RunUMAP(no_leuk_obj_anno,
reduction = "pca",
dims = 1:35)
no_leuk_obj_anno <- FindNeighbors(no_leuk_obj_anno,
reduction = "pca",
dims = 1:35)
no_leuk_obj_anno <- FindClusters(no_leuk_obj_anno,
resolution = 0.3)
qs_save(no_leuk_obj_anno, standard_path)
cat("Standard integration complete — saved\n")
} else {
no_leuk_obj_anno <- qs_read(standard_path)
cat("Loaded existing standard object\n")
}## Loaded existing standard objectcat("Cells:", ncol(no_leuk_obj_anno), "\n")## Cells: 17699Step 7 — Assess Mouse_ID Effects
Before deciding whether to apply Harmony correction, the UMAP is inspected for Mouse_ID clustering. If individual mice intermix well across cell type clusters no correction is needed. If discrete Mouse_ID islands are visible — particularly in shared cell types like macrophages or T cells — Harmony correction by Mouse_ID is warranted.
p_tissue <- DimPlot(no_leuk_obj_anno,
group.by = "Tissue",
cols = tissue_cols,
alpha = 0.5) +
ggtitle("by tissue")
p_mouse <- DimPlot(no_leuk_obj_anno,
group.by = "Mouse_ID",
cols = mouse_cols,
alpha = 0.5) +
ggtitle("by Mouse_ID")
p_tissue + p_mouse
n_types <- length(unique(no_leuk_obj_anno$cell_annotation))
anno_cols <- setNames(
colorRampPalette(c("#1a3a6b","#e67e22","#27ae60","#8b0000",
"#9b59b6","#2980b9","#e74c3c","#1abc9c",
"#f39c12","#7f8c8d","#d35400","#2ecc71",
"#8e44ad","#16a085","#c0392b"))(n_types),
sort(unique(no_leuk_obj_anno$cell_annotation))
)
p <- DimPlot(no_leuk_obj_anno,
group.by = "cell_annotation",
cols = anno_cols,
label = FALSE,
pt.size = 0.3) +
ggtitle("Immune cells — standard UMAP") +
theme(legend.position = "right",
legend.text = element_text(size = 7)) +
guides(colour = guide_legend(override.aes = list(size = 3)))
LabelClusters(p,
id = "cell_annotation",
fontface = "bold",
size = 3,
repel = TRUE,
box = FALSE,
force = 15,
max.overlaps = 30)
Step 8 — Optional Harmony Correction by Mouse_ID
Harmony correction by Mouse_ID is applied only if individual mouse effects are clearly driving UMAP structure assessed in Step 7. Unlike the leukaemia pipeline where correction for transplant source removes a known technical confound (donor genetics), correction of endogenous immune cells by Mouse_ID removes genuine inter-individual biological variation alongside any technical effects and should only be applied if Mouse_ID effects are clearly visible and distorting cell type separation.
# Set apply_harmony <- TRUE if Mouse_ID effects are clearly visible in Step 7
apply_harmony <- FALSE # UPDATE AFTER INSPECTING STEP 7 UMAP
harmony_path <- "/exports/eddie/scratch/aduguid3/harmony_clustering/11_immune_harmony.qs"
if (apply_harmony) {
if (!file.exists(harmony_path)) {
no_leuk_obj_anno@reductions[["umap"]] <- NULL
no_leuk_obj_anno@reductions[["harmony"]] <- NULL
no_leuk_obj_anno <- harmony::RunHarmony(no_leuk_obj_anno,
group.by.vars = "Mouse_ID",
reduction = "pca",
reduction.save = "harmony",
max.iter.harmony = 20,
theta = 2,
verbose = TRUE)
no_leuk_obj_anno <- RunUMAP(no_leuk_obj_anno,
reduction = "harmony",
dims = 1:35)
no_leuk_obj_anno <- FindNeighbors(no_leuk_obj_anno,
reduction = "harmony",
dims = 1:35)
no_leuk_obj_anno <- FindClusters(no_leuk_obj_anno,
resolution = 0.3)
qs_save(no_leuk_obj_anno, harmony_path)
cat("Harmony (Mouse_ID) integration complete — saved\n")
} else {
no_leuk_obj_anno <- qs_read(harmony_path)
cat("Loaded existing Harmony object\n")
}
# Reassess Mouse_ID mixing after correction
p_mouse_post <- DimPlot(no_leuk_obj_anno,
group.by = "Mouse_ID",
cols = mouse_cols,
alpha = 0.5) +
ggtitle("by Mouse_ID — post Harmony")
p_tissue_post <- DimPlot(no_leuk_obj_anno,
group.by = "Tissue",
cols = tissue_cols,
alpha = 0.5) +
ggtitle("by tissue — post Harmony")
p_tissue_post + p_mouse_post
} else {
cat("Harmony correction not applied — using standard PCA/UMAP\n")
cat("Set apply_harmony <- TRUE in this chunk if Mouse_ID effects\n")
cat("are clearly visible in the Step 7 UMAP\n")
}## Harmony correction not applied — using standard PCA/UMAP
## Set apply_harmony <- TRUE in this chunk if Mouse_ID effects
## are clearly visible in the Step 7 UMAPStep 9 — Cell Type Composition by Tissue
comp <- no_leuk_obj_anno@meta.data %>%
group_by(Tissue, cell_annotation) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(Tissue) %>%
mutate(prop = n / sum(n))
ggplot(comp, aes(x = Tissue, y = prop, fill = cell_annotation)) +
geom_bar(stat = "identity") +
scale_fill_manual(values = anno_cols) +
scale_y_continuous(labels = scales::percent) +
theme_classic(base_size = 12) +
labs(x = NULL, y = "Proportion", fill = NULL,
title = "Immune cell composition by tissue") +
theme(legend.position = "right",
legend.text = element_text(size = 8))
comp %>%
mutate(prop = round(prop, 3)) %>%
pivot_wider(names_from = Tissue,
values_from = c(n, prop),
values_fill = 0) %>%
arrange(desc(n_CNS)) %>%
kable(caption = "Immune cell composition by tissue — counts and proportions")| cell_annotation | n_BM | n_CNS | prop_BM | prop_CNS |
|---|---|---|---|---|
| Microglia-like macrophages | 64 | 1476 | 0.005 | 0.292 |
| Macrophages | 1127 | 904 | 0.089 | 0.179 |
| Stroma | 719 | 723 | 0.057 | 0.143 |
| Non-classical monocytes | 356 | 250 | 0.028 | 0.049 |
| Pre-B cells | 2893 | 205 | 0.229 | 0.041 |
| Effector CD4+ T cells | 631 | 183 | 0.050 | 0.036 |
| γδ-T cells (unclassified) | 219 | 169 | 0.017 | 0.033 |
| Vδ4 (other Vγ) | 78 | 144 | 0.006 | 0.028 |
| γδ-T cells (Vδ unknown) | 210 | 141 | 0.017 | 0.028 |
| CD8_Naive | 1639 | 129 | 0.130 | 0.025 |
| Vγ6Vδ4 | 2 | 119 | 0.000 | 0.024 |
| CD8+ NKT-like cells | 173 | 115 | 0.014 | 0.023 |
| CD8_Tex | 183 | 101 | 0.014 | 0.020 |
| Natural killer cells | 887 | 92 | 0.070 | 0.018 |
| Plasma B cells | 754 | 90 | 0.060 | 0.018 |
| Neutrophils | 1103 | 42 | 0.087 | 0.008 |
| Myeloid Dendritic cells | 76 | 41 | 0.006 | 0.008 |
| T_eff | 32 | 34 | 0.003 | 0.007 |
| Basophils | 523 | 32 | 0.041 | 0.006 |
| CD8_Teff (unclassified) | 49 | 31 | 0.004 | 0.006 |
| CD8_Tpex | 113 | 30 | 0.009 | 0.006 |
| Progenitor cells | 709 | 5 | 0.056 | 0.001 |
| Naive B cells | 94 | 3 | 0.007 | 0.001 |
| Vδ5 | 5 | 1 | 0.000 | 0.000 |
Step 10 — γδ T Cell Annotations on UMAP
gdt_labels <- c("Vγ6Vδ4", "Vδ4 (other Vγ)", "Vδ5",
"γδ-T cells (Vδ unknown)",
"γδ-T cells (unclassified)")
gdt_cols <- c(
"Vγ6Vδ4" = "#e74c3c",
"Vδ4 (other Vγ)" = "#e67e22",
"Vδ5" = "#f39c12",
"γδ-T cells (Vδ unknown)" = "#9b59b6",
"γδ-T cells (unclassified)" = "#bdc3c7",
"Other" = "grey90"
)
no_leuk_obj_anno$gdt_highlight <- ifelse(
no_leuk_obj_anno$cell_annotation %in% gdt_labels,
no_leuk_obj_anno$cell_annotation,
"Other"
)
DimPlot(no_leuk_obj_anno,
group.by = "gdt_highlight",
cols = gdt_cols,
order = gdt_labels,
pt.size = 0.5) +
ggtitle("γδ T cells — V gene populations") +
theme(legend.position = "bottom",
legend.text = element_text(size = 8))
no_leuk_obj_anno$vg6vd4_highlight <- ifelse(
no_leuk_obj_anno$cell_annotation == "Vγ6Vδ4",
"Vγ6Vδ4", "Other"
)
DimPlot(no_leuk_obj_anno,
group.by = "vg6vd4_highlight",
cols = c("Vγ6Vδ4" = "#e74c3c", "Other" = "grey90"),
order = "Vγ6Vδ4",
split.by = "Tissue",
pt.size = 0.8) +
ggtitle("Vγ6Vδ4 T cells — BM vs CNS") +
theme(legend.position = "bottom")
no_leuk_obj_anno@meta.data %>%
filter(cell_annotation %in% gdt_labels) %>%
group_by(cell_annotation, Tissue) %>%
summarise(n = n(), .groups = "drop") %>%
pivot_wider(names_from = Tissue, values_from = n, values_fill = 0) %>%
mutate(total = BM + CNS,
pct_CNS = round(CNS / total * 100, 1)) %>%
arrange(desc(pct_CNS)) %>%
kable(caption = "γδ T cell populations by tissue")| cell_annotation | BM | CNS | total | pct_CNS |
|---|---|---|---|---|
| Vγ6Vδ4 | 2 | 119 | 121 | 98.3 |
| Vδ4 (other Vγ) | 78 | 144 | 222 | 64.9 |
| γδ-T cells (unclassified) | 219 | 169 | 388 | 43.6 |
| γδ-T cells (Vδ unknown) | 210 | 141 | 351 | 40.2 |
| Vδ5 | 5 | 1 | 6 | 16.7 |
gdt_scored <- no_leuk_obj_anno@meta.data %>%
filter(!is.na(AntiTumorScore1))
cat("Cells with bias scores (barcode-matched):", nrow(gdt_scored), "\n")## Cells with bias scores (barcode-matched): 357if (nrow(gdt_scored) > 0) {
gdt_scored %>%
pivot_longer(c(AntiTumorScore1, ProTumorScore1),
names_to = "score_type",
values_to = "score") %>%
mutate(score_type = recode(score_type,
"AntiTumorScore1" = "Anti-tumor",
"ProTumorScore1" = "Pro-tumor")) %>%
ggplot(aes(x = Tissue, y = score, fill = Tissue)) +
geom_violin(trim = FALSE, alpha = 0.8) +
geom_boxplot(width = 0.1, fill = "white", outlier.size = 0.2) +
scale_fill_manual(values = tissue_cols) +
facet_wrap(~score_type) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "grey40") +
theme_classic(base_size = 12) +
labs(y = "Module score",
title = "γδ T cell functional bias scores by tissue") +
theme(legend.position = "none")
}
Step 11 — CD8 T Cell Annotations on UMAP
cd8_labels <- c("CD8_Tex","CD8_Tpex","T_eff",
"CD8_Teff (unclassified)","CD8_Naive")
cd8_cols <- c(
"CD8_Tex" = "#8b0000",
"CD8_Tpex" = "#e67e22",
"T_eff" = "#27ae60",
"CD8_Teff (unclassified)"= "#f39c12",
"CD8_Naive" = "#2980b9",
"Other" = "grey90"
)
no_leuk_obj_anno$cd8_highlight <- ifelse(
no_leuk_obj_anno$cell_annotation %in% cd8_labels,
no_leuk_obj_anno$cell_annotation,
"Other"
)
DimPlot(no_leuk_obj_anno,
group.by = "cd8_highlight",
cols = cd8_cols,
order = cd8_labels,
pt.size = 0.5) +
ggtitle("CD8 T cell subtypes") +
theme(legend.position = "bottom",
legend.text = element_text(size = 8))
no_leuk_obj_anno@meta.data %>%
filter(cell_annotation %in% cd8_labels) %>%
group_by(cell_annotation, Tissue) %>%
summarise(n = n(), .groups = "drop") %>%
pivot_wider(names_from = Tissue, values_from = n, values_fill = 0) %>%
mutate(total = BM + CNS,
pct_CNS = round(CNS / total * 100, 1)) %>%
arrange(desc(pct_CNS)) %>%
kable(caption = "CD8 T cell subtypes by tissue")| cell_annotation | BM | CNS | total | pct_CNS |
|---|---|---|---|---|
| T_eff | 32 | 34 | 66 | 51.5 |
| CD8_Teff (unclassified) | 49 | 31 | 80 | 38.8 |
| CD8_Tex | 183 | 101 | 284 | 35.6 |
| CD8_Tpex | 113 | 30 | 143 | 21.0 |
| CD8_Naive | 1639 | 129 | 1768 | 7.3 |
Step 12 — IL-17 Ligand and Receptor Expression Across Immune Cell Types
Leukaemia cells express the non-canonical IL-17 receptor heterodimer Il17ra/Il17rb at ~47% but produce no IL-17 ligands themselves (script 06). Here all IL-17 ligands and receptors are assessed across immune cell types to identify candidate ligand sources. Expression of Il17b or Il25 in CNS-enriched populations — particularly Vγ6Vδ4 T cells or Microglia-like macrophages — would identify a candidate extrinsic signal for the non-canonical Il17ra/Il17rb pathway in leukaemia cells.
il17_ligands <- c("Il17a","Il17b","Il17c","Il17d","Il17f","Il25")
il17_receptors <- c("Il17ra","Il17rb","Il17rc","Il17rd")
il17_detected <- c(il17_ligands, il17_receptors)[
c(il17_ligands, il17_receptors) %in% rownames(no_leuk_obj_anno)]
cat("IL-17 genes detected:", paste(il17_detected, collapse = ", "), "\n")## IL-17 genes detected: Il17a, Il17b, Il17c, Il17d, Il17f, Il25, Il17ra, Il17rb, Il17rc, Il17rdIdents(no_leuk_obj_anno) <- "cell_annotation"
DotPlot(no_leuk_obj_anno,
features = il17_detected,
cols = c("lightgrey", "firebrick"),
dot.scale = 7) +
RotatedAxis() +
ggtitle("IL-17 ligands and receptors — immune cell types") +
theme(axis.text.x = element_text(size = 9),
axis.text.y = element_text(size = 8))
gdt_obj <- subset(no_leuk_obj_anno,
subset = cell_annotation %in% gdt_labels)
Idents(gdt_obj) <- "cell_annotation"
DotPlot(gdt_obj,
features = il17_detected,
cols = c("lightgrey", "#e74c3c"),
dot.scale = 8) +
RotatedAxis() +
ggtitle("IL-17 expression — γδ T cell populations") +
theme(axis.text.x = element_text(size = 9))
il17l_detected <- il17_ligands[il17_ligands %in% rownames(no_leuk_obj_anno)]
DotPlot(no_leuk_obj_anno,
features = il17l_detected,
cols = c("lightgrey", "firebrick"),
split.by = "Tissue",
dot.scale = 5) +
RotatedAxis() +
ggtitle("IL-17 ligands — BM vs CNS across immune cell types") +
theme(axis.text.x = element_text(size = 8),
axis.text.y = element_text(size = 7))
Step 13 — Save Object
output_path <- "/exports/eddie/scratch/aduguid3/harmony_clustering/11_immune_annotated.qs"
qs_save(no_leuk_obj_anno, output_path)
cat("Saved:", output_path, "\n")## Saved: /exports/eddie/scratch/aduguid3/harmony_clustering/11_immune_annotated.qsexpected_cols <- c("Tissue", "Mouse_ID",
"cell_annotation", "annotations",
"scType_classification",
"gdt_annotation", "gdt_bias", "gdt_class",
"AntiTumorScore1", "ProTumorScore1",
"cd8_annotation",
"macro_refined", "Microglia_sig1", "BMMacro_sig1",
"gdt_highlight", "vg6vd4_highlight", "cd8_highlight")
cat("\nKey metadata columns:\n")##
## Key metadata columns:for (col in expected_cols) {
cat(" ", col, ":", col %in% colnames(no_leuk_obj_anno@meta.data), "\n")
}## Tissue : TRUE
## Mouse_ID : TRUE
## cell_annotation : TRUE
## annotations : TRUE
## scType_classification : TRUE
## gdt_annotation : TRUE
## gdt_bias : TRUE
## gdt_class : TRUE
## AntiTumorScore1 : TRUE
## ProTumorScore1 : TRUE
## cd8_annotation : TRUE
## macro_refined : TRUE
## Microglia_sig1 : TRUE
## BMMacro_sig1 : TRUE
## gdt_highlight : TRUE
## vg6vd4_highlight : TRUE
## cd8_highlight : TRUEcat("\nFinal cell count:", ncol(no_leuk_obj_anno), "\n")##
## Final cell count: 17699print(table(no_leuk_obj_anno$Tissue))##
## BM CNS
## 12639 5060Session Information
sessionInfo()## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Rocky Linux 9.5 (Blue Onyx)
##
## Matrix products: default
## BLAS/LAPACK: /opt/intel/oneapi/mkl/2024.0/lib/libmkl_gf_lp64.so.2; LAPACK version 3.10.1
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/London
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggrepel_0.9.6 harmony_1.2.4 Rcpp_1.0.12 scales_1.4.0
## [5] tibble_3.3.1 tidyr_1.3.1 viridis_0.6.5 viridisLite_0.4.2
## [9] qs2_0.1.7 knitr_1.51 dplyr_1.1.4 patchwork_1.3.2
## [13] ggplot2_4.0.2 Seurat_5.4.0 SeuratObject_5.3.0 sp_2.1-4
##
## loaded via a namespace (and not attached):
## [1] deldir_2.0-4 pbapply_1.7-4 gridExtra_2.3
## [4] rlang_1.1.7 magrittr_2.0.3 RcppAnnoy_0.0.23
## [7] otel_0.2.0 spatstat.geom_3.7-0 matrixStats_1.5.0
## [10] ggridges_0.5.6 compiler_4.4.1 png_0.1-8
## [13] vctrs_0.6.5 reshape2_1.4.4 stringr_1.5.1
## [16] pkgconfig_2.0.3 fastmap_1.2.0 labeling_0.4.3
## [19] utf8_1.2.4 promises_1.5.0 rmarkdown_2.30
## [22] purrr_1.0.2 xfun_0.56 cachem_1.1.0
## [25] jsonlite_2.0.0 goftest_1.2-3 later_1.4.6
## [28] spatstat.utils_3.2-1 irlba_2.3.7 parallel_4.4.1
## [31] cluster_2.1.6 R6_2.6.1 ica_1.0-3
## [34] spatstat.data_3.1-9 bslib_0.7.0 stringi_1.8.4
## [37] RColorBrewer_1.1-3 reticulate_1.45.0 spatstat.univar_3.1-6
## [40] parallelly_1.37.1 lmtest_0.9-40 jquerylib_0.1.4
## [43] scattermore_1.2 tensor_1.5.1 future.apply_1.11.2
## [46] zoo_1.8-15 sctransform_0.4.3 httpuv_1.6.16
## [49] Matrix_1.7-0 splines_4.4.1 igraph_2.2.2
## [52] tidyselect_1.2.1 abind_1.4-5 rstudioapi_0.16.0
## [55] dichromat_2.0-0.1 yaml_2.3.12 stringfish_0.18.0
## [58] spatstat.random_3.4-4 codetools_0.2-20 miniUI_0.1.2
## [61] spatstat.explore_3.7-0 listenv_0.9.1 lattice_0.22-6
## [64] plyr_1.8.9 withr_3.0.2 shiny_1.12.1
## [67] S7_0.2.1 ROCR_1.0-12 evaluate_1.0.5
## [70] Rtsne_0.17 future_1.33.2 fastDummies_1.7.5
## [73] survival_3.6-4 RcppParallel_5.1.8 polyclip_1.10-7
## [76] fitdistrplus_1.2-6 pillar_1.9.0 KernSmooth_2.23-24
## [79] plotly_4.12.0 generics_0.1.3 RcppHNSW_0.6.0
## [82] globals_0.16.3 xtable_1.8-4 glue_1.8.0
## [85] lazyeval_0.2.2 tools_4.4.1 data.table_1.15.4
## [88] RSpectra_0.16-2 RANN_2.6.2 dotCall64_1.2
## [91] cowplot_1.2.0 grid_4.4.1 nlme_3.1-164
## [94] cli_3.6.5 spatstat.sparse_3.1-0 spam_2.11-3
## [97] fansi_1.0.6 uwot_0.2.4 gtable_0.3.6
## [100] sass_0.4.9 digest_0.6.35 progressr_0.18.0
## [103] htmlwidgets_1.6.4 farver_2.1.2 htmltools_0.5.8.1
## [106] lifecycle_1.0.5 httr_1.4.7 mime_0.12
## [109] MASS_7.3-60.2