---
title: "Infer spatially resolved temporal dynamics of cell-cell communication at cellular resolution"
author: "Jialin Liu, Yichen Wang"
date: "2024-01-03"
output: 
  html_document:
    toc: true
    toc_float:
      toc_collapsed: false
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Infer temporal cell-cell communication change in spatial at single-cell level with RNA velocity information}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

### Please run the [CytoSignal workflow](https://htmlpreview.github.io/?https://github.com/welch-lab/cytosignal/blob/master/doc/cytosignal.html) before this vignette.

## Introduction

CytoSignal detects cell-cell signaling from spatial transcriptomic data at single-cell resolution. CytoSignal performs a nonparametric statistical test to identify which cells within a tissue have significant activity for a particular signaling interaction. CytoSignal considers multi-component interactions and separately models interactions mediated by diffusion vs. contact-dependent molecules. 

As an extension to CytoSignal, we have also developed VeloCytoSignal to infer temporal dynamic changes in signaling strength using spatial data from only a single time point. VeloCytoSignal predicts the instantaneous rate of change for a signaling interaction at a specific position within a tissue by combining RNA velocity from the ligand and receptor genes to calculate the time derivative (dS/dt) of the CytoSignal LRscore (S). 

## Example data

In this tutorial, we demonstrate how to preform VeloCytoSignal on a mouse embryonic E15 brain data captured by Slide-seq V2, a spatial transcriptomic protocol at single-cell resolution. This data was originally produced in the work of [Stickels R.R. et al., 2020](https://doi.org/10.1038/s41587-020-0739-1), and is now publicly available on Single-Cell Portal via [SCP815](https://singlecell.broadinstitute.org/single_cell/study/SCP815). 

To calculate RNA velocity estimates, we utilized our previously published tool, [VeloVAE](https://github.com/welch-lab/VeloVAE). VeloVAE uses gene-shared latent times to calculate much more accurate velocity estimates for individual genes than approaches that infers by single-gene. However, VeloCytoSignal is widely applicable on most kinds of RNA velocity approaches, as long as it provides the spliced ($s$) and unspliced ($u$) RNA velocity estimates in the gene-by-cell matrix format.

Here we provide the processed CytoSignal object and pre-calculated RNA velocity estimates.

- CytoSignal object with regular CytoSignal analysis performed: [SCP815_190921_19_CytoSignal.rds](https://figshare.com/ndownloader/files/44865757)
- $s$ matrix stands for the spliced RNA velocity estimates: [SCP815_190921_19_velo_spliced_matrix.rds](https://figshare.com/ndownloader/files/44865748)
- $u$ matrix stands for the unspliced RNA velocity estimates: [SCP815_190921_19_velo_unspliced_matrix.rds](https://figshare.com/ndownloader/files/44865751)

Download the files listed above from the provided link into a desired local directory. Next, load the downloaded data into the R environment using the following codes. We assume that users have the files placed at the current working directory (as can be shown by `getwd()`). Alternatively, users can also specify the path to the files.

```{r download, include=FALSE}
if (!file.exists("SCP815_190921_19_CytoSignal.rds")) {
    download.file("https://figshare.com/ndownloader/files/44865757",
                  "SCP815_190921_19_CytoSignal.rds")
}
if (!file.exists("SCP815_190921_19_velo_spliced_matrix.rds")) {
    download.file("https://figshare.com/ndownloader/files/44865748",
                  "SCP815_190921_19_velo_spliced_matrix.rds")
}
if (!file.exists("SCP815_190921_19_velo_unspliced_matrix.rds")) {
    download.file("https://figshare.com/ndownloader/files/44865751",
                  "SCP815_190921_19_velo_unspliced_matrix.rds")
}
```

```{r load}
library(cytosignal)

vcs <- readRDS("SCP815_190921_19_CytoSignal.rds")
velo.s <- readRDS("SCP815_190921_19_velo_spliced_matrix.rds")
velo.u <- readRDS("SCP815_190921_19_velo_unspliced_matrix.rds")
```

Add the velocity matrices to the CytoSignal object using the `addVelo()` method. The method takes the CytoSignal object and the spliced and unspliced velocity matrices as input.

```{r addVelo}
vcs <- addVelo(vcs, velo.s = velo.s, velo.u = velo.u)
```

## Imputation

Similar to the CytoSignal workflow, we impute the velocity estimates for the ligands and receptors from the spatial neighbors of each location using function `imputeVeloLR`. We use the same strategy as we did for calculating LRscore in the CytoSignal workflow, but with the RNA velocity estimates. Hence, there is no need to do the neighbor finding here because it is already done during the CytoSignal workflow (with `findNN()`.

```{r impute}
vcs <- imputeVeloLR(vcs)
```

## Calculating velocity estimates for signaling interactions

Next, we call `inferIntrVelo()` to calculate velocity estimates for each signaling interaction. Because the CytoSignal test statistic is simply $S = LR$, we can calculate the time derivative of $S$ using the product rule of differential calculus: $dS/dt = RdL/dt + L dR/dt$. 

Similarly to the CytoSignal workflow, the velocity estimates for diffusible and contact-dependent interactions are treated differently. By default, they are derived using `"GauEps-Raw"` and `"DT-Raw"` strategies, respectively. When the sparsity of the data is relatively high, the receptor expressions can be smoothed with Delaunay Triangulation with setting `recep.smooth = TRUE` and hence triggering `"GauEps-DT"` and `"DT-DT"` mode. Please consult the CytoSignal tutorial for explanation on these tags.

```{r veloLR}
vcs <- inferIntrVelo(vcs)
```

To access the exact value for the LRvelo estimates, users can explore the object with the following command with replacing `<nn.type>` by strings `"GauEps-Raw"` for diffusion-dependent interactions, or `"DT-Raw"` for contact-dependent interactions. The returned result is a sparse matrix with barcodes as rows and interaction ID as columns and values indicating the velocity estimates for each interaction in each cell.

```{r veloLRRes, eval=FALSE}
vcs@lrvelo[[<nn.type>]]@intr.velo
```

## Visualization

VeloCytoSignal makes cellular-resolution, spatially-resolved signaling inference. We developed several new visualizations for plotting each inferred signaling interaction.

Most of the functions take two arguments, `intr` and `slot.use`, for specifying individual interactions. Users can first show a list of significant interactions that is available for visualization with `showIntr()`. `intr` should then be the unique ID of available interaction(s) and `slot.use` should be a selection as explained above.

```{r showIntr, results=TRUE}
# Show significant diffusible interaction IDs associated with interaction names
allIntrs <- showIntr(vcs, slot.use = "GauEps-Raw", return.name = TRUE)
print(head(allIntrs))
```

Setting `return.name = TRUE` displays both interaction unique IDs and the interaction names for easier interpretation. Interaction names are shown as a "ligand-receptor" gene symbol. Note that the list of available interactions varies for different LRscore inferences. For listing significant interactions with diffusible ligands, most of the time we have `slot.use = "GauEps-Raw"` if users did the LRscore inference with default setting. Similarly, for finding the significant contact dependent interactions, we need to set `slot.use = "DT-Raw"` if everything above went with default. 

In the plotting functions we show next, users need to specify the unique ID(s) of interests and the LRscore slot, where the interaction is found, at the same time. This is for ensuring that the correct result is captured when multiple inferences with parameter tweaks are performed for the same type of interactions.

### Visualizing velocity estimates of interaction signaling in 3D

We developed a way (`plotVelo`) to show the velocity estimates of an interaction at cellular resolution in 3D, with the directions of the arrows indicating whether the signaling activity is currently increasing (pointing upward) or decreasing (pointing downward) at each spatial position.

```{r plotVelo}
intr.use <- names(allIntrs)[1]
plotVelo(vcs, intr.use, slot.use = "GauEps-Raw", pt.size = 0.2, arrow.width = .1)
```

### Combination plots of 3D velocity, gene expression, LRscore and cluster annotation

We provide a function named `plotSignif2` for making a general combination figure. Specify `edge = TRUE` and `velo = TRUE` to include the 3D edge plot and the velocity estimates besides the gene expression and LRscore. 

```{r signif, fig.width=9, fig.height=10}
plotSignif2(vcs, intr = intr.use, slot.use = "GauEps-Raw", return.plot = TRUE, 
            edge = TRUE, velo = TRUE, pt.size = 0.3)
```

**Recommended:** When plotting a large number of interactions, it’s recommended to use numeric index with `intr` and turn to `return.plot = FALSE` (by default). In this way, plots of all input interactions at high resolution will be stored on the device instead of on the screen.

The codes below show the top 5 significant interactions ranked by SPARK-X for both diffusion-dependent and contact-dependent interactions.

```{r plotsignif2, eval=FALSE}
# For diffusion dependent interactions
## Choose a level of significance metric
signif.use <- "result.spx"
lrscore.slot <- "GauEps-Raw"
plot_dir <- paste0("path_to_result/Velo-diff-dep_", signif.use, "_", lrscore.slot, "/")
plotSignif2(vcs, intr = 1:5, slot.use = lrscore.slot, signif.use = signif.use, plot_dir = plot_dir)

# For contact dependent interactions
## Choose a level of significance metric
signif.use <- "result.spx"
lrscore.slot <- "DT-Raw"
plot_dir <- paste0("path_to_result/Velo-cont-dep_", signif.use, "_", lrscore.slot, "/")
plotSignif2(vcs, intr = 1:5, slot.use = lrscore.slot, signif.use = signif.use, plot_dir = plot_dir)
```


```{r cleanup, include=FALSE}
unlink("csSignifPlot", recursive = TRUE)
```
