Coercion of GeoMxSet to Seurat and SpatialExperiment Objects

Overview

This tutorial demonstrates how to coerce GeoMxSet objects into Seurat or SpatialExperiment objects and the subsequent analyses. For more examples of what analyses are available in these objects, look at these Seurat or SpatialExperiment vignettes.

Data Processing

Data Processing should occur in GeomxTools. Due to the unique nature of the regions of interest (ROIs), it is recommended to use the preproccesing steps available in GeomxTools rather than the single-cell made preprocessing available in Seurat.

library(GeomxTools)
library(Seurat)
library(SpatialDecon)
library(patchwork)
datadir <- system.file("extdata", "DSP_NGS_Example_Data",
                       package="GeomxTools")
DCCFiles <- dir(datadir, pattern=".dcc$", full.names=TRUE)
PKCFiles <- unzip(zipfile = file.path(datadir,  "/pkcs.zip"))
SampleAnnotationFile <- file.path(datadir, "annotations.xlsx")

demoData <-
  suppressWarnings(readNanoStringGeoMxSet(dccFiles = DCCFiles,
                                          pkcFiles = PKCFiles,
                                          phenoDataFile = SampleAnnotationFile,
                                          phenoDataSheet = "CW005",
                                          phenoDataDccColName = "Sample_ID",
                                          protocolDataColNames = c("aoi",
                                                                   "cell_line",
                                                                   "roi_rep",
                                                                   "pool_rep",
                                                                   "slide_rep")))

After reading in the object, we will do a couple of QC steps.

  1. Shift all 0 counts by 1
  2. Flag low quality ROIs
  3. Flag low quality probes
  4. Remove low quality ROIs and probes
demoData <- shiftCountsOne(demoData, useDALogic=TRUE)
demoData <- setSegmentQCFlags(demoData, qcCutoffs = list(percentSaturation = 45))
demoData <- setBioProbeQCFlags(demoData)

# low sequenced ROIs
lowSaturation <- which(protocolData(demoData)[["QCFlags"]]$LowSaturation)

# probes that are considered outliers 
lowQCprobes <- which(featureData(demoData)[["QCFlags"]]$LowProbeRatio | 
                       featureData(demoData)[["QCFlags"]]$GlobalGrubbsOutlier)

# remove low quality ROIs and probes
passedQC <- demoData[-lowQCprobes, -lowSaturation]

dim(demoData)
## Features  Samples 
##     8707       88
dim(passedQC)
## Features  Samples 
##     8698       83

Objects must be aggregated to Target level data before coercing. This changes the row (gene) information to be the gene name rather than the probe ID.

featureType(passedQC)
## [1] "Probe"
data.frame(assayData(passedQC)[["exprs"]][seq_len(3), seq_len(3)])
##            DSP.1001250002642.A02.dcc DSP.1001250002642.A03.dcc
## RTS0039454                       294                       239
## RTS0039455                       270                       281
## RTS0039456                       255                       238
##            DSP.1001250002642.A04.dcc
## RTS0039454                         6
## RTS0039455                         6
## RTS0039456                         3
target_demoData <- aggregateCounts(passedQC)

featureType(target_demoData)
## [1] "Target"
data.frame(assayData(target_demoData)[["exprs"]][seq_len(3), seq_len(3)])
##       DSP.1001250002642.A02.dcc DSP.1001250002642.A03.dcc
## ACTA2                328.286182                323.490808
## FOXA2                  4.919019                  4.919019
## NANOG                  2.954177                  4.128918
##       DSP.1001250002642.A04.dcc
## ACTA2                  6.081111
## FOXA2                  6.942503
## NANOG                  8.359554

It is recommended to normalize using a GeoMx specific model before coercing. The normalized data is now in the assayData slot called “q_norm”.

norm_target_demoData <- normalize(target_demoData, norm_method="quant",
                                  desiredQuantile = .75, toElt = "q_norm")

assayDataElementNames(norm_target_demoData)
## [1] "exprs"  "q_norm"
data.frame(assayData(norm_target_demoData)[["q_norm"]][seq_len(3), seq_len(3)])
##       DSP.1001250002642.A02.dcc DSP.1001250002642.A03.dcc
## ACTA2                349.571598                344.257297
## FOXA2                  5.237958                  5.234795
## NANOG                  3.145720                  4.393974
##       DSP.1001250002642.A04.dcc
## ACTA2                  3.968122
## FOXA2                  4.530208
## NANOG                  5.454880

Seurat

Seurat Coercion

The three errors that can occur when trying to coerce to Seurat are:

  1. object must be on the target level
  2. object should be normalized, if you want raw data you can set forceRaw to TRUE
  3. normalized count matrix name must be valid
as.Seurat(demoData)
## Error in as.Seurat.NanoStringGeoMxSet(demoData): Data must be on Target level before converting to a Seurat Object
as.Seurat(target_demoData, normData = "exprs")
## Error in as.Seurat.NanoStringGeoMxSet(target_demoData, normData = "exprs"): It is NOT recommended to use Seurat's normalization for GeoMx data. 
##              Normalize using GeomxTools::normalize() or set forceRaw to TRUE if you want to continue with Raw data
as.Seurat(norm_target_demoData, normData = "exprs_norm")
## Error in as.Seurat.NanoStringGeoMxSet(norm_target_demoData, normData = "exprs_norm"): The normData name "exprs_norm" is not a valid assay name. Valid names are: exprs, q_norm

After coercing to a Seurat object all of the metadata is still accessible.

