---

title: "ROI Analysis"

output:

  html_document:

    toc: yes

    toc_float:

      collapsed: no

      smooth_scroll: no

---

<style>

.title{

  display: none;

}

body {

  text-align: justify

}

.center {

  display: block;

  margin-left: auto;

  margin-right: auto;

}

table, th, td {

  border-collapse: collapse;

  align-self: center;

  padding-right: 10px;

  padding-left: 10px;

}

</style>

```{css, echo=FALSE}

.watch-out {

  color: black;

}

```

```{r setup, include=FALSE}

knitr::opts_chunk$set(highlight = TRUE, echo = TRUE)

```

<br>

# ROI Analysis

VoltRon is capable of analyzing readouts from distinct spatial technologies including **segmentation (ROI)-based transciptomics assays** that capture large polygonic regions on tissue sections. VoltRon recognizes such readouts including ones from commercially available tools and allows users to implement a workflow similar to ones conducted on bulk RNA-Seq datasets. In this tutorial, we will analyze morphological images and gene expression profiles provided by the readouts of the [Nanostring's GeoMx Digital Spatial Profiler](https://nanostring.com/products/geomx-digital-spatial-profiler/geomx-dsp-overview/) platform, a high-plex spatial profiling technology which produces segmentation-based protein and RNA assays. 

In this use case, **eight tissue micro array (TMA) tissue sections** were fitted into the scan area of the slide loaded on the GeoMx DSP instrument. These sections were cut from **control and COVID-19 lung tissues** of donors categorized based on disease durations (acute and prolonged). See [GSE190732](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190732) for more information.  

You can download these tutorial files from here:

<table>

  <tr>

    <th>File</th>

    <th>Link</th>

  </tr>

  <tr>

    <td>Counts</td>

    <td><a href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/ROIanalysis/GeoMx/DCC-20230427.zip">DDC files</a></td>

  </tr>

  <tr>

    <td>GeoMx Human Whole Transcriptome Atlas</td>

    <td><a href="https://nanostring.com/wp-content/uploads/Hs_R_NGS_WTA_v1.0.pkc_.zip">Human RNA WTA for NGS</a></td>

  </tr>

  <tr>

    <td>Segment Summary</td>

    <td><a href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/ROIanalysis/GeoMx/segmentSummary.csv"> ROI Metadata file </a></td>

  </tr>

  <tr>

    <td>Morphology Image</td>

    <td><a href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/ROIanalysis/GeoMx/Lu1A1B5umtrueexp.tif"> Image file </a></td>

  </tr>

  <tr>

    <td>OME.TIFF Image </td>

    <td><a href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/ROIanalysis/GeoMx/Lu1A1B5umtrueexp.ome.tiff"> OME.TIFF file </a></td>

  </tr>

    <tr>

    <td>OME.TIFF Image (XML) </td>

    <td><a href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/ROIanalysis/GeoMx/Lu1A1B5umtrueexp.ome.tiff.xml" download target="_blank"> OME.TIFF (XML) file </a></td>

  </tr>

</table>

<p></p>

We now import the GeoMx data, and start analyzing 87 user selected segments (i.e. region of interests, **ROI**) to check spatial localization of signals. The **importGeoMx** function requires:

* The path to the Digital Count Conversion file, **dcc.path**, and Probe Kit Configuration file, **pkc.file**, which are both provided as the output of the [GeoMx NGS Pipeline](https://emea.illumina.com/products/by-type/informatics-products/basespace-sequence-hub/apps/nanostring-geomxr-ngs-pipeline.html). 

* The path the to the metadata file, **summarySegment**, and the specific excel sheet that the metadata is found, **summarySegmentSheetName**, the path to the main morphology **image** and the original **ome.tiff (xml)** file, all of which are provided and imported from the DSP Control Center. Please see [GeoMx DSP Data Analysis User Manual](https://nanostring.com/support-documents/geomx-dsp-data-analysis-user-manual/) for more information.

```{r eval = FALSE, class.source="watch-out"}

library(VoltRon)

library(xlsx)

GeoMxR1 <- importGeoMx(dcc.path = "DCC-20230427/", 

                       pkc.file = "Hs_R_NGS_WTA_v1.0.pkc",

                       summarySegment = "segmentSummary.csv",

                       image = "Lu1A1B5umtrueexp.tif",

                       ome.tiff = "Lu1A1B5umtrueexp.ome.tiff.xml",

                       sample_name = "GeoMxR1")

```

We can import the GeoMx ROI segments from the **Lu1A1B5umtrueexp.ome.tiff** image file directly by replacing the .xml file with the .ome.tiff file in the **ome.tiff** argument. Note that you need to call the **RBioFormats** library. If you are getting a **java error** when running importGeoMx, try increasing the maximum heap size by supplying the **-Xmx** parameter. Run the code below before rerunning **importGeoMx** again.

```{r eval = FALSE, class.source="watch-out"}

options(java.parameters = "-Xmx4g")

library(RBioFormats)

```

<br>

## Quality Control

Once the GeoMx data is imported, we can start off with examining key quality control measures and statistics on each segment to investigate each individual ROI such as sequencing saturation and the number of cells (nuclei) within each segment. VoltRon also provides the total number of unique transcripts per ROI and stores in the metadata. 

```{r eval = FALSE, class.source="watch-out"}

library(ggplot2)

vrBarPlot(GeoMxR1, 

          features = c("Count", "Nuclei.count", "Sequencing.saturation"),

          x.label = "ROI.name", group.by = "ROI.type") + 

  theme(axis.text.x = element_text(size = 3))

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_barplot.png" class="center">

For measuring the quality of individual ROIs, we can add a new metadata column, called **CountPerNuclei**, to check the average quality of cells per ROI. It seems some number of ROIs with low counts per nuclei also have low sequencing saturation. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1$CountPerNuclei <- GeoMxR1$Count/GeoMxR1$Nuclei.count

vrBarPlot(GeoMxR1, 

          features = c("Count", "Nuclei.count", 

                       "Sequencing.saturation", "CountPerNuclei"),

          x.label = "ROI.name", group.by = "ROI.type", ncol = 2) + 

  theme(axis.text.x = element_text(size = 5))

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_barplot2.png" class="center">

<br>

## Processing

We can now filter ROIs based on our earlier observation of them having low count per nuclei where some also have low sequencing saturation. 

```{r eval = FALSE, class.source="watch-out"}

# Filter for count per nuclei

GeoMxR1 <- subset(GeoMxR1, subset = CountPerNuclei > 500)

```

We then filter genes with low counts by extracting the count matrix and putting aside all genes whose maximum count across all 87 ROIs are less than 10.

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1_data <- vrData(GeoMxR1, norm = FALSE)

feature_ind <- apply(GeoMxR1_data, 1, function(x) max(x) > 10)

selected_features <- vrFeatures(GeoMxR1)[feature_ind]

GeoMxR1_lessfeatures <- subset(GeoMxR1, features = selected_features)

```

VoltRon is capable of normalizing data provided by a diverse set of spatial technologies, including the quantile normalization method suggested by the [GeoMx DSP Data Analysis User Manual](https://nanostring.com/support-documents/geomx-dsp-data-analysis-user-manual/) which scales the ROI profiles to the third quartile followed by the geometric mean of all third quartiles multipled to the scaled profile. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1 <- normalizeData(GeoMxR1, method = "Q3Norm")

```

<br>

## Interactive Subsetting

Spatially informed genomic technologies heavily depend on image manipulation as images provide an additional set of information. Hence, VoltRon incorporates several interactive built-in utilities. One such functionality allows manipulating images of VoltRon assays where users can interactively choose subsets of images. However, we first resize the morphology image by providing the width of the new image (thus height will be reduced to preserve the aspect ratio). 

```{r eval = FALSE, class.source="watch-out"}

# resizing the image

# GeoMxR1 <- resizeImage(GeoMxR1, size = 4000)

```

VoltRon provides **a mini Shiny app** for subsetting spatial points of a VoltRon object by using the image as a reference. This app is particularly useful when multiple tissue sections were fitted to a scan area of a slide, such as the one from GeoMx DSP instrument. We use **interactive = TRUE** option in the subset function to call the mini Shiny app, and select bounding boxes of each newly created subset. **Please continue until all eight sections are selected and subsetted**. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1_subset <- subset(GeoMxR1, interactive = TRUE)

```

<img width="85%" height="85%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/interactivesubset.gif" class="center">

We can now merge the list of subsets, or samples, each associated with one of eight sections. You can provide a list of names for the newly subsetted samples. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1_subset_list <- GeoMxR1_subset$subsets

GeoMxR1 <- merge(GeoMxR1_subset_list[[1]], GeoMxR1_subset_list[-1])

GeoMxR1

```

```

VoltRon Object 

prolonged case 4: 

  Layers: Section1 

acute case 3: 

  Layers: Section1 

control case 2: 

  Layers: Section1 

acute case 1: 

  Layers: Section1 

acute case 2: 

  Layers: Section1 

... 

There are 8 samples in total 

Assays: GeoMx(Main) 

```

<br>

You may also save the selected image subsets and reproduce the interactive subsetting operation for later use. 

```{r eval = FALSE, class.source="watch-out"}

samples <- c("prolonged case 4", "acute case 3", "control case 2", 

             "acute case 1", "acute case 2", "prolonged case 5", 

             "prolonged case 3", "control case 1")

subset_info_list <- GeoMxR1_subset$subset_info_list[[1]]

GeoMxR1_subset_list <- list()

for(i in 1:length(subset_info_list)){

  GeoMxR1_subset_list[[i]]  <- subset(GeoMxR1, image = subset_info_list[i])

  GeoMxR1_subset_list[[i]] <- samples[i]

}

GeoMxR1 <- merge(GeoMxR1_subset_list[[1]], GeoMxR1_subset_list[-1])

```

<br>

## Visualization

We will now select sections of interests from the VoltRon object, and visualize ROIs and features for each of these sections.

The function **vrSpatialPlot** plots categorical attributes associated with ROIs. In this case, we will visualize types of ROIs that were labelled and annotated during ROI selection. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1_subset <- subset(GeoMxR1, sample = c("prolonged case 4","acute case 3"))

vrSpatialPlot(GeoMxR1_subset, group.by = "ROI.type", ncol = 3, alpha = 0.6)

```

<img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_spatialplot.png" class="center">

The function **vrSpatialFeaturePlot** detects the number of assays within each VoltRon object and visualizes each feature per each spatial image. A grid of images are visualized either the number of assays or the number of features are larger than 1. Here, we visualize the normalized expression of COL1A1 and C1S, two fibrotic markers, across ROIs of two prolonged covid cases. One may observe that the fibrotic regions of prolonged case 5 have considerably more COL1A1 and C1S compared to prolonged case 4. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1_subset <- subset(GeoMxR1, sample = c("prolonged case 4","prolonged case 5"))

vrSpatialFeaturePlot(GeoMxR1_subset, features = c("COL1A1", "C1S"), group.by = "ROI.name", 

                     log = TRUE, label = TRUE, keep.scale = "feature", title.size = 15)

```

<img width="85%" height="85%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_featureplot.png" class="center">

<br>

## Dimensionality Reduction

We can now process the normalized and demultiplexed samples to map ROIs across all sections onto lower dimensional spaces. The functions **getFeatures** and **getPCA** select features (i.e. genes) of interest from the data matrix across all samples and reduce it to a selected number of principal components. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1 <- getFeatures(GeoMxR1)

GeoMxR1 <- getPCA(GeoMxR1, dims = 30)

```

The function **vrEmbeddingPlot** can be used to visualize embedding spaces (pca, umap, etc.) for any spatial point supported by VoltRon, hence cells, spots and ROI are all visualized using the same set of functions. Here we generate a new metadata column that represents the **disease durations (control, acute and prolonged case)**, then map gene profiles to the first two principal components.

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1$Condition <- gsub(" [0-9]+$", "", GeoMxR1$Sample)

vrEmbeddingPlot(GeoMxR1, group.by = c("Condition"), embedding = "pca", pt.size = 3)

```

<img width="70%" height="70%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_embeddingplotsingle.png" class="center">

VoltRon provides additional dimensionality reduction techniques such as **UMAP**. 

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1 <- getUMAP(GeoMxR1)

```

Gene expression profiles of ROIs associated with prolonged case sections seem to show some heterogeneity. We now color segments by section (or replicate, **Sample**) to observe the sources of variability. Three replicates of prolonged cases exhibit three different clusters of ROIs.

```{r eval = FALSE, class.source="watch-out"}

vrEmbeddingPlot(GeoMxR1, group.by = c("Condition"), embedding = "pca", pt.size = 3)

vrEmbeddingPlot(GeoMxR1, group.by = c("ROI.type"), embedding = "pca", pt.size = 3)

vrEmbeddingPlot(GeoMxR1, group.by = c("ROI.type"), embedding = "umap", pt.size = 3)

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_embeddingplot.png" class="center">

## Differential Expression Analysis

VoltRon provides wrapping functions for calling tools and methods from popular differential expression analysis package such as [DESeq2](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8). We utilize **DESeq2** to find differentially expressed genes across each pair of ROI types conditions.

```{r eval = FALSE, class.source="watch-out"}

# get DE for all conditions

library(DESeq2)

library(dplyr)

DEresults <- getDiffExp(GeoMxR1, group.by = "ROI.type",

                        group.base = "vessel", method = "DESeq2")

DEresults_sig <- DEresults %>% filter(!is.na(padj)) %>% 

  filter(padj < 0.05, abs(log2FoldChange) > 1)

head(DEresults_sig)

```

<div><pre><code style="font-size: 11.7px;">         baseMean log2FoldChange     lfcSE     stat       pvalue         padj     gene             comparison

ACTA2    33.48395       1.508701 0.3458464 4.362343 1.286768e-05 4.902586e-03    ACTA2 ROI.type_vessel_epithelium

ADAMTS1  22.29160       1.152556 0.2272085 5.072680 3.922515e-07 4.109815e-04  ADAMTS1 ROI.type_vessel_epithelium

C11orf96 27.48924       1.142085 0.3041057 3.755554 1.729585e-04 2.588819e-02 C11orf96 ROI.type_vessel_epithelium

CNN1     16.87670       1.112662 0.2680222 4.151381 3.304757e-05 9.766004e-03     CNN1 ROI.type_vessel_epithelium

CRYAB    21.85960       1.264747 0.2173272 5.819552 5.900570e-09 2.472929e-05    CRYAB ROI.type_vessel_epithelium

FLNA     44.50957       1.270025 0.3243115 3.916066 9.000548e-05 1.985331e-02     FLNA ROI.type_vessel_epithelium</code></pre></div>

<br>

The **vrHeatmapPlot** takes a set of features for any type of spatial point (cells, spots and ROIs) and visualizes scaled data per each feature. We select **highlight.some = TRUE** to annotate features which could be large in size where you can also select the number of such highlighted genes with **n_highlight**. There seems to be two groups of fibrotic regions that most likely associated with two prolonged case samples. 

```{r eval = FALSE, class.source="watch-out"}

# get DE for all conditions

library(ComplexHeatmap)

vrHeatmapPlot(GeoMxR1, features = unique(DEresults_sig$gene), group.by = "ROI.type", 

              highlight.some = TRUE, n_highlight = 40)

```

<img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_heatmap.png" class="center">

<br>

In order to get a deeper understanding of differences between fibrotic regions across two prolonged case replicates. We first subset the GeoMx data to only sections with fibrotic ROIs. 

```{r eval = FALSE, class.source="watch-out"}

fibrotic_ROI <- vrSpatialPoints(GeoMxR1)[GeoMxR1$ROI.type == "fibrotic"]

GeoMxR1_subset <- subset(GeoMxR1, spatialpoints = fibrotic_ROI)

vrSpatialPlot(GeoMxR1_subset, group.by = "ROI.type", ncol = 3, alpha = 0.4)

```

<img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_spatialplot_fibrotic.png" class="center">

The **getDiffExp** function is capable of testing differential expression across all metadata columns, including the **Samples** column that captures labels of different tissue sections. Macrophage markers such as CD68 and CD163 are differentially expressed in fibrotic regions of case 5 compared to case 4, including FN1, a profibrotic gene whose expression can be found on macrophages.  

```{r eval = FALSE, class.source="watch-out"}

DEresults <- getDiffExp(GeoMxR1_subset, group.by = "Sample", 

                        group.base = "prolonged case 5", method = "DESeq2")

DEresults_sig <- DEresults %>% filter(!is.na(padj)) %>%

  filter(padj < 0.05, abs(log2FoldChange) > 1)

DEresults_sig <- DEresults_sig[order(DEresults_sig$log2FoldChange, decreasing = TRUE),]

head(DEresults_sig)

```

<div><pre><code style="font-size: 10.7px;">       baseMean log2FoldChange     lfcSE      stat       pvalue         padj   gene                               comparison

COL3A1 708.5599       6.635411 0.5198805 12.763338 2.626978e-37 1.596809e-33 COL3A1 Sample_prolonged case 5_prolonged case 4

COL1A2 836.0790       5.237228 0.4407380 11.882861 1.453071e-32 4.416246e-29 COL1A2 Sample_prolonged case 5_prolonged case 4

COL1A1 460.2184       5.175153 0.5237868  9.880267 5.069785e-23 3.081669e-20 COL1A1 Sample_prolonged case 5_prolonged case 4

FN1    278.6594       5.083687 0.3717299 13.675754 1.417301e-42 1.723013e-38    FN1 Sample_prolonged case 5_prolonged case 4

HBB    202.7693       4.944228 0.4884175 10.122955 4.370193e-24 3.794888e-21    HBB Sample_prolonged case 5_prolonged case 4

A2M    466.4328       4.925236 0.4542682 10.842133 2.173435e-27 3.774636e-24    A2M Sample_prolonged case 5_prolonged case 4</code></pre></div>

<br>

<!-- Markers of each individual tissue section for each disease duration is shown on the Heatmap.  -->

<!-- ```{r eval = FALSE, class.source="watch-out"} -->

<!-- # get DE for all conditions -->

<!-- vrHeatmapPlot(GeoMxR1, features = unique(DEresults_sig$gene),  -->

<!--               group.by = "Sample", highlight.some = TRUE) -->

<!-- ``` -->

<!-- <img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_heatmap2.png" class="center"> -->

<!-- <br> -->

## ROI Deconvolution

VoltRon supports multiple bulk RNA deconvolution algorithms to analyze the cellular composition of both ROIs and spots. In order to integrate the scRNA data and the GeoMx data sets within the VoltRon objects, we will use the [MuSiC](https://xuranw.github.io/MuSiC/articles/MuSiC.html) package. We will use a human lung scRNA dataset (GSE198864) as reference, which is found [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/ROIanalysis/GSE198864_lung.rds). 

```{r eval = FALSE, class.source="watch-out"}

set.seed(1)

library(Seurat)

library(SingleCellExperiment)

seu <- readRDS("GSE198864_lung.rds")

scRNAlung <- seu[,sample(1:ncol(seu), 10000, replace = FALSE)]

# Visualize clusters

DimPlot(scRNAlung, reduction = "umap", label = T, group.by = "Clusters") + NoLegend()

```

<img width="60%" height="60%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_deconvolution_singlecell.png" class="center">

We utilize the **getDeconvolution** function to call wrapper functions for deconvolution algorithms (see ). For all layers with GoeMx assays, an additional assay within the same layer with **\_decon** postfix will be created. The **sc.object** argument can either be a **Seurat** or **SingleCellExperiment** object. 

```{r eval = FALSE, class.source="watch-out"}

library(MuSiC)

GeoMxR1 <- getDeconvolution(GeoMxR1,

                            sc.object = scRNAlung, sc.assay = "RNA", 

                            sc.cluster = "Clusters", sc.samples = "orig.ident")

```

```

VoltRon Object 

prolonged case 4: 

  Layers: Section1 

acute case 3: 

  Layers: Section1 

control case 2: 

  Layers: Section1 

acute case 1: 

  Layers: Section1 

acute case 2: 

  Layers: Section1 

... 

There are 8 samples in total 

Assays: GeoMx(Main) 

Features: RNA(Main) NegProbe Decon 

```

We can now visualize cell type compositions of each ROI. Before running **vrProportionPlot** function, we need to set the main feature type as **Decon**. One may see the increased proportion of cells NK cells and T cells in immune ROIs. 

```{r eval = FALSE, class.source="watch-out"}

vrMainFeatureType(GeoMxR1) <- "Decon"

vrProportionPlot(GeoMxR1, x.label = "ROI.name")

```

<img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_deconvolution.png" class="center">

Here, we can subset the GeoMx object further to dive deep into the cellular proportions of each fibrotic region of prolonged cases. Comparing prolonged case 5 against case 4, we see an increase in the cellular population of the stromal cluster defined in [Mothes et. al 2023](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9922044/) (that are both vascular and airway smooth muscle cells), and an increase abundance of macrophages which may indicate that these macrophages are profibrotic being within and close to fibrotic regions with increased gene expression of FN1. 

```{r eval = FALSE, class.source="watch-out"}

# subsetting fibrotic regions

spatialpoints <- vrSpatialPoints(GeoMxR1)[GeoMxR1$ROI.type == "fibrotic"]

GeoMxR1_subset <- subset(GeoMxR1, spatialpoints = spatialpoints)

# Proportion plot

vrProportionPlot(GeoMxR1_subset, x.label = "ROI.name", split.method = "facet_grid", 

                 split.by = "Sample")

# barplot 

vrBarPlot(GeoMxR1_subset, features = c("stromal", "macrophages"), group.by = "Sample",

          x.label = "ROI.name", split.by = "Sample")

```

<img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_heatmap_fibrotic.png" class="center">

<img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_barplot_prop_fibrotic.png" class="center">

## H&E Image Integration

Questions may arise if in fact these ROIs are fibrotic even though they were initially annotated as fibrotic regions. VoltRon provides advanced image registration workflows which we can use to H&E images of from the same TMA blocks where GeoMx slides were established.

We first import the H&E image of the prolonged case 4 TMA section using the **importImageData** function. This will import the H&E image as a standalone VoltRon object. For more information on image-based VoltRon objects, see the [Downstream analysis of Pixels](pixelanalysis.html) section.

We will use the H&E image of TMA section taken from the same block as the Prolonged case 4. You can download the image from [here](https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/ROIanalysis/GeoMx/prolonged_case4_HE.tif).

```{r eval = FALSE, class.source="watch-out"}

vrHEimage <- importImageData("prolonged_case4_HE.tif", channel_names = "H&E")

vrImages(vrHEimage, scale.perc = 40)

```

<img width="50%" height="50%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/prolonged_case4_HE.png" class="center">

<br>

We can now register the GeoMx slide with the newly imported H&E based VoltRon object. Since two images are associated with TMA sections that are not adjacent, we have to rely on the localization of vessels visible in both images for alignment. VoltRon allows manipulating multiple channels of an image object two choose the best background image for manual landmark selection. For more information on both manual and automated registration of VoltRon objects, see [here](registration.html). 

VoltRon provides multiple registration methods to align images. Here, after running the **registerSpatialData** function, we choose the **Homography + Non-rigid (TPS)** method which utilizes a perspective transformation on the H&E image followed by a thin plate spline (TPS) alignment. The perspective transformation performs a global alignment between the H&E image and the main GeoMx image (here the scan image with combined antibody channels), and the TPS method allows correct local deformations between the perspective transformed H&E image and the GeoMx image for a more accurate

```{r eval = FALSE, class.source="watch-out"}

GeoMxR1_subset <- subset(GeoMxR1, sample = c("prolonged case 4"))

GeoMxReg <- registerSpatialData(reference_spatdata = GeoMxR1_subset, 

                                query_spatdata = vrHEimage)

```

<img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_manualregistration.png" class="center">

The result of the registration is a list of registered VoltRon objects which we can use for parsing the registered image as well. In this case, the second element of the resulting list would be the registered H&E based VoltRon object. 

```{r eval = FALSE, class.source="watch-out"}

vrHEimage_reg <- GeoMxReg$registered_spat[[2]]

```

We can check the additional coordinate system now available for the registered H&E image. One is the coordinate system with the original image, and the other with the one that is registered.

```{r eval = FALSE, class.source="watch-out"}

vrSpatialNames(vrHEimage_reg)

```

```

[1] "main"     "main_reg"

```

We can also exchange images where the H&E image now is registered to the perspective of the GeoMx channels, and can be defined as a new channel in the original GeoMx object.

```{r eval = FALSE, class.source="watch-out"}

vrImages(GeoMxR1_subset[["Assay1"]], name = "main", channel = "H&E") <- vrImages(vrHEimage_reg, name = "main_reg", channel = "H&E")

```

We can now observe the new channels (H&E) available for the GeoMx assay using **vrImageChannelNames**. H&E is saved as a separate channel along with the originally available antibody channels of the original GeoMx experiment.

```{r eval = FALSE, class.source="watch-out"}

vrImageChannelNames(GeoMxR1_subset)

```

<div><pre><code style="font-size: 12px;">       Assay    Layer           Sample Spatial                                               Channels

Assay1 GeoMx Section1 prolonged case 4    main scanimage,DNA,PanCK,CD45,Alpha Smooth Muscle Actin,H&E</code></pre></div>

We can now visualize the ROIs and their annotations where the registered H&E visible in the background. We define the spatial key **main** and the channel name **H&E**. 

```{r eval = FALSE, class.source="watch-out"}

vrSpatialPlot(GeoMxR1_subset, group.by = "ROI.type", alpha = 0.7, 

              spatial = "main", channel = "H&E")

```

<img width="70%" height="70%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_spatialplot_withHE.png" class="center">

Interactive Visualization can also be used to zoom in to ROIs and investigate the pathology associated with GeoMx ROIs labeled as fibrotic. 

```{r eval = FALSE, class.source="watch-out"}

vrSpatialPlot(GeoMxR1_subset, group.by = "ROI.type", alpha = 0.7, 

              spatial = "main", channel = "H&E", 

              interactive = TRUE)

```

<img width="60%" height="60%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/roi_spatialplot_withHE_interactive.png" class="center">