Spatial analysis of 10X example Visium dataset
Lambda Moses
2024-11-23
Source:vignettes/visium_10x_spatial.Rmd
visium_10x_spatial.RmdIntroduction
In a more introductory vignette, we performed basic non-spatial analyses on a mouse olfactory bulb Visium dataset from the 10X website. In this vignette, we perform spatial analyses in histological space as well as in gene expression space.
Here we load the packages used in this vignette:
library(Voyager)
library(SpatialFeatureExperiment)
library(SingleCellExperiment)
library(ggplot2)
library(scater)
library(scuttle)
library(scran)
library(stringr)
library(patchwork)
library(bluster)
library(rjson)
library(EBImage)
library(terra)
library(rlang)
library(sf)
library(rmapshaper)
library(dplyr)
library(BiocParallel)
library(BiocNeighbors)
library(reticulate)
theme_set(theme_bw())
# Specify Python version to use gget
PY_PATH <- system("which python", intern = TRUE)
use_python(PY_PATH)
py_config()<<<<<<< HEAD
Here we download the data from the 10X website. This is the unfiltered gene count matrix:
if (!file.exists("visium_ob.tar.gz"))
download.file("https://cf.10xgenomics.com/samples/spatial-exp/2.0.0/Visium_Mouse_Olfactory_Bulb/Visium_Mouse_Olfactory_Bulb_raw_feature_bc_matrix.tar.gz",
destfile = "visium_ob.tar.gz")This is the spatial information:
if (!file.exists("visium_ob_spatial.tar.gz"))
download.file("https://cf.10xgenomics.com/samples/spatial-exp/2.0.0/Visium_Mouse_Olfactory_Bulb/Visium_Mouse_Olfactory_Bulb_spatial.tar.gz",
destfile = "visium_ob_spatial.tar.gz")Decompress the downloaded content:
if (!dir.exists("outs")) {
dir.create("outs")
system("tar -xvf visium_ob.tar.gz -C outs")
system("tar -xvf visium_ob_spatial.tar.gz -C outs")
}Contents of the outs directory as from Space Ranger is
explained in the
introductory vignette.
Here we read the data into R as an SFE object.
(sfe <- read10xVisiumSFE(samples = ".", type = "sparse", data = "raw"))
#> >>> 10X Visium data will be loaded: outs
#> class: SpatialFeatureExperiment
#> dim: 32285 4992
#> metadata(0):
#> assays(1): counts
#> rownames(32285): ENSMUSG00000051951 ENSMUSG00000089699 ...
#> ENSMUSG00000095019 ENSMUSG00000095041
#> rowData names(8): symbol Feature.Type ...
#> Median.Normalized.Average.Counts_sample01
#> Barcodes.Detected.per.Feature_sample01
#> colnames(4992): AAACAACGAATAGTTC-1 AAACAAGTATCTCCCA-1 ...
#> TTGTTTGTATTACACG-1 TTGTTTGTGTAAATTC-1
#> colData names(4): in_tissue array_row array_col sample_id
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : pxl_col_in_fullres pxl_row_in_fullres
#> imgData names(4): sample_id image_id data scaleFactor
#>
#> unit: full_res_image_pixel
#> Geometries:
#> colGeometries: spotPoly (POLYGON)
#>
#> Graphs:
#> sample01:Here we add QC metrics, already plotted in the introductory vignette.
is_mt <- str_detect(rowData(sfe)$symbol, "^mt-")
sfe <- addPerCellQCMetrics(sfe, subsets = list(mito = is_mt))Tissue segmentation
While Space Ranger can automatically detect which spots are in tissue
and the Loupe browser can be used to manually annotate which spots are
in tissue, it may be interesting to get the tissue outline polygon, so
we would know how much each spot overlaps with the tissue and plot the
outline. The tissue boundary polygon can be manually annotated with
QuPath, which saves the polygon as a GeoJSON and can be directly read
into R with st_read().
Or we can segment the tissue computationally. R generally isn’t great
for image processing, but there are some packages that can perform the
segmentation, such as EBImage,
which is based on its own in house C and C++ code, and imager, which
is based on CImg.
Here we don’t have the full resolution image. We will perform tissue
segmentation on the high resolution downsampled image and then scale it
to make the coordinates of the tissue boundary match those of the spots.
The EBImage package is used here. Compared to OpenCV,
EBImage is slow on the full resolution image, but should be
fine here for the downsized image.
When rendered as a static webpage, the image is static, but when run interactively, this image will be shown in an interactive widget where you can zoom and pan.
Here we show the RGB channels separately
hist(img)
The tissue can be discerned with thresholding. The tall peak on the right is the background. The much lower peaks from around 0.6 to 0.85 must be the tissue. To capture the faint bluish region, the blue channel is used for thresholding. The threshold here is chosen based on the histogram and experimenting with nearby values.
mask <- img2[,,3] < 0.87
display(mask)
Then we use an opening operation (erosion followed by dilation) to denoise
There are some small holes in the tissue, which can be removed by a closing operation (dilation followed by erosion):
There are some larger holes in the tissue mask, which may be real holes or faint regions with few nuclei missed by thresholding. They might not be large enough to affect which Visium spots intersect the tissue.
Now the main piece of tissue is clear. It must be the object with the
largest area. However, there are two small pieces that should belong to
the tissue at the top left. The debris and fiducials can be removed by
setting all pixels in the mask outside the bounding box of the main
piece to 0. Here we assign a different value to each contiguous object
with bwlabel(), and use
computeFeatures.shape() to find the area among other shape
features (e.g. perimeter) of each object.
mask_label <- bwlabel(mask_close)
fts <- computeFeatures.shape(mask_label)
head(fts)
#> s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
#> 1 39 25 3.428773 1.3542219 1.4176036 5.762777
#> 2 20 14 2.032665 0.3439068 1.5000000 2.500000
#> 3 8 8 1.144123 0.4370160 0.7071068 1.581139
#> 4 14 10 1.689175 0.2160726 1.5811388 2.121320
#> 5 15 12 1.716761 0.4684015 1.0000000 2.236068
#> 6 9 8 1.207107 0.2071068 1.0000000 1.414214
summary(fts[,"s.area"])
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 8.0 14.0 51.0 595.9 345.0 496326.0
max_ind <- which.max(fts[,"s.area"])
inds <- which(as.array(mask_label) == max_ind, arr.ind = TRUE)
head(inds)
#> row col
#> [1,] 1168 562
#> [2,] 1169 562
#> [3,] 1170 562
#> [4,] 1158 563
#> [5,] 1159 563
#> [6,] 1160 563
row_inds <- c(seq_len(min(inds[,1])-1), seq(max(inds[,1])+1, nrow(mask_label), by = 1))
col_inds <- c(seq_len(min(inds[,2])-1), seq(max(inds[,2])+1, nrow(mask_label), by = 1))
mask_label[row_inds, ] <- 0
mask_label[,col_inds] <- 0
display(mask_label)
Then remove the small pieces that are debris.
unique(as.vector(mask_label))
#> [1] 0 421 425 429 430 438 450 458 461 469 473 483 487 505 523 633 640 642 651
#> [20] 678 739 741 757 762 775 778 789 791 797 805 810 813 820 821 822 826 831 838
#> [39] 839 840 843 845 848 849 861 862 863
fts2 <- fts[unique(as.vector(mask_label))[-1],]
fts2 <- fts2[order(fts2[,"s.area"], decreasing = TRUE),]
plot(fts2[,1][-1], type = "l", ylab = "Area")
head(fts2, 10)
#> s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
#> 421 496326 3151 395.118732 68.6493949 234.1605637 485.715835
#> 450 217 55 7.840627 1.2883202 5.0010247 10.458197
#> 849 211 63 7.961248 2.0228566 3.8569142 12.189753
#> 797 182 56 7.547362 2.3368359 3.0772370 11.839919
#> 461 136 54 7.186020 3.2751628 0.9255555 12.479805
#> 741 92 56 6.365661 2.8382829 1.3273276 11.653219
#> 840 69 33 4.503526 1.6417370 1.6026264 7.076974
#> 862 63 37 4.854424 2.4445530 0.6361407 8.837838
#> 839 45 25 3.305562 0.7074306 1.9320455 4.526897
#> 775 32 22 2.887407 1.1755375 0.5300865 4.543636Object number 797 is a piece of debris at the bottom left. The other pieces with area over 100 pixels are tissue. Since debris really is small bits of tissue, so the boundary between debris and tissue can be blurry. Here the two are distinguished by morphology on the H&E image and proximity to the main tissue.
#display(mask_label == 797)Here we remove the debris from the mask
mask_label[mask_label %in% c(797, as.numeric(rownames(fts2)[fts2[,1] < 100]))] <- 0Since most holes in the mask are faint regions of the tissue missed by thresholding, the holes will be filled
mask_label <- fillHull(mask_label)
display(paintObjects(mask_label, img, col=c("red", "yellow"), opac=c(1, 0.3)))
This segmentation process took a lot of manual oversight, in choosing the threshold, choosing kernel size and shape in the opening and closing operations, deciding whether to fill the holes, and deciding what is debris and what is tissue.
Convert tissue mask to polygon
Now we have the tissue mask, which we will convert to polygon. While
OpenCV can directly perform the conversion, as there isn’t a
comprehensive R wrapper of OpenCV, this conversion is more convoluted in
R. We first convert the Image object to a raster as
implemented in terra, the core R package for geospatial
raster data. Then terra can convert the raster to polygon.
As this image is downsized, the polygon will look quite pixelated. To
mitigate this pixelation and save memory, the ms_simplify()
function is used to simplify the polygon, only keeping a small
proportion of all vertices. The st_simplify() function in
sf can also simplify the polygons, but we can’t specify
what proportion of vertices to keep.
raster2polygon <- function(seg, keep = 0.2) {
r <- rast(as.array(seg), extent = ext(0, nrow(seg), 0, ncol(seg))) |>
trans() |> flip()
r[r < 1] <- NA
contours <- st_as_sf(as.polygons(r, dissolve = TRUE))
simplified <- ms_simplify(contours, keep = keep)
list(full = contours,
simplified = simplified)
}
tb <- raster2polygon(mask_label)Before adding the geometry to the SFE object, it needs to be scaled to match the coordinates of the spots
scale_factors <- fromJSON(file = "outs/spatial/scalefactors_json.json")
tb$simplified$geometry <- tb$simplified$geometry / scale_factors$tissue_hires_scalef
tissueBoundary(sfe) <- tb$simplified
plotSpatialFeature(sfe, "sum", annotGeometryName = "tissueBoundary",
annot_fixed = list(fill = NA, color = "black"),
image_id = "lowres") +
theme_void()
The mouse olfactory bulb is conventionally plotted horizontally. The entire SFE object can be transposed in histologial space to make the olfactory bulb horizontal.
sfe <- SpatialFeatureExperiment::transpose(sfe)
plotSpatialFeature(sfe, "sum", annotGeometryName = "tissueBoundary",
annot_fixed = list(fill = NA, color = "black"),
image_id = "lowres")
Then we can use geometric operations to find which spots intersect tissue, which spots are covered by tissue, and how much of each spot intersects tissue.
# Which spots intersect tissue
sfe$int_tissue <- annotPred(sfe, colGeometryName = "spotPoly",
annotGeometryName = "tissueBoundary",
pred = st_intersects)
sfe$cov_tissue <- annotPred(sfe, colGeometryName = "spotPoly",
annotGeometryName = "tissueBoundary",
pred = st_covered_by)Discrepancies between Space Ranger’s annotation and the annotation based on tissue segmentation here:
sfe$diff_sr <- case_when(sfe$in_tissue == sfe$int_tissue ~ "same",
sfe$in_tissue & !sfe$int_tissue ~ "Space Ranger",
sfe$int_tissue & !sfe$in_tissue ~ "segmentation") |>
factor(levels = c("Space Ranger", "same", "segmentation"))
plotSpatialFeature(sfe, "diff_sr",
annotGeometryName = "tissueBoundary",
annot_fixed = list(fill = NA, size = 0.5, color = "black")) +
scale_fill_brewer(type = "div", palette = 4)
#> Scale for fill is already present.
#> Adding another scale for fill, which will replace the existing scale.
Spots at the margin can intersect the tissue without being covered by it.
sfe$diff_int_cov <- sfe$int_tissue != sfe$cov_tissue
plotSpatialFeature(sfe, "diff_int_cov",
annotGeometryName = "tissueBoundary",
annot_fixed = list(fill = NA, size = 0.5, color = "black"))
We can also get the geometries of the intersections between the tissue and the Visium spots, and then calculate what percentage of each spot is in tissue. However, this percentage may not be very useful if the tissue segmentation is subject to error. This percentage may be more useful for pathologist annotated histological regions or objects such as nuclei and myofibers.
spot_ints <- annotOp(sfe, colGeometryName = "spotPoly",
annotGeometryName = "tissueBoundary", op = st_intersection)
sfe$pct_tissue <- st_area(spot_ints) / st_area(spotPoly(sfe)) * 100For spots that intersect tissue, does total counts relate to percentage of the spot in tissue?
sfe_tissue <- sfe[,sfe$int_tissue]
plotColData(sfe_tissue, x = "pct_tissue", y = "sum", color_by = "diff_int_cov")
Spots that are not fully covered by tissue have lower total UMI counts, which can be due to both that they are not fully in tissue and the cell types with lower total counts in the histological region near the edge, as some spots fully covered by tissue also have low UMI counts.
Spatial autocorrelation of QC metrics
colGraph(sfe_tissue, "visium") <- findVisiumGraph(sfe_tissue)
qc_features <- c("sum", "detected", "subsets_mito_percent")
sfe_tissue <- colDataUnivariate(sfe_tissue, "moran.mc", qc_features, nsim = 200)
plotMoranMC(sfe_tissue, qc_features)
sfe_tissue <- colDataUnivariate(sfe_tissue, "sp.correlogram", qc_features,
order = 8)
plotCorrelogram(sfe_tissue, qc_features)
sfe_tissue <- colDataUnivariate(sfe_tissue, "localmoran", qc_features)
plotLocalResult(sfe_tissue, "localmoran", qc_features, ncol = 2,
colGeometryName = "spotPoly", divergent = TRUE,
diverge_center = 0, image_id = "lowres", maxcell = 5e4)
sfe_tissue <- colDataUnivariate(sfe_tissue, "LOSH", qc_features)
plotLocalResult(sfe_tissue, "LOSH", qc_features, ncol = 2,
colGeometryName = "spotPoly", image_id = "lowres", maxcell = 5e4)
sfe_tissue <- colDataUnivariate(sfe_tissue, "moran.plot", qc_features)
moranPlot(sfe_tissue, "subsets_mito_percent")
Spatial autocorrelation of gene expression
Normalize data with the scran method, and find highly
variable genes
#clusters <- quickCluster(sfe_tissue)
#sfe_tissue <- computeSumFactors(sfe_tissue, clusters=clusters)
#sfe_tissue <- sfe_tissue[, sizeFactors(sfe_tissue) > 0]
sfe_tissue <- logNormCounts(sfe_tissue)
dec <- modelGeneVar(sfe_tissue)
hvgs <- getTopHVGs(dec, n = 2000)Find Moran’s I for all highly variable genes:
sfe_tissue <- runMoransI(sfe_tissue, features = hvgs, BPPARAM = MulticoreParam(2))
plotRowDataHistogram(sfe_tissue, "moran_sample01")
#> Warning: Removed 30285 rows containing non-finite outside the scale range
#> (`stat_bin()`).
The vast majority of genes have positive Moran’s I. Here we’ll find the genes with the highest Moran’s I:
top_moran <- rownames(sfe_tissue)[order(rowData(sfe_tissue)$moran_sample01,
decreasing = TRUE)[1:9]]We can use the gget info module from the gget package to get additional information on these genes, such as their descriptions, synonyms, transcripts and more from a collection of reference databases including Ensembl, UniProt and NCBI Here, we are showing their gene descriptions from NCBI: <<<<<<< HEAD
>>>>>>> documentation
gget_info <- gget$info(top_moran)
rownames(gget_info) <- gget_info$primary_gene_name
select(gget_info, ncbi_description)Plot the genes with the highest Moran’s I:
plotSpatialFeature(sfe_tissue, top_moran, ncol = 3, image_id = "lowres",
maxcell = 5e4, swap_rownames = "symbol")
Here global Moran’s I seems to be more about tissue structure.
Some genes have negative Moran’s I that might not be statistically significant:
neg_moran <- rownames(sfe_tissue)[order(rowData(sfe_tissue)$moran_sample01,
decreasing = FALSE)[1:9]]<<<<<<< HEAD
>>>>>>> documentation
# Display NCBI descriptions for these genes
gget_info_neg <- gget$info(neg_moran)
rownames(gget_info_neg) <- gget_info_neg$primary_gene_name
select(gget_info_neg, ncbi_description)
plotSpatialFeature(sfe_tissue, neg_moran, ncol = 3, swap_rownames = "symbol",
image_id = "lowres", maxcell = 5e4)
sfe_tissue <- runUnivariate(sfe_tissue, "moran.mc", neg_moran,
colGraphName = "visium", nsim = 200, alternative = "less")
plotMoranMC(sfe_tissue, neg_moran, swap_rownames = "symbol")
rowData(sfe_tissue)[neg_moran, c("moran_sample01", "moran.mc_p.value_sample01")]
#> DataFrame with 9 rows and 2 columns
#> moran_sample01 moran.mc_p.value_sample01
#> <numeric> <numeric>
#> ENSMUSG00000041426 -0.0531915 0.00497512
#> ENSMUSG00000048277 -0.0451179 0.00497512
#> ENSMUSG00000021236 -0.0445148 0.01492537
#> ENSMUSG00000025241 -0.0419121 0.01492537
#> ENSMUSG00000026126 -0.0399917 0.00497512
#> ENSMUSG00000034342 -0.0393964 0.00995025
#> ENSMUSG00000047205 -0.0381599 0.00995025
#> ENSMUSG00000031431 -0.0369456 0.01492537
#> ENSMUSG00000037552 -0.0368969 0.01990050As there are 2000 highly variable genes and 2000 tests, these would no longer be significant after correcting for multiple testing.
Does global Moran’s I relate to gene expression level?
sfe_tissue <- addPerFeatureQCMetrics(sfe_tissue)
names(rowData(sfe_tissue))
#> [1] "symbol"
#> [2] "Feature.Type"
#> [3] "I_sample01"
#> [4] "P.value_sample01"
#> [5] "Adjusted.p.value_sample01"
#> [6] "Feature.Counts.in.Spots.Under.Tissue_sample01"
#> [7] "Median.Normalized.Average.Counts_sample01"
#> [8] "Barcodes.Detected.per.Feature_sample01"
#> [9] "subsets_mito"
#> [10] "moran_sample01"
#> [11] "K_sample01"
#> [12] "moran.mc_statistic_sample01"
#> [13] "moran.mc_parameter_sample01"
#> [14] "moran.mc_p.value_sample01"
#> [15] "moran.mc_alternative_sample01"
#> [16] "moran.mc_method_sample01"
#> [17] "moran.mc_res_sample01"
#> [18] "mean"
#> [19] "detected"
plotRowData(sfe_tissue, x = "mean", y = "moran_sample01") +
scale_x_log10() +
annotation_logticks(sides = "b") +
geom_density2d()
#> Warning in scale_x_log10(): log-10 transformation introduced infinite values.
#> log-10 transformation introduced infinite values.
#> Warning: Removed 30285 rows containing non-finite outside the scale range
#> (`stat_density2d()`).
#> Warning: Removed 30285 rows containing missing values or values outside the scale range
#> (`geom_point()`).
Genes that are more highly expressed overall tend to have higher Moran’s I.
Apply spatial analysis methods to gene expression space
Spatial statistics that require a spatial neighborhood graph can also be applied to the k nearest neighbor graph not in histological space but in gene expression space. This is done in more depth in this vignette.
sfe_tissue <- runPCA(sfe_tissue, ncomponents = 30, subset_row = hvgs,
scale = TRUE) # scale as in Seurat
foo <- findKNN(reducedDim(sfe_tissue, "PCA")[,1:10], k=10, BNPARAM=AnnoyParam())
# Split by row
foo_nb <- asplit(foo$index, 1)
dmat <- 1/foo$distance
# Row normalize the weights
dmat <- sweep(dmat, 1, rowSums(dmat), FUN = "/")
glist <- asplit(dmat, 1)
# Sort based on index
ord <- lapply(foo_nb, order)
foo_nb <- lapply(seq_along(foo_nb), function(i) foo_nb[[i]][ord[[i]]])
class(foo_nb) <- "nb"
glist <- lapply(seq_along(glist), function(i) glist[[i]][ord[[i]]])
listw <- list(style = "W",
neighbours = foo_nb,
weights = glist)
class(listw) <- "listw"
attr(listw, "region.id") <- colnames(sfe_tissue)
colGraph(sfe_tissue, "knn10") <- listw
sfe_tissue <- runMoransI(sfe_tissue, features = hvgs, BPPARAM = MulticoreParam(2),
colGraphName = "knn10", name = "moran_ns")Here we store the results in “moran_ns”, not to be confused with spatial Moran’s I results.
These are the genes that tend to be more similar to their neighbors in the 10 nearest neighbor graph in PCA space for gene expression rather than in histological space:
top_moran2 <- rownames(sfe_tissue)[order(rowData(sfe_tissue)$moran_ns_sample01,
decreasing = TRUE)[1:9]]<<<<<<< HEAD
>>>>>>> documentation
# Display NCBI descriptions for these genes
gget_info2 <- gget$info(top_moran2)
rownames(gget_info2) <- gget_info2$primary_gene_name
select(gget_info2, ncbi_description)
plotSpatialFeature(sfe_tissue, top_moran2, ncol = 3, swap_rownames = "symbol",
image_id = "lowres", maxcell = 5e4)
Although this Moran’s I was not computed in histological space, these genes with the highest Moran’s I in PCA space also show spatial structure, as different cell types reside in different spatial regions.
Session info
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] reticulate_1.40.0 BiocNeighbors_2.0.0
#> [3] BiocParallel_1.40.0 dplyr_1.1.4
#> [5] rmapshaper_0.5.0 sf_1.0-19
#> [7] rlang_1.1.4 terra_1.7-83
#> [9] EBImage_4.48.0 rjson_0.2.23
#> [11] bluster_1.16.0 patchwork_1.3.0
#> [13] stringr_1.5.1 scran_1.34.0
#> [15] scater_1.34.0 scuttle_1.16.0
#> [17] ggplot2_3.5.1 SingleCellExperiment_1.28.1
#> [19] SummarizedExperiment_1.36.0 Biobase_2.66.0
#> [21] GenomicRanges_1.58.0 GenomeInfoDb_1.42.0
#> [23] IRanges_2.40.0 S4Vectors_0.44.0
#> [25] BiocGenerics_0.52.0 MatrixGenerics_1.18.0
#> [27] matrixStats_1.4.1 Voyager_1.8.1
#> [29] SpatialFeatureExperiment_1.9.4
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.4.2 bitops_1.0-9
#> [3] tibble_3.2.1 R.oo_1.27.0
#> [5] lifecycle_1.0.4 edgeR_4.4.0
#> [7] lattice_0.22-6 MASS_7.3-61
#> [9] magrittr_2.0.3 limma_3.62.1
#> [11] sass_0.4.9 rmarkdown_2.29
#> [13] jquerylib_0.1.4 yaml_2.3.10
#> [15] metapod_1.14.0 sp_2.1-4
#> [17] cowplot_1.1.3 DBI_1.2.3
#> [19] RColorBrewer_1.1-3 multcomp_1.4-26
#> [21] abind_1.4-8 spatialreg_1.3-5
#> [23] zlibbioc_1.52.0 R.utils_2.12.3
#> [25] RCurl_1.98-1.16 TH.data_1.1-2
#> [27] sandwich_3.1-1 GenomeInfoDbData_1.2.13
#> [29] ggrepel_0.9.6 irlba_2.3.5.1
#> [31] units_0.8-5 RSpectra_0.16-2
#> [33] dqrng_0.4.1 pkgdown_2.1.1
#> [35] DelayedMatrixStats_1.28.0 codetools_0.2-20
#> [37] DropletUtils_1.26.0 DelayedArray_0.32.0
#> [39] tidyselect_1.2.1 UCSC.utils_1.2.0
#> [41] memuse_4.2-3 farver_2.1.2
#> [43] ScaledMatrix_1.14.0 viridis_0.6.5
#> [45] jsonlite_1.8.9 geojsonsf_2.0.3
#> [47] e1071_1.7-16 survival_3.7-0
#> [49] systemfonts_1.1.0 dbscan_1.2-0
#> [51] tools_4.4.2 ggnewscale_0.5.0
#> [53] ragg_1.3.3 Rcpp_1.0.13-1
#> [55] glue_1.8.0 gridExtra_2.3
#> [57] SparseArray_1.6.0 mgcv_1.9-1
#> [59] xfun_0.49 HDF5Array_1.34.0
#> [61] withr_3.0.2 fastmap_1.2.0
#> [63] boot_1.3-31 rhdf5filters_1.18.0
#> [65] fansi_1.0.6 spData_2.3.3
#> [67] digest_0.6.37 rsvd_1.0.5
#> [69] R6_2.5.1 textshaping_0.4.0
#> [71] colorspace_2.1-1 wk_0.9.4
#> [73] LearnBayes_2.15.1 jpeg_0.1-10
#> [75] R.methodsS3_1.8.2 utf8_1.2.4
#> [77] generics_0.1.3 data.table_1.16.2
#> [79] class_7.3-22 httr_1.4.7
#> [81] htmlwidgets_1.6.4 S4Arrays_1.6.0
#> [83] spdep_1.3-6 pkgconfig_2.0.3
#> [85] scico_1.5.0 gtable_0.3.6
#> [87] XVector_0.46.0 htmltools_0.5.8.1
#> [89] fftwtools_0.9-11 scales_1.3.0
#> [91] png_0.1-8 SpatialExperiment_1.16.0
#> [93] knitr_1.49 coda_0.19-4.1
#> [95] nlme_3.1-166 curl_6.0.1
#> [97] proxy_0.4-27 cachem_1.1.0
#> [99] zoo_1.8-12 rhdf5_2.50.0
#> [101] KernSmooth_2.23-24 parallel_4.4.2
#> [103] vipor_0.4.7 desc_1.4.3
#> [105] s2_1.1.7 pillar_1.9.0
#> [107] grid_4.4.2 vctrs_0.6.5
#> [109] BiocSingular_1.22.0 beachmat_2.22.0
#> [111] sfheaders_0.4.4 cluster_2.1.6
#> [113] beeswarm_0.4.0 evaluate_1.0.1
#> [115] isoband_0.2.7 zeallot_0.1.0
#> [117] magick_2.8.5 mvtnorm_1.3-2
#> [119] cli_3.6.3 locfit_1.5-9.10
#> [121] compiler_4.4.2 crayon_1.5.3
#> [123] labeling_0.4.3 classInt_0.4-10
#> [125] fs_1.6.5 ggbeeswarm_0.7.2
#> [127] stringi_1.8.4 viridisLite_0.4.2
#> [129] deldir_2.0-4 munsell_0.5.1
#> [131] tiff_0.1-12 V8_6.0.0
#> [133] Matrix_1.7-1 sparseMatrixStats_1.18.0
#> [135] Rhdf5lib_1.28.0 statmod_1.5.0
#> [137] igraph_2.1.1 bslib_0.8.0