demoSeurat <- as.Seurat(x = norm_target_demoData, normData = "q_norm")

demoSeurat # overall data object
## An object of class Seurat 
## 1821 features across 83 samples within 1 assay 
## Active assay: GeoMx (1821 features, 0 variable features)
##  2 layers present: counts, data
head(demoSeurat, 3) # most important ROI metadata
##                           orig.ident nCount_GeoMx nFeature_GeoMx
## DSP-1001250002642-A02.dcc      GeoMx     63524.55           1821
## DSP-1001250002642-A03.dcc      GeoMx     62357.01           1821
## DSP-1001250002642-A04.dcc      GeoMx     82370.45           1821
##                                              slide name
## DSP-1001250002642-A02.dcc 6panel-old-slide1 (PTL-10891)
## DSP-1001250002642-A03.dcc 6panel-old-slide1 (PTL-10891)
## DSP-1001250002642-A04.dcc 6panel-old-slide1 (PTL-10891)
##                                          scan name
## DSP-1001250002642-A02.dcc cw005 (PTL-10891) Slide1
## DSP-1001250002642-A03.dcc cw005 (PTL-10891) Slide1
## DSP-1001250002642-A04.dcc cw005 (PTL-10891) Slide1
##                                                                                        panel
## DSP-1001250002642-A02.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
## DSP-1001250002642-A03.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
## DSP-1001250002642-A04.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
##                           roi           segment     area
## DSP-1001250002642-A02.dcc   1 Geometric Segment 31318.73
## DSP-1001250002642-A03.dcc   2 Geometric Segment 31318.73
## DSP-1001250002642-A04.dcc   3 Geometric Segment 31318.73
##                           NegGeoMean_Six-gene_test_v1_v1.1
## DSP-1001250002642-A02.dcc                         1.487738
## DSP-1001250002642-A03.dcc                         2.518775
## DSP-1001250002642-A04.dcc                         2.847315
##                           NegGeoMean_VnV_GeoMx_Hs_CTA_v1.2
## DSP-1001250002642-A02.dcc                         3.722751
## DSP-1001250002642-A03.dcc                         3.068217
## DSP-1001250002642-A04.dcc                         3.556275
##                           NegGeoSD_Six-gene_test_v1_v1.1
## DSP-1001250002642-A02.dcc                       1.560397
## DSP-1001250002642-A03.dcc                       1.820611
## DSP-1001250002642-A04.dcc                       1.654831
##                           NegGeoSD_VnV_GeoMx_Hs_CTA_v1.2 q_norm_qFactors
## DSP-1001250002642-A02.dcc                       1.796952       0.9391100
## DSP-1001250002642-A03.dcc                       1.806070       0.9396774
## DSP-1001250002642-A04.dcc                       1.762066       1.5324910
##                                        SampleID                       aoi
## DSP-1001250002642-A02.dcc DSP-1001250002642-A02 Geometric Segment-aoi-001
## DSP-1001250002642-A03.dcc DSP-1001250002642-A03 Geometric Segment-aoi-001
## DSP-1001250002642-A04.dcc DSP-1001250002642-A04 Geometric Segment-aoi-001
##                           cell_line roi_rep pool_rep slide_rep
## DSP-1001250002642-A02.dcc    HS578T       1        1         1
## DSP-1001250002642-A03.dcc    HS578T       2        1         1
## DSP-1001250002642-A04.dcc       HEL       1        1         1
demoSeurat@misc[1:8] # experiment data
## $PKCFileName
##            VnV_GeoMx_Hs_CTA_v1.2            Six-gene_test_v1_v1.1 
## "VnV Cancer Transcriptome Atlas"           "Six gene test custom" 
## 
## $PKCModule
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##    "VnV_GeoMx_Hs_CTA"    "Six-gene_test_v1" 
## 
## $PKCFileVersion
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                   1.2                   1.1 
## 
## $PKCFileDate
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##              "200518"              "200707" 
## 
## $AnalyteType
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                 "RNA"                 "RNA" 
## 
## $MinArea
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                 16000                 16000 
## 
## $MinNuclei
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                   200                   200 
## 
## $shiftedByOne
## [1] TRUE
head(demoSeurat@misc$sequencingMetrics) # sequencing metrics
##                           FileVersion SoftwareVersion       Date      Plate_ID
## DSP-1001250002642-A02.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A03.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A04.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A05.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A06.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A07.dcc         0.1           1.0.0 2020-07-14 1001250002642
##                           Well            SeqSetId    Raw Trimmed Stitched
## DSP-1001250002642-A02.dcc  A02 VH00121:3:AAAG2YWM5 646250  646250   616150
## DSP-1001250002642-A03.dcc  A03 VH00121:3:AAAG2YWM5 629241  629241   603243
## DSP-1001250002642-A04.dcc  A04 VH00121:3:AAAG2YWM5 831083  831083   798188
## DSP-1001250002642-A05.dcc  A05 VH00121:3:AAAG2YWM5 884485  884485   849060
## DSP-1001250002642-A06.dcc  A06 VH00121:3:AAAG2YWM5 781936  781936   751930
## DSP-1001250002642-A07.dcc  A07 VH00121:3:AAAG2YWM5 703034  703034   674815
##                           Aligned umiQ30 rtsQ30 DeduplicatedReads
## DSP-1001250002642-A02.dcc  610390 0.9785 0.9804            312060
## DSP-1001250002642-A03.dcc  597280 0.9784 0.9811            305528
## DSP-1001250002642-A04.dcc  791804 0.9785 0.9801            394981
## DSP-1001250002642-A05.dcc  842133 0.9796 0.9814            424162
## DSP-1001250002642-A06.dcc  744669 0.9779 0.9803            355121
## DSP-1001250002642-A07.dcc  668726 0.9776 0.9797            341008
##                                              NTC_ID NTC Trimmed (%)
## DSP-1001250002642-A02.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A03.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A04.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A05.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A06.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A07.dcc DSP-1001250002642-A01.dcc   7         100
##                           Stitched (%) Aligned (%) Saturated (%)
## DSP-1001250002642-A02.dcc     95.34236    94.45106      48.87531
## DSP-1001250002642-A03.dcc     95.86836    94.92071      48.84677
## DSP-1001250002642-A04.dcc     96.04191    95.27376      50.11632
## DSP-1001250002642-A05.dcc     95.99484    95.21168      49.63242
## DSP-1001250002642-A06.dcc     96.16260    95.23401      52.31156
## DSP-1001250002642-A07.dcc     95.98611    95.12001      49.00632
head(demoSeurat@misc$QCMetrics$QCFlags) # QC metrics
##                           LowReads LowTrimmed LowStitched LowAligned
## DSP-1001250002642-A02.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A03.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A04.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A05.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A06.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A07.dcc    FALSE      FALSE       FALSE      FALSE
##                           LowSaturation LowNegatives HighNTC LowArea
## DSP-1001250002642-A02.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A03.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A04.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A05.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A06.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A07.dcc         FALSE         TRUE   FALSE   FALSE
head(demoSeurat@assays$GeoMx@meta.data) # gene metadata
##   TargetName                Module  CodeClass        GeneID SystematicName
## 1      ACTA2 VnV_GeoMx_Hs_CTA_v1.2 Endogenous            59          ACTA2
## 2      FOXA2 VnV_GeoMx_Hs_CTA_v1.2 Endogenous          3170          FOXA2
## 3      NANOG VnV_GeoMx_Hs_CTA_v1.2 Endogenous 79923, 388112 NANOG, NANOGP8
## 4       TRAC VnV_GeoMx_Hs_CTA_v1.2 Endogenous          <NA>           TRAC
## 5    TRBC1/2 VnV_GeoMx_Hs_CTA_v1.2 Endogenous          <NA>          TRBC1
## 6       TRDC VnV_GeoMx_Hs_CTA_v1.2 Endogenous          <NA>           TRDC
##   Negative
## 1    FALSE
## 2    FALSE
## 3    FALSE
## 4    FALSE
## 5    FALSE
## 6    FALSE

