Xenium breast cancer dataset
Lambda Moses
2024-05-17
Source:vignettes/vig5_xenium.Rmd
vig5_xenium.RmdIntroduction
Xenium is a new technology from 10X genomics for single cell resolution smFISH based spatial transcriptomics. The first Xenium dataset is for formalin fixed paraffin embedded (FFPE) human breast tumor, reported in (Janesick et al. 2022) and downloaded from the 10X website.
The gene count matrix was downloaded as an HDF5 file and read into R
as a SingleCellExperiment (SCE) object with
DropletUtils::read10xCounts(). The gene count matrix is
originally a DelayedArray, so the data is not all loaded
into memory. For now, the matrix has been converted into an in memory
dgCMatrix. However, for the next release, we would like to
write another vignette on on disk analyses. The challenge is
representing sf data frames on disk, perhaps with sedona
and SQLDataFrame.
The cell metadata (including centroid coordinates) and cell
segmentation polygons were downloaded as parquet files, a
more compact way to store columnar data than CSV, and read into R as
data frames with read_parquet in the arrow package. The cell polygons
were converted into sf data frame with
SpatialFeatureExperiment::df2sf(). Then the SCE object was
converted into SpatialFeatureExperiment (SFE) and the
polygon geometry was added to the SFE object, which is in the
SFEData package.
Here we load the packages used in this vignette.
library(Voyager)
library(SFEData)
library(SingleCellExperiment)
library(SpatialExperiment)
library(SpatialFeatureExperiment)
library(ggplot2)
library(stringr)
library(scater)
library(scuttle)
library(BiocParallel)
library(BiocSingular)
library(bluster)
library(scran)
library(patchwork)
theme_set(theme_bw())
(sfe <- JanesickBreastData(dataset = "rep2"))
#> see ?SFEData and browseVignettes('SFEData') for documentation
#> downloading 1 resources
#> retrieving 1 resource
#> loading from cache
#> class: SpatialFeatureExperiment
#> dim: 541 118708
#> metadata(1): Samples
#> assays(1): counts
#> rownames(541): ABCC11 ACTA2 ... BLANK_0497 BLANK_0499
#> rowData names(6): ID Symbol ... vars cv2
#> colnames: NULL
#> colData names(10): Sample Barcode ... nCounts nGenes
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x_centroid y_centroid
#> imgData names(1): sample_id
#>
#> unit:
#> Geometries:
#> colGeometries: centroids (POINT), cellSeg (POLYGON), nucSeg (GEOMETRY)
#>
#> Graphs:
#> sample01:There are 118708 cells in this dataset, a little more than in the CosMX dataset.
The SFE object doesn’t have column names (i.e. cell IDs). Here we assign cell IDs.
This is what the tissue, with the cell outlines, looks like
plotGeometry(sfe, "cellSeg")
Plot cell density in space
plotCellBin2D(sfe, hex = TRUE)
Quality control
Cells
Some QC metrics are precomputed and are stored in
colData
names(colData(sfe))
#> [1] "Sample" "Barcode"
#> [3] "transcript_counts" "control_probe_counts"
#> [5] "control_codeword_counts" "cell_area"
#> [7] "nucleus_area" "sample_id"
#> [9] "nCounts" "nGenes"Since there’re more cells, it would be better to plot the tissue larger, so we’ll plot the histogram of QC metrics and the spatial plots separately, unlike in the CosMx vignette.
n_panel <- 313
colData(sfe)$nCounts_normed <- sfe$nCounts/n_panel
colData(sfe)$nGenes_normed <- sfe$nGenes/n_panelHere we divided nCounts by the total number of genes probed, so this histogram is comparable to those from other smFISH-based datasets.
plotColDataHistogram(sfe, c("nCounts_normed", "nGenes_normed"))
Compared to the FFPE CosMX non-small cell lung cancer dataset, more transcripts per gene on average and a larger proportion of all genes are detected in this dataset, which is also FFPE. However, this should be interpreted with care, since these two datasets are from different tissues and have different gene panels, so this may or may not indicate that Xenium has better detection efficiency than CosMX.
plotSpatialFeature(sfe, "nCounts", colGeometryName = "cellSeg")
There seem to be FOV artifacts. However, the cell ID and FOV information were unavailable so we cannot examine them.
plotSpatialFeature(sfe, "nGenes", colGeometryName = "cellSeg")
A standard examination is to look at the relationship between nCounts and nGenes:
plotColData(sfe, x="nCounts", y="nGenes", bins = 100)
There appear to be two branches.
Here we plot the distribution of cell area
plotColDataHistogram(sfe, c("cell_area", "nucleus_area"), scales = "free_y")
That should be in pixels. There’s a very long tail. The nuclei are much smaller than the cells.
How is cell area distributed in space?
plotSpatialFeature(sfe, "cell_area", colGeometryName = "cellSeg")
Cells in the sparse region tend to be larger than those in the dense region. This may be biological or an artifact of the cell segmentation algorithm or both.
Here the nuclei segmentations are plotted instead of cell segmentation. The nuclei are much smaller to the extent that they are difficult to see.
plotSpatialFeature(sfe, "nucleus_area", colGeometryName = "nucSeg")
There’s an outlier near the right edge of the section, throwing off the dynamic range of the plot. Upon inspection of the H&E image, the outlier is a bit of tissue debris that doesn’t look like a cell. But we can still that cells in the dense, gland like regions tend to have larger nuclei. This may be biological, or that nuclei are so densely packed in those regions that they are more likely to be undersegmented, i.e. when multiple nuclei are counted as one by the nuclei segmentation program, or both.
These observations motivate an examination of the relationship between cell area and nuclei area:
plotColData(sfe, x="cell_area", y="nucleus_area", bins = 100)
Again, there are two branches, probably related to cell density and cell type. The nucleus outlier also has large cell area, though it is not as much an outlier in cell area. However, it is a spatial outlier as it’s unusually large compared to its neighbors (scroll up two plots back).
Next we calculate the proportion of cell in this z-plane taken up by the nucleus, and examine the distribution:
colData(sfe)$prop_nuc <- sfe$nucleus_area / sfe$cell_area
plotColDataHistogram(sfe, "prop_nuc")
This distribution could have been generated from two peaks that were combined. From the histogram, there do not seem to be cells without nuclei or segmentation artifacts where the nucleus is larger than the cell. However, there are so many cells in this dataset and it is possible that just a few cells would not be visible on this histogram. We double check:
# No nucleus
sum(sfe$nucleus_area < 1)
#> [1] 0
# Nucleus larger than cell
sum(sfe$nucleus_area > sfe$cell_area)
#> [1] 0So there are no cells without nuclei or nuclei larger than their cells. Here we plot the nuclei proportion in space:
plotSpatialFeature(sfe, "prop_nuc", colGeometryName = "cellSeg")
Cells in some histological regions have larger proportions occupied by the nuclei. It is interesting to check, controlling for cell type, how cell area, nucleus area, and the proportion of cell occupied by nucleus relate to gene expression. However, a problem in performing such an analysis is that cell segmentation is only available for one z-plane here and these areas also relate to where this z-plane intersects each cell.
Below we plot a 2D histogram to better show the density of points on this plot:
plotColData(sfe, x="cell_area", y="prop_nuc")
Smaller cells tend to have higher proportion occupied by the nucleus. This can be related to cell type, or it could be a limitation in how small the nuclei can be in this tissue.
We also examine the relationship between nucleus area and the proportion of cell occupied by the nucleus:
plotColData(sfe, x="nucleus_area", y="prop_nuc", bins = 100)
The outlier is obvious. There are more cells with both small nuclei and low proportion of area occupied by the nucleus.
Negative controls
Since there are only a few hundred genes plus negative control probes, all row names of the SFE object can be printed out to find what the negative control probes are called.
rownames(sfe)
#> [1] "ABCC11" "ACTA2"
#> [3] "ACTG2" "ADAM9"
#> [5] "ADGRE5" "ADH1B"
#> [7] "ADIPOQ" "AGR3"
#> [9] "AHSP" "AIF1"
#> [11] "AKR1C1" "AKR1C3"
#> [13] "ALDH1A3" "ANGPT2"
#> [15] "ANKRD28" "ANKRD29"
#> [17] "ANKRD30A" "APOBEC3A"
#> [19] "APOBEC3B" "APOC1"
#> [21] "AQP1" "AQP3"
#> [23] "AR" "AVPR1A"
#> [25] "BACE2" "BANK1"
#> [27] "BASP1" "BTNL9"
#> [29] "C15orf48" "C1QA"
#> [31] "C1QC" "C2orf42"
#> [33] "C5orf46" "C6orf132"
#> [35] "CAV1" "CAVIN2"
#> [37] "CCDC6" "CCDC80"
#> [39] "CCL20" "CCL5"
#> [41] "CCL8" "CCND1"
#> [43] "CCPG1" "CCR7"
#> [45] "CD14" "CD163"
#> [47] "CD19" "CD1C"
#> [49] "CD247" "CD27"
#> [51] "CD274" "CD3D"
#> [53] "CD3E" "CD3G"
#> [55] "CD4" "CD68"
#> [57] "CD69" "CD79A"
#> [59] "CD79B" "CD80"
#> [61] "CD83" "CD86"
#> [63] "CD8A" "CD8B"
#> [65] "CD9" "CD93"
#> [67] "CDC42EP1" "CDH1"
#> [69] "CEACAM6" "CEACAM8"
#> [71] "CENPF" "CLCA2"
#> [73] "CLDN4" "CLDN5"
#> [75] "CLEC14A" "CLEC9A"
#> [77] "CLECL1" "CLIC6"
#> [79] "CPA3" "CRHBP"
#> [81] "CRISPLD2" "CSF3"
#> [83] "CTH" "CTLA4"
#> [85] "CTSG" "CTTN"
#> [87] "CX3CR1" "CXCL12"
#> [89] "CXCL16" "CXCL5"
#> [91] "CXCR4" "CYP1A1"
#> [93] "CYTIP" "DAPK3"
#> [95] "DERL3" "DMKN"
#> [97] "DNAAF1" "DNTTIP1"
#> [99] "DPT" "DSC2"
#> [101] "DSP" "DST"
#> [103] "DUSP2" "DUSP5"
#> [105] "EDN1" "EDNRB"
#> [107] "EGFL7" "EGFR"
#> [109] "EIF4EBP1" "ELF3"
#> [111] "ELF5" "ENAH"
#> [113] "EPCAM" "ERBB2"
#> [115] "ERN1" "ESM1"
#> [117] "ESR1" "FAM107B"
#> [119] "FAM49A" "FASN"
#> [121] "FBLIM1" "FBLN1"
#> [123] "FCER1A" "FCER1G"
#> [125] "FCGR3A" "FGL2"
#> [127] "FLNB" "FOXA1"
#> [129] "FOXC2" "FOXP3"
#> [131] "FSTL3" "GATA3"
#> [133] "GJB2" "GLIPR1"
#> [135] "GNLY" "GPR183"
#> [137] "GZMA" "GZMB"
#> [139] "GZMK" "HAVCR2"
#> [141] "HDC" "HMGA1"
#> [143] "HOOK2" "HOXD8"
#> [145] "HOXD9" "HPX"
#> [147] "IGF1" "IGSF6"
#> [149] "IL2RA" "IL2RG"
#> [151] "IL3RA" "IL7R"
#> [153] "ITGAM" "ITGAX"
#> [155] "ITM2C" "JUP"
#> [157] "KARS" "KDR"
#> [159] "KIT" "KLF5"
#> [161] "KLRB1" "KLRC1"
#> [163] "KLRD1" "KLRF1"
#> [165] "KRT14" "KRT15"
#> [167] "KRT16" "KRT23"
#> [169] "KRT5" "KRT6B"
#> [171] "KRT7" "KRT8"
#> [173] "LAG3" "LARS"
#> [175] "LDHB" "LEP"
#> [177] "LGALSL" "LIF"
#> [179] "LILRA4" "LPL"
#> [181] "LPXN" "LRRC15"
#> [183] "LTB" "LUM"
#> [185] "LY86" "LYPD3"
#> [187] "LYZ" "MAP3K8"
#> [189] "MDM2" "MEDAG"
#> [191] "MKI67" "MLPH"
#> [193] "MMP1" "MMP12"
#> [195] "MMP2" "MMRN2"
#> [197] "MNDA" "MPO"
#> [199] "MRC1" "MS4A1"
#> [201] "MUC6" "MYBPC1"
#> [203] "MYH11" "MYLK"
#> [205] "MYO5B" "MZB1"
#> [207] "NARS" "NCAM1"
#> [209] "NDUFA4L2" "NKG7"
#> [211] "NOSTRIN" "NPM3"
#> [213] "OCIAD2" "OPRPN"
#> [215] "OXTR" "PCLAF"
#> [217] "PCOLCE" "PDCD1"
#> [219] "PDCD1LG2" "PDE4A"
#> [221] "PDGFRA" "PDGFRB"
#> [223] "PDK4" "PECAM1"
#> [225] "PELI1" "PGR"
#> [227] "PIGR" "PIM1"
#> [229] "PLD4" "POLR2J3"
#> [231] "POSTN" "PPARG"
#> [233] "PRDM1" "PRF1"
#> [235] "PTGDS" "PTN"
#> [237] "PTPRC" "PTRHD1"
#> [239] "QARS" "RAB30"
#> [241] "RAMP2" "RAPGEF3"
#> [243] "REXO4" "RHOH"
#> [245] "RORC" "RTKN2"
#> [247] "RUNX1" "S100A14"
#> [249] "S100A4" "S100A8"
#> [251] "SCD" "SCGB2A1"
#> [253] "SDC4" "SEC11C"
#> [255] "SEC24A" "SELL"
#> [257] "SERHL2" "SERPINA3"
#> [259] "SERPINB9" "SFRP1"
#> [261] "SFRP4" "SH3YL1"
#> [263] "SLAMF1" "SLAMF7"
#> [265] "SLC25A37" "SLC4A1"
#> [267] "SLC5A6" "SMAP2"
#> [269] "SMS" "SNAI1"
#> [271] "SOX17" "SOX18"
#> [273] "SPIB" "SQLE"
#> [275] "SRPK1" "SSTR2"
#> [277] "STC1" "SVIL"
#> [279] "TAC1" "TACSTD2"
#> [281] "TCEAL7" "TCF15"
#> [283] "TCF4" "TCF7"
#> [285] "TCIM" "TCL1A"
#> [287] "TENT5C" "TFAP2A"
#> [289] "THAP2" "TIFA"
#> [291] "TIGIT" "TIMP4"
#> [293] "TMEM147" "TNFRSF17"
#> [295] "TOMM7" "TOP2A"
#> [297] "TPD52" "TPSAB1"
#> [299] "TRAC" "TRAF4"
#> [301] "TRAPPC3" "TRIB1"
#> [303] "TUBA4A" "TUBB2B"
#> [305] "TYROBP" "UCP1"
#> [307] "USP53" "VOPP1"
#> [309] "VWF" "WARS"
#> [311] "ZEB1" "ZEB2"
#> [313] "ZNF562" "NegControlProbe_00042"
#> [315] "NegControlProbe_00041" "NegControlProbe_00039"
#> [317] "NegControlProbe_00035" "NegControlProbe_00034"
#> [319] "NegControlProbe_00033" "NegControlProbe_00031"
#> [321] "NegControlProbe_00025" "NegControlProbe_00024"
#> [323] "NegControlProbe_00022" "NegControlProbe_00019"
#> [325] "NegControlProbe_00017" "NegControlProbe_00016"
#> [327] "NegControlProbe_00014" "NegControlProbe_00013"
#> [329] "NegControlProbe_00012" "NegControlProbe_00009"
#> [331] "NegControlProbe_00004" "NegControlProbe_00003"
#> [333] "NegControlProbe_00002" "antisense_PROKR2"
#> [335] "antisense_ULK3" "antisense_SCRIB"
#> [337] "antisense_TRMU" "antisense_MYLIP"
#> [339] "antisense_LGI3" "antisense_BCL2L15"
#> [341] "antisense_ADCY4" "NegControlCodeword_0500"
#> [343] "NegControlCodeword_0501" "NegControlCodeword_0502"
#> [345] "NegControlCodeword_0503" "NegControlCodeword_0504"
#> [347] "NegControlCodeword_0505" "NegControlCodeword_0506"
#> [349] "NegControlCodeword_0507" "NegControlCodeword_0508"
#> [351] "NegControlCodeword_0509" "NegControlCodeword_0510"
#> [353] "NegControlCodeword_0511" "NegControlCodeword_0512"
#> [355] "NegControlCodeword_0513" "NegControlCodeword_0514"
#> [357] "NegControlCodeword_0515" "NegControlCodeword_0516"
#> [359] "NegControlCodeword_0517" "NegControlCodeword_0518"
#> [361] "NegControlCodeword_0519" "NegControlCodeword_0520"
#> [363] "NegControlCodeword_0521" "NegControlCodeword_0522"
#> [365] "NegControlCodeword_0523" "NegControlCodeword_0524"
#> [367] "NegControlCodeword_0525" "NegControlCodeword_0526"
#> [369] "NegControlCodeword_0527" "NegControlCodeword_0528"
#> [371] "NegControlCodeword_0529" "NegControlCodeword_0530"
#> [373] "NegControlCodeword_0531" "NegControlCodeword_0532"
#> [375] "NegControlCodeword_0533" "NegControlCodeword_0534"
#> [377] "NegControlCodeword_0535" "NegControlCodeword_0536"
#> [379] "NegControlCodeword_0537" "NegControlCodeword_0538"
#> [381] "NegControlCodeword_0539" "NegControlCodeword_0540"
#> [383] "BLANK_0006" "BLANK_0013"
#> [385] "BLANK_0037" "BLANK_0069"
#> [387] "BLANK_0072" "BLANK_0087"
#> [389] "BLANK_0110" "BLANK_0114"
#> [391] "BLANK_0120" "BLANK_0147"
#> [393] "BLANK_0180" "BLANK_0186"
#> [395] "BLANK_0272" "BLANK_0278"
#> [397] "BLANK_0319" "BLANK_0321"
#> [399] "BLANK_0337" "BLANK_0350"
#> [401] "BLANK_0351" "BLANK_0352"
#> [403] "BLANK_0353" "BLANK_0354"
#> [405] "BLANK_0355" "BLANK_0356"
#> [407] "BLANK_0357" "BLANK_0358"
#> [409] "BLANK_0359" "BLANK_0360"
#> [411] "BLANK_0361" "BLANK_0362"
#> [413] "BLANK_0363" "BLANK_0364"
#> [415] "BLANK_0365" "BLANK_0366"
#> [417] "BLANK_0367" "BLANK_0368"
#> [419] "BLANK_0369" "BLANK_0370"
#> [421] "BLANK_0371" "BLANK_0372"
#> [423] "BLANK_0373" "BLANK_0374"
#> [425] "BLANK_0375" "BLANK_0376"
#> [427] "BLANK_0377" "BLANK_0378"
#> [429] "BLANK_0379" "BLANK_0380"
#> [431] "BLANK_0381" "BLANK_0382"
#> [433] "BLANK_0383" "BLANK_0384"
#> [435] "BLANK_0385" "BLANK_0386"
#> [437] "BLANK_0387" "BLANK_0388"
#> [439] "BLANK_0389" "BLANK_0390"
#> [441] "BLANK_0391" "BLANK_0392"
#> [443] "BLANK_0393" "BLANK_0394"
#> [445] "BLANK_0395" "BLANK_0396"
#> [447] "BLANK_0397" "BLANK_0398"
#> [449] "BLANK_0399" "BLANK_0400"
#> [451] "BLANK_0401" "BLANK_0402"
#> [453] "BLANK_0403" "BLANK_0404"
#> [455] "BLANK_0405" "BLANK_0406"
#> [457] "BLANK_0407" "BLANK_0408"
#> [459] "BLANK_0409" "BLANK_0410"
#> [461] "BLANK_0411" "BLANK_0412"
#> [463] "BLANK_0413" "BLANK_0414"
#> [465] "BLANK_0415" "BLANK_0416"
#> [467] "BLANK_0417" "BLANK_0418"
#> [469] "BLANK_0419" "BLANK_0420"
#> [471] "BLANK_0421" "BLANK_0422"
#> [473] "BLANK_0423" "BLANK_0424"
#> [475] "BLANK_0425" "BLANK_0426"
#> [477] "BLANK_0427" "BLANK_0428"
#> [479] "BLANK_0429" "BLANK_0430"
#> [481] "BLANK_0431" "BLANK_0432"
#> [483] "BLANK_0433" "BLANK_0434"
#> [485] "BLANK_0435" "BLANK_0436"
#> [487] "BLANK_0437" "BLANK_0438"
#> [489] "BLANK_0439" "BLANK_0440"
#> [491] "BLANK_0441" "BLANK_0442"
#> [493] "BLANK_0443" "BLANK_0444"
#> [495] "BLANK_0445" "BLANK_0446"
#> [497] "BLANK_0447" "BLANK_0448"
#> [499] "BLANK_0449" "BLANK_0450"
#> [501] "BLANK_0451" "BLANK_0452"
#> [503] "BLANK_0453" "BLANK_0454"
#> [505] "BLANK_0455" "BLANK_0456"
#> [507] "BLANK_0457" "BLANK_0458"
#> [509] "BLANK_0459" "BLANK_0460"
#> [511] "BLANK_0461" "BLANK_0462"
#> [513] "BLANK_0463" "BLANK_0464"
#> [515] "BLANK_0465" "BLANK_0466"
#> [517] "BLANK_0467" "BLANK_0468"
#> [519] "BLANK_0469" "BLANK_0470"
#> [521] "BLANK_0471" "BLANK_0472"
#> [523] "BLANK_0473" "BLANK_0474"
#> [525] "BLANK_0475" "BLANK_0476"
#> [527] "BLANK_0477" "BLANK_0478"
#> [529] "BLANK_0479" "BLANK_0480"
#> [531] "BLANK_0481" "BLANK_0482"
#> [533] "BLANK_0483" "BLANK_0484"
#> [535] "BLANK_0485" "BLANK_0486"
#> [537] "BLANK_0487" "BLANK_0488"
#> [539] "BLANK_0489" "BLANK_0497"
#> [541] "BLANK_0499"According to the Xenium paper (Janesick et al. 2022), there are 3 types of controls:
- probe controls to assess non-specific binding to RNA,
- decoding controls to assess misassigned genes, and
- genomic DNA (gDNA) controls to ensure the signal is from RNA.
The paper does not explain in detail how those control probes were designed, nor explain what the blank probes are. But the blank probes can be used as a negative control.
is_blank <- str_detect(rownames(sfe), "^BLANK_")
sum(is_blank)
#> [1] 159This should be number 1, the probe control
is_neg <- str_detect(rownames(sfe), "^NegControlProbe")
sum(is_neg)
#> [1] 20This should be number 2, the decoding control
is_neg2 <- str_detect(rownames(sfe), "^NegControlCodeword")
sum(is_neg2)
#> [1] 41This must be number 3, gDNA control
is_anti <- str_detect(rownames(sfe), "^antisense")
sum(is_anti)
#> [1] 8Also make an indicator of whether a feature is any sort of negative control
is_any_neg <- is_blank | is_neg | is_neg2 | is_antiThe addPerCellQCMetrics() function in the
scuttle package can conveniently add transcript counts,
proportion of total counts, and number of features detected for any
subset of features to the SCE object. Here we do this for the SFE
object, as SFE inherits from SCE.
sfe <- addPerCellQCMetrics(sfe, subsets = list(blank = is_blank,
negProbe = is_neg,
negCodeword = is_neg2,
anti = is_anti,
any_neg = is_any_neg))
names(colData(sfe))
#> [1] "Sample" "Barcode"
#> [3] "transcript_counts" "control_probe_counts"
#> [5] "control_codeword_counts" "cell_area"
#> [7] "nucleus_area" "sample_id"
#> [9] "nCounts" "nGenes"
#> [11] "nCounts_normed" "nGenes_normed"
#> [13] "prop_nuc" "sum"
#> [15] "detected" "subsets_blank_sum"
#> [17] "subsets_blank_detected" "subsets_blank_percent"
#> [19] "subsets_negProbe_sum" "subsets_negProbe_detected"
#> [21] "subsets_negProbe_percent" "subsets_negCodeword_sum"
#> [23] "subsets_negCodeword_detected" "subsets_negCodeword_percent"
#> [25] "subsets_anti_sum" "subsets_anti_detected"
#> [27] "subsets_anti_percent" "subsets_any_neg_sum"
#> [29] "subsets_any_neg_detected" "subsets_any_neg_percent"
#> [31] "total"Next we plot the proportion of transcript counts coming from any negative control.
cols_use <- names(colData(sfe))[str_detect(names(colData(sfe)), "_percent$")]
plotColDataHistogram(sfe, cols_use, bins = 100, ncol = 3)
#> Warning: Removed 285 rows containing non-finite outside the scale range
#> (`stat_bin()`).
The histogram is dominated by the bin at zero and there are some extreme outliers too few to be seen but evident from the scale of the x axis. We also plot the histogram only for cells with at least 1 count from a negative control. The NA’s come from cells that got segmented but have no transcripts detected.
plotColDataHistogram(sfe, cols_use, bins = 100, ncol = 3) +
scale_x_log10() +
annotation_logticks(sides = "b")
#> Warning in scale_x_log10(): log-10 transformation introduced
#> infinite values.
#> Warning: Removed 565577 rows containing non-finite outside the scale range
#> (`stat_bin()`).
The vast majority of these cells have less than 1% of transcript counts from negative controls, but there are outliers with up to 50%.
Next we plot the distribution of the number of negative control counts per cell:
cols_use2 <- names(colData(sfe))[str_detect(names(colData(sfe)), "_detected$")]
plotColDataHistogram(sfe, cols_use2, bins = 20, ncol = 3) +
# Avoid decimal breaks on x axis unless there're too few breaks
scale_x_continuous(breaks = scales::breaks_extended(Q = c(1,2,5)))
The counts are low, mostly zero, but there are outliers with up to 10 counts of all types aggregated. Then the outlier with 50% of counts from negative controls must have very low total real transcript counts to begin with.
The scuttle package can detect outliers, but by default
it assigns anything above zero as an outlier, since that is over 3
median absolute deviations (MADs) away from the median, which is 0, and
the MAD is 0 since the vast majority of cells don’t have any negative
control count. But it makes sense to allow a small proportion of
negative controls. Here we use the distribution just for cells with at
least 1 negative control count to find outliers. This distribution has a
very long tail and some definite outliers.
The code below extracts the outliers, based only on cells with at least one negative control count
get_neg_ctrl_outliers <- function(col, sfe) {
inds <- colData(sfe)$nCounts > 0 & colData(sfe)[[col]] > 0
df <- colData(sfe)[inds,]
outlier_inds <- isOutlier(df[[col]], type = "higher")
outliers <- rownames(df)[outlier_inds]
col2 <- str_remove(col, "^subsets_")
col2 <- str_remove(col2, "_percent$")
new_colname <- paste("is", col2, "outlier", sep = "_")
colData(sfe)[[new_colname]] <- colnames(sfe) %in% outliers
sfe
}
cols_use <- names(colData(sfe))[str_detect(names(colData(sfe)), "_percent$")]
for (n in cols_use) {
sfe <- get_neg_ctrl_outliers(n, sfe)
}
names(colData(sfe))
#> [1] "Sample" "Barcode"
#> [3] "transcript_counts" "control_probe_counts"
#> [5] "control_codeword_counts" "cell_area"
#> [7] "nucleus_area" "sample_id"
#> [9] "nCounts" "nGenes"
#> [11] "nCounts_normed" "nGenes_normed"
#> [13] "prop_nuc" "sum"
#> [15] "detected" "subsets_blank_sum"
#> [17] "subsets_blank_detected" "subsets_blank_percent"
#> [19] "subsets_negProbe_sum" "subsets_negProbe_detected"
#> [21] "subsets_negProbe_percent" "subsets_negCodeword_sum"
#> [23] "subsets_negCodeword_detected" "subsets_negCodeword_percent"
#> [25] "subsets_anti_sum" "subsets_anti_detected"
#> [27] "subsets_anti_percent" "subsets_any_neg_sum"
#> [29] "subsets_any_neg_detected" "subsets_any_neg_percent"
#> [31] "total" "is_blank_outlier"
#> [33] "is_negProbe_outlier" "is_negCodeword_outlier"
#> [35] "is_anti_outlier" "is_any_neg_outlier"Below we examine where the outliers are located in space:
plotSpatialFeature(sfe, "is_blank_outlier", colGeometryName = "cellSeg")
We find that the outliers are difficult to see:
plotColData(sfe, y = "is_blank_outlier", x = "cell_area",
point_fun = function(...) list()) 
The analysis reveals that the outliers seem to be smaller. Outliers for negative probe controls and negative codeword controls are also hard to see on the plot, so their plots are skipped here. But the top left region in the tissue tends to have more counts from antisense controls.
plotSpatialFeature(sfe, "is_anti_outlier", colGeometryName = "cellSeg")
Now that we have identified the outliers, we can remove them along with empty cells before proceeding to further analysis:
inds_keep <- sfe$nCounts > 0 & sfe$nucleus_area < 400 & !sfe$is_anti_outlier &
!sfe$is_blank_outlier & !sfe$is_negCodeword_outlier & !sfe$is_negProbe_outlier
(sfe <- sfe[,inds_keep])
#> class: SpatialFeatureExperiment
#> dim: 541 117503
#> metadata(1): Samples
#> assays(1): counts
#> rownames(541): ABCC11 ACTA2 ... BLANK_0497 BLANK_0499
#> rowData names(6): ID Symbol ... vars cv2
#> colnames(117503): 1 2 ... 118707 118708
#> colData names(36): Sample Barcode ... is_anti_outlier
#> is_any_neg_outlier
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x_centroid y_centroid
#> imgData names(1): sample_id
#>
#> unit:
#> Geometries:
#> colGeometries: centroids (POINT), cellSeg (POLYGON), nucSeg (GEOMETRY)
#>
#> Graphs:
#> sample01:Over 1000 cells were removed.
Next we check how many negative control features are detected per cell:
plotColDataHistogram(sfe, cols_use2, bins = 20, ncol = 3) +
# Avoid decimal breaks on x axis unless there're too few breaks
scale_x_continuous(breaks = scales::breaks_extended(3, Q = c(1,2,5)))
There are at most 3 counts per cell per type. For the non-outliers, each type is at most around 1%, so this data looks good.
Genes
Here we look at the mean and variance of each gene
Real genes generally have higher mean expression across cells than negative controls.
rowData(sfe)$is_neg <- is_any_neg
plotRowData(sfe, x = "means", y = "is_neg") +
scale_y_log10() +
annotation_logticks(sides = "b")
Here the real genes and negative controls are plotted in different colors
plotRowData(sfe, x="means", y="vars", color_by = "is_neg") +
geom_abline(slope = 1, intercept = 0, color = "red") +
scale_x_log10() + scale_y_log10() +
annotation_logticks() +
coord_equal() +
labs(color = "Negative control")
The red line \(y = x\) is expected if the data follows a Poisson distribution. Negative controls and real genes form mostly separate clusters. Negative controls stick close to the line, while real genes are overdispersed. Unlike in the CosMX dataset, the negative controls don’t seem overdispersed.
Spatial autocorrelation of QC metrics
There’s a sparse and a dense region. This poses the question of what type of neighborhood graph to use, e.g. it is conceivable that cells in the sparse region should just be singletons. Furthermore, it is unclear what the length scale of their influence might be. It might depend on the cell type and how contact and secreted signals are used in the cell type, and length scale of the influence. If k nearest neighbors are used, then the neighbors in the dense region are much closer together than those in the sparse region. If distance based neighbors are used, then cells in the dense region will have more neighbors than cells in the sparse region, and the sparse region can break into multiple compartments if the distance cutoff is not long enough.
For the purpose of demonstration, we use k nearest neighbors with \(k = 5\), with inverse distance weighting. Note that using more neighbors leads to longer computation time of spatial autocorrelation metrics.
system.time(
colGraph(sfe, "knn5") <- findSpatialNeighbors(sfe, method = "knearneigh",
dist_type = "idw", k = 5,
style = "W")
)
#> user system elapsed
#> 7.035 0.042 7.086
sfe <- colDataMoransI(sfe, c("nCounts", "nGenes", "cell_area", "nucleus_area"),
colGraphName = "knn5")
colFeatureData(sfe)[c("nCounts", "nGenes", "cell_area", "nucleus_area"),]
#> DataFrame with 4 rows and 2 columns
#> moran_sample01 K_sample01
#> <numeric> <numeric>
#> nCounts 0.422387 5.55509
#> nGenes 0.401395 3.13694
#> cell_area 0.628837 7.57098
#> nucleus_area 0.377248 6.88331Global Moran’s I indicatse positive spatial autocorrelation. As the strength of spatial autocorrelation can vary spatially, we also run local Moran’s I.
sfe <- colDataUnivariate(sfe, type = "localmoran",
features = c("nCounts", "nGenes", "cell_area",
"nucleus_area"),
colGraphName = "knn5", BPPARAM = MulticoreParam(2))The pointsize argument adjusts the point size in
scattermore. The default is 0, meaning single pixels, but
since the cells in the sparse region are hard to see that way, we
increase pointsize. We would still plot the polygons in
larger single panel plots, but use scattermore in
multi-panel plots where the polygons in each panel are invisible anyway
due to the small size to save some time.
plotLocalResult(sfe, "localmoran",
features = c("nCounts", "nGenes", "cell_area", "nucleus_area"),
colGeometryName = "centroids", scattermore = TRUE,
divergent = TRUE, diverge_center = 0, pointsize = 1)
Interestingly, nCounts is more homogeneous in the interior of the dense region, while nGenes is more homogeneous by the edge of the dense region. As expected, cell area is more homogeneous in the sparse region. However, the nucleus area is more homogeneous in the interior of the dense region.
Moran plot for nCounts
sfe <- colDataUnivariate(sfe, "moran.plot", "nCounts", colGraphName = "knn5")
p1 <- moranPlot(sfe, "nCounts", binned = TRUE, plot_influential = FALSE)
p2 <- moranPlot(sfe, "nCounts", binned = TRUE)
p1 / p2 + plot_layout(guides = "collect")
There are no obvious clusters here. In the lower panel, the 2D histogram of influential points is plotted in red.
Moran’s I
By default, for gene expression, the log normalized counts are used in spatial autocorrelation metrics, so before running Moran’s I, we normalize the data.
sfe <- logNormCounts(sfe)Use more cores if available to speed this up.
system.time(
sfe <- runMoransI(sfe, colGraphName = "knn5", BPPARAM = MulticoreParam(2))
)
#> user system elapsed
#> 33.566 3.288 19.238
rowData(sfe)$is_neg <- is_any_neg
plotRowData(sfe, x = "moran_sample01", y = "is_neg")
As expected, generally the negative controls are tightly clustered around 0, while the real genes have positive Moran’s I, which means there is generally no technical artifact spatial trend. No significantly negative Moran’s I is observed. Why is negative spatial autocorrelation so rare in gene expression?
What are the two negative controls with a sizable Moran’s I?
ord <- order(rowData(sfe)$moran_sample01[is_any_neg], decreasing = TRUE)[1:2]
top_neg <- rownames(sfe)[is_any_neg][ord]
plotSpatialFeature(sfe, top_neg, colGeometryName = "centroids",
scattermore = TRUE, pointsize = 1)
There is somewhat a spatial trend for that antisense probe, with more detected in the upper left. However, this might not significantly affect other results since there are at most 2 counts and at most about 1% of all counts in each cell. The negative control codeword has at most 1 count per cell and the cells with this negative control detected seem to be few and far between.
These are the most detected negative controls, and the most detected one is also the one with the highest Moran’s I among negative controls. However, the other negative control with higher Moran’s I is not among the most detected.
head(sort(rowData(sfe)$means[is_any_neg], decreasing = TRUE), 15)
#> antisense_PROKR2 antisense_SCRIB antisense_BCL2L15
#> 0.0192761036 0.0131741317 0.0066806805
#> antisense_TRMU antisense_MYLIP antisense_ULK3
#> 0.0042041480 0.0030807724 0.0028169494
#> BLANK_0485 antisense_ADCY4 antisense_LGI3
#> 0.0023403658 0.0019829281 0.0017871884
#> BLANK_0430 NegControlProbe_00035 NegControlProbe_00012
#> 0.0015063445 0.0010978443 0.0010382714
#> NegControlProbe_00033 NegControlProbe_00014 BLANK_0120
#> 0.0009446567 0.0009361463 0.0009276359What are the genes with the highest Moran’s I?
top_moran <- rownames(sfe)[order(rowData(sfe)$moran_sample01, decreasing = TRUE)[1:6]]
plotSpatialFeature(sfe, top_moran, colGeometryName = "centroids",
scattermore = TRUE, ncol = 2, pointsize = 0.5)
They all highlight the same histological regions, as in the CosMX vignette. How does Moran’s I relate to gene expression level?
plotRowData(sfe, x = "means", y = "moran_sample01")
Very highly expressed genes have higher Moran’s I, but there are some less expressed genes with higher Moran’s I as well.
Non-spatial dimension reduction and clustering
Here we run non-spatial PCA as for scRNA-seq data
set.seed(29)
sfe <- runPCA(sfe, ncomponents = 30, scale = TRUE, BSPARAM = IrlbaParam())
ElbowPlot(sfe, ndims = 30)
plotDimLoadings(sfe, dims = 1:6)
spatialReducedDim(sfe, "PCA", 6, colGeometryName = "centroids", divergent = TRUE,
diverge_center = 0, ncol = 2, scattermore = TRUE, pointsize = 0.5)
While spatial region is not explicitly used, the PC’s highlight spatial regions due to spatial autocorrelation in gene expression and histological regions with different cell types.
Non-spatial clustering and locating the clusters in space
colData(sfe)$cluster <- clusterRows(reducedDim(sfe, "PCA")[,1:15],
BLUSPARAM = SNNGraphParam(
cluster.fun = "leiden",
cluster.args = list(
resolution_parameter = 0.5,
objective_function = "modularity")))
#> Warning in (function (to_check, X, clust_centers, clust_info, dtype, nn, :
#> detected tied distances to neighbors, see ?'BiocNeighbors-ties'Now the scater can also rasterize the plots with lots of
points with the rasterise argument, but with a different
mechanism from scattermore that requires more system
dependencies.
plotPCA(sfe, ncomponents = 4, colour_by = "cluster", rasterise = FALSE)
Plot the location of the clusters in space
plotSpatialFeature(sfe, "cluster", colGeometryName = "cellSeg")
Differential expression
Cluster marker genes are found with Wilcoxon rank sum test as commonly done for scRNA-seq.
markers <- findMarkers(sfe, groups = colData(sfe)$cluster,
test.type = "wilcox", pval.type = "all", direction = "up")It’s already sorted by p-values:
markers[[6]]
#> DataFrame with 541 rows and 16 columns
#> p.value FDR summary.AUC AUC.1 AUC.2 AUC.3
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> TENT5C 1.63326e-296 8.83591e-294 0.967126 0.925332 0.970773 0.973994
#> SEC11C 7.21967e-230 1.95292e-227 0.910578 0.899842 0.926434 0.943356
#> MZB1 1.53225e-208 2.76316e-206 0.890901 0.954750 0.952539 0.953965
#> PRDM1 1.85488e-171 2.50872e-169 0.851328 0.902618 0.844973 0.845018
#> SLAMF7 1.36571e-121 1.47770e-119 0.797638 0.973023 0.928639 0.964871
#> ... ... ... ... ... ... ...
#> TRAF4 1 1 0.190188 0.190188 0.498645 0.463686
#> USP53 1 1 0.237375 0.237375 0.515006 0.439300
#> VWF 1 1 0.148324 0.521143 0.492783 0.472879
#> ZEB2 1 1 0.227203 0.764394 0.227203 0.317134
#> ZNF562 1 1 0.286342 0.286342 0.426846 0.479020
#> AUC.4 AUC.5 AUC.7 AUC.8 AUC.9 AUC.10 AUC.11
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> TENT5C 0.978959 0.957859 0.975948 0.951552 0.970808 0.957569 0.967970
#> SEC11C 0.955726 0.912545 0.930013 0.929966 0.935488 0.925262 0.924159
#> MZB1 0.957403 0.947673 0.955833 0.909592 0.954841 0.954570 0.946052
#> PRDM1 0.886670 0.720187 0.838240 0.903536 0.876730 0.880933 0.863019
#> SLAMF7 0.970315 0.931636 0.966127 0.973165 0.967341 0.927913 0.951534
#> ... ... ... ... ... ... ... ...
#> TRAF4 0.498581 0.510160 0.476258 0.267629 0.424212 0.485015 0.482499
#> USP53 0.367034 0.532702 0.506453 0.317515 0.461329 0.516505 0.492849
#> VWF 0.436092 0.501213 0.148324 0.548782 0.541546 0.545200 0.438237
#> ZEB2 0.296810 0.491759 0.278638 0.784944 0.645498 0.568365 0.301803
#> ZNF562 0.468102 0.525932 0.449265 0.297023 0.407262 0.512553 0.480637
#> AUC.12 AUC.13 AUC.14
#> <numeric> <numeric> <numeric>
#> TENT5C 0.953956 0.967126 0.979535
#> SEC11C 0.902004 0.910578 0.926748
#> MZB1 0.939414 0.890901 0.952151
#> PRDM1 0.834607 0.851328 0.899233
#> SLAMF7 0.944967 0.797638 0.971807
#> ... ... ... ...
#> TRAF4 0.475006 0.352237 0.389822
#> USP53 0.552792 0.545264 0.333889
#> VWF 0.507755 0.508495 0.557926
#> ZEB2 0.368603 0.251563 0.771073
#> ZNF562 0.525134 0.543883 0.375170The code below extracts the significant markers for each cluster:
genes_use <- vapply(markers, function(x) rownames(x)[1], FUN.VALUE = character(1))
plotExpression(sfe, genes_use, x = "cluster", point_fun = function(...) list())
This allows for plotting more top marker genes in a heatmap:
genes_use2 <- unique(unlist(lapply(markers, function(x) rownames(x)[1:5])))
plotGroupedHeatmap(sfe, genes_use2, group = "cluster", colour = scales::viridis_pal()(100))
Local spatial statistics of marker genes
First we plot those genes in space as a reference
plotSpatialFeature(sfe, genes_use, colGeometryName = "centroids", ncol = 3,
pointsize = 0.3, scattermore = TRUE)
Global Moran’s I of these marker genes is shown below:
setNames(rowData(sfe)[genes_use, "moran_sample01"], genes_use)
#> FOXA1 FGL2 LUM ADIPOQ CD3E TENT5C CD93 GATA3
#> 0.7421765 0.2604219 0.6812312 0.5493112 0.4015241 0.2806543 0.3250982 0.6558350
#> MYLK APOC1 CPA3 MS4A1 LILRA4 KRT15
#> 0.4653625 0.2696177 0.1904912 0.2144728 0.1092981 0.5425155All these marker genes have positive spatial autocorrelation, but some stronger than others.
Local Moran’s I of these marker genes is shown below:
sfe <- runUnivariate(sfe, "localmoran", features = genes_use, colGraphName = "knn5",
BPPARAM = MulticoreParam(2))
plotLocalResult(sfe, "localmoran", features = genes_use,
colGeometryName = "centroids", ncol = 3, divergent = TRUE,
diverge_center = 0, scattermore = TRUE, pointsize = 0.3)
It seems that some histological regions tend to be more spatially homogenous in gene expression than others. The epithelial region tends to be more homogenous. For some genes, regions with higher expression also have higher local Moran’s I, such as FOXA1 and GATA3, while for some genes, this is not the case, such as FGL2 and LUM.
Finally, we assess local spatial heteroscdasticity (LOSH) for these marker genes to find local heterogeneity:
sfe <- runUnivariate(sfe, "LOSH", features = genes_use, colGraphName = "knn5",
BPPARAM = MulticoreParam(2))
plotLocalResult(sfe, "LOSH", features = genes_use,
colGeometryName = "centroids", ncol = 3, scattermore = TRUE,
pointsize = 0.3)
Again, just like in the CosMX dataset, LOSH is higher where the gene is more highly expressed in some (e.g. CD3E, LUM, TENT5C) but not all cases (e.g. FOXA1, GATA3). This may be due to spatial distribution of different cell types.
Session info
sessionInfo()
#> R version 4.4.0 (2024-04-24)
#> Platform: x86_64-apple-darwin20
#> Running under: macOS Ventura 13.6.6
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: UTC
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] patchwork_1.2.0 scran_1.32.0
#> [3] bluster_1.14.0 BiocSingular_1.20.0
#> [5] BiocParallel_1.38.0 scater_1.32.0
#> [7] scuttle_1.14.0 stringr_1.5.1
#> [9] ggplot2_3.5.1 SpatialExperiment_1.14.0
#> [11] SingleCellExperiment_1.26.0 SummarizedExperiment_1.34.0
#> [13] Biobase_2.64.0 GenomicRanges_1.56.0
#> [15] GenomeInfoDb_1.40.0 IRanges_2.38.0
#> [17] S4Vectors_0.42.0 BiocGenerics_0.50.0
#> [19] MatrixGenerics_1.16.0 matrixStats_1.3.0
#> [21] SFEData_1.6.0 Voyager_1.6.0
#> [23] SpatialFeatureExperiment_1.6.1
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.4.0 bitops_1.0-7
#> [3] filelock_1.0.3 tibble_3.2.1
#> [5] R.oo_1.26.0 lifecycle_1.0.4
#> [7] sf_1.0-16 edgeR_4.2.0
#> [9] lattice_0.22-6 magrittr_2.0.3
#> [11] limma_3.60.0 sass_0.4.9
#> [13] rmarkdown_2.27 jquerylib_0.1.4
#> [15] yaml_2.3.8 metapod_1.12.0
#> [17] sp_2.1-4 cowplot_1.1.3
#> [19] RColorBrewer_1.1-3 DBI_1.2.2
#> [21] abind_1.4-5 zlibbioc_1.50.0
#> [23] purrr_1.0.2 R.utils_2.12.3
#> [25] RCurl_1.98-1.14 rappdirs_0.3.3
#> [27] GenomeInfoDbData_1.2.12 ggrepel_0.9.5
#> [29] irlba_2.3.5.1 terra_1.7-71
#> [31] pheatmap_1.0.12 units_0.8-5
#> [33] RSpectra_0.16-1 dqrng_0.4.0
#> [35] pkgdown_2.0.9 DelayedMatrixStats_1.26.0
#> [37] codetools_0.2-20 DropletUtils_1.24.0
#> [39] DelayedArray_0.30.1 tidyselect_1.2.1
#> [41] farver_2.1.2 UCSC.utils_1.0.0
#> [43] memuse_4.2-3 ScaledMatrix_1.12.0
#> [45] viridis_0.6.5 BiocFileCache_2.12.0
#> [47] jsonlite_1.8.8 BiocNeighbors_1.22.0
#> [49] e1071_1.7-14 systemfonts_1.1.0
#> [51] tools_4.4.0 ggnewscale_0.4.10
#> [53] ragg_1.3.2 Rcpp_1.0.12
#> [55] glue_1.7.0 gridExtra_2.3
#> [57] SparseArray_1.4.3 mgcv_1.9-1
#> [59] xfun_0.44 EBImage_4.46.0
#> [61] dplyr_1.1.4 HDF5Array_1.32.0
#> [63] withr_3.0.0 BiocManager_1.30.23
#> [65] fastmap_1.2.0 boot_1.3-30
#> [67] rhdf5filters_1.16.0 fansi_1.0.6
#> [69] spData_2.3.0 digest_0.6.35
#> [71] rsvd_1.0.5 R6_2.5.1
#> [73] mime_0.12 textshaping_0.3.7
#> [75] colorspace_2.1-0 wk_0.9.1
#> [77] scattermore_1.2 jpeg_0.1-10
#> [79] RSQLite_2.3.6 R.methodsS3_1.8.2
#> [81] hexbin_1.28.3 utf8_1.2.4
#> [83] generics_0.1.3 data.table_1.15.4
#> [85] class_7.3-22 httr_1.4.7
#> [87] htmlwidgets_1.6.4 S4Arrays_1.4.0
#> [89] spdep_1.3-3 pkgconfig_2.0.3
#> [91] scico_1.5.0 gtable_0.3.5
#> [93] blob_1.2.4 XVector_0.44.0
#> [95] htmltools_0.5.8.1 fftwtools_0.9-11
#> [97] scales_1.3.0 png_0.1-8
#> [99] knitr_1.46 rjson_0.2.21
#> [101] nlme_3.1-164 curl_5.2.1
#> [103] proxy_0.4-27 cachem_1.1.0
#> [105] rhdf5_2.48.0 BiocVersion_3.19.1
#> [107] KernSmooth_2.23-24 parallel_4.4.0
#> [109] vipor_0.4.7 AnnotationDbi_1.66.0
#> [111] desc_1.4.3 s2_1.1.6
#> [113] pillar_1.9.0 grid_4.4.0
#> [115] vctrs_0.6.5 dbplyr_2.5.0
#> [117] beachmat_2.20.0 sfheaders_0.4.4
#> [119] cluster_2.1.6 beeswarm_0.4.0
#> [121] evaluate_0.23 zeallot_0.1.0
#> [123] magick_2.8.3 cli_3.6.2
#> [125] locfit_1.5-9.9 compiler_4.4.0
#> [127] rlang_1.1.3 crayon_1.5.2
#> [129] labeling_0.4.3 classInt_0.4-10
#> [131] fs_1.6.4 ggbeeswarm_0.7.2
#> [133] stringi_1.8.4 viridisLite_0.4.2
#> [135] deldir_2.0-4 munsell_0.5.1
#> [137] Biostrings_2.72.0 tiff_0.1-12
#> [139] Matrix_1.7-0 ExperimentHub_2.12.0
#> [141] sparseMatrixStats_1.16.0 bit64_4.0.5
#> [143] Rhdf5lib_1.26.0 KEGGREST_1.44.0
#> [145] statmod_1.5.0 highr_0.10
#> [147] AnnotationHub_3.12.0 igraph_2.0.3
#> [149] memoise_2.0.1 bslib_0.7.0
#> [151] bit_4.0.5