Multi-omics
Multi-omics Analysis
VoltRon is a end-to-end spatial multi-omics platform that enables both automated and manual data alignment workflows to synchronize coordinate frameworks across spatial omics assays. However, upon harmonization of the spatial observations, VoltRon can be used to further transfer and integrate information across omic datasets captured on same and adjacent tissue sections. These assays are not just limited to datasets with single cell and spot resolutions but also to assays solely composed of molecule and image tile observations.
Below, we demonstrate data integration capabilities of VoltRon across distinct modalities using several use cases. In the Spatial Data Alignment, we have highlighted several scenarios to align spatial assays within tissue blocks. However in this tutorial, we demonstrate analysis starting with data alignment and further downstream data integration.
Xenium COVID-19 Lung Data
Xenium in situ platform is compatible with H&E and immunofluorescence imaging , thus researchers can stain whole slides post-Xenium to capture additional information on the morphology of the tissue and integrate it with the Xenium readouts. In this usecase, we would like to integrate the post-Xenium H&E image of tissue micro array sections composed of control and COVID-19 lung biopsy punches. Our objective is to asses the localization of SARS-COV-2 virus molecules by integrating pathology image annotations.
VoltRon is capable of analyzing molecule-level subcellular datasets independent of single cells, and specifically those that may also found outside of these cells. The Xenium run was performed with custom Xenium probes designed to hybridize distinct open reading frames of SARS-COV-2 virus molecules. These subcellular entities which then can be detected both within and outside of cells which allows to understand proliferation mechanics of the virus across the tissue. By aligning the image and the associated pathology annotations of the postXenium H&E of this slide, we can transfer labels of annotations to Xenium cells and molecules, thus inferring the abundance of viral RNA further.
We analyse readouts of eight tissue micro array (TMA) tissue sections that were fitted into the scan area of a slide loaded on a Xenium in situ instrument. These sections were cut from control and COVID-19 lung tissues of donors categorized based on disease durations (acute and prolonged). You can download the standard Xenium output folder here. For more information on the TMA blocks, see GSE190732 for more information.
Single Cells and Molecules
Across these eight TMA section, we investigate a section of acute case which is originated from a lung with extreme number of detected open reading frames of virus molecules. For convenience, we load a VoltRon object where cells are already annotated. You can also find the RDS file here.
vr2_merged_acute1 <- readRDS(file = "acutecase1_annotated.rds")Lets visualize both the UMAP and Spatial plot of the annotated cells.
vrEmbeddingPlot(vr2_merged_acute1, assay = "Xenium", embedding = "umap",
group.by = "CellType", label = TRUE)
vrSpatialPlot(vr2_merged_acute1, assay = "Xenium", group.by = "CellType",
plot.segments = TRUE)
|
|
We incorporate two open reading frames (ORFs), named S2_N and S2_orf1ab, which represent unreplicated and replicated virus molecules, respectively. We can visualize again tile the count of these virus expressions by specifically selecting these two ORFs.
vrSpatialPlot(vr2_merged_acute1, assay = "Xenium_mol", group.by = "gene",
group.ids = c("S2_N", "S2_orf1ab"), n.tile = 500)
We can even zoom into specific locations at the tissue to investigate virus particles more closely by droping the n.tile argument and calling interactive visualization.
vrSpatialPlot(vr2_merged_acute1, assay = "Xenium_mol", group.by = "gene",
group.ids = c("S2_N", "S2_orf1ab"), interactive = TRUE, pt.size = 0.1)
In the following sections, we will integrate pathological images and annotations with Xenium datasets to further understand the spatial localization of virus ORF types over the tissue.
Hot Spot Analysis
VoltRon platform allows users to find hot spots of several types of spatial entities, for spots, cells, and even molecules. We first learn spatial neighborhoods from molecules of interests, in this case, the extracellular virus particles and their ORFs.
vr2_merged_acute1VoltRon Object
acute case 1:
Layers: Section1
Assays: Xenium(Main) Xenium_mol We switch to the molecule assay of the VoltRon object, and select virus ORFs. We also look for other virus ORFs that are 50 pixel distance from each virus molecule to pin point neighborhoods that are rich in virus.
vr2_merged_acute1 <- getSpatialNeighbors(vr2_merged_acute1, assay = "Xenium_mol",
group.by = "gene", group.ids = c("S2_N", "S2_orf1ab"),
method = "radius", radius = 50)We can now observe the new spatial neighborhood graph saved in the VoltRon object with name radius.
vrGraphNames(vr2_merged_acute1)[1] "radius"We now select the feature type and graph name to run an hot spot analysis which will label each molecule if their neighborhood are dense in other virus molecules.
vr2_merged_acute1 <- getHotSpotAnalysis(vr2_merged_acute1, assay = "Xenium_mol",
features = "gene", graph.type = "radius")
vrSpatialPlot(vr2_merged_acute1, assay = "Xenium_mol",
group.by = "gene_hotspot_flag", group.ids = c("cold", "hot"),
alpha = 1, colors = list(cold = "blue", hot = "red"),
background.color = "white")
Automated H&E Registration
VoltRon can analyze and also integrate information from distinct spatial data types such as images, annotations (as regions of interests, i.e. ROIs) and molecules independently. Using such advanced utilities, we can make use of histological information and generate new metadata level information for molecule datasets.
We will first import both histological images and manual annotations using the importImageData function which accepts both images to generate tile/pixel level datasets but also allows one to import either a list of segments or GeoJSON objects for create ROI-level datasets as separate assays in a single VoltRon layer. The .geojson file was generated using QuPath on the same section H&E image of one Xenium section with the acute COVID-19 case. We also have to flip the coordinates of ROI annotations also for they were directly imported from QuPath which incorporates a reverse coordinate system on the y-axis.
Once imported, the resulting VoltRon object will have two assays in a single layer, one for tile dataset of the H&E image and the other for ROI based annotations of again the same image. You can download the H&E image from here, and download the json file from here.
# get image
imgdata <- importImageData("acutecase1_HE.jpg",
segments = "acutecase1_membrane.geojson",
sample_name = "acute case 1 (HE)")
imgdata <- flipCoordinates(imgdata, assay = "ROIAnnotation")
imgdataVoltRon Object
acute case 1 (HE):
Layers: Section1
Assays: ImageData(Main) ROIAnnotation
Features: main(Main) By visualizing these annotations interactively, we can see that the ROIs point to the hyaline membranes. Here, we likely find extra-cellular SARS-COV-2 molecules mostly outside of any single cell.
vrSpatialPlot(imgdata, assay = "ROIAnnotation", group.by = "Sample",
alpha = 0.7, interactive = TRUE)
At a first glance, although the Xenium (DAPI) and H&E images are associated with the same TMA core, they were captured in a different perspective; that is, one image is almost the 90 degree rotated version of the other. We will account for this rotation when we (automatically) align the Xenium data with the H&E image.
vr2_merged_acute1 <- modulateImage(vr2_merged_acute1, brightness = 300, channel = "DAPI")
vrImages(vr2_merged_acute1, assay = "Assay7", scale.perc = 20)
vrImages(imgdata, assay = "Assay1", scale.perc = 20)
We now register the H&E image and annotations to the DAPI image of Xenium section using the registerSpatialData function. We select FLANN (with “Homography” method) automated registration mode, negate the DAPI image of the Xenium slide, turn 90 degrees to the left and set the scale parameter of both images to width = 1859. See Spatial Data Alignment tutorial for more information.
xen_reg <- registerSpatialData(object_list = list(vr2_merged_acute1, imgdata))
Once registered, we can isolate the registered H&E data and use it further analysis. We can transfer images and their channels across assays of one or multiple VoltRon objects. Here, we save the registered image of the H&E data as an additional channel of Xenium section of the acuse case 1 sample with molecule data.
imgdata_reg <- xen_reg$registered_spat[[2]]
vrImages(vr2_merged_acute1[["Assay7"]], name = "main", channel = "H&E") <-
vrImages(imgdata_reg, assay = "Assay1", name = "main_reg")
vrImages(vr2_merged_acute1[["Assay8"]], name = "main", channel = "H&E") <-
vrImages(imgdata_reg, assay = "Assay1", name = "main_reg")We can now observe the new channels available for the both molecule and cell-level assays of Xenium data.
vrImageChannelNames(vr2_merged_acute1) Assay Layer Sample Spatial Channels
Assay7 Xenium Section1 acute case 1 main DAPI,H&E
Assay8 Xenium_mol Section1 acute case 1 main DAPI,H&EYou can also add the VoltRon object of H&E data as an additional assay of the Xenium section such that one layer includes cell, molecules, ROI Annotations and images in the same time. Specifically, we add the ROI annotation to the Xenium VoltRon object using the addAssay function where we choose the destination sample/block and the layer of the assay.
vr2_merged_acute1 <- addAssay(vr2_merged_acute1,
assay = imgdata_reg[["Assay2"]],
metadata = Metadata(imgdata_reg, assay = "ROIAnnotation"),
assay_name = "ROIAnnotation",
sample = "acute case 1", layer = "Section1")
vr2_merged_acute1VoltRon Object
acute case 1:
Layers: Section1
Assays: ROIAnnotation(Main) Xenium_mol Xenium Interactive Visualization
Once the H&E image is registered and transfered to the Xenium data, we can convert the VoltRon object into an Anndata object (h5ad file) and use TissUUmaps tool for interactive visualization.
# convert VoltRon object to h5ad
as.AnnData(vr2_merged_acute1, assay = "Xenium", file = "vr2_merged_acute1.h5ad",
flip_coordinates = TRUE, name = "main", channel = "H&E")To run TissUUmaps please follow installation instructions here, then you can simply drag and drop both the h5ad file and png file to the application.
Label Transfer
Now we can transfer ROI annotations as additional metadata features of the molecule assay. We will refer the new metadata column as “Region” which will indicate if the molecule is within any annotated ROI in the same layer.
vrMainAssay(vr2_merged_acute1) <- "ROIAnnotation"
vr2_merged_acute1$Region <- vrSpatialPoints(vr2_merged_acute1)
# set the spatial coordinate system of ROI Annotations assay
vrMainSpatial(vr2_merged_acute1[["Assay9"]]) <- "main_reg"
# transfer ROI annotations to molecules
vr2_merged_acute1 <- transferData(object = vr2_merged_acute1, from = "Assay9", to = "Assay8",
features = "Region")
# Metadata of molecules
Metadata(vr2_merged_acute1, assay = "Xenium_mol") id assay_id overlaps_nucleus gene qv Assay Layer Sample Region
<char> <char> <int> <char> <num> <char> <char> <char> <char>
1: 281651070371256_cb791e Assay8 0 ENAH 40.00000 Xenium_mol Section1 acute case 1 undefined
2: 281651070371258_cb791e Assay8 1 CD274 40.00000 Xenium_mol Section1 acute case 1 undefined
3: 281651070372515_cb791e Assay8 0 CD163 40.00000 Xenium_mol Section1 acute case 1 undefined
4: 281651070374059_cb791e Assay8 1 CTSL 33.97290 Xenium_mol Section1 acute case 1 undefined
5: 281651070374411_cb791e Assay8 0 C1S 40.00000 Xenium_mol Section1 acute case 1 undefined
---
785787: 281874408669584_cb791e Assay8 0 SFRP2 40.00000 Xenium_mol Section1 acute case 1 undefined
785788: 281874408671410_cb791e Assay8 0 S2_N 40.00000 Xenium_mol Section1 acute case 1 undefined
785789: 281874408672438_cb791e Assay8 0 S100A8 40.00000 Xenium_mol Section1 acute case 1 undefined
785790: 281874408673446_cb791e Assay8 0 TIMP1 29.86642 Xenium_mol Section1 acute case 1 undefined
785791: 281874408673882_cb791e Assay8 0 S2_N 40.00000 Xenium_mol Section1 acute case 1 undefinedThis annotations can be accessed from the default molecule level metadata of the VoltRon object.
vrMainAssay(vr2_merged_acute1) <- "Xenium_mol"
head(table(vr2_merged_acute1$Region)) ROI1_Assay9 ROI10_Assay9 ROI100_Assay9 ROI101_Assay9 ROI102_Assay9 ROI103_Assay9
583 624 784 357 215 200 Now we will grab these annotations from molecule metadata and calculate the ratio of N to orf1ab frequency of SARS-COV-2 particles across all annotated molecules.
library(dplyr)
s2_summary_hyaline <-
Metadata(vr2_merged_acute1, assay = "Xenium_mol") %>%
filter(gene %in% c("S2_N", "S2_orf1ab"),
Region != "undefined") %>%
summarise(S2_N = sum(gene == "S2_N"),
S2_orf1ab = sum(gene == "S2_orf1ab"),
ratio = sum(gene == "S2_N")/sum(gene == "S2_orf1ab")) %>%
as.matrix()
s2_summary_hyaline S2_N S2_orf1ab ratio
[1,] 50977 33532 1.520249Manual annotation
To compare the proportion of S2_N and S2_orf1ab molecules in hyaline membranes with tissue sites of possible infection. We focus on another tissue niche where a large accumulation of cells with high S2_N and S2_orf1ab counts.
By visualizing and zooming on a spatial plot with annotated cells, we can detect a site of cells with high infection which are also accompanied by a group of T cells.
vrSpatialPlot(vr2_merged_acute1, assay = "Xenium", group.by = "CellType",
group.ids = c("H.I. Cells", "T cells"), plot.segments = TRUE,
alpha = 0.6, spatial = "main", channel = "H&E", interactive = TRUE)
Now we use the annotateSpatialData function to create an additional annotation of Xenium molecules directly. By selecting this region of infection we can directly annotate the ‘Xenium_mol’ assay and the metadata of molecules. We use the H&E image as the background again, and generate a new molecule-level metadata column called “Infected”.
vr2_merged_acute1 <- annotateSpatialData(vr2_merged_acute1, assay = "Xenium_mol",
label = "Region", use.image = TRUE,
channel = "H&E", annotation_assay = "ROIAnnotation1")
Visualize annotations
Before visualizing the annotations and molecule localizations in the same time, lets make some changes to the metadata. We basically wanna change the annotation of all ROIs originated from the GeoJson file to have the label Hyaline Membrane.
# update molecule metadata
vrMainAssay(vr2_merged_acute1) <- "Xenium_mol"
vr2_merged_acute1$Region <- gsub("ROI[0-9]+", "Hyaline Membrane",
vr2_merged_acute1$Region)
# update ROI metadata
vrMainAssay(vr2_merged_acute1) <- "ROIAnnotation"
vr2_merged_acute1$Region <- gsub("ROI[0-9]+", "Hyaline Membrane",
vr2_merged_acute1$Region)Now we can visualize two assays together. We use the addSpatialLayer function to overlay molecule locations with both the Hyaline Membrane and the infected region annotations.
vrSpatialPlot(vr2_merged_acute1, assay = "Xenium_mol", group.by = "gene",
group.ids = c("S2_N", "S2_orf1ab"), n.tile = 500) |>
addSpatialLayer(vr2_merged_acute1, assay = "ROIAnnotation",
group.by = "Region", alpha = 0.3,
colors = list(`Hyaline Membrane` = "blue", `Infected` = "yellow"))
Comparison of annotations
By using the newly annotated infection-associated virus molecules, we can do a comparison of infected regions versus the hyaline membranes.
library(dplyr)
s2_summary_infected <-
Metadata(vr2_merged_acute1, assay = "Xenium_mol") %>%
filter(gene %in% c("S2_N", "S2_orf1ab"),
Region == "Infected") %>%
summarise(S2_N = sum(gene == "S2_N"),
S2_orf1ab = sum(gene == "S2_orf1ab"),
ratio = sum(gene == "S2_N")/sum(gene == "S2_orf1ab")) %>%
as.matrix()
s2_summary_infected S2_N S2_orf1ab ratio
[1,] 6376 2372 2.688027Comparison of both the ratio between the infected region and the hyaline membranes show a considerable difference of the population of S2_N and S2_orf1ab molecules.
S2_table <- matrix(c(s2_summary_hyaline[,1:2],
s2_summary_infected[,1:2]),
dimnames = list(Region = c("Hyaline", "Infected"),
S2 = c("N", "orf1ab")),
ncol = 2, byrow = TRUE)
S2_table S2
Region N orf1ab
Hyaline 50977 33532
Infected 6376 2372A quick test of independance on this contingency table show a significant difference of ratios across Hyaline and Infected regions.
fisher.test(S2_table, alternative = "two.sided")
Fisher's Exact Test for Count Data
data: S2_table
p-value < 2.2e-16
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.5382305 0.5941683
sample estimates:
odds ratio
0.5655696 Xenium Breast Cancer Data
VoltRon allows analyzing biological images, such as H&E stainings of tissue sections, independant to omics assays on either the same or adjacent tissue section. VoltRon defines image tiles, which are collections of image pixels of the same size, and generate latent variables that are analyzed downstream. Unsupervised clustering of these tiles partition regions over the image with respect to morphological features. Here, these regions might be associated to tumor, immune or strama niches can then be transfered to cells in the same tissue section learn collections of cells with similar features.
We will use the H&E images generated from the same sections as the Xenium In Situ platform, and cluster tiles of the H&E to later transfer groups of tiles to Xenium cells.
You can download the Xenium readout and the H&E image of the same tissue section from the 10x Genomics website (specifically, import In Situ Replicate 1 and Supplemental: Post-Xenium H&E image (TIFF)).
Import Spatial Data
We import both Xenium and post-Xenium H&E images into two separate VoltRon objects and overlay H&E images.
Xen_R1 <- importXenium("Xenium_R1/outs",
sample_name = "XeniumR1",
resolution_level = 3,
overwrite_resolution = TRUE)
Xen_R1Xen_R1_image <- importImageData("Xenium_FFPE_Human_Breast_Cancer_Rep1_he_image.tif",
sample_name = "XeniumR1",
channel_names = "H&E",
tile.size = 128)
Xen_R1_image Clustering
Using VoltRon we can also integrate tile-level spatial entities and image tiles with existing image foundation models such as DinoV2 or UNI. We can parse tiles from VoltRon, and get embedding values from a Foundation model to conduct downstream clustering using PCA and kmeans clustering. See Image tiles/pixels analysis tutorial for more information. You can also download the DinoV2 embeddings from here.
embed_data <- readRDS("img_128_processed.rds")
ind <- sapply(vrSpatialPoints(Xen_R1_image), function(x) strsplit(x, split = "_")[[1]][1])
embed_data <- embed_data[ind,]
rownames(embed_data) <- vrSpatialPoints(Xen_R1_image)
vrEmbeddings(Xen_R1_image, type = "dinov2", overwrite = TRUE) <- embed_dataXen_R1_image <- getPCA(Xen_R1_image, dims = 50, data.type = "dinov2")
Xen_R1_image <- getClusters(Xen_R1_image, method = "kmeans", label = "Cluster_kmeans", data.type = "pca", nclus = 10)vrSpatialPlot(Xen_R1_image, group.by = "Cluster_kmeans", alpha = 1)
Automated H&E Registration
In order to transfer the cluster labels of the tiles that partition regions, we need to first syncronize the spatial coordinate systems of the H&E image and the Xenium cells. Here we use the H&E image as the reference image and warp the coordinate system of the Xenium cells using the DAPI channel of the Xenium assay. See Spatial Data Alignment tutorial for more information on image registeration and alignment.
xen_reg <- registerSpatialData(object_list = list(Xen_R1_image, Xen_R1))
We can now parse the registered Xenium data from the result of the registerSpatialData function, and add this assay to the same section as the H&E image.
Xen_R1_reg <- xen_reg$registered_spat[[2]]
Xen_R1_merged <- Xen_R1_image
Xen_R1_merged <- addAssay(Xen_R1_merged,
metadata = Metadata(Xen_R1_reg),
assay = Xen_R1_reg[["Assay1"]],
assay_name = "Xenium",
sample = "XeniumR1")Label Transfer
We can view the combined dataset with the H&E tile data and Xenium assay that are in the same tissue section.
SampleMetadata(Xen_R1_merged) Assay Layer Sample
Assay1 ImageData Section1 XeniumR1
Assay2 Xenium Section1 XeniumR1We would like to transfer the labels of the tiles to cells whose centroids are within a single tile. Thus we transfer Cluster_kmeans from Assay1 (H&E tiles) to Assay2 (Xenium cells). We can also transfer the H&E image to the Xenium assay the second (registered) coordinate system of Xenium data is sync with the H&E image.
# transfer label
Xen_R1_merged <- transferData(Xen_R1_merged,
from = "Assay1",
to = "Assay2",
features = "Cluster_kmeans",
expand = FALSE)
# transfer image
vrImages(Xen_R1_merged[["Assay2"]], name = "main_reg", channel = "H&E") <-
vrImages(Xen_R1_merged, assay = "Assay1")You can visualize the transferred tile labels in two different spatial coordinate systems and with different images (DAPI, H&E etc.) using spatial and channel arguements, respectively.
vrMainAssay(Xen_R1_merged) <- "Xenium"
vrSpatialPlot(Xen_R1_merged, group.by = "Cluster_kmeans", spatial = "main")
vrSpatialPlot(Xen_R1_merged, group.by = "Cluster_kmeans", spatial = "main_reg",
channel = "H&E")
|
|
Marker Analysis
Now that we have labelled cells whose partitioning is determined by the morphological (image foundayion model) features, we can perform marker analysis to determine signatures associated with groups. We convert the cells assays of the VoltRon object to a Seurat object prior to marker analysis, find markers associated with groups and then parse significant markers to visualize later. See Conversion tutorial for more information on VoltRon’s interoperability features.
library(Seurat)
Xen_data_seu <- VoltRon::as.Seurat(Xen_R1_merged, cell.assay = "Xenium", type = "image")
Xen_data_seu <- NormalizeData(Xen_data_seu, scale.factor = 100)
Idents(Xen_data_seu) <- "Cluster_kmeans"
markers_image <- FindAllMarkers(Xen_data_seu)
markers_image_top <- markers_image %>%
group_by(cluster) %>%
dplyr::filter(avg_log2FC > 1, p_val_adj < 0.05, pct.1 > 0.5)
head(as.data.frame(markers_image_top)) p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
1 0 1.679807 0.787 0.335 0 8 MLPH
2 0 1.406676 0.869 0.430 0 8 KRT8
3 0 1.442903 0.782 0.359 0 8 CDH1
4 0 1.693801 0.808 0.385 0 8 FASN
5 0 1.526783 0.874 0.458 0 8 EPCAM
6 0 1.493334 0.670 0.259 0 8 MYO5BHeatmap visualization of these clusters reveal markers of tumor, stroma and immune niches.
Xen_R1_merged <- normalizeData(Xen_R1_merged, sizefactor = 1000)
vrHeatmapPlot(Xen_R1_merged, group.by = "Cluster_kmeans",
features = unique(markers_image_top$gene),
show_row_names = TRUE)
Inspecting the heatmap and the spatial plot simultanuously, we can see that some tiles are associated with the H&E background where others are representing tumor, immune and stroma niches.
annotation <- Xen_R1_image$Cluster_kmeans
annotation[annotation %in% c(6,8)] <- "Tumor Niche"
annotation[annotation %in% c(1)] <- "Immune Niche"
annotation[annotation %in% c(2,10,9)] <- "Stroma"
annotation[annotation %in% c(3,5,7,4)] <- "Background"
Xen_R1_image$annotation <- annotationUsers can also visualize tile and cell simultaneously and their features.
vrSpatialPlot(Xen_R1_image_subset, group.by = "Cluster_kmeans", alpha = 0.7)
vrSpatialPlot(Xen_R1_image_subset, assay = "Xenium", group.by = "annotation",
alpha = 0.7, plot.segments = TRUE, channel = "H&E")
vrSpatialFeaturePlot(Xen_R1_image_subset, assay = "Xenium",
features = c("PTPRC", "EPCAM", "MMP2"),
alpha = 0.7, plot.segments = TRUE, ncol = 2)