Infer spatially resolved temporal dynamics of cell-cell communication at cellular resolution

Please run the CytoSignal workflow 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, and is now publicly available on Single-Cell Portal via SCP815.

To calculate RNA velocity estimates, we utilized our previously published tool, 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
• s matrix stands for the spliced RNA velocity estimates: SCP815_190921_19_velo_spliced_matrix.rds
• u matrix stands for the unspliced RNA velocity estimates: SCP815_190921_19_velo_unspliced_matrix.rds

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.

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.

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().

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 + LdR/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.

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.

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.

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

##           CPI-CS0AB19226A           CPI-CS0AA769F5B           CPI-CS0D238C22B 
## "SEM3A-PlexinA4_complex1" "SEM3A-PlexinA3_complex1" "SEM3A-PlexinA2_complex1"

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.

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

plotSignif2(vcs, intr = intr.use, slot.use = "GauEps-Raw", return.plot = TRUE, 
            edge = TRUE, velo = TRUE, pt.size = 0.3)

## $`SEM3A-PlexinA4_complex1`

# 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)