#!/usr/bin/env python
"""
Analyze Slide-seqV2 data
========================
This tutorial shows how to apply Squidpy for the analysis of Slide-seqV2 data.
The data used here was obtained from :cite:`Stickels2020-rf`.
We provide a pre-processed subset of the data, in :class:`anndata.AnnData` format.
We would like to thank @tudaga for providing cell-type level annotation.
For details on how it was pre-processed, please refer to the original paper.
.. seealso::
See :ref:`sphx_glr_auto_tutorials_tutorial_visium_hne.py` and
:ref:`sphx_glr_auto_tutorials_tutorial_seqfish.py` for additional analysis examples.
Import packages & data
----------------------
To run the notebook locally, create a conda environment as *conda env create -f environment.yml* using this
`environment.yml `_.
"""
import squidpy as sq
print(f"squidpy=={sq.__version__}")
# load the pre-processed dataset
adata = sq.datasets.slideseqv2()
adata
###############################################################################
# First, let's visualize cluster annotation in spatial context
# with :func:`squidpy.pl.spatial_scatter`.
sq.pl.spatial_scatter(adata, color="cluster", size=1, shape=None)
###############################################################################
# Neighborhood enrichment analysis
# --------------------------------
# Similar to other spatial data, we can investigate spatial organization of clusters
# in a quantitative way, by computing a neighborhood enrichment score.
# You can compute such score with the following function: :func:`squidpy.gr.nhood_enrichment`.
# In short, it's an enrichment score on spatial proximity of clusters:
# if spots belonging to two different clusters are often close to each other,
# then they will have a high score and can be defined as being *enriched*.
# On the other hand, if they are far apart, the score will be low
# and they can be defined as *depleted*.
# This score is based on a permutation-based test, and you can set
# the number of permutations with the `n_perms` argument (default is 1000).
#
# Since the function works on a connectivity matrix, we need to compute that as well.
# This can be done with :func:`squidpy.gr.spatial_neighbors`.
# Please see :ref:`sphx_glr_auto_examples_graph_compute_spatial_neighbors.py` and
# :ref:`sphx_glr_auto_examples_graph_compute_nhood_enrichment.py` for more details
# of how these functions works.
#
# Finally, we'll directly visualize the results with :func:`squidpy.pl.nhood_enrichment`.
# We'll add a dendrogram to the heatmap computed with linkage method *ward*.
sq.gr.spatial_neighbors(adata, coord_type="generic")
sq.gr.nhood_enrichment(adata, cluster_key="cluster")
sq.pl.nhood_enrichment(adata, cluster_key="cluster", method="single", cmap="inferno", vmin=-50, vmax=100)
###############################################################################
# Interestingly, there seems to be an enrichment between the *Endothelial_Tip*,
# the *Ependymal* cells. Another putative enrichment is between the *Oligodendrocytes*
# and *Polydendrocytes* cells. We can visualize the spatial organization of such clusters.
# For this, we'll use :func:`squidpy.pl.spatial_scatter` again.
sq.pl.spatial_scatter(
adata,
shape=None,
color="cluster",
groups=["Endothelial_Tip", "Ependymal", "Oligodendrocytes", "Polydendrocytes"],
size=3,
)
###############################################################################
# Ripley's statistics
# -------------------
# In addition to the neighbor enrichment score, we can further investigate spatial
# organization of cell types in tissue by means of the Ripley's statistics.
# Ripley's statistics allow analyst to evaluate whether a discrete annotation (e.g. cell-type)
# appears to be clustered, dispersed or randomly distributed on the area of interest.
# In Squidpy, we implement three closely related Ripley's statistics, that can be
# easily computed with :func:`squidpy.gr.ripley`. Here, we'll showcase the Ripley's L statistic,
# which is a variance-stabilized version of the Ripley's K statistics.
# We'll visualize the results with :func:`squidpy.pl.ripley`.
# Check :ref:`sphx_glr_auto_examples_graph_compute_ripley.py` for more details.
mode = "L"
sq.gr.ripley(adata, cluster_key="cluster", mode=mode, max_dist=500)
sq.pl.ripley(adata, cluster_key="cluster", mode=mode)
###############################################################################
# The plot highlight how some cell-types have a more clustered pattern,
# like *Astrocytes* and *CA11_CA2_CA3_Subiculum* cells, whereas other have a more
# dispersed pattern, like *Mural* cells. To confirm such interpretation, we can
# selectively visualize again their spatial organization.
sq.pl.spatial_scatter(adata, color="cluster", groups=["Mural", "CA1_CA2_CA3_Subiculum"], size=3, shape=None)
###############################################################################
# Ligand-receptor interaction analysis
# ------------------------------------
# The analysis showed above has provided us with quantitative information on
# cellular organization and communication at the tissue level.
# We might be interested in getting a list of potential candidates that might be driving
# such cellular communication.
# This naturally translates in doing a ligand-receptor interaction analysis.
# In Squidpy, we provide a fast re-implementation the popular method CellPhoneDB :cite:`cellphonedb`
# (`code `_ )
# and extended its database of annotated ligand-receptor interaction pairs with
# the popular database *Omnipath* :cite:`omnipath`.
# You can run the analysis for all clusters pairs, and all genes (in seconds,
# without leaving this notebook), with :func:`squidpy.gr.ligrec`.
#
# Let's perform the analysis and visualize the result for three clusters of
# interest: *Polydendrocytes* and *Oligodendrocytes*.
# For the visualization, we will filter out annotations
# with low-expressed genes (with the ``means_range`` argument)
# and decreasing the threshold
# for the adjusted p-value (with the ``alpha`` argument)
# Check :ref:`sphx_glr_auto_examples_graph_compute_ligrec.py` for more details.
sq.gr.ligrec(
adata,
n_perms=100,
cluster_key="cluster",
clusters=["Polydendrocytes", "Oligodendrocytes"],
)
sq.pl.ligrec(
adata,
cluster_key="cluster",
source_groups="Oligodendrocytes",
target_groups=["Polydendrocytes"],
pvalue_threshold=0.05,
swap_axes=True,
)
###############################################################################
# The dotplot visualization provides an interesting set of candidate interactions
# that could be involved in the tissue organization of the cell types of interest.
# It should be noted that this method is a pure re-implementation of the original
# permutation-based test, and therefore retains all its caveats
# and should be interpreted accordingly.
###############################################################################
# Spatially variable genes with spatial autocorrelation statistics
# ----------------------------------------------------------------
# Lastly, with Squidpy we can investigate spatial variability of gene expression.
# :func:`squidpy.gr.spatial_autocorr` conveniently wraps two
# spatial autocorrelation statistics: *Moran's I* and *Geary's C**.
# They provide a score on the degree of spatial variability of gene expression.
# The statistic as well as the p-value are computed for each gene, and FDR correction
# is performed. For the purpose of this tutorial, let's compute the *Moran's I* score.
# See :ref:`sphx_glr_auto_examples_graph_compute_moran.py` for more details.
sq.gr.spatial_autocorr(adata, mode="moran")
adata.uns["moranI"].head(10)
###############################################################################
# The results are stored in `adata.uns["moranI"]` and we can visualize selected genes
# with :func:`squidpy.pl.spatial_scatter`.
sq.pl.spatial_scatter(
adata,
shape=None,
color=["Ttr", "Plp1", "Mbp", "Hpca", "Enpp2"],
size=0.1,
)