All Seurat functionality is available after coercing. Outputs might differ if the ident value is set or not.

VlnPlot(demoSeurat, features = "nCount_GeoMx", pt.size = 0.1)

demoSeurat <- as.Seurat(norm_target_demoData, normData = "q_norm", ident = "cell_line")
VlnPlot(demoSeurat, features = "nCount_GeoMx", pt.size = 0.1)

Simple GeoMx data workflow

Here is an example of a typical dimensional reduction workflow.

demoSeurat <- FindVariableFeatures(demoSeurat)
demoSeurat <- ScaleData(demoSeurat)
demoSeurat <- RunPCA(demoSeurat, assay = "GeoMx", verbose = FALSE)
demoSeurat <- FindNeighbors(demoSeurat, reduction = "pca", dims = seq_len(30))
demoSeurat <- FindClusters(demoSeurat, verbose = FALSE)
demoSeurat <- RunUMAP(demoSeurat, reduction = "pca", dims = seq_len(30))

DimPlot(demoSeurat, reduction = "umap", label = TRUE, group.by = "cell_line")

In-depth GeoMx data workflow

Here is a work through of a more indepth DSP dataset. This is a non-small cell lung cancer (nsclc) tissue sample that has an ROI strategy to simulate a visium dataset (55 um circles evenly spaced apart). It was segmented on tumor and non-tumor.

data("nsclc", package = "SpatialDecon")
nsclc
## NanoStringGeoMxSet (storageMode: lockedEnvironment)
## assayData: 1700 features, 199 samples 
##   element names: exprs, exprs_norm 
## protocolData
##   sampleNames: ROI01Tumor ROI01TME ... ROI100TME (199 total)
##   varLabels: Mask.type Raw ... hkFactors (17 total)
##   varMetadata: labelDescription
## phenoData
##   sampleNames: ROI01Tumor ROI01TME ... ROI100TME (199 total)
##   varLabels: Sample_ID Tissue ... istumor (10 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: ABCF1 ABL1 ... LAG3 (1700 total)
##   fvarLabels: TargetName HUGOSymbol ... Negative (9 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation: kiloplex with cell type spike-in [legacy panel] 
## signature: none
## feature: Target
## analyte: RNA
dim(nsclc)
## Features  Samples 
##     1700      199
data.frame(exprs(nsclc)[seq_len(5), seq_len(5)])
##        ROI01Tumor ROI01TME ROI02Tumor ROI02TME ROI03Tumor
## ABCF1          55       26         47       30        102
## ABL1           21       22         27       18         47
## ACVR1B         89       30         57       29        122
## ACVR1C          9        7          4        8         14
## ACVR2A         14       15          9       12         22
head(pData(nsclc))
##                                      Sample_ID Tissue Slide.name   ROI AOI.name
## ROI01Tumor ICP20th-L11-ICPKilo-ROI01-Tumor-A02    L11    ICPKilo ROI01    Tumor
## ROI01TME     ICP20th-L11-ICPKilo-ROI01-TME-A03    L11    ICPKilo ROI01      TME
## ROI02Tumor ICP20th-L11-ICPKilo-ROI02-Tumor-A04    L11    ICPKilo ROI02    Tumor
## ROI02TME     ICP20th-L11-ICPKilo-ROI02-TME-A05    L11    ICPKilo ROI02      TME
## ROI03Tumor ICP20th-L11-ICPKilo-ROI03-Tumor-A06    L11    ICPKilo ROI03    Tumor
## ROI03TME     ICP20th-L11-ICPKilo-ROI03-TME-A07    L11    ICPKilo ROI03      TME
##            AOI.annotation    x    y nuclei istumor
## ROI01Tumor          PanCK    0 8000    572    TRUE
## ROI01TME              TME    0 8000    733   FALSE
## ROI02Tumor          PanCK  600 8000    307    TRUE
## ROI02TME              TME  600 8000    697   FALSE
## ROI03Tumor          PanCK 1200 8000    583    TRUE
## ROI03TME              TME 1200 8000    484   FALSE

