Inferring a spatial code of cell-cell interactions across a whole animal body

Abstract

Cell-cell interactions shape cellular function and ultimately organismal phenotype. Interacting cells can sense their mutual distance using combinations of ligand-receptor pairs, suggesting the existence of a spatial code, i.e., signals encoding spatial properties of cellular organization. However, this code driving and sustaining the spatial organization of cells remains to be elucidated. Here we present a computational framework to infer the spatial code underlying cell-cell interactions from the transcriptomes of the cell types across the whole body of a multicellular organism. As core of this framework, we introduce our tool cell2cell, which uses the coexpression of ligand-receptor pairs to compute the potential for intercellular interactions, and we test it across the Caenorhabditis elegans' body. Leveraging a 3D atlas of C. elegans' cells, we also implement a genetic algorithm to identify the ligand-receptor pairs most informative of the spatial organization of cells across the whole body. Validating the spatial code extracted with this strategy, the resulting intercellular distances are negatively correlated with the inferred cell-cell interactions. Furthermore, for selected cell-cell and ligand-receptor pairs, we experimentally confirm the communicatory behavior inferred with cell2cell and the genetic algorithm. Thus, our framework helps identify a code that predicts the spatial organization of cells across a whole-animal body.

Fig 1. Calculation of the modified Bray-Curtis CCI score.

(A) To represent the overall interaction potential between cell A and cell B, our CCI score is computed from two vectors representing the ligands and receptors independently expressed in each cell. If only the ligands from one cell and the cognate receptors on the other are considered (“Cell A to Cell B” half or “Cell B to Cell A” half, independently), the score would be a directed score for representing the interaction (one cell is the sender and the other is the receiver). However, our score is undirected by considering both ligands and receptors of each cell to build the vector (both halves simultaneously, indicated with the yellow rectangle on the left). Thus, the vector of each cell is built with both directed halves of molecule production (e.g., top half possess ligands of cell A while the bottom half considers its receptors, generating a unique vector with both the ligands and the receptors of cell A). (B) Toy examples for computing our score for the interaction of Cell A and Cell B. Here, both possible directions of interaction are represented to show that they result in the same (undirected) score.

Fig 2. Cell-cell interactions and communication in C. elegans.

(A) Heatmap of CCI scores obtained for each pair of cell types using the curated list of LR pairs. An agglomerative hierarchical clustering was performed on a dissimilarity-like metric by taking the complement (1-score) of CCI scores, disregarding autocrine interactions. Cell types are colored by their lineages as indicated in the legend. Lineages and colors were assigned previously [16]. (B) UMAP visualization of CCIs. Dots represent pairs of interacting cells and they were projected based on their Jaccard distances, which were computed from the LR pairs expressed in the directed interactions between cells (one cell is producing the ligands and the other the receptors). Dots are colored by either the sender cell (left) or the receiver cell (right), depending on their lineages as indicated in the legend of (A). A readable version of the data used for this projection is available in S2 Table, where names of LR pairs and their communication scores are specified for each cell pair. Another UMAP visualization based on a more appropriate similarity metric is available in S1 Fig, which uses the Rand index that accounts for both active and inactive LR pairs. Using the Rand index still represents the same behavior of sender cells driving similarities. (C) Receiver operating characteristic (ROC) curves of random forest models for classifying cell-cell pairs from their CCI scores computed with different approaches as indicated in the legend. These models predict the intercellular distance range (short-, mid-, or long-range distance, as defined in the N1 Fig in S1 Text). For each classifier, the mean (solid line) ± standard deviation (transparent area) of the ROCs were computed with 3-fold stratified cross validations. The area under the curve (AUC) for the ROC curves is shown in the legend, detailing the mean ± standard deviation across all distance-range classifications. Separate evaluations for the distance ranges are provided in S3 Fig.

Fig 3. Changes in the relative abundances of signaling functions across initial GA-LR pairs.

