Applying GCNG to spatial transcriptomics data

While a number of recent methods have been suggested for spatial profiling of single cell RNA-Seq data [7–11], we decided to focus on the two methods that currently provide expression levels for the most number of genes in such experiments. The first is seqFISH+ [16]. We tested two datasets that used seqFISH+. The first contained information on the expression of 10,000 genes in 913 cells in the mouse cortex and the second profiled 2050 cells in mouse olfactory bulb (OB) in seven separate fields of view. The second method we used is MERFISH [17] for which we analyzed a dataset consisting of 10,050 genes in 1368 cells. Unlike the seqFISH+ data that profiled the expression in-vivo, the MERFISH data is from in vitro cultured cells and so does not include a diverse set of cell types. Still, as the authors of the MERFISH paper observed, even within this population, there are differences in spatial expression and so the data can be used to study extracellular gene-gene interaction. We normalized the expression data such that expression levels for all genes in each cell sum to the same value as was previously done [16]. See the “Methods” section and Additional file 1 for complete details on both datasets.

GCN requires labeled data for supervised training. While the exact set of signaling interactions between cells in the spatial data we studied is unknown, we used as true interactions a curated list of interacting ligands and receptors [20]. Ligands are proteins that are secreted by cells and they then interact with membrane receptor proteins on the cell itself or on neighboring to activate signaling pathways within the receiving cell [21]. See the “Methods” section for complete details on the positive and negative pairs used for training.

For evaluation, GCNG adopted a tenfold cross-validation. Train and test sets were completely separated to avoid any information leakage (Fig. 1b). See the “Methods” section and Additional file 1 for details.

GCNG correctly infers ligand-receptor interactions between cells

We first evaluated GCNG’s ability to predict ligand-receptor interactions. For this, we used two datasets. The first is seqFISH+ mouse cortex tissue which contains the expression of 10,000 genes in 913 cells. Our labeled set consisted of 1056 known interactions between 309 ligands and 481 receptors. The second is a MERFISH dataset with 10,050 genes from 1368 cells, and 841 known interactions between 270 ligands and 376 receptors.

We enforced a strict separation between the training and test sets in the 10-fold cross-validation (CV) (Fig. 1b). Negative pairs were also ligand-receptor and were randomly selected from non-interacting training (test) data genes. We also used the 10-fold CV to select hyper-parameters to determine the neighborhood for each cell (Methods). We compared three possible GCNG models: diagonal GCNG that only uses a diagonal matrix to represent the graph so that only autocrine interactions are possible, exocrine GCNG where only exocrine interaction between cells are allowed, and autocrine plus (+) GCNG that allows for both autocrine and exocrine interactions. To evaluate the performance of GCNG, we compared it to a number of prior methods that were recently used to predict genes involved in extracellular gene interactions from spatial expression data. These include computing the spatial Pearson correlation (PC) between ligand and receptors in neighboring cells and Giotto [12] which first calculates a similarity score for all pairs of genes in all pairs of neighboring cell types and we then rank pairs based on their average score.

We also compared GCNG to two alternative methods that do not use spatial information at all to determine the contribution of neighborhood data. These included Pearson’s correlation between the expression of ligand and receptors within each cell [22] and our diagonal GCNG method with only autocrine interactions. Finally, we compared GCNG performance on the real data to results when applied to permutation of both, the neighborhood information for each cell and the set of interacting ligand-receptors used for training and testing. We also tested additional variants of GCNG including variants that utilized cell type information (encoded as a node attribute), edge weight (using the distance between cells), and variants using other GNN architectures including EdgeconditionConv [23] and graph attention [24] (Additional file 1: Fig. S1).

Results are presented in Fig. 2. As can be seen, GCNG achieved the best results in both datasets (Fig. 2a–d). Specifically, for seqFISH+ cortex data, autocrine+, diagonal, and exocrine GCNG reached mean (median) AUROC/AUPRC of 0.65/0.73 (0.99/1.0), 0.59/0.70 (0.99/1.0), and 0.60/0.69 (0.99/1.0), respectively. In contrast, for this data spatial PC, Giotto and single cell PC all performed much worse with mean (median) AUROC/AUPRC of 0.54/0.65 (0.75/0.79), 0.45/0.58 (0.25/0.33), and 0.48/0.60 (0.38/0.38), respectively. For MERFISH data, autocrine+, diagonal, and exocrine GCNG reached mean (median) AUROC/AUPRC of 0.69/0.76 (0.99/1.0), 0.60/0.69 (0.99/1.0), and 0.61/0.71 (0.75/0.79), respectively, again improving on the other methods we compared to. See also Additional file 1: Figs. S2&3 for detailed performance values. Overall, for both datasets, GCNG achieves a relative improvement of at least 20% for mean AUROC/AUPRC when compared to prior methods. In addition, the fact that autocrine+ GCNG outperformed diagonal GCNG for both datasets confirms the importance of spatial information for this task.

