Load packages
# we would like to demonstrate the hybrid strategy using Seurat workflow
library(Seurat)
library(cowplot)
# packages for cell hashing demultiplexing
library(HTOreader)
library(cellhashR) # this package integrated the following methods: multiseq, htodemux, gmm_demux, bff_raw, bff_cluster and dropletutilsStep 1 Load data from 10X folder
We load 8pool-CA dataset from the 10X folder using Seurat
# Load in the RNA UMI matrix for doublet negative (DN) dataset and antigen
# specific (AS) dataset
DN.data <- Read10X(data.dir = "../data/8pool-Carr-RNA/filtered_feature_bc_matrix/")split different modalities: Cell hashing (HTO), antigen probes (Probe) and CITE-seq surface protein panel (ADT) In 8pool DN and AS datasets, we have applied 8 cell hashing tags, 18 antigen probes, and a large CITE-seq panel
DN.probe <- DN.data[["Antibody Capture"]][1:18, ] # expression of antigen probes, Carrier data set doesn't have this measurement
DN.hto <- DN.data[["Antibody Capture"]][19:26, ] # expression of cell hashing tags
DN.ADT <- DN.data[["Antibody Capture"]][27:163, ] # expression of CITE-seq surface protein panelStep 2 Create Seurat object and attached multi-modality data
# create Seurat object
DN <- CreateSeuratObject(counts = DN.data$`Gene Expression`, project = "Carr")
# create multi-modality assay and attach them to Seurat object
DN[["HTO"]] <- CreateAssayObject(counts = DN.hto)
DN[["ADT"]] <- CreateAssayObject(counts = DN.ADT)
# percent of MT genes
DN[["percent.mt"]] <- PercentageFeatureSet(object = DN, pattern = "^MT-")
# percent of IG genes (for B cell study)
DN[["percent.ig"]] <- PercentageFeatureSet(object = DN, pattern = "^IG[HKL]")Step 3 Pre-process
Filter out unwanted cells (optional)
DN <- subset(DN, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt <
10)Normalize HTO and ADT data
DN <- NormalizeData(DN, assay = "HTO", normalization.method = "CLR")
DN <- NormalizeData(DN, assay = "ADT", normalization.method = "CLR")Step 4 Cell hashing demultiplexing
1) Demultiplexing of cell hashing using existing methods integrated in cellhashR: multiseq, gmm_demux, bff_raw, bff_cluster, dropletutils
load data (we hide the plots of this step)
# parse cell hashing output, providing a barcode whitelist. (please refer to
# cellhashR tutorial for more information:
# https://github.com/BimberLab/cellhashR?tab=readme-ov-file#example)
barcodeData <- ProcessCountMatrix(rawCountData = "../data/8pool-Carr-RNA/filtered_feature_bc_matrix/",
barcodeWhitelist = c("hashtag1", "hashtag2", "hashtag3", "hashtag4", "hashtag5",
"hashtag6", "hashtag7", "hashtag8"), datatypeName = "Antibody Capture")rownames(barcodeData) <- c("S282", "S283", "S284", "S344", "S397", "S417", "S421",
"S423")QC (we hide the plots of this step)
# Create QC plots of barcode normalization
PlotNormalizationQC(barcodeData)Calling cell hashing using different methods (we hide the plots of this step)
# Generate the final cell hashing calls
calls <- GenerateCellHashingCalls(barcodeMatrix = barcodeData, methods = c("multiseq",
"gmm_demux", "bff_raw", "bff_cluster", "dropletutils"))# Inspect negative cells:
SummarizeCellsByClassification(calls = calls, barcodeMatrix = barcodeData)
2) Demultiplexing of cell hashing using Seurat (htodemux)
DN <- NormalizeData(DN, assay = "HTO", normalization.method = "CLR")
DN <- HTODemux(DN, assay = "HTO", positive.quantile = 0.99)3) Demultiplexing of cell hashing using HTOreader
DN <- HTOClassification(DN, assay = "HTO", method = "log")
plots <- PlotHTO(DN, assay = "HTO", method = "log")
plot_grid(plotlist = plots, nrow = 2)
4)integrate HTOreader, Seurat results with cellhashR results
rownames(calls) <- calls$cellbarcode
DN_barcode <- names(DN@active.ident)
DN_barcode <- sub("-1", "", DN_barcode)
DN_hq_barcode <- intersect(DN_barcode, calls$cellbarcode)
calls_hq <- calls[DN_hq_barcode, ]
HTOid.hq <- DN$HTOid
names(HTOid.hq) <- sub("-1", "", names(HTOid.hq))
calls_hq$HTOid <- HTOid.hq[rownames(calls_hq)]
Seuratid.hq <- DN$hash.ID
names(Seuratid.hq) <- sub("-1", "", names(Seuratid.hq))
calls_hq$seurat <- as.character(Seuratid.hq[rownames(calls_hq)])res_table <- matrix(data = c(table(calls_hq$multiseq), table(calls_hq$gmm_demux),
table(calls_hq$bff_raw), table(calls_hq$bff_cluster), table(calls_hq$dropletutils),
table(calls_hq$HTOid), table(calls_hq$seurat)), nrow = 7, ncol = 10, byrow = TRUE,
dimnames = list(c("multiseq", "gmm_demux", "bff_raw", "bff_cluster", "dropletutils",
"HTOreader", "htodemux(Seurat)"), names(table(calls_hq$multiseq))))
# singlet rate
singlet_count <- rowSums(res_table[, which(colnames(res_table) %in% rownames(barcodeData))])
singlet_rate <- singlet_count/rowSums(res_table)
final_table <- cbind(res_table, singlet_count, singlet_rate)
colnames(final_table) <- c(colnames(res_table), "singlet", "singlet rate")Cell hashing Demultiplexing results:
final_table## Doublet Negative S282 S283 S284 S344 S397 S417 S421 S423
## multiseq 2103 160 312 1738 1546 1548 1866 1196 780 1635
## gmm_demux 2631 27 310 1696 1498 1483 1706 1165 730 1638
## bff_raw 1673 579 0 1786 1515 1585 2094 1179 794 1679
## bff_cluster 456 404 0 2011 1713 1771 2380 1376 890 1883
## dropletutils 1801 599 267 1769 1463 1512 1953 1184 791 1545
## HTOreader 1893 170 317 1771 1502 1560 2055 1174 780 1662
## htodemux(Seurat) 2759 125 302 1698 1389 1476 1635 1148 727 1625
## singlet singlet rate
## multiseq 10621 0.8243558
## gmm_demux 10226 0.7936976
## bff_raw 10632 0.8252096
## bff_cluster 12024 0.9332505
## dropletutils 10484 0.8137224
## HTOreader 10821 0.8398789
## htodemux(Seurat) 10000 0.7761565We found two BFF models failed to identify S282. The best singlet rate is 83.9%, average rate is ~80% (exlcuding BFF models due to failure of S282)
Step 4, Apply hybrid strategy: combine cell hashing demultiplexing with SNP-based method
In order to apply hybrid demultiplexing, we need to add SNP-based labels into Seurat object:
# add soupercell labels into Seurat object
souporcell.res <- read.csv("./8pool/clusters.tsv", sep = "\t")
souporcell.res$identity <- souporcell.res$status
souporcell.res$identity[which(souporcell.res$status == "singlet")] <- paste0(souporcell.res$status[which(souporcell.res$status ==
"singlet")], souporcell.res$assignment[which(souporcell.res$status == "singlet")])
rownames(souporcell.res) <- souporcell.res$barcode
DN$souporcell <- souporcell.res[names(DN@active.ident), "identity"]Then we reveal cell labels using hybrid demultiplexing (“HybridDemultiplexing” function in HTOreader package). In this case, we use HTOreader result for cell hashing demultiplexing
DN <- HybridDemultiplexing(object = DN, cellhashing_label = "HTOid", genotype_label = "souporcell",
hto_names = c("S282", "S283", "S284", "S344", "S397", "S417", "S421", "S423"))## We matched:
## S282 with singlet6
## S283 with singlet3
## S284 with singlet2
## S344 with singlet1
## S397 with singlet0
## S417 with singlet4
## S421 with singlet7
## S423 with singlet5
## The convergence score between cell hashing demultiplexing and SNP-based demultiplexing is: 0.993251476239573## Hybrid demultiplexing done! Results in hybridIDtable(DN$hybridID)##
## doublet S282 S283 S284 S344 S397 S417
## 1173 278 1859 1629 1613 2243 1252
## S421 S423 unassigned
## 857 1749 233BTW, The Hybrid strategy is highly flexiable to cell hashing demultiplexing method and quality, it works with all cell hashing methods. For example:
# Hybrid demultiplexing using HTOdemux (Seurat)
DN <- HybridDemultiplexing(object = DN, cellhashing_label = "hash.ID", genotype_label = "souporcell",
hto_names = c("S282", "S283", "S284", "S344", "S397", "S417", "S421", "S423"))## We matched:
## S282 with singlet6
## S283 with singlet3
## S284 with singlet2
## S344 with singlet1
## S397 with singlet0
## S417 with singlet4
## S421 with singlet7
## S423 with singlet5
## The convergence score between cell hashing demultiplexing and SNP-based demultiplexing is: 0.993932038834951## Hybrid demultiplexing done! Results in hybridIDtable(DN$hybridID)##
## doublet S282 S283 S284 S344 S397 S417
## 1213 279 1859 1628 1614 2242 1256
## S421 S423 unassigned
## 865 1749 181Furthermore, The Hybrid strategy even works with two BFF models that fail to identify S282. HybridDemultiplexing correctly linked singlet6 with S282.
rownames(calls) <- paste0(calls$cellbarcode, "-1")
DN$bff_cluster <- calls[names(DN@active.ident), "bff_cluster"]
DN <- HybridDemultiplexing(object = DN, cellhashing_label = "bff_cluster", genotype_label = "souporcell",
hto_names = c("S282", "S283", "S284", "S344", "S397", "S417", "S421", "S423"))## The missing HTO S282 is assigned to singlet6## We matched:
## S283 with singlet3
## S284 with singlet2
## S344 with singlet1
## S397 with singlet0
## S417 with singlet4
## S421 with singlet7
## S423 with singlet5
## S282 with singlet6
## The convergence score between cell hashing demultiplexing and SNP-based demultiplexing is: 0.985078676071622## Hybrid demultiplexing done! Results in hybridIDWe observe that results from hybrid runs using various cell hashing methods may exhibit minor variations. This is due to our adoption of a conservative strategy, which categorizes all Case 2/3 inconsistency as unassigned to guarantee the accuracy and precision of the outcomes. For detailed information, please refer to the Methods section of our paper.
table(DN$hybridID)##
## doublet S282 S283 S284 S344 S397 S417
## 311 235 1846 1615 1599 2222 1240
## S421 S423 unassigned
## 893 1737 1188Step 5, hybrid demultiplexing for multiple datasets from the same group of donors (optional)
In this case, we have 8pool-AS, which shares the same group of human donors but doesn’t have cell hashing for some reason. Using Hybrid demultiplexing, we’re still able to handle it. We load the 8pool-AS data first:
AS.data <- Read10X(data.dir = "../data/8pool-AS-RNA/filtered_feature_bc_matrix/")
AS.probe <- AS.data[["Antibody Capture"]][1:18, ] # expression of antigen probes
AS.hto <- AS.data[["Antibody Capture"]][19:26, ] # expression of cell hashing tags, Antigen specific data set doesn't have this measurement, they have 0 counts
AS.ADT <- AS.data[["Antibody Capture"]][27:163, ] # expression of CITE-seq surface protein panel
AS <- CreateSeuratObject(counts = AS.data$`Gene Expression`, project = "AS")
AS[["ADT"]] <- CreateAssayObject(counts = AS.ADT)
AS[["Probe"]] <- CreateAssayObject(counts = AS.probe)
AS[["percent.mt"]] <- PercentageFeatureSet(object = AS, pattern = "^MT-")
AS[["percent.ig"]] <- PercentageFeatureSet(object = AS, pattern = "^IG[HKL]")
AS <- subset(AS, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt <
10)
AS <- NormalizeData(AS, assay = "Probe", normalization.method = "CLR")
AS <- NormalizeData(AS, assay = "ADT", normalization.method = "CLR")We ran souporcell on the entire 8pool dataset (AS + DN) so that the genotypes are across two datasets. We annotate souporcell genotype clusters to 8pool-AS dataset
souporcell.res <- read.csv("./8pool/clusters.tsv", sep = "\t")
souporcell.res$identity <- souporcell.res$status
souporcell.res$identity[which(souporcell.res$status == "singlet")] <- paste0(souporcell.res$status[which(souporcell.res$status ==
"singlet")], souporcell.res$assignment[which(souporcell.res$status == "singlet")])
rownames(souporcell.res) <- souporcell.res$barcode
AS$souporcell <- souporcell.res[names(AS@active.ident), "identity"]We get high singlet rate for 8pool-AS using hybrid demultiplexing and saved reagents of cell hashing on 8pool-AS. Beside the economy, there is a more critical reason for not using cell hashing in 8pool-AS dataset: we aim to measure the binding affinity between antigen specific B cells and our antigen probes. We don’t want any disruptions from the cell hashing due to they share the same channel(Antibody Capture)
Singlet rate of hybrid on 8pool-AS:
(1 - length(which(AS$souporcell %in% c("doublet", "unassigned")))/length(AS$souporcell)) *
100## [1] 84.76845table(AS$souporcell)##
## doublet singlet0 singlet1 singlet2 singlet3 singlet4 singlet5
## 416 391 868 71 130 447 28
## singlet6 singlet7 unassigned
## 118 232 5Supplementary, hybrid demultiplexing for dataset that has high doublet rate
We have illustrated the effectiveness of hybrid demultiplexing on “high quality” data, where cell hashing and SNP-based methods are highly consistent with each other. However, there are instances where the two modalities do not align well. In such cases, our hybrid demultiplexing approach adeptly identifies significant misalignment issues via the convergence score and issues a warning to the user. Discrepancies between the two modalities often suggest a high doublet rate or a complex genetic background, necessitating careful examination by users with relevant biological expertise. To exemplify this, we use the dataset labeled ‘9pool-CA’:
DN.9pool.data <- Read10X(data.dir = "../data/9pool-CA-RNA/filtered_feature_bc_matrix/")## 10X data contains more than one type and is being returned as a list containing matrices of each type.DN.9pool.hto <- DN.9pool.data[["Antibody Capture"]][138:146, ]
DN.9pool.data[["Antibody Capture"]] <- DN.9pool.data[["Antibody Capture"]][1:137,
]
DN.9pool <- CreateSeuratObject(counts = DN.9pool.data$`Gene Expression`, project = "Carr")
DN.9pool[["HTO"]] <- CreateAssayObject(counts = DN.9pool.hto)
DN.9pool$subject <- "CA"
DN.9pool <- RenameCells(object = DN.9pool, add.cell.id = "CA")
DN.9pool[["percent.mt"]] <- PercentageFeatureSet(object = DN.9pool, pattern = "^MT-")
DN.9pool[["percent.ig"]] <- PercentageFeatureSet(object = DN.9pool, pattern = "^IG[HKL]")
DN.9pool <- subset(DN.9pool, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 &
percent.mt < 10)
DN.9pool <- NormalizeData(DN.9pool, assay = "HTO", normalization.method = "CLR")## Normalizing across featuresCell hashing demultiplexing using HTOreader
DN.9pool <- HTOClassification(DN.9pool, assay = "HTO", method = "log")## Normalization method is: log
## Cutoff for S289 is 3.35181680088473
## Cutoff for S297 is 4.30147252317663
## Cutoff for S299 is 3.97183902805626
## Cutoff for S302 is 3.83673446295452
## Cutoff for S354 is 3.90049907157204
## Cutoff for S365 is 2.71044460199379
## Cutoff for S374 is 2.86249365243517
## Cutoff for S403 is 2.65461369188085
## Cutoff for S409 is 3.21699969089809Load souporcell results. There are three runs, 1) using all cells, 2) only remove cells that expressed both B and T cell markers, 3) remove all doublet cells that identified by cell hashing
# souporcell results using all cells
sampleID <- read.table(file = "9pool/clusters.csv", header = TRUE, sep = ",", row.names = 1)
soupID <- as.character(DN.9pool@meta.data[["orig.ident"]])
i <- 1
for (barcode in names(DN.9pool@active.ident)) {
barcode1 <- gsub("^[A-Za-z]+_", "", barcode)
if (barcode1 %in% rownames(sampleID)) {
soupID[i] <- sampleID[barcode1, "identity"]
} else {
soupID[i] <- "unassigned"
}
i <- i + 1
}
DN.9pool$soupID1 <- soupID
# souporcell results on the dataset, after excluding cells that expressed both
# B and T cell markers
sampleID <- read.table(file = "9pool/clusters1.csv", header = TRUE, sep = ",", row.names = 1)
soupID <- as.character(DN.9pool@meta.data[["orig.ident"]])
i <- 1
for (barcode in names(DN.9pool@active.ident)) {
barcode1 <- gsub("^[A-Za-z]+_", "", barcode)
if (barcode1 %in% rownames(sampleID)) {
soupID[i] <- sampleID[barcode1, "identity"]
} else {
soupID[i] <- "unassigned"
}
i <- i + 1
}
DN.9pool$soupID2 <- soupID
# souporcell results on the dataset, , after excluding all doublet cells that
# identified by cell hashing
sampleID <- read.table(file = "9pool/clusters2.csv", header = TRUE, sep = ",", row.names = 1)
soupID <- as.character(DN.9pool@meta.data[["orig.ident"]])
i <- 1
for (barcode in names(DN.9pool@active.ident)) {
barcode1 <- gsub("^[A-Za-z]+_", "", barcode)
if (barcode1 %in% rownames(sampleID)) {
soupID[i] <- sampleID[barcode1, "identity"]
} else {
soupID[i] <- "unassigned"
}
i <- i + 1
}
DN.9pool$soupID3 <- soupIDWe then run HybridDemultiplexing for using all three genotype labels. For soupID1, a low convergence score was identified, and a warning was poped. HybridID was set to the cell hashing labels cause the genotype labels are not accurate.
DN.9pool <- HybridDemultiplexing(object = DN.9pool, cellhashing_label = "HTOid",
genotype_label = "soupID1", hto_names = DN.9pool.hto@Dimnames[[1]])## We matched:
## S289 with singlet8
## S297 with singlet1
## S299 with singlet4
## S302 with singlet2
## S354 with singlet6
## S365 with singlet7
## S374 with singlet3
## S403 with singlet5
## S409 with singlet0
## The convergence score between cell hashing demultiplexing and SNP-based demultiplexing is: 0.106513222331048## Warning in HybridDemultiplexing(object = DN.9pool, cellhashing_label = "HTOid",
## : Corrolation between your cell hashing demultiplexing and SNP-based
## demultiplexing is very poor. The convergence score is 0.106513222331048 , which
## is much lower than normal. In this case, we can not trust your SNP-based
## demultiplexing results. We set your cell hashing ID to the hybrid labels.
## However, a low convergance score usually indicates high doublet rates or even
## messed genetic background among different samples. We highly recommand the
## users to double check the data.For soupID2, the convergence score was improved, but still too low. HybridID was set to the cell hashing labels cause the genotype labels are not accurate.
DN.9pool <- HybridDemultiplexing(object = DN.9pool, cellhashing_label = "HTOid",
genotype_label = "soupID2", hto_names = DN.9pool.hto@Dimnames[[1]])## We matched:
## S289 with singlet1
## S297 with singlet3
## S299 with singlet5
## S302 with singlet2
## S354 with singlet4
## S365 with singlet6
## S374 with singlet8
## S403 with singlet7
## S409 with singlet0
## The convergence score between cell hashing demultiplexing and SNP-based demultiplexing is: 0.552840009051822## Warning in HybridDemultiplexing(object = DN.9pool, cellhashing_label = "HTOid",
## : Corrolation between your cell hashing demultiplexing and SNP-based
## demultiplexing is very poor. The convergence score is 0.552840009051822 , which
## is much lower than normal. In this case, we can not trust your SNP-based
## demultiplexing results. We set your cell hashing ID to the hybrid labels.
## However, a low convergance score usually indicates high doublet rates or even
## messed genetic background among different samples. We highly recommand the
## users to double check the data.For soupID3, the convergence score was highly improved. HybridID was determined by both genotype labels and cell hashing labels.
DN.9pool <- HybridDemultiplexing(object = DN.9pool, cellhashing_label = "HTOid",
genotype_label = "soupID3", hto_names = DN.9pool.hto@Dimnames[[1]])## We matched:
## S289 with singlet7
## S297 with singlet5
## S299 with singlet3
## S302 with singlet0
## S354 with singlet6
## S365 with singlet1
## S374 with singlet4
## S403 with singlet8
## S409 with singlet2
## The convergence score between cell hashing demultiplexing and SNP-based demultiplexing is: 0.956576825616217## Hybrid demultiplexing done! Results in hybridIDSession info
sessionInfo()## R version 4.2.2 (2022-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cellhashR_1.0.3 HTOreader_0.1.0 cowplot_1.1.1 SeuratObject_4.1.3
## [5] Seurat_4.3.0.1
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.16 colorspace_2.1-0 deldir_1.0-9
## [4] ellipsis_0.3.2 modeltools_0.2-23 ggridges_0.5.4
## [7] spatstat.data_3.0-1 rstudioapi_0.15.0 leiden_0.4.3
## [10] listenv_0.9.0 hash_2.2.6.2 ggrepel_0.9.3
## [13] flexmix_2.3-19 fansi_1.0.4 R.methodsS3_1.8.2
## [16] codetools_0.2-18 splines_4.2.2 cachem_1.0.8
## [19] knitr_1.43 polyclip_1.10-4 jsonlite_1.8.7
## [22] ica_1.0-3 cluster_2.1.4 R.oo_1.25.0
## [25] png_0.1-8 uwot_0.1.16 spatstat.sparse_3.0-2
## [28] shiny_1.7.4.1 sctransform_0.3.5 BiocManager_1.30.19
## [31] compiler_4.2.2 httr_1.4.6 Matrix_1.5-3
## [34] fastmap_1.1.1 lazyeval_0.2.2 cli_3.6.1
## [37] formatR_1.14 later_1.3.1 htmltools_0.5.5
## [40] tools_4.2.2 igraph_1.5.1 gtable_0.3.3
## [43] glue_1.6.2 RANN_2.6.1 reshape2_1.4.4
## [46] dplyr_1.1.2 Rcpp_1.0.11 scattermore_1.2
## [49] jquerylib_0.1.4 vctrs_0.6.3 nlme_3.1-161
## [52] spatstat.explore_3.2-1 progressr_0.13.0 lmtest_0.9-40
## [55] spatstat.random_3.1-5 xfun_0.39 stringr_1.5.0
## [58] globals_0.16.2 mime_0.12 miniUI_0.1.1.1
## [61] lifecycle_1.0.3 irlba_2.3.5.1 goftest_1.2-3
## [64] future_1.33.0 MASS_7.3-58.1 zoo_1.8-12
## [67] scales_1.2.1 spatstat.utils_3.0-3 promises_1.2.0.1
## [70] parallel_4.2.2 RColorBrewer_1.1-3 yaml_2.3.7
## [73] gridExtra_2.3 reticulate_1.30 pbapply_1.7-2
## [76] ggplot2_3.4.2 sass_0.4.6 stringi_1.7.12
## [79] rlang_1.1.1 pkgconfig_2.0.3 matrixStats_1.0.0
## [82] evaluate_0.21 lattice_0.20-45 tensor_1.5
## [85] ROCR_1.0-11 purrr_1.0.1 patchwork_1.1.2
## [88] htmlwidgets_1.6.2 tidyselect_1.2.0 parallelly_1.36.0
## [91] RcppAnnoy_0.0.21 bookdown_0.34 plyr_1.8.8
## [94] magrittr_2.0.3 R6_2.5.1 generics_0.1.3
## [97] pillar_1.9.0 fitdistrplus_1.1-11 prettydoc_0.4.1
## [100] abind_1.4-5 survival_3.5-0 sp_2.0-0
## [103] nnet_7.3-18 tibble_3.2.1 future.apply_1.11.0
## [106] KernSmooth_2.23-20 utf8_1.2.3 spatstat.geom_3.2-2
## [109] plotly_4.10.2 rmarkdown_2.23 grid_4.2.2
## [112] data.table_1.14.8 rmdformats_1.0.4 digest_0.6.33
## [115] xtable_1.8-4 tidyr_1.3.0 httpuv_1.6.11
## [118] R.utils_2.12.2 stats4_4.2.2 munsell_0.5.0
## [121] viridisLite_0.4.2 bslib_0.5.0