Boxplots summarizing the changes of the relative abundances for each of the signaling functions that LR pairs are associated with (y-axis). Changes were computed from the fold change (FC) between the relative abundance in each of the 100 runs of the genetic algorithm (GA) with respect to the corresponding relative abundance in the complete list of LR pairs (S1 Table), and shown as the log10(FC+1) transformation (x-axis). Here, relative abundance is the number of LR pairs involved in a given pathway with respect to the total number of LR pairs in the list of GA-LR pairs. A two-tailed Wilcoxon’s test was performed to evaluate the significance of the changes. An adjusted P-value is reported to the right of each boxplot (FDR < 1%). All GA runs are shown in each boxplot (gray dots); dashed-gray lines indicate a change of at least 50% either decreasing (left line, FC = 0.5) or increasing (right line, FC = 1.5), while the dashed-red line indicates the value of no change (FC = 1).

Fig 4. CCI analyses based on LR pairs associated with intercellular distances.

(A) Heatmap of CCI scores obtained for each pair of cells using the consensus GA-LR pairs. An agglomerative hierarchical clustering was performed on a dissimilarity-like metric by taking the complement of CCI scores (1-score), excluding autocrine interactions. Cell types are colored by their lineage groups as indicated. (B) Heatmaps representing the presence or absence of ligands (left) and receptors (right) after expression thresholding (>10 TPM) in sender and receiver cells, respectively. Lines at the center connect ligands with their cognate receptors according to the GA-selected interactions. Cell types are colored as in (A).

Fig 5. Spatial enrichment and depletion of communication pathways.

(A) To study the anteroposterior use of communication pathways, the body of C. elegans was divided into three sections along the anteroposterior axis (top) and cell-type composition of each section (bottom) given a previously published 3D atlas. The mid-body section is defined by the presence of the intestine cells, and the head and tail are the anterior and posterior sections to it, respectively. Cells in the 3D atlas (top) are colored according to the cell types as delineated in barplots (y-axis, bottom). (B) Enrichment/depletion (FDR < 1%) of ligand-receptor pairs (y-axis) in each of the three sections (x-axis), calculated from their usage across all pairs of cells of each section. Communication pathways are also colored by their annotated functions (left column). (C) Circos plots for representing the importance of cell-cell communication occurring at different ranges of distance. A Fisher exact test was performed to find enriched/depleted LR pairs among all pairs of cells for a given proximity. The ranges of distances were defined as explained in Fig N1C. Nodes represent ligands or receptors and edges connect those ligands and receptors that interact in the GA-LR pairs (S3 Table). The color of the nodes represent whether they are ligands or receptors and the color of the edges indicate the negative value of the logarithmic transformation on the Benjamini-Hochberg adjusted P-values, according to the colored bar at the bottom. Interactions that resulted significantly enriched or depleted (FDR < 1%) are equivalent to the color assigned to a value of 2.0 or bigger.

Fig 6. Validation of the spatial expression of specific GA-LR pairs.

Single-molecule Fluorescent In Situ Hybridization of genes encoding three GA-LR pairs in C. elegans L2 larvae. (A) Intestinal/rectal cells expressing arg-1 (magenta) and non-seam hypodermal cells (arrow) expressing lin-12 (green) are adjacent (see rectangle in the merge channel and S1 Movie). (B) Non-seam hypodermal cells expressing let-756 (magenta) and amphid sheath cells (arrows) expressing ver-1 (green). Amphid sheath cells are surrounded by hypodermal cells (see rectangle in the merge channel and S2 Movie). (C) Seam cells in the tail (arrows) expressing lin-17 (magenta) and non-seam hypodermal cells in the tail (arrowheads) expressing lin-44 (green). The two genes are expressed in proximal cells (see ends of rectangle in the merge channel and S3 Movie). In all cases a DAPI staining was performed to distinguish cell nuclei. Scale bar = 10μm.