Data exploration

As example dataset, we use the one by Vicari et al. (2024), which includes replicated measurements of spatial transcriptomics and metabolomics coprofiled on mouse brain sections. A subset of the data is available in the sbivar package and can be loaded as:

data(Vicari)
names(Vicari)

## [1] "TranscriptCoords"   "TranscriptOutcomes" "MetaboliteCoords"  
## [4] "MetaboliteOutcomes"

This object consists of four lists of length six: TranscriptCoords and MetaboliteCoords for the coordinates, and TranscriptOutcomes and MetaboliteOutcomes for the feature measurements. Only the 5 most abundant features per modality are included in the package, for real datasets we also recommend to filter on the most abundant features, or apply the method to the entire dataset. The sbivar package assumes that the data have been pre-aligned, hence it is always a good idea to check the alignment visually:

par(mfrow = c(2, 3))
plotCoordsMulti(Vicari$TranscriptCoords, Vicari$MetaboliteCoords, cex = 0.1)

par(mfrow = c(1, 1))

The alignment seems successful, notice that some tissue areas are only covered by one modality though.

Single-image analysis

The Vicari data consist of six images which can be analysed jointly, but for didactical purposes we first analyse a single image here, the sample “V11T16-085_A1”.

singleSample <- "V11T16-085_A1"
singleStxCoords <- Vicari$TranscriptCoords[[singleSample]]
singleStx <- Vicari$TranscriptOutcomes[[singleSample]]
singleMetCoords <- Vicari$MetaboliteCoords[[singleSample]]
singleMet <- Vicari$MetaboliteOutcomes[[singleSample]]

We analyse this single sample using bivariate Moran’s I. We choose to normalize both modalities to relative abundances (see ?normMat).

moransIRes <- sbivar(singleStx, singleMet, singleStxCoords, singleMetCoords,
    method = "Moran's I", normX = "rel", normY = "rel"
)

## Starting sbivar analysis of a single image on 2 computing cores

## Testing significance of bivariate Moran's I for 25 feature pairs

## Calculating bivariate Moran's I statistics ...

## Fitting variograms for first modality (5 features) ...

## Fitting variograms for second modality (5 features) ...

## Calculating variances of bivariate Moran's I statistics ...

Have a look at the results, which are sorted by increasing p-value:

head(moransIRes$result)

##   Modality_X Modality_Y     Ixy_5e.06     Ixy_2e.04      Ixy_0.02 SE.Ixy._5e.06
## 1      Sparc    Taurine -1.275865e-04 -8.789649e-05 -4.676381e-05  3.830270e-05
## 2      Sparc       GABA  1.243888e-04  1.093865e-04  6.397515e-05  4.263167e-05
## 3       Fth1  Histidine -9.796108e-05 -7.611385e-05 -2.593648e-05  4.979809e-05
## 4      Sparc   Dopamine  8.103384e-05  5.156486e-05  1.695850e-05  4.482058e-05
## 5       Gnas    Taurine -3.533777e-05 -1.832929e-05 -6.049097e-06  2.143787e-05
## 6    mt.Atp6  Histidine  8.383458e-05  6.505115e-05  2.256473e-05  5.304026e-05
##   SE.Ixy._2e.04 SE.Ixy._0.02        pVal       pAdj
## 1  3.603835e-05 3.095988e-05 0.002437642 0.06094105
## 2  4.060163e-05 3.502142e-05 0.006819546 0.08524432
## 3  4.810448e-05 4.309704e-05 0.102890944 0.76609491
## 4  4.270522e-05 3.607202e-05 0.168602858 0.76609491
## 5  1.570091e-05 1.304224e-05 0.218341722 0.76609491
## 6  5.152161e-05 4.674568e-05 0.222927571 0.76609491

Plot the most significantly spatially associated gene-metabolite pair :

plotTopPair(moransIRes,
    X = singleStx, Y = singleMet, Cx = singleStxCoords,
    Ey = singleMetCoords
)

We see a clear negative relationship.

In addition, we employ generalized additive models (GAMs), with a negative binomial outcome distribution for the transcriptome data, and a gamma outcome distribution for the metabolome data. A log-link is used in both cases. For non-normal data, such as spatial transcriptomics and metabolomics data, this approach is usually preferable to the bivariate Moran’s I, as this latter one fails to account for heteroskedasticity,

FamiliesVicari <- list("X" = mgcv::nb(), "Y" = Gamma(lin = "log"))
gamRes <- sbivar(singleStx, singleMet, singleStxCoords, singleMetCoords,
    method = "GAMs", families = FamiliesVicari
)

## Starting sbivar analysis of a single image on 2 computing cores

## Fitting GAMs for first modality (5 features) ...

## Fitting GAMs for second modality (5 features) ...

## Performing 25 pairwise tests on fitted GAMs ...

