<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta name="generator" content="pandoc" />
<meta http-equiv="X-UA-Compatible" content="IE=EDGE" />
<title>Cell/Spot Analysis</title>
<script src="site_libs/header-attrs-2.29/header-attrs.js"></script>
<script src="site_libs/jquery-3.6.0/jquery-3.6.0.min.js"></script>
<meta name="viewport" content="width=device-width, initial-scale=1" />
<link href="site_libs/bootstrap-3.3.5/css/flatly.min.css" rel="stylesheet" />
<script src="site_libs/bootstrap-3.3.5/js/bootstrap.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/html5shiv.min.js"></script>
<script src="site_libs/bootstrap-3.3.5/shim/respond.min.js"></script>
<style>h1 {font-size: 34px;}
h1.title {font-size: 38px;}
h2 {font-size: 30px;}
h3 {font-size: 24px;}
h4 {font-size: 18px;}
h5 {font-size: 16px;}
h6 {font-size: 12px;}
code {color: inherit; background-color: rgba(0, 0, 0, 0.04);}
pre:not([class]) { background-color: white }</style>
<script src="site_libs/jqueryui-1.13.2/jquery-ui.min.js"></script>
<link href="site_libs/tocify-1.9.1/jquery.tocify.css" rel="stylesheet" />
<script src="site_libs/tocify-1.9.1/jquery.tocify.js"></script>
<script src="site_libs/navigation-1.1/tabsets.js"></script>
<link href="site_libs/highlightjs-9.12.0/textmate.css" rel="stylesheet" />
<script src="site_libs/highlightjs-9.12.0/highlight.js"></script>
<link href="site_libs/font-awesome-6.5.2/css/all.min.css" rel="stylesheet" />
<link href="site_libs/font-awesome-6.5.2/css/v4-shims.min.css" rel="stylesheet" />
<style type="text/css">
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
span.underline{text-decoration: underline;}
div.column{display: inline-block; vertical-align: top; width: 50%;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
</style>
<style type="text/css">code{white-space: pre;}</style>
<script type="text/javascript">
if (window.hljs) {
hljs.configure({languages: []});
hljs.initHighlightingOnLoad();
if (document.readyState && document.readyState === "complete") {
window.setTimeout(function() { hljs.initHighlighting(); }, 0);
}
}
</script>
<style type = "text/css">
.main-container {
max-width: 940px;
margin-left: auto;
margin-right: auto;
}
img {
max-width:100%;
}
.tabbed-pane {
padding-top: 12px;
}
.html-widget {
margin-bottom: 20px;
}
button.code-folding-btn:focus {
outline: none;
}
summary {
display: list-item;
}
details > summary > p:only-child {
display: inline;
}
pre code {
padding: 0;
}
</style>
<style type="text/css">
.dropdown-submenu {
position: relative;
}
.dropdown-submenu>.dropdown-menu {
top: 0;
left: 100%;
margin-top: -6px;
margin-left: -1px;
border-radius: 0 6px 6px 6px;
}
.dropdown-submenu:hover>.dropdown-menu {
display: block;
}
.dropdown-submenu>a:after {
display: block;
content: " ";
float: right;
width: 0;
height: 0;
border-color: transparent;
border-style: solid;
border-width: 5px 0 5px 5px;
border-left-color: #cccccc;
margin-top: 5px;
margin-right: -10px;
}
.dropdown-submenu:hover>a:after {
border-left-color: #adb5bd;
}
.dropdown-submenu.pull-left {
float: none;
}
.dropdown-submenu.pull-left>.dropdown-menu {
left: -100%;
margin-left: 10px;
border-radius: 6px 0 6px 6px;
}
</style>
<script type="text/javascript">
// manage active state of menu based on current page
$(document).ready(function () {
// active menu anchor
href = window.location.pathname
href = href.substr(href.lastIndexOf('/') + 1)
if (href === "")
href = "index.html";
var menuAnchor = $('a[href="' + href + '"]');
// mark the anchor link active (and if it's in a dropdown, also mark that active)
var dropdown = menuAnchor.closest('li.dropdown');
if (window.bootstrap) { // Bootstrap 4+
menuAnchor.addClass('active');
dropdown.find('> .dropdown-toggle').addClass('active');
} else { // Bootstrap 3
menuAnchor.parent().addClass('active');
dropdown.addClass('active');
}
// Navbar adjustments
var navHeight = $(".navbar").first().height() + 15;
var style = document.createElement('style');
var pt = "padding-top: " + navHeight + "px; ";
var mt = "margin-top: -" + navHeight + "px; ";
var css = "";
// offset scroll position for anchor links (for fixed navbar)
for (var i = 1; i <= 6; i++) {
css += ".section h" + i + "{ " + pt + mt + "}\n";
}
style.innerHTML = "body {" + pt + "padding-bottom: 40px; }\n" + css;
document.head.appendChild(style);
});
</script>
<!-- tabsets -->
<style type="text/css">
.tabset-dropdown > .nav-tabs {
display: inline-table;
max-height: 500px;
min-height: 44px;
overflow-y: auto;
border: 1px solid #ddd;
border-radius: 4px;
}
.tabset-dropdown > .nav-tabs > li.active:before, .tabset-dropdown > .nav-tabs.nav-tabs-open:before {
content: "\e259";
font-family: 'Glyphicons Halflings';
display: inline-block;
padding: 10px;
border-right: 1px solid #ddd;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before {
content: "\e258";
font-family: 'Glyphicons Halflings';
border: none;
}
.tabset-dropdown > .nav-tabs > li.active {
display: block;
}
.tabset-dropdown > .nav-tabs > li > a,
.tabset-dropdown > .nav-tabs > li > a:focus,
.tabset-dropdown > .nav-tabs > li > a:hover {
border: none;
display: inline-block;
border-radius: 4px;
background-color: transparent;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li {
display: block;
float: none;
}
.tabset-dropdown > .nav-tabs > li {
display: none;
}
</style>
<!-- code folding -->
<style type="text/css">
#TOC {
margin: 25px 0px 20px 0px;
}
@media (max-width: 768px) {
#TOC {
position: relative;
width: 100%;
}
}
@media print {
.toc-content {
/* see https://github.com/w3c/csswg-drafts/issues/4434 */
float: right;
}
}
.toc-content {
padding-left: 30px;
padding-right: 40px;
}
div.main-container {
max-width: 1200px;
}
div.tocify {
width: 20%;
max-width: 260px;
max-height: 85%;
}
@media (min-width: 768px) and (max-width: 991px) {
div.tocify {
width: 25%;
}
}
@media (max-width: 767px) {
div.tocify {
width: 100%;
max-width: none;
}
}
.tocify ul, .tocify li {
line-height: 20px;
}
.tocify-subheader .tocify-item {
font-size: 0.90em;
}
.tocify .list-group-item {
border-radius: 0px;
}
.tocify-subheader {
display: inline;
}
.tocify-subheader .tocify-item {
font-size: 0.95em;
}
</style>
</head>
<body>
<div class="container-fluid main-container">
<!-- setup 3col/9col grid for toc_float and main content -->
<div class="row">
<div class="col-xs-12 col-sm-4 col-md-3">
<div id="TOC" class="tocify">
</div>
</div>
<div class="toc-content col-xs-12 col-sm-8 col-md-9">
<div class="navbar navbar-default navbar-fixed-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-bs-toggle="collapse" data-target="#navbar" data-bs-target="#navbar">
<span class="icon-bar"></span>
<span class="icon-bar"></span>
<span class="icon-bar"></span>
</button>
<a class="navbar-brand" href="index.html">VoltRon</a>
</div>
<div id="navbar" class="navbar-collapse collapse">
<ul class="nav navbar-nav">
<li>
<a href="tutorials.html">Explore</a>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
Vignette
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li class="dropdown-submenu">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">Spatial Data Integration</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="registration.html">Spatial Data Alignment</a>
</li>
<li>
<a href="multiomic.html">Multi-omic Integration</a>
</li>
<li>
<a href="nicheclustering.html">Niche Clustering</a>
</li>
</ul>
</li>
<li class="dropdown-submenu">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">Downstream Analysis</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="roianalysis.html">ROI Analysis</a>
</li>
<li>
<a href="spotanalysis.html">Cell/Spot Analysis</a>
</li>
<li>
<a href="moleculeanalysis.html">Molecule Analysis</a>
</li>
<li>
<a href="pixelanalysis.html">Pixels (Image Only) Analysis</a>
</li>
</ul>
</li>
<li class="dropdown-submenu">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">Utilities</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="interactive.html">Interactive Utilities</a>
</li>
<li>
<a href="importingdata.html">Importing Spatial Data</a>
</li>
<li>
<a href="voltronobjects.html">Working with VoltRon Objects</a>
</li>
<li>
<a href="conversion.html">Converting VoltRon Objects</a>
</li>
<li>
<a href="ondisk.html">OnDisk-based Analysis Utilities</a>
</li>
</ul>
</li>
</ul>
</li>
</ul>
<ul class="nav navbar-nav navbar-right">
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
<span class="fa fa-envelope-o"></span>
Contact
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="https://bioinformatics.mdc-berlin.de">Altuna Lab/BIMSB Bioinfo</a>
</li>
<li>
<a href="https://www.mdc-berlin.de/landthaler">Landthaler Lab/BIMSB</a>
</li>
</ul>
</li>
<li class="dropdown">
<a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" data-bs-toggle="dropdown" aria-expanded="false">
<span class="fa fa-github"></span>
GitHub
<span class="caret"></span>
</a>
<ul class="dropdown-menu" role="menu">
<li>
<a href="https://github.com/BIMSBbioinfo/VoltRon">VoltRon</a>
</li>
<li>
<a href="https://github.com/BIMSBbioinfo">BIMSB Bioinfo</a>
</li>
</ul>
</li>
</ul>
</div><!--/.nav-collapse -->
</div><!--/.container -->
</div><!--/.navbar -->
<div id="header">
<h1 class="title toc-ignore">Cell/Spot Analysis</h1>
</div>
<style>
.title{
display: none;
}
body {
text-align: justify
}
.center {
display: block;
margin-left: auto;
margin-right: auto;
}
</style>
<style type="text/css">
.watch-out {
color: black;
}
</style>
<p><br></p>
<div id="xenium-data-analysis" class="section level1">
<h1>Xenium Data Analysis</h1>
<p>VoltRon is an end-to-end spatial omic analysis package which also
supports investigating spatial points in single cell resolution. VoltRon
includes essential built-in functions capable of
<strong>filtering</strong>, <strong>processing</strong> and
<strong>clustering</strong> as well as <strong>visualizing</strong>
spatial datasets with a goal of cell type discovery and annotation.</p>
<p>In this use case, we analyse readouts of the experiments conducted on
example tissue sections analysed by the <a
href="https://www.10xgenomics.com/platforms/xenium">Xenium In Situ</a>
platform. Two tissue sections of 5 <span
class="math inline">\(\mu\)</span>m tickness are derived from a single
formalin-fixed, paraffin-embedded (FFPE) breast cancer tissue block.
More information on the spatial datasets and the study can be also be
found on the <a
href="https://www.biorxiv.org/content/10.1101/2022.10.06.510405v1">BioArxiv
preprint</a>.</p>
<p>You can import these readouts from the <a
href="https://www.10xgenomics.com/products/xenium-in-situ/preview-dataset-human-breast">10x
Genomics website</a> (specifically, import <strong>In Situ Replicate
1/2</strong>). Alternatively, you can <strong>download a zipped
collection of Xenium readouts</strong> from <a
href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/SpatialDataAlignment/Xenium_vs_Visium/10X_Xenium_Visium.zip">here</a>.</p>
<p><br></p>
<div id="building-voltron-objects" class="section level2">
<h2>Building VoltRon objects</h2>
<p>VoltRon includes built-in functions for converting readouts of Xenium
experiments into VoltRon objects. The <strong>importXenium</strong>
function locates all readout documents under the output folder of the
Xenium experiment, and forms a VoltRon object. We will import both
Xenium replicates separately, and merge them after some image
manipulation.</p>
<pre class="r watch-out"><code>library(VoltRon)
Xen_R1 <- importXenium("Xenium_R1/outs", sample_name = "XeniumR1", import_molecules = TRUE)
Xen_R2 <- importXenium("Xenium_R2/outs", sample_name = "XeniumR2", import_molecules = TRUE)</code></pre>
<p>Before moving on to the downstream analysis of the imaging-based
data, we can inspect both Xenium images. We use the
<strong>vrImages</strong> function to call and visualize reference
images of all VoltRon objects. Observe that the DAPI image of the second
Xenium replicate is dim, hence we might need to increase the
brightness.</p>
<pre class="r watch-out"><code>vrImages(Xen_R1)
vrImages(Xen_R2)</code></pre>
<table>
<tbody>
<tr style="vertical-align: center">
<td style="width:50%; vertical-align: center">
<img src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/xeniumr1.png" class="center">
</td>
<td style="width:50%; vertical-align: center">
<img src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/xeniumr2.png" class="center">
</td>
</tr>
</tbody>
</table>
<p><br></p>
<p>We can adjust the brightness of the second Xenium replicate using the
<strong>modulateImage</strong> function where we can change the
brightness and saturation of the reference image of this VoltRon object.
This functionality is optional for VoltRon objects and should be used
when images require further adjustments.</p>
<pre class="r watch-out"><code>Xen_R2 <- modulateImage(Xen_R2, brightness = 800)
vrImages(Xen_R2)</code></pre>
<p><img width="40%" height="40%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/xeniumr2_new.png" class="center"></p>
<p><br></p>
<p>Once both VoltRon objects are created and images are well-tuned, we
can merge these two into a single VoltRon object.</p>
<pre class="r watch-out"><code>Xen_list <- list(Xen_R1, Xen_R2)
Xen_data <- merge(Xen_list[[1]], Xen_list[-1])</code></pre>
<pre><code>VoltRon Object
XeniumR1:
Layers: Section1
XeniumR2:
Layers: Section1
Assays: Xenium(Main) </code></pre>
<p><br></p>
</div>
<div id="spatial-visualization" class="section level2">
<h2>Spatial Visualization</h2>
<p>With <strong>vrSpatialPlot</strong>, we can visualize Xenium
experiments in both cellular and subcellular context. Since we have not
yet started analyzing raw counts of cells, we can first visualize some
transcripts of interest. We first visualize mRNAs of ACTA2, a marker for
smooth muscle cell actin, and TCF7, an early exhausted t cell marker. We
can interactively select a subset of interest within the tissue section
and visualize the localization of these transcripts. Here we subset a
ductal carcinoma niche, and visualize visualize mRNAs of
<strong>(i)</strong> ACTA2, a marker for smooth muscle cell actin, and
<strong>(ii)</strong> TCF7, an early exhausted t cell marker.</p>
<pre class="r watch-out"><code>Xen_R1_subsetinfo <- subset(Xen_R1, interactive = TRUE)
Xen_R1_subset <- Xen_R1_subsetinfo$subsets[[1]]
vrSpatialPlot(Xen_R1_subset, assay = "Xenium_mol", group.by = "gene",
group.id = c("ACTA2", "KRT15", "TACSTD2", "CEACAM6"), pt.size = 0.2, legend.pt.size = 5)</code></pre>
<p><img width="70%" height="70%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_transcripts_visualize.png" class="center"></p>
<p>We can also visualize count data of cells in the Xenium replicates.
The behaviour of <strong>vrSpatialFeaturePlot</strong> (and most
plotting functions in VoltRon) depend on the number of assays associated
with the assay type (e.g. Xenium is both cell and subcellular type).
Here, we have two assays, and we visualize two features, hence the
resulting plot would include four panels. Prior to spatial
visualization, we can normalize the counts to correct for count depth of
cells by <strong>(i)</strong> dividing counts with total counts in each
cell, <strong>(ii)</strong> multiply with some constant (default:
10000), and followed by <strong>(iii)</strong> log transformation of the
counts.</p>
<pre class="r watch-out"><code>Xen_data <- normalizeData(Xen_data, sizefactor = 1000)
vrSpatialFeaturePlot(Xen_data, features = c("ACTA2", "TCF7"), alpha = 1, pt.size = 0.7)</code></pre>
<p><img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatialfeature_xenium.png" class="center"></p>
</div>
<div id="processing-and-embedding" class="section level2">
<h2>Processing and Embedding</h2>
<p>Some number of cells in both Xenium replicates might have extremely
low counts. Although cells are detected at these locations, the low
total counts of cells would make it challenging for phenotyping and
clustering these cells. Hence, we remove such cells from the VoltRon
objects.</p>
<pre class="r watch-out"><code>Xen_data <- subset(Xen_data, Count > 5)</code></pre>
<p>VoltRon is capable of reducing dimensionality of datasets using both
PCA and UMAP which we gonna use to build profile-specific neighborhood
graphs and partition the data into cell types.</p>
<pre class="r watch-out"><code>Xen_data <- getPCA(Xen_data, dims = 20)
Xen_data <- getUMAP(Xen_data, dims = 1:20)</code></pre>
<p>We can also visualize the normalized expression of these features on
embedding spaces (e.g. UMAP) using
<strong>vrEmbeddingFeaturePlot</strong> function.</p>
<pre class="r watch-out"><code>vrEmbeddingFeaturePlot(Xen_data, features = c("LRRC15", "TCF7"), embedding = "umap",
pt.size = 0.4)</code></pre>
<p><img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_featureplot_xenium.png" class="center"></p>
<p><br></p>
</div>
<div id="clustering" class="section level2">
<h2>Clustering</h2>
<p>Next, we build neighborhood graphs with the <strong>shared nearest
neighbors (SNN)</strong> of cells which are constructed from
dimensionally reduced gene expression profiles. The function
<strong>getProfileNeighbors</strong> also has an option of building
<strong>k-nearest neighbors (kNN)</strong> graphs.</p>
<pre class="r watch-out"><code>Xen_data <- getProfileNeighbors(Xen_data, dims = 1:20, method = "SNN")
vrGraphNames(Xen_data)</code></pre>
<pre><code>[1] "SNN"</code></pre>
<p>We can later conduct a clustering of cells using the <strong>leiden’s
method</strong> from the igraph package, which is utilized with the
<strong>getClusters</strong> function.</p>
<pre class="r watch-out"><code>Xen_data <- getClusters(Xen_data, resolution = 1.0, label = "Clusters", graph = "SNN")</code></pre>
<p>Now we can label each cell with the associated clustering index and
take a look at the clustering accuracy on the embedding space, and we
can also visualize these clusters on a spatial context.</p>
<pre class="r watch-out"><code>vrEmbeddingPlot(Xen_data, group.by = "Clusters", embedding = "umap",
pt.size = 0.4, label = TRUE)</code></pre>
<p><img width="60%" height="60%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_embedplot_xenium.png" class="center"></p>
<pre class="r watch-out"><code>vrSpatialPlot(Xen_data, group.by = "Clusters", pt.size = 0.18, background.color = "black")</code></pre>
<p><img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatial_xenium.png" class="center"></p>
<p><br></p>
</div>
<div id="annotation" class="section level2">
<h2>Annotation</h2>
<p>We can annotate each of these clusters according to their positive
markers across 313 features. One can use the
<strong>FindAllMarkers</strong> from the <a
href="https://satijalab.org/seurat/">Seurat</a> package to pinpoint
these markers by first utilizing the <strong>as.Seurat</strong> function
first on the Xenium assays of the VoltRon object.</p>
<p>For more information on conversion to other packages, please visit
the <a href="conversion.html">Converting VoltRon Objects</a>.</p>
<p>Let us create a new metadata feature from the
<strong>Clusters</strong> column, called <strong>CellType</strong>, we
can insert this new metadata column directly to the object.</p>
<pre class="r watch-out"><code>clusters <- factor(Xen_data$Clusters, levels = sort(unique(Xen_data$Clusters)))
levels(clusters) <- c("DCIS_1",
"DCIS_2",
"CD4_TCells",
"Adipocytes",
"PLD4+_LILRA4+_CD4+_Cells",
"ACTA2_myoepithelial",
"IT_2",
"Macrophages",
"MastCells",
"Bcells",
"StromalCells",
"CD8_TCells",
"CD8_TCells",
"EndothelialCells",
"StromalCells",
"MyelomaCells",
"IT_1",
"IT_2",
"ACTA2_myoepithelial",
"DCIS_2",
"IT_3",
"KRT15_myoepithelial")
Xen_data$CellType <- as.character(clusters)</code></pre>
<p><strong>vrSpatialPlot</strong> function can visualize multiple types
of metadata columns, and users can change the location of the legends as
well.</p>
<pre class="r watch-out"><code>vrSpatialPlot(Xen_data, group.by = "CellType", pt.size = 0.13, background.color = "black",
legend.loc = "top", n.tile = 500)</code></pre>
<p><img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatial_xenium_annotated.png" class="center"></p>
<p><br></p>
</div>
</div>
<div id="visium-data-analysis" class="section level1">
<h1>Visium Data Analysis</h1>
<p>Spot-based spatial transcriptomic assays capture spatially-resolved
gene expression profiles that are somewhat closer to single cell
resolution. However, each spot still include a few number of cells that
are likely from a combination of cell types within the tissue of origin.
VoltRon analyzes spot level spatial data sets and even allows selecting
a highly variable subset of features to cluster spots into meaningful
groups of in situ spots for detecting niches of interests</p>
<div id="import-st-data" class="section level2">
<h2>Import ST Data</h2>
<p>For this tutorial we will analyze spot-based transcriptomic assays
from Mouse Brain generated by the <a
href="https://www.10xgenomics.com/products/spatial-gene-expression">Visium</a>
instrument.</p>
<p>You can find and download readouts of all four Visium sections <a
href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/Visium/MouseBrainSerialSections.zip">here</a>.
The <strong>Mouse Brain Serial Section 1</strong> datasets can be
downloaded from <a
href="https://www.10xgenomics.com/resources/datasets?menu%5Bproducts.name%5D=Spatial%20Gene%20Expression&query=&page=1&configure%5BhitsPerPage%5D=50&configure%5BmaxValuesPerFacet%5D=1000">here</a>
(specifically, please filter for <strong>Species=Mouse</strong>,
<strong>AnatomicalEntity=brain</strong>, <strong>Chemistry=v1</strong>
and <strong>PipelineVersion=v1.1.0</strong>).</p>
<p>We will now import each of four samples separately and merge them
into one VoltRon object. There are four sections in total given two
serial anterior and serial posterior sections, hence we have <strong>two
tissue blocks each having two layers</strong>.</p>
<pre class="r watch-out"><code>library(VoltRon)
Ant_Sec1 <- importVisium("Sagittal_Anterior/Section1/", sample_name = "Anterior1")
Pos_Sec1 <- importVisium("Sagittal_Posterior/Section1/", sample_name = "Posterior1")
# merge datasets
MBrain_Sec <- merge(Ant_Sec1, Pos_Sec1, samples = c("Anterior", "Posterior"))
MBrain_Sec</code></pre>
<pre><code>VoltRon Object
Anterior:
Layers: Section1
Posterior:
Layers: Section1
Assays: Visium(Main) </code></pre>
<p>VoltRon maps metadata features on the spatial images, multiple
features can be provided for all assays/layers associated with the main
assay (Visium).</p>
<pre class="r watch-out"><code>vrSpatialFeaturePlot(MBrain_Sec, features = "Count", crop = TRUE, alpha = 1, ncol = 2)</code></pre>
<p><img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_firstplot.png" class="center"></p>
<p><br></p>
</div>
<div id="feature-selection" class="section level2">
<h2>Feature Selection</h2>
<p>VoltRon captures the nearly full transcriptome of the Visium data
which then can be filtered from a list of features ranked by their
variance and importance. We use the <strong>variance stabilization
transformation (vst)</strong> on each individual assay using the
<strong>getFeatures</strong> function and combine these ranked list to
capture features important for all assay of the Visium data later with
<strong>getVariableFeatures</strong> function.</p>
<pre class="r watch-out"><code>head(vrFeatures(MBrain_Sec))</code></pre>
<pre><code>[1] "Xkr4" "Gm1992" "Gm19938" "Gm37381" "Rp1" "Sox17" </code></pre>
<pre class="r watch-out"><code>length(vrFeatures(MBrain_Sec))</code></pre>
<pre><code>[1] 33502</code></pre>
<pre class="r watch-out"><code>MBrain_Sec <- normalizeData(MBrain_Sec)
MBrain_Sec <- getFeatures(MBrain_Sec, n = 3000)
head(vrFeatureData(MBrain_Sec))</code></pre>
<pre><code> mean var adj_var rank
Xkr4 0.0248608534 0.0249941807 0.02800216 14114
Gm1992 0.0000000000 0.0000000000 0.00000000 0
Gm19938 0.0285714286 0.0322197476 0.03224908 13889
Gm37381 0.0000000000 0.0000000000 0.00000000 0
Rp1 0.0003710575 0.0003710575 0.00000000 0
Sox17 0.1907235622 0.2219629135 0.23715920 10304</code></pre>
<pre class="r watch-out"><code>selected_features <- getVariableFeatures(MBrain_Sec)
head(selected_features, 20)</code></pre>
<pre><code>[1] "Bc1" "mt-Co1" "mt-Co3" "mt-Atp6" "mt-Co2" "mt-Cytb" "mt-Nd4" "mt-Nd1" "mt-Nd2"
[2] "Fth1" "Hbb-bs" "Cst3" "Gapdh" "Tmsb4x" "Mbp" "Rplp1" "Ttr" "Ppia"
[3] "Ckb" "mt-Nd3" </code></pre>
</div>
<div id="embedding" class="section level2">
<h2>Embedding</h2>
<p>Now we can learn and visualize PCA and UMAP embeddings on this
smaller number of selected features</p>
<pre class="r watch-out"><code>MBrain_Sec <- getPCA(MBrain_Sec, features = selected_features, dims = 30)
MBrain_Sec <- getUMAP(MBrain_Sec, dims = 1:30)
vrEmbeddingPlot(MBrain_Sec, embedding = "umap")</code></pre>
<p><img width="65%" height="65%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_umap.png" class="center"></p>
<p><br></p>
</div>
<div id="clustering-1" class="section level2">
<h2>Clustering</h2>
<pre class="r watch-out"><code>MBrain_Sec <- getProfileNeighbors(MBrain_Sec, dims = 1:30, k = 10, method = "SNN")
vrGraphNames(MBrain_Sec)</code></pre>
<pre><code>[1] "SNN"</code></pre>
<pre class="r watch-out"><code>MBrain_Sec <- getClusters(MBrain_Sec, resolution = 0.5, label = "Clusters", graph = "SNN")
vrEmbeddingPlot(MBrain_Sec, embedding = "umap", group.by = "Clusters")</code></pre>
<p><img width="65%" height="65%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_umap_clusters.png" class="center"></p>
<p><br></p>
<pre class="r watch-out"><code>vrSpatialPlot(MBrain_Sec, group.by = "Clusters")</code></pre>
<p><img width="65%" height="65%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_visium_spatial_clusters1.png" class="center"></p>
<p><br></p>
</div>
</div>
<div id="melc-data-analysis" class="section level1">
<h1>MELC Data Analysis</h1>
<p>VoltRon also provides support for imaging based proteomics assays. In
this next use case, we analyze cells characterized by
<strong>multi-epitope ligand cartography (MELC)</strong> with a panel of
44 parameters. We use the already segmented cells on which expression of
<strong>43 protein features</strong> (excluding DAPI) were mapped to
these cells.</p>
<p>We use the segmented cells over microscopy images collected from
<strong>control</strong> and <strong>COVID-19</strong> lung tissues of
donors categorized based on disease durations (<strong>control</strong>,
<strong>acute</strong>, <strong>chronic</strong> and
<strong>prolonged</strong>). Each image is associated with one of few
field of views (FOVs) from a single tissue section of a donor. See <a
href="https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190732">GSE190732</a>
for more information. You can download the <strong>IFdata.csv</strong>
file and the folder with the <strong>DAPI</strong> images <a
href="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Cellanalysis/MELC/GSE190732.zip">here</a>.</p>
<p>We import the <strong>protein intensities</strong>,
<strong>metadata</strong> and <strong>coordinates</strong> associated
with segmented cells across FOVs of samples.</p>
<pre class="r watch-out"><code>library(VoltRon)
IFdata <- read.csv("IFdata.csv")
data <- IFdata[,c(2:43)]
metadata <- IFdata[,c("disease_state", "object_id", "cluster", "Clusters",
"SourceID", "Sample", "FOV", "Section")]
coordinates <- as.matrix(IFdata[,c("posX","posY")], rownames.force = TRUE)</code></pre>
<p><br></p>
<div id="importing-melc-data" class="section level2">
<h2>Importing MELC data</h2>
<p>Before analyzing MELC assays across FOVs, we should <strong>build a
VoltRon object</strong> for each individual FOV/Section by using the
<strong>formVoltron</strong> function. We then merge these sections to
respective tissue blocks by defining their samples of origins. We can
also define <strong>assay names</strong>, <strong>assay types</strong>
and <strong>sample (i.e. block) names</strong> of these objects.</p>
<pre class="r watch-out"><code>library(dplyr)
library(magick)
vr_list <- list()
sample_metadata <- metadata %>% select(Sample, FOV, Section) %>% distinct()
for(i in 1:nrow(sample_metadata)){
vrassay <- sample_metadata[i,]
cells <- rownames(metadata)[metadata$Section == vrassay$Section]
image <- image_read(paste0("DAPI/", vrassay$Sample, "/DAPI_", vrassay$FOV, ".tif"))
vr_list[[vrassay$Section]] <- formVoltRon(data = t(data[cells,]),
metadata = metadata[cells,],
image = image,
coords = coordinates[cells,],
main.assay = "MELC",
assay.type = "cell",
sample_name = vrassay$Section)
}</code></pre>
<p>Before moving forward with merging FOVs, we should <strong>flip
coordinates</strong> of cells and perhaps also then
<strong>resize</strong> these images. The main reason for this
coordinate flipping is that the y-axis of most digital images are of the
opposite direction to the commonly used coordinate spaces.</p>
<pre class="r watch-out"><code>for(i in 1:nrow(sample_metadata)){
vrassay <- sample_metadata[i,]
vr_list[[vrassay$Section]] <- flipCoordinates(vr_list[[vrassay$Section]])
vr_list[[vrassay$Section]] <- resizeImage(vr_list[[vrassay$Section]], size = 600)
}</code></pre>
<p>Finally, we merge these assays into one VoltRon object. The
<strong>samples</strong> arguement in the merge function determines
which assays are layers of a single tissue sample/block.</p>
<pre class="r watch-out"><code>vr_merged <- merge(vr_list[[1]], vr_list[-1], samples = sample_metadata$Sample)
vr_merged </code></pre>
<pre><code>VoltRon Object
control_case_3:
Layers: Section1 Section2
control_case_2:
Layers: Section1 Section2
control_case_1:
Layers: Section1 Section2 Section3
acute_case_3:
Layers: Section1 Section2
acute_case_1:
Layers: Section1 Section2
...
There are 13 samples in total
Assays: MELC(Main) </code></pre>
<p><br></p>
<p>The prolonged case 4 has two fields of views (FOVs). By subsetting on
the sample of a prolonged case, we can visualize only these two
sections, and visualize the protein expression of CD31 and
Pancytokeratin which are markers of endothelial and epithelial
cells.</p>
<pre class="r watch-out"><code>vr_subset <- subset(vr_merged, samples = "prolonged_case_4")
g1 <- vrSpatialFeaturePlot(vr_subset, features = c("CD31", "Pancytokeratin"), alpha = 1,
pt.size = 0.7, background.color = "black")</code></pre>
<p><img width="90%" height="90%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_spatialfeature.png" class="center"></p>
<p><br></p>
</div>
<div id="dimensionality-reduction" class="section level2">
<h2>Dimensionality Reduction</h2>
<p>We can utilize dimensional reduction of the available protein markers
using the getPCA and getUMAP functions, but now with relatively lower
numbers of principal components which are enough to capture the
information across 44 features.</p>
<pre class="r watch-out"><code>vr_merged <- getPCA(vr_merged, dims = 10)
vr_merged <- getUMAP(vr_merged, dims = 1:10)
vrEmbeddingFeaturePlot(vr_merged, features = c("CD31", "Pancytokeratin"), embedding = "umap")</code></pre>
<p><img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_embedding.png" class="center"></p>
<p><br></p>
</div>
<div id="clustering-2" class="section level2">
<h2>Clustering</h2>
<p>Now we can visualize the clusters across these sections and perhaps
also check for clusters that may reside in only specific disease
conditions.</p>
<pre class="r watch-out"><code># SNN graph and clusters
vr_merged <- getProfileNeighbors(vr_merged, dims = 1:10, k = 10, method = "SNN")
vrGraphNames(vr_merged)</code></pre>
<pre><code>[1] "SNN"</code></pre>
<pre class="r watch-out"><code>vr_merged <- getClusters(vr_merged, resolution = 0.8, label = "MELC_Clusters", graph = "SNN")
# install patchwork package
if (!requireNamespace("patchwork", quietly = TRUE))
install.packages("patchwork")
library(patchwork)
# visualize conditions and clusters
vr_merged$Condition <- gsub("_[0-9]$", "", vr_merged$Sample)
g1 <- vrEmbeddingPlot(vr_merged, group.by = c("Condition"), embedding = "umap")
g2 <- vrEmbeddingPlot(vr_merged, group.by = c("MELC_Clusters"), embedding = "umap",
label = TRUE)
g1 | g2</code></pre>
<p><img width="100%" height="100%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_embeddingclusters.png" class="center"></p>
<p><br></p>
</div>
<div id="visualization-of-markers" class="section level2">
<h2>Visualization of Markers</h2>
<p>VoltRon provides both violin plots (<strong>vrViolinPlot</strong>)
and heatmaps (<strong>vrHeatmapPlot</strong>) to further investigate the
enrichment of markers across newly clustered datasets.
<strong>Note:</strong> the vrHeatmapPlot function would require you to
have the <strong>ComplexHeatmap</strong> package in your namespace.</p>
<pre class="r watch-out"><code># install patchwork package
if (!requireNamespace("ComplexHeatmap", quietly = TRUE))
BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)
# Visualize Markers
vrHeatmapPlot(vr_merged, features = vrFeatures(vr_merged),
group.by = "MELC_Clusters", show_row_names = TRUE)</code></pre>
<p><img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_heatmapclusters.png" class="center"></p>
<p><br></p>
<pre class="r watch-out"><code>vrViolinPlot(vr_merged, features = c("CD3", "SMA", "Pancytokeratin", "CCR2"),
group.by = "MELC_Clusters", ncol = 2)</code></pre>
<p><img width="80%" height="80%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_violinclusters.png" class="center"></p>
<p><br></p>
</div>
<div id="neighborhood-analysis" class="section level2">
<h2>Neighborhood Analysis</h2>
<p>We use the <strong>vrNeighbourhoodEnrichment</strong> function to
detect cell type pairs that co occur within each others’ neighborhoods.
First, we establish <strong>spatial neighborhood graphs</strong> that
determine the neighbors of each cell on tissue sections.</p>
<p><a
href="https://en.wikipedia.org/wiki/Delaunay_triangulation">Delaunay
tesselations</a> or graphs are commonly used to determine neighbors of
spatial entities. The function <strong>getSpatialNeighbors</strong>
builds a delaunay graph of all assays of a certain type and detects
neighbors of cells in a VoltRon object.</p>
<pre class="r watch-out"><code>vr_merged <- getSpatialNeighbors(vr_merged, method = "delaunay")</code></pre>
<p>The graph <strong>delaunay</strong>, which we will use for
spatially-aware neightborhood analysis, is now the second graph
available in the VoltRon object along with <strong>SNN</strong>.</p>
<pre class="r watch-out"><code>vrGraphNames(vr_merged)</code></pre>
<pre><code>[1] "SNN" "delaunay"</code></pre>
<p>Once neighbors are founds, we can apply a <strong>permutation
test</strong> that compares the number of cell type occurances with an
expected number of these occurances under multiple permutations of
labels in the tissue (fixed coordinates but cells are randomly
labelled). A similar approach is used to by several spatial analysis
frameworks and packages (<a
href="https://www.nature.com/articles/nmeth.4391">Schapiro et. al
2017</a>, <a
href="https://www.nature.com/articles/s41592-021-01358-2">Palla et. al
2022</a>).</p>
<p>Here, we will use the original cell type labels annotated by <a
href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9922044/">Mothes et.
al 2023</a>.</p>
<pre class="r watch-out"><code>neighborhood_results <- vrNeighbourhoodEnrichment(vr_merged, group.by = "Clusters")</code></pre>
<p>The neighborhood analysis provides the results of:</p>
<ul>
<li>the <strong>association</strong> tests (whether cell types are
within each other’s neighborhood)</li>
<li>the <strong>segregation</strong> tests (whether cell types are
clustered separately)</li>
</ul>
<p>between all cell type pairs across each layers and assay.</p>
<p>The number of each cell in a pair in each section is reported to
assess the impact of the results of the test (i.e. low number of
abundance in one cell type may indicate low impact).</p>
<pre class="r watch-out"><code>head(neighborhood_results)</code></pre>
<div>
<pre><code style="font-size: 10px;"> from_value to_value p_assoc p_segreg p_assoc_adj p_segreg_adj n_from n_to AssayID Assay Layer Sample
Assay1.1 CD163+ macs CD163+ macs 0.0000000 1.00000000 0.0000 1.00000000 41 41 Assay1 MELC Section1 control_case_3
Assay1.2 CD163+ macs CD4+ T cells 0.9380000 0.03300000 0.9980 0.09762866 41 48 Assay1 MELC Section1 control_case_3
Assay1.3 CD163+ macs CD8+ Tcells 0.8779011 0.04339051 0.9980 0.09762866 41 11 Assay1 MELC Section1 control_case_3
Assay1.4 CD163+ macs NK cells 0.8190000 0.08700000 0.9980 0.15660000 41 15 Assay1 MELC Section1 control_case_3
Assay1.5 CD163+ macs endothelia 0.1230000 0.85100000 0.5535 0.95737500 41 139 Assay1 MELC Section1 control_case_3
Assay1.6 CD163+ macs epithelia 0.9320000 0.03600000 0.9980 0.09762866 41 39 Assay1 MELC Section1 control_case_3</code></pre>
</div>
<p><br></p>
<pre class="r watch-out"><code>vrNeighbourhoodEnrichmentPlot(neighborhood_results, assay = "Assay1", type = "assoc")</code></pre>
<p><img width="70%" height="70%" src="https://bimsbstatic.mdc-berlin.de/landthaler/VoltRon/Package/images/cellspot_neighenrichment.png" class="center"></p>
</div>
</div>
</div>
</div>
</div>
<script>
// add bootstrap table styles to pandoc tables
function bootstrapStylePandocTables() {
$('tr.odd').parent('tbody').parent('table').addClass('table table-condensed');
}
$(document).ready(function () {
bootstrapStylePandocTables();
});
</script>
<!-- tabsets -->
<script>
$(document).ready(function () {
window.buildTabsets("TOC");
});
$(document).ready(function () {
$('.tabset-dropdown > .nav-tabs > li').click(function () {
$(this).parent().toggleClass('nav-tabs-open');
});
});
</script>
<!-- code folding -->
<script>
$(document).ready(function () {
// temporarily add toc-ignore selector to headers for the consistency with Pandoc
$('.unlisted.unnumbered').addClass('toc-ignore')
// move toc-ignore selectors from section div to header
$('div.section.toc-ignore')
.removeClass('toc-ignore')
.children('h1,h2,h3,h4,h5').addClass('toc-ignore');
// establish options
var options = {
selectors: "h1,h2,h3",
theme: "bootstrap3",
context: '.toc-content',
hashGenerator: function (text) {
return text.replace(/[.\\/?&!#<>]/g, '').replace(/\s/g, '_');
},
ignoreSelector: ".toc-ignore",
scrollTo: 0
};
options.showAndHide = false;
options.smoothScroll = false;
// tocify
var toc = $("#TOC").tocify(options).data("toc-tocify");
});
</script>
<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
(function () {
var script = document.createElement("script");
script.type = "text/javascript";
script.src = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
document.getElementsByTagName("head")[0].appendChild(script);
})();
</script>
</body>
</html>
Datasets
Models
