Skip to contents

Introduction

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.

(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.

colnames(sfe) <- seq_len(ncol(sfe))

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_panel

Here 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] 0

So 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:

  1. probe controls to assess non-specific binding to RNA,
  2. decoding controls to assess misassigned genes, and
  3. 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] 159

This should be number 1, the probe control

is_neg <- str_detect(rownames(sfe), "^NegControlProbe")
sum(is_neg)
#> [1] 20

This should be number 2, the decoding control

is_neg2 <- str_detect(rownames(sfe), "^NegControlCodeword")
sum(is_neg2)
#> [1] 41

This must be number 3, gDNA control

is_anti <- str_detect(rownames(sfe), "^antisense")
sum(is_anti)
#> [1] 8

Also make an indicator of whether a feature is any sort of negative control

is_any_neg <- is_blank | is_neg | is_neg2 | is_anti

The 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

rowData(sfe)$means <- rowMeans(counts(sfe))
rowData(sfe)$vars <- rowVars(counts(sfe))

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.88331

Global 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.0009276359

What 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.375170

The 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.5425155

All 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

References

Janesick, Amanda, Robert Shelansky, Andrew Gottscho, Florian Wagner, Morgane Rouault, Ghezal Beliakoff, Michelli Faria de Oliveira, et al. 2022. “High Resolution Mapping of the Breast Cancer Tumor Microenvironment Using Integrated Single Cell, Spatial and in Situ Analysis of FFPE Tissue.” bioRxiv.