Plot the most significant pair. Normalized data are shown, while the original analysis is on raw data, but we scale the size of the spots by the sample sums to reflect their higher weight in the analysis:

plotTopPair(gamRes,
    X = singleStx, Y = singleMet, Cx = singleStxCoords, scaleBySampleSums = TRUE,
    Ey = singleMetCoords, normX = "rel", normY = "rel"
)

An other feature is found most significant by GAMs, now with a positive relationship. For GAMs, we can also plot the corresponding spline fit, with the contributions to the correlation 0.291:

plotGAMsTopResults(gamRes,
    X = singleStx, Y = singleMet, Cx = singleStxCoords,
    Ey = singleMetCoords
)

Another available analysis for single images are bivariate Gaussian processes (GPs). Yet fitting GPs and calculating the score tests statistics can be computation and memory intensive! Hence the code is shown, but not run:

gpRes <- sbivar(singleStx, singleMet, singleStxCoords, singleMetCoords, method = "GPs", normX = "rel", normY = "rel")
head(gpRes$result)

The analysis results can be written to an external spreadsheet using the following one-liner, e.g. for the GAM results.

writeSbivarToXlsx(gamRes, file = tmpFile <- tempfile(fileext = ".xlsx"))

Multi-image analysis

Next, we analyse the six images jointly. Let’s look at their names:

names(Vicari$TranscriptOutcomes)

## [1] "V11L12-109_A1" "V11L12-109_B1" "V11L12-109_C1" "V11T16-085_A1"
## [5] "V11T16-085_B1" "V11T16-085_C1"

We see that the first three and the last three are replicates from different mice (V11L12-109 and V11T16-085). We construct a variable identifying the mouse, consisting of the first 10 characters of the names:

mouse <- substr(names(Vicari$TranscriptOutcomes), 1, 10)

For the multi-image case, we can also use bivariate Moran’s I as measure of spatial association, here with a nearest-neighbour (nn) weighting scheme for illustration (if the required RANN package is available). For computational reasons, we limit calculations to the first 500 observations:

VicariMultiTest <- lapply(Vicari, function(x) lapply(x, function(y) y[1:500, ]))
multiMoranRes <- sbivar(VicariMultiTest$TranscriptOutcomes, VicariMultiTest$MetaboliteOutcomes,
    normX = "rel", normY = "rel",
    VicariMultiTest$TranscriptCoords, VicariMultiTest$MetaboliteCoords,
    method = "Moran", wo = if (require(RANN)) "nn" else "Gauss"
)

## Loading required package: RANN

## Starting sbivar analysis (Moran's I) of 6 images on 2 computing cores

## Image 1 of 6

## Image 2 of 6

## Image 3 of 6

## Image 4 of 6

## Image 5 of 6

## Image 6 of 6

Note that for this analysis, we only calculate the bivariate Moran’s I statistics and their maxima, but not their variances, as we will quantify variability across replicates. Next we plug the calculated Moran’s I values into a linear model, with random effects for the individual mice, and adding a bogus fixed variable for illustrative purposes:

designDf <- data.frame("mouse" = mouse, "bogus" = factor(c(1, 2, 1, 2, 1, 2)))
multiMoranLmms <- fitLinModels(multiMoranRes, designDf, Formula = ~ bogus + (1 | mouse))

## Fitting 25 mixed effects models on 2 cores

Extract the results and show the results for the intercept, i.e. the overall effect

multiMoranLmmsRes <- extractResultsMulti(multiMoranLmms, designDf)
head(multiMoranLmmsRes$result$Intercept)

##   Modality_X Modality_Y      numNN_4      numNN_8     numNN_24       pVal
## 1       Gnas Tocopherol  0.008676723  0.010162118  0.013522596 0.01834492
## 2    mt.Atp6 Tocopherol  0.010649459  0.011038879  0.002456757 0.04915305
## 3       Fth1    Taurine  0.008847759  0.008113868 -0.007303785 0.06798428
## 4    mt.Atp6    Taurine -0.005607521 -0.003873033  0.011118747 0.14156282
## 5       Pcp4    Taurine -0.009708262 -0.011691945 -0.007060228 0.16938645
## 6       Fth1 Tocopherol -0.012294981 -0.013746067 -0.006134808 0.32288705
##        pAdj
## 1 0.4586231
## 2 0.5665357
## 3 0.5665357
## 4 0.8469322
## 5 0.8469322
## 6 0.9880843

No features are significantly associated after multiplicity correction. This is in part caused by different experimental procedures used for the different sections, i.e. these are not real replicates yet, as technology is still being developed. For illustration, we plot the feature pair with the smallest p-values nevertheless:

plotTopPair(multiMoranLmmsRes,
    Xl = Vicari$TranscriptOutcomes, Yl = Vicari$MetaboliteOutcomes,
    Cxl = Vicari$TranscriptCoords, Eyl = Vicari$MetaboliteCoords, size = 0.4
)