Open in a separate window

Fig. 2

GCNG correctly infers extracellular ligand-receptor interaction. a, b AUROC and AUPRC curves for seqFISH+ using autocrine+ GCNG model. Here, each gray line represents results for one ligand (a total of 91 curves), red line represents the median curve, and the light green part represents the region between 40 and 60 quantile. Median and mean of area under the curves are shown on top of each panel. c, d AUROC and AUPRC curves for MERFISH when using the autocrine+ GCNG model (73 curves). e, f Overall comparison in terms of AUROC and AUPRC for single cell Pearson correlation, Giotto, spatial Pearson correlation, diagonal GCNG, autocrine+ GCNG, and exocrine GCNG models. For Giotto on MERFISH (f), we set its all AUROC and AUPRC as 0.5 since the data only has one cell type

To test if the interactions identified are likely active in the tissue tested, we further compared the performance when using the real interaction data to running GCNG on a permuted training and test dataset (in which we permute the set of interacting ligand-receptor pairs). We observed a large drop in performance when using the randomized data (autocrine plus GCNG in terms of mean AUROC for MERFISH (seqFISH+), real vs. random: 0.69 vs. 0.55 (0.65 vs. 0.52); exocrine GCNG vs. random: 0.58 vs. 0.50 (0.60 vs. 0.43). See Additional file 1: Fig. S1 for detailed comparison results on both datasets.

Analysis of co-expression patterns identified by GCNG

To further explore the predictions of our GCNG and to map them back to the original spatial representation, we looked at some of the top correctly predicted pairs. For each such pair (a ligand and receptor predicted to interact), we projected their expression on the cell distribution figures. Figure 3 presents such projection for two local ligands (col4a1 and lamc1) with their positive and negative receptor partners in seqFISH+ cortex data (Fig. 3a, b for col4a1, and c, d for lamc1) Since here the genes are fixed while cells need to be selected (in contrast to common cases where for each cell highly expressed genes are selected), for a gene, cells are defined to “highly expressing” it if the expression of the gene in that cell is in the top 100 expression levels for that gene among all cells. For the positive col4a1-cd93 pair, cells highly expressing col4a1 and cd93 are both concentrated in the 1st and 5th fields, which have the most cells highly expressing the ligand or receptor genes (see Additional file 1: Fig. S4A and B for plots of all fields of view). In contrast, for the negative col4a1-hrh3 pair, cells highly expressing hrh3 do not seem to reside next to cells expressing col4a1. Similar pattern comparison can also be observed for ligand lamc1 with positive (itgb1) and negative receptor (lyve1) (Fig. 3c, d (see Additional file 1: Fig. S4C and D for all fields of view)). The ability of GCNG to predict such interactions based on a subset of the data highlights the usefulness of this approach compared to global analysis methods including PC. Cell type plots (Additional file 1: Fig. S5) indicate that correctly predicted pairs can be found in both, neighboring cells from the same type and cells from different types. These results indicate that the GCNG method can generalize well and can be used to correctly identify several different types of interactions.

Open in a separate window

Fig. 3

Spatial expression patterns for selected ligand-receptor pairs. a, b Spatial expression distribution of correctly predicted positive (cd93) and negative (hrh3) receptors for ligand col4a1. Two of the 7 fields of view (FOV) profiled are shown (see Additional file 1: Fig. S4 for all FOVs). As can be seen, the correctly predicted pairs are indeed much better spatially correlated than the negative pair. c, d Spatial expression distribution of correctly predicted positive (itgb1) and negative (lyve1) receptors for ligand lamc1. Cells highly expressing lamc1 and itgb1 are both concentrated in the 5th and 7th fields as shown

Inferring causal interactions

While correlation-based methods can be used to identify gene co-expression interactions and networks [25, 26], these methods cannot be used to infer causality since their outcome is symmetric. Causality information may be trivial for ligand receptors, since the direction for such pair is known. However, for other interacting genes across cells, the direction is often not clear. Thus, a method that can infer both interactions and directionality may be beneficial for studying spatial transcriptomics data. Unlike prior unsupervised methods, our supervised framework can be trained to identify causal interactions based on the pair-wise spatial expression pattern if training data exists, inspired by a recent work [27]. We thus trained a GCNG on a subset of known causal pairs (ligand and receptors) and then used it to predict directionality for other pairs. To generate train and test data for this, for each known ligand-receptor (L_a,R_b) gene pair, we introduce a negative sample (R_b,L_a) with label 0. The same tenfold cross-validation strategy is used to evaluate GCNG’s performance here. Results are presented in Fig. 4. As can be seen for seqFISH+ cortex and MERFISH datasets, GCNG performs well on this task with mean (median) AUROC/AUPRC of 0.636/0.728 (0.99/1.0) and 0.642/0.734 (0.99/1.0), respectively. Thus, for top predicted pairs, the direction predictions of GCNG can be used to further assign causality.

Open in a separate window

Fig. 4

GCNG can infer the direction of extracellular ligand-receptor interactions. a, b AUROC and AUPRC of GCNG for the direction prediction task using the seqFISH+ dataset. c, d AUROC and AUPRC of GCNG for the direction prediction task using the MRFISH dataset

Graph representation for spatial transcriptomics data

To determine the neighbors of each cell, we calculated the Euclidean distance in the image coordinates for all cells and used a distance threshold to select neighbors. The threshold value was selected using 10-fold cross-validation (see “Train and test strategy” section below for details), which for the 2D images, we used seemed to represent the number of neighbors that were in physical contact with the cell (Additional file 1: Tabs. S1&S2). Consider the seqFISH+ cortex data as an example. Given the set of neighbors, we constructed an adjacency matrix size of 913 × 913 (where 913 is the number of cells in the seqFISH+ dataset), which we term A. In other words, for the symmetric network, A, A_ij=A_ji=1 if i and j are neighbors and 0 otherwise. Using the adjacency matrix A, the normalized (symmetric) Laplacian matrix L_N is defined as [48]:

where $D_{\mathit{ii}}=\sum _j{A}_{\mathit{ij}}$, and I is the identity matrix. L_N can be used as a graph matrix. We also tried other graph representation matrices in GCNG. One such method which we adopted first normalizes the adjacency matrix to get A_N,

where $D_{\mathit{ii}}=\sum _j{A}_{\mathit{ij}}$,

and then computes the normalized Laplacian using the normalized adjacency matrix, L_NN,

where $D_{N_{\mathit{ii}}}=\sum _j{A}_{N_{\mathit{ij}}}$.

Finally, we also tested the following formulation of a graph matrix:

where $D_{\mathit{ii}}^{\prime }=\sum _j{A}_{\mathit{ij}}^{\prime }$, and A^′=A+I.

In this paper, we use the graph matrix L_NN in the exocrine GCNG, and L^′ in the autocrine+ GCNG model, and diagonal matrix for diagonal GCNG. See Additional file 1 for a detailed discussion about graph representations.

Labeled data

GCN requires labeled data for supervised training. While the exact set of signaling interactions between cells in the cortex data we studied is unknown, we used as true interactions a curated list of interacting ligands and receptors [20], consisting of 708 ligands, and 691 receptors with 2557 known interactions. Of these, 309 ligands and 481 (involved in 1056 known interactions) are profiled by the seqFISH+ datasets we studied, and 270 ligands, 376 (841 known interactions) are included in the MERFISH dataset. These were used for training and testing. See “Train and test strategy” section below for details on how we define and generate positive and negative examples during training.

GCNG network architecture

To construct GCNG, we used the python packages of “spektral,” “Keras,” and “Tensorflow.” See Fig. 1a for the architecture of GCNG. GCNG uses two types of input: one encodes the graph structure, L_NN of size of 913× 913 as discussed above, while the second encodes node-specific values, and is generated as the normalized expression of a pair of candidate genes using a matrix of dimension of 913 × 2. GCNG consists of two 32-channel graph convolutional layers, one flatten layer, one 512-dimension dense layer, and one sigmoid function output layer for classification. The graph convolutional layer is defined as:

where X is the expression matrix with dimension of 913 × 2 and L is the 913 × 913 graph matrix. W is a weight matrix of filters (also termed the convolution kernel) with a size of 2 × 32, where 2 corresponds to the two-dimension gene expression of each node, and 32 represents 32 filters or feature maps. b is the bias vector term with a size of 1 × 32. The “elu” (exponential linear unit) function is defined as:

where α=1 by default. Here, Z represents the embedding vectors of all cell nodes with a size of 913 × 32. Note that we are using two convolutional layers here allowing the method to learn indirect graph relationships as well. Since the impact of secreted proteins can be larger than just direct neighbors, such an approach allows the method to infer interactions that may be missed by only considering direct neighbors.

The first graph convolutional layer combines the two inputs and converts them to embedding vectors for cell nodes of dimension of 913 × 32. The second graph convolutional layer combines the embedding vector of each cell with the one learned for its direct neighbors. The flatten layer then converts the matrices generated by the second layer to a vector using ReLU activation function. Finally, a dense layer with one-dimensional output is used to predict the interaction probability based on the sigmoid activation function. The activation functions used by the different layers are defined below.

And the objective function for the entire GCNG model is as follows:

where i represents the ith sample, y_i represents the label for the ith sample, and Θ represents all parameters that need to be optimized in GCNG.

We tested one, two, and three graph convolutional layer networks and determined that two layers network led to the best performance. In addition, we also tried GCNG with cell type information as node attribute using one-hot encoding, distance value as edge attributes, and other GNN architectures including “EdgeconditionConv” model [23] and graph attention model [24], see detailed results for these in Additional file 1: Fig. S1.