When coercing, we can add the coordinate columns allowing for spatial graphing using Seurat.

nsclcSeurat <- as.Seurat(nsclc, normData = "exprs_norm", ident = "AOI.annotation", 
                         coordinates = c("x", "y"))

nsclcSeurat
## An object of class Seurat 
## 1700 features across 199 samples within 1 assay 
## Active assay: GeoMx (1700 features, 0 variable features)
##  2 layers present: counts, data
##  1 image present: image
VlnPlot(nsclcSeurat, features = "nCount_GeoMx", pt.size = 0.1)

nsclcSeurat <- FindVariableFeatures(nsclcSeurat)
nsclcSeurat <- ScaleData(nsclcSeurat)
nsclcSeurat <- RunPCA(nsclcSeurat, assay = "GeoMx", verbose = FALSE)
nsclcSeurat <- FindNeighbors(nsclcSeurat, reduction = "pca", dims = seq_len(30))
nsclcSeurat <- FindClusters(nsclcSeurat, verbose = FALSE)
nsclcSeurat <- RunUMAP(nsclcSeurat, reduction = "pca", dims = seq_len(30))

DimPlot(nsclcSeurat, reduction = "umap", label = TRUE, group.by = "AOI.name")

Spatial Graphing

Because this dataset is segmented, we need to separate the tumor and TME sections before using the spatial graphing. These Seurat functions were created for Visium data, so they can only plot the same sized circles.

Here we are showing the gene counts in each ROI separated by segment.

tumor <- suppressMessages(SpatialFeaturePlot(nsclcSeurat[,nsclcSeurat$AOI.name == "Tumor"], 
                                             features = "nCount_GeoMx", pt.size.factor = 12) + 
  labs(title = "Tumor") + 
  theme(legend.position = "none") + 
  scale_fill_continuous(type = "viridis",
                        limits = c(min(nsclcSeurat$nCount_GeoMx), 
                                   max(nsclcSeurat$nCount_GeoMx))))

TME <- suppressMessages(SpatialFeaturePlot(nsclcSeurat[,nsclcSeurat$AOI.name == "TME"], 
                                           features = "nCount_GeoMx", pt.size.factor = 12) + 
  labs(title = "TME") + 
  theme(legend.position = "right") +
  scale_fill_continuous(type = "viridis", 
                        limits = c(min(nsclcSeurat$nCount_GeoMx),
                                   max(nsclcSeurat$nCount_GeoMx))))

wrap_plots(tumor, TME)

Here we show the count for A2M

tumor <- suppressMessages(SpatialFeaturePlot(nsclcSeurat[,nsclcSeurat$AOI.name == "Tumor"], 
                                             features = "A2M", pt.size.factor = 12) + 
  labs(title = "Tumor") + 
  theme(legend.position = "none") + 
  scale_fill_continuous(type = "viridis",
                        limits = c(min(nsclcSeurat@assays$GeoMx$data["A2M",]), 
                                   max(nsclcSeurat@assays$GeoMx$data["A2M",]))))

TME <- suppressMessages(SpatialFeaturePlot(nsclcSeurat[,nsclcSeurat$AOI.name == "TME"], 
                                           features = "A2M", pt.size.factor = 12) + 
  labs(title = "TME") + 
  theme(legend.position = "right") +
  scale_fill_continuous(type = "viridis", 
                        limits = c(min(nsclcSeurat@assays$GeoMx$data["A2M",]),
                                   max(nsclcSeurat@assays$GeoMx$data["A2M",]))))

wrap_plots(tumor, TME)

Using the FindMarkers built in function from Seurat, we can determine the most differentially expressed genes in Tumor and TME

Idents(nsclcSeurat) <- nsclcSeurat$AOI.name
de_genes <- FindMarkers(nsclcSeurat, ident.1 = "Tumor", ident.2 = "TME")

de_genes <- de_genes[order(abs(de_genes$avg_log2FC), decreasing = TRUE),]
de_genes <- de_genes[is.finite(de_genes$avg_log2FC) & de_genes$p_val < 1e-25,]



for(i in rownames(de_genes)[1:2]){
  print(data.frame(de_genes[i,]))
  
  tumor <- suppressMessages(SpatialFeaturePlot(nsclcSeurat[,nsclcSeurat$AOI.name == "Tumor"], 
                                               features = i, pt.size.factor = 12) + 
  labs(title = "Tumor") + 
  theme(legend.position = "none") + 
  scale_fill_continuous(type = "viridis",
                        limits = c(min(nsclcSeurat@assays$GeoMx$data[i,]), 
                                   max(nsclcSeurat@assays$GeoMx$data[i,]))))

  TME <- suppressMessages(SpatialFeaturePlot(nsclcSeurat[,nsclcSeurat$AOI.name == "TME"], 
                                             features = i, pt.size.factor = 12) + 
    labs(title = "TME") + 
    theme(legend.position = "right") +
    scale_fill_continuous(type = "viridis", 
                          limits = c(min(nsclcSeurat@assays$GeoMx$data[i,]),
                                     max(nsclcSeurat@assays$GeoMx$data[i,]))))
  
  print(wrap_plots(tumor, TME))
}
##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## CEACAM6 1.756187e-31    754.903     1     1 2.985517e-28

##            p_val avg_log2FC pct.1 pct.2    p_val_adj
## LYZ 9.276639e-32  -615.2034     1     1 1.577029e-28

SpatialExperiment

SpatialExperiment is an S4 class inheriting from SingleCellExperiment. It is meant as a data storage object rather than an analysis suite like Seurat. Because of this, this section won’t have the fancy analysis outputs like the Seurat section had but will show where in the object all the pieces are stored.

library(SpatialExperiment)

SpatialExperiment Coercion

The three errors that can occur when trying to coerce to SpatialExperiment are:

  1. object must be on the target level
  2. object should be normalized, if you want raw data you can set forceRaw to TRUE
  3. normalized count matrix name must be valid
as.SpatialExperiment(demoData)
## Error in as.SpatialExperiment.NanoStringGeoMxSet(demoData): Data must be on Target level before converting to a SpatialExperiment Object
as.SpatialExperiment(target_demoData, normData = "exprs")
## Error in as.SpatialExperiment.NanoStringGeoMxSet(target_demoData, normData = "exprs"): It is NOT recommended to use Seurat's normalization for GeoMx data. 
##              Normalize using GeomxTools::normalize() or set forceRaw to TRUE if you want to continue with Raw data
as.SpatialExperiment(norm_target_demoData, normData = "exprs_norm")
## Error in as.SpatialExperiment.NanoStringGeoMxSet(norm_target_demoData, : The normData name "exprs_norm" is not a valid assay name. Valid names are: exprs, q_norm

After coercing to a SpatialExperiment object all of the metadata is still accessible.

demoSPE <- as.SpatialExperiment(norm_target_demoData, normData = "q_norm")

demoSPE # overall data object
## class: SpatialExperiment 
## dim: 1821 83 
## metadata(11): PKCFileName PKCModule ... sequencingMetrics QCMetrics
## assays(1): GeoMx
## rownames(1821): ACTA2 FOXA2 ... C1orf43 SNRPD3
## rowData names(6): TargetName Module ... SystematicName Negative
## colnames(83): DSP-1001250002642-A02.dcc DSP-1001250002642-A03.dcc ...
##   DSP-1001250002642-H04.dcc DSP-1001250002642-H05.dcc
## colData names(18): slide name scan name ... slide_rep sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(0) :
## imgData names(0):
data.frame(head(colData(demoSPE))) # most important ROI metadata
##                                              slide.name
## DSP-1001250002642-A02.dcc 6panel-old-slide1 (PTL-10891)
## DSP-1001250002642-A03.dcc 6panel-old-slide1 (PTL-10891)
## DSP-1001250002642-A04.dcc 6panel-old-slide1 (PTL-10891)
## DSP-1001250002642-A05.dcc 6panel-old-slide1 (PTL-10891)
## DSP-1001250002642-A06.dcc 6panel-old-slide1 (PTL-10891)
## DSP-1001250002642-A07.dcc 6panel-old-slide1 (PTL-10891)
##                                          scan.name
## DSP-1001250002642-A02.dcc cw005 (PTL-10891) Slide1
## DSP-1001250002642-A03.dcc cw005 (PTL-10891) Slide1
## DSP-1001250002642-A04.dcc cw005 (PTL-10891) Slide1
## DSP-1001250002642-A05.dcc cw005 (PTL-10891) Slide1
## DSP-1001250002642-A06.dcc cw005 (PTL-10891) Slide1
## DSP-1001250002642-A07.dcc cw005 (PTL-10891) Slide1
##                                                                                        panel
## DSP-1001250002642-A02.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
## DSP-1001250002642-A03.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
## DSP-1001250002642-A04.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
## DSP-1001250002642-A05.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
## DSP-1001250002642-A06.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
## DSP-1001250002642-A07.dcc (v1.2) VnV Cancer Transcriptome Atlas, (v1.0) Six gene test custom
##                           roi           segment     area
## DSP-1001250002642-A02.dcc   1 Geometric Segment 31318.73
## DSP-1001250002642-A03.dcc   2 Geometric Segment 31318.73
## DSP-1001250002642-A04.dcc   3 Geometric Segment 31318.73
## DSP-1001250002642-A05.dcc   4 Geometric Segment 31318.73
## DSP-1001250002642-A06.dcc   5 Geometric Segment 31318.73
## DSP-1001250002642-A07.dcc   6 Geometric Segment 31318.73
##                           NegGeoMean_Six.gene_test_v1_v1.1
## DSP-1001250002642-A02.dcc                         1.487738
## DSP-1001250002642-A03.dcc                         2.518775
## DSP-1001250002642-A04.dcc                         2.847315
## DSP-1001250002642-A05.dcc                         2.632148
## DSP-1001250002642-A06.dcc                         2.275970
## DSP-1001250002642-A07.dcc                         2.059767
##                           NegGeoMean_VnV_GeoMx_Hs_CTA_v1.2
## DSP-1001250002642-A02.dcc                         3.722751
## DSP-1001250002642-A03.dcc                         3.068217
## DSP-1001250002642-A04.dcc                         3.556275
## DSP-1001250002642-A05.dcc                         3.785600
## DSP-1001250002642-A06.dcc                         4.064107
## DSP-1001250002642-A07.dcc                         4.153701
##                           NegGeoSD_Six.gene_test_v1_v1.1
## DSP-1001250002642-A02.dcc                       1.560397
## DSP-1001250002642-A03.dcc                       1.820611
## DSP-1001250002642-A04.dcc                       1.654831
## DSP-1001250002642-A05.dcc                       2.042221
## DSP-1001250002642-A06.dcc                       1.812577
## DSP-1001250002642-A07.dcc                       1.952628
##                           NegGeoSD_VnV_GeoMx_Hs_CTA_v1.2 q_norm_qFactors
## DSP-1001250002642-A02.dcc                       1.796952       0.9391100
## DSP-1001250002642-A03.dcc                       1.806070       0.9396774
## DSP-1001250002642-A04.dcc                       1.762066       1.5324910
## DSP-1001250002642-A05.dcc                       1.793823       1.6725916
## DSP-1001250002642-A06.dcc                       1.839165       1.2351225
## DSP-1001250002642-A07.dcc                       1.626391       1.2229991
##                                        SampleID                       aoi
## DSP-1001250002642-A02.dcc DSP-1001250002642-A02 Geometric Segment-aoi-001
## DSP-1001250002642-A03.dcc DSP-1001250002642-A03 Geometric Segment-aoi-001
## DSP-1001250002642-A04.dcc DSP-1001250002642-A04 Geometric Segment-aoi-001
## DSP-1001250002642-A05.dcc DSP-1001250002642-A05 Geometric Segment-aoi-001
## DSP-1001250002642-A06.dcc DSP-1001250002642-A06 Geometric Segment-aoi-001
## DSP-1001250002642-A07.dcc DSP-1001250002642-A07 Geometric Segment-aoi-001
##                           cell_line roi_rep pool_rep slide_rep sample_id
## DSP-1001250002642-A02.dcc    HS578T       1        1         1  sample01
## DSP-1001250002642-A03.dcc    HS578T       2        1         1  sample01
## DSP-1001250002642-A04.dcc       HEL       1        1         1  sample01
## DSP-1001250002642-A05.dcc       HEL       2        1         1  sample01
## DSP-1001250002642-A06.dcc    U118MG       1        1         1  sample01
## DSP-1001250002642-A07.dcc    U118MG       2        1         1  sample01
demoSPE@metadata[1:8] # experiment data
## $PKCFileName
##            VnV_GeoMx_Hs_CTA_v1.2            Six-gene_test_v1_v1.1 
## "VnV Cancer Transcriptome Atlas"           "Six gene test custom" 
## 
## $PKCModule
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##    "VnV_GeoMx_Hs_CTA"    "Six-gene_test_v1" 
## 
## $PKCFileVersion
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                   1.2                   1.1 
## 
## $PKCFileDate
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##              "200518"              "200707" 
## 
## $AnalyteType
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                 "RNA"                 "RNA" 
## 
## $MinArea
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                 16000                 16000 
## 
## $MinNuclei
## VnV_GeoMx_Hs_CTA_v1.2 Six-gene_test_v1_v1.1 
##                   200                   200 
## 
## $shiftedByOne
## [1] TRUE
head(demoSPE@metadata$sequencingMetrics) # sequencing metrics
##                           FileVersion SoftwareVersion       Date      Plate_ID
## DSP-1001250002642-A02.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A03.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A04.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A05.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A06.dcc         0.1           1.0.0 2020-07-14 1001250002642
## DSP-1001250002642-A07.dcc         0.1           1.0.0 2020-07-14 1001250002642
##                           Well            SeqSetId    Raw Trimmed Stitched
## DSP-1001250002642-A02.dcc  A02 VH00121:3:AAAG2YWM5 646250  646250   616150
## DSP-1001250002642-A03.dcc  A03 VH00121:3:AAAG2YWM5 629241  629241   603243
## DSP-1001250002642-A04.dcc  A04 VH00121:3:AAAG2YWM5 831083  831083   798188
## DSP-1001250002642-A05.dcc  A05 VH00121:3:AAAG2YWM5 884485  884485   849060
## DSP-1001250002642-A06.dcc  A06 VH00121:3:AAAG2YWM5 781936  781936   751930
## DSP-1001250002642-A07.dcc  A07 VH00121:3:AAAG2YWM5 703034  703034   674815
##                           Aligned umiQ30 rtsQ30 DeduplicatedReads
## DSP-1001250002642-A02.dcc  610390 0.9785 0.9804            312060
## DSP-1001250002642-A03.dcc  597280 0.9784 0.9811            305528
## DSP-1001250002642-A04.dcc  791804 0.9785 0.9801            394981
## DSP-1001250002642-A05.dcc  842133 0.9796 0.9814            424162
## DSP-1001250002642-A06.dcc  744669 0.9779 0.9803            355121
## DSP-1001250002642-A07.dcc  668726 0.9776 0.9797            341008
##                                              NTC_ID NTC Trimmed (%)
## DSP-1001250002642-A02.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A03.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A04.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A05.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A06.dcc DSP-1001250002642-A01.dcc   7         100
## DSP-1001250002642-A07.dcc DSP-1001250002642-A01.dcc   7         100
##                           Stitched (%) Aligned (%) Saturated (%)
## DSP-1001250002642-A02.dcc     95.34236    94.45106      48.87531
## DSP-1001250002642-A03.dcc     95.86836    94.92071      48.84677
## DSP-1001250002642-A04.dcc     96.04191    95.27376      50.11632
## DSP-1001250002642-A05.dcc     95.99484    95.21168      49.63242
## DSP-1001250002642-A06.dcc     96.16260    95.23401      52.31156
## DSP-1001250002642-A07.dcc     95.98611    95.12001      49.00632
head(demoSPE@metadata$QCMetrics$QCFlags) # QC metrics
##                           LowReads LowTrimmed LowStitched LowAligned
## DSP-1001250002642-A02.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A03.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A04.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A05.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A06.dcc    FALSE      FALSE       FALSE      FALSE
## DSP-1001250002642-A07.dcc    FALSE      FALSE       FALSE      FALSE
##                           LowSaturation LowNegatives HighNTC LowArea
## DSP-1001250002642-A02.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A03.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A04.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A05.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A06.dcc         FALSE         TRUE   FALSE   FALSE
## DSP-1001250002642-A07.dcc         FALSE         TRUE   FALSE   FALSE
data.frame(head(rowData(demoSPE))) # gene metadata
##         TargetName                Module  CodeClass        GeneID
## ACTA2        ACTA2 VnV_GeoMx_Hs_CTA_v1.2 Endogenous            59
## FOXA2        FOXA2 VnV_GeoMx_Hs_CTA_v1.2 Endogenous          3170
## NANOG        NANOG VnV_GeoMx_Hs_CTA_v1.2 Endogenous 79923, 388112
## TRAC          TRAC VnV_GeoMx_Hs_CTA_v1.2 Endogenous          <NA>
## TRBC1/2    TRBC1/2 VnV_GeoMx_Hs_CTA_v1.2 Endogenous          <NA>
## TRDC          TRDC VnV_GeoMx_Hs_CTA_v1.2 Endogenous          <NA>
##         SystematicName Negative
## ACTA2            ACTA2    FALSE
## FOXA2            FOXA2    FALSE
## NANOG   NANOG, NANOGP8    FALSE
## TRAC              TRAC    FALSE
## TRBC1/2          TRBC1    FALSE
## TRDC              TRDC    FALSE

When coercing, we can add the coordinate columns and they will be appended to the correct location in SpatialExperiment

nsclcSPE <- as.SpatialExperiment(nsclc, normData = "exprs_norm", 
                                 coordinates = c("x", "y"))

nsclcSPE
## class: SpatialExperiment 
## dim: 1700 199 
## metadata(1): sequencingMetrics
## assays(1): GeoMx
## rownames(1700): ABCF1 ABL1 ... TNFSF4 LAG3
## rowData names(9): TargetName HUGOSymbol ... GlobalOutliers Negative
## colnames(199): ROI01Tumor ROI01TME ... ROI100Tumor ROI100TME
## colData names(20): Sample_ID Tissue ... hkFactors sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(2) : x y
## imgData names(0):
data.frame(head(spatialCoords(nsclcSPE)))
##               x    y
## ROI01Tumor    0 8000
## ROI01TME      0 8000
## ROI02Tumor  600 8000
## ROI02TME    600 8000
## ROI03Tumor 1200 8000
## ROI03TME   1200 8000

With the coordinates and the metadata, we can create spatial graphing figures similar to Seurat’s. To get the same orientation as Seurat we must swap the axes and flip the y.

figureData <- as.data.frame(cbind(colData(nsclcSPE), spatialCoords(nsclcSPE)))

figureData <- cbind(figureData, A2M=as.numeric(nsclcSPE@assays@data$GeoMx["A2M",]))

tumor <- ggplot(figureData[figureData$AOI.name == "Tumor",], aes(x=y, y=-x, color = A2M))+
  geom_point(size = 6)+
  scale_color_continuous(type = "viridis",
                        limits = c(min(figureData$A2M), 
                                   max(figureData$A2M)))+
  theme(legend.position = "none", panel.grid = element_blank(),
        panel.background = element_rect(fill = "white"),
        axis.title = element_blank(), axis.text = element_blank(), 
        axis.ticks = element_blank(), axis.line = element_blank())+
  labs(title = "Tumor")


TME <- ggplot(figureData[figureData$AOI.name == "TME",], aes(x=y, y=-x, color = A2M))+
  geom_point(size = 6)+
  scale_color_continuous(type = "viridis",
                        limits = c(min(figureData$A2M), 
                                   max(figureData$A2M))) +
  theme(panel.grid = element_blank(), 
        panel.background = element_rect(fill = "white"), axis.title = element_blank(), 
        axis.text = element_blank(), axis.ticks = element_blank(), axis.line = element_blank())+
  labs(title = "TME")

wrap_plots(tumor, TME)

Image Overlays

The free-handed nature of Region of Interest (ROI) selection in a GeoMx experiment makes visualization on top of the image difficult in packages designed for different data. We created SpatialOmicsOverlay specifically to visualize and analyze these types of ROIs in a GeoMx experiment and the immunofluorescent-guided segmentation process.

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] SpatialExperiment_1.15.1    SingleCellExperiment_1.27.2
##  [3] SummarizedExperiment_1.35.1 GenomicRanges_1.57.1       
##  [5] GenomeInfoDb_1.41.1         IRanges_2.39.2             
##  [7] MatrixGenerics_1.17.0       matrixStats_1.3.0          
##  [9] patchwork_1.2.0             SpatialDecon_1.15.0        
## [11] Seurat_5.1.0                SeuratObject_5.0.2         
## [13] sp_2.1-4                    ggiraph_0.8.10             
## [15] EnvStats_3.0.0              GeomxTools_3.9.0           
## [17] NanoStringNCTools_1.13.0    ggplot2_3.5.1              
## [19] S4Vectors_0.43.2            Biobase_2.65.1             
## [21] BiocGenerics_0.51.0         rmarkdown_2.28             
## 
## loaded via a namespace (and not attached):
##   [1] RcppAnnoy_0.0.22        splines_4.4.1           later_1.3.2            
##   [4] R.oo_1.26.0             tibble_3.2.1            cellranger_1.1.0       
##   [7] polyclip_1.10-7         fastDummies_1.7.4       lifecycle_1.0.4        
##  [10] globals_0.16.3          lattice_0.22-6          MASS_7.3-61            
##  [13] magrittr_2.0.3          plotly_4.10.4           sass_0.4.9             
##  [16] jquerylib_0.1.4         yaml_2.3.10             httpuv_1.6.15          
##  [19] sctransform_0.4.1       spam_2.10-0             spatstat.sparse_3.1-0  
##  [22] reticulate_1.38.0       cowplot_1.1.3           pbapply_1.7-2          
##  [25] buildtools_1.0.0        minqa_1.2.8             RColorBrewer_1.1-3     
##  [28] abind_1.4-5             zlibbioc_1.51.1         R.cache_0.16.0         
##  [31] Rtsne_0.17              R.utils_2.12.3          purrr_1.0.2            
##  [34] GenomeInfoDbData_1.2.12 ggrepel_0.9.5           irlba_2.3.5.1          
##  [37] listenv_0.9.1           spatstat.utils_3.1-0    maketools_1.3.0        
##  [40] pheatmap_1.0.12         goftest_1.2-3           RSpectra_0.16-2        
##  [43] spatstat.random_3.3-1   fitdistrplus_1.2-1      parallelly_1.38.0      
##  [46] DelayedArray_0.31.11    leiden_0.4.3.1          codetools_0.2-20       
##  [49] tidyselect_1.2.1        farver_2.1.2            UCSC.utils_1.1.0       
##  [52] lme4_1.1-35.5           spatstat.explore_3.3-2  jsonlite_1.8.8         
##  [55] progressr_0.14.0        ggridges_0.5.6          survival_3.7-0         
##  [58] systemfonts_1.1.0       tools_4.4.1             ica_1.0-3              
##  [61] Rcpp_1.0.13             glue_1.7.0              SparseArray_1.5.31     
##  [64] gridExtra_2.3           xfun_0.47               ggthemes_5.1.0         
##  [67] dplyr_1.1.4             withr_3.0.1             numDeriv_2016.8-1.1    
##  [70] fastmap_1.2.0           GGally_2.2.1            repmis_0.5             
##  [73] boot_1.3-31             fansi_1.0.6             digest_0.6.37          
##  [76] R6_2.5.1                mime_0.12               colorspace_2.1-1       
##  [79] scattermore_1.2         tensor_1.5              spatstat.data_3.1-2    
##  [82] R.methodsS3_1.8.2       utf8_1.2.4              tidyr_1.3.1            
##  [85] generics_0.1.3          data.table_1.16.0       S4Arrays_1.5.7         
##  [88] httr_1.4.7              htmlwidgets_1.6.4       ggstats_0.6.0          
##  [91] uwot_0.2.2              pkgconfig_2.0.3         gtable_0.3.5           
##  [94] lmtest_0.9-40           XVector_0.45.0          sys_3.4.2              
##  [97] htmltools_0.5.8.1       dotCall64_1.1-1         scales_1.3.0           
## [100] png_0.1-8               logNormReg_0.5-0        spatstat.univar_3.0-0  
## [103] knitr_1.48              reshape2_1.4.4          rjson_0.2.22           
## [106] uuid_1.2-1              nlme_3.1-166            nloptr_2.1.1           
## [109] cachem_1.1.0            zoo_1.8-12              stringr_1.5.1          
## [112] KernSmooth_2.23-24      parallel_4.4.1          miniUI_0.1.1.1         
## [115] vipor_0.4.7             pillar_1.9.0            grid_4.4.1             
## [118] vctrs_0.6.5             RANN_2.6.2              promises_1.3.0         
## [121] xtable_1.8-4            cluster_2.1.6           beeswarm_0.4.0         
## [124] evaluate_0.24.0         magick_2.8.4            cli_3.6.3              
## [127] compiler_4.4.1          rlang_1.1.4             crayon_1.5.3           
## [130] future.apply_1.11.2     labeling_0.4.3          plyr_1.8.9             
## [133] ggbeeswarm_0.7.2        stringi_1.8.4           viridisLite_0.4.2      
## [136] deldir_2.0-4            lmerTest_3.1-3          munsell_0.5.1          
## [139] Biostrings_2.73.1       lazyeval_0.2.2          spatstat.geom_3.3-2    
## [142] Matrix_1.7-0            RcppHNSW_0.6.0          future_1.34.0          
## [145] shiny_1.9.1             highr_0.11              ROCR_1.0-11            
## [148] igraph_2.0.3            bslib_0.8.0             readxl_1.4.3