Again, the results of the multi-image analysis can be written to an external spreadsheet. Setting sigLevel to 1 means all feature pairs are written to the file:

writeSbivarToXlsx(multiMoranLmmsRes, file = "myfile2.xlsx", sigLevel = 1)

SpatialExperiment and MultiAssayExperiment

The SpatialExperiment and MultiAssayExperiment classes are used in Bioconductor to store spatial (multi) omics experiments. The sbivar package provides functionality for doing analysis on these objects directly, for the general case of disjoint coordinate sets between the modalities.

library(SpatialExperiment, quietly = TRUE)
StxSE <- SpatialExperiment(
    assays = list("transcripts" = t(singleStx)),
    spatialCoords = singleStxCoords
)
MetSE <- SpatialExperiment(
    assays = list("metabolites" = t(singleMet)),
    spatialCoords = singleMetCoords
)

StxSE

## class: SpatialExperiment 
## dim: 5 2460 
## metadata(0):
## assays(1): transcripts
## rownames(5): Fth1 mt.Atp6 Sparc Gnas Pcp4
## rowData names(0):
## colnames(2460): AAACAACGAATAGTTC-1 AAACAAGTATCTCCCA-1 ...
##   TTGTTGTGTGTCAAGA-1 TTGTTTGTGTAAATTC-1
## colData names(1): sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(2) : x y
## imgData names(0):

resSpatExp <- sbivar(StxSE, MetSE,
    assayX = "transcripts", assayY = "metabolites",
    families = FamiliesVicari, method = "GAMs", verbose = FALSE
)

plotTopPair(resSpatExp, X = StxSE, Y = MetSE, topRank = 2, size = 0.4)

library(MultiAssayExperiment)
MAE <- MultiAssayExperiment(experiments = list(
    "Stx" = StxSE,
    "Msi" = MetSE
))
MAE

## A MultiAssayExperiment object of 2 listed
##  experiments with user-defined names and respective classes.
##  Containing an ExperimentList class object of length 2:
##  [1] Stx: SpatialExperiment with 5 rows and 2460 columns
##  [2] Msi: SpatialExperiment with 5 rows and 3278 columns
## Functionality:
##  experiments() - obtain the ExperimentList instance
##  colData() - the primary/phenotype DataFrame
##  sampleMap() - the sample coordination DataFrame
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment
##  *Format() - convert into a long or wide DataFrame
##  assays() - convert ExperimentList to a SimpleList of matrices
##  exportClass() - save data to flat files

resMAE <- sbivar(MAE,
    assayX = "transcripts", assayY = "metabolites",
    experimentX = "Stx", experimentY = "Msi",
    families = FamiliesVicari, method = "GAMs", normX = "rel", normY = "rel"
)

StxSElist <- lapply(selfName(names(Vicari$TranscriptOutcomes)), function(nam) {
    SpatialExperiment(
        assays = list("transcripts" = t(Vicari$TranscriptOutcomes[[nam]])),
        spatialCoords = Vicari$TranscriptCoords[[nam]]
    )
})
MetSElist <- lapply(selfName(names(Vicari$MetaboliteOutcomes)), function(nam) {
    SpatialExperiment(
        assays = list("metabolites" = t(Vicari$MetaboliteOutcomes[[nam]])),
        spatialCoords = Vicari$MetaboliteCoords[[nam]]
    )
})
multiMoranResSE <- sbivar(
    X = StxSElist, Y = MetSElist, assayX = "transcripts",
    assayY = "metabolites", method = "GAM", verbose = FALSE
)
multiMoranLmmsSE <- fitLinModels(multiMoranResSE, designDf, Formula = ~ (1 | mouse))

## Fitting 25 mixed effects models on 2 cores

multiMoranLmmsResSE <- extractResultsMulti(multiMoranLmmsSE, designDf)
head(multiMoranLmmsResSE$result$Intercept)

##   Modality_X Modality_Y   Estimate        SE      pVal     pAdj
## 1    mt.Atp6 Tocopherol -0.3930271 0.2157481 0.1281296 0.960651
## 2       Fth1   Dopamine -0.3677403 0.2399068 0.1858886 0.960651
## 3    mt.Atp6  Histidine -0.4259019 0.2378129 0.3241990 0.960651
## 4    mt.Atp6    Taurine -0.4457374 0.2640538 0.3404712 0.960651
## 5      Sparc Tocopherol  0.4037679 0.2700480 0.3752822 0.960651
## 6    mt.Atp6       GABA -0.3386076 0.3152997 0.4773179 0.960651

SEsingleStx <- SpatialExperiment(
    assays = list("transcripts" = t(Reduce(Vicari$TranscriptOutcomes, f = rbind))),
    spatialCoords = Reduce(Vicari$TranscriptCoords, f = rbind),
    colData = data.frame("sample_id" = rep(names(Vicari$TranscriptOutcomes),
        times = vapply(Vicari$TranscriptOutcomes, FUN.VALUE = 0L, nrow)
    ))
)
SEsingleMet <- SpatialExperiment(
    assays = list("metabolites" = t(Reduce(Vicari$MetaboliteOutcomes, f = rbind))),
    spatialCoords = Reduce(Vicari$MetaboliteCoords, f = rbind),
    colData = data.frame("sample_id" = rep(names(Vicari$MetaboliteOutcomes),
        times = vapply(Vicari$MetaboliteOutcomes, FUN.VALUE = 0L, nrow)
    ))
)

SEsingleStx

## class: SpatialExperiment 
## dim: 5 16064 
## metadata(0):
## assays(1): transcripts
## rownames(5): Fth1 mt.Atp6 Sparc Gnas Pcp4
## rowData names(0):
## colnames(16064): AAACAAGTATCTCCCA-1 AAACATTTCCCGGATT-1 ...
##   TTGTTCTAGATACGCT-1 TTGTTGTGTGTCAAGA-1
## colData names(1): sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
## spatialCoords names(2) : x y
## imgData names(0):

sbivarMultiSingleSE <- sbivar(SEsingleStx, SEsingleMet,
    sample_id_x = "sample_id",
    assayX = "transcripts", assayY = "metabolites",
    method = "GAMs", normX = "rel", normY = "rel"
)
multiMoranLmmsSingleSE <- fitLinModels(sbivarMultiSingleSE, designDf,
    Formula = ~ (1 | mouse)
)
multiMoranLmmsResSingleSE <- extractResultsMulti(multiMoranLmmsSingleSE, designDf)

AnnData

AnnData is python data format often used to store spatial experiment results, which can also be read into R. An AnnData object is e.g. created by reading an ‘.h5ad’ with the functionality of the anndata R-package. Here we adapted one a file from the study by Godfrey et al. (2025), keeping its structure but subsetting it to only necessary components and the 10 most abundant features per modality. For illustration, this file is included in the sbivar package and can be read in (if the anndata and zellkonverter packages are installed) as:

if (reqadzk <- (require(anndata) && require(zellkonverter))) {
    annDataObj <- read_h5ad(system.file("h5adData", "LC_170_subset.h5ad", package = "sbivar"))
}

We suggest converting it to a SingleCellExperiment object from the eponymous R-package using the zellkonverter package:

if (reqadzk) {
    sceStx <- AnnData2SCE(annDataObj, X_name = "X")
    reducedDims(sceStx) <- annDataObj$obsm
}

## For native R and reading and writing of H5AD files, an R <AnnData> object, and
## conversion to <SingleCellExperiment> or <Seurat> objects, check out the
## anndataR package:
## ℹ Install it from Bioconductor with `BiocManager::install("anndataR")`
## ℹ See more at <https://bioconductor.org/packages/anndataR/>
## This message is displayed once per session.

## Warning: The passed object is a 'AnnDataR6' object, conversion is likely to be less
## reliable

Note that a SingleCellExperiment object can only contain one modality at the time, here the transcriptomics. From there on, it is only a small step to a SpatialExperiment object as accepted by sbivar. The only tricky part are the spatial coordinates, which have no fixed slot in the AnnData object. Here, they have been stored in the “obsm” slot of the AnnData object, which means they have ended up in the “reducedDims” slot of the SingleCellExperiment object. In addition, we create a SpatialExperiment object for the metabolomics data, which has oddly been stored in the unstructured (“uns”) slot of the AnnData object. To keep build times of this vignette low, we again only look at a subset of 10 features

if (reqadzk) {
    library(SpatialExperiment)
    speStx <- SpatialExperiment(list("Stx" = assay(sceStx, "X")[1:10, ]),
        colData = colData(sceStx),
        spatialCoords = reducedDims(sceStx)$spatial
    )
    # Now for the metabolomics data too
    speMet <- SpatialExperiment(list("msi" = t(annDataObj$uns$msi)[1:10, ]),
        colData = colData(sceStx),
        spatialCoords = reducedDims(sceStx)$spatial
    )
    rownames(speMet) <- annDataObj$uns[["mz_features"]][1:10]
    # Need to set column names manually
}

Note that in this example, the samples of the two modalities have been pre-matched (the “spatialCoords” slot is the same), which is in general not recommended. The above code chunks should not be applied thoughtlessly, as the AnnData structure is not very strict on the data types it includes. The analysis is then run and plotted as:

if (reqadzk) {
    GodfreyRes <- sbivar(speStx, speMet,
        assayX = "Stx", assayY = "msi",
        normX = "rel", normY = "rel", verbose = FALSE
    )
    plotTopPair(GodfreyRes, speStx, speMet)
}