SpaCcLink: exploring downstream signaling regulations with graph attention network for systematic inference of spatial cell–cell communication

Abstract

Background

Cellular communication is vital for the proper functioning of multicellular organisms. A comprehensive analysis of cellular communication demands the consideration not only of the binding between ligands and receptors but also of a series of downstream signal transduction reactions within cells. Thanks to the advancements in spatial transcriptomics technology, we are now able to better decipher the process of cellular communication within the cellular microenvironment. Nevertheless, the majority of existing spatial cell–cell communication algorithms fail to take into account the downstream signals within cells.

Results

In this study, we put forward SpaCcLink, a cell–cell communication analysis method that takes into account the downstream influence of individual receptors within cells and systematically investigates the spatial patterns of communication as well as downstream signal networks. Analyses conducted on real datasets derived from humans and mice have demonstrated that SpaCcLink can help in identifying more relevant ligands and receptors, thereby enabling us to systematically decode the downstream genes and signaling pathways that are influenced by cell–cell communication. Comparisons with other methods suggest that SpaCcLink can identify downstream genes that are more closely associated with biological processes and can also discover reliable ligand-receptor relationships.

Conclusions

By means of SpaCcLink, a more profound and all-encompassing comprehension of the mechanisms underlying cellular communication can be achieved, which in turn promotes and deepens our understanding of the intricate complexity within organisms.

Background

In multicellular organisms, cell–cell communication is of fundamental importance for sustaining the harmonious and coordinated operation of tissues, organs, or systems. Cellular communication empowers cells to identify one another, interchange information, and collaborate in a synchronized manner to fulfill a diverse array of biological tasks, such as cell differentiation, proliferation, development, and the orchestration of immune responses [1, 2]. Thus, the dissection of cell–cell interactions holds utmost significance in comprehensively grasping the structure, functionality, and overall “landscape” of biological systems. It offers crucial revelations into the elaborate network of intercellular communication and coordination and also plays a pivotal part in deciphering disease mechanisms and guiding therapeutic strategies.

Cell–cell interactions transpire via a multiplicity of mechanisms, encompassing autocrine, paracrine, juxtacrine, and endocrine signaling [3, 4]. Autocrine and juxtacrine interactions entail short-range effects, whereas paracrine and endocrine interactions can exert their influence over more extended distances. Within tissues, signaling cells can attach to receptors on recipient cells by means of specific ligands, thereby instigating a sequence of signal transduction reactions. These signal transduction reactions are capable of either activating or suppressing particular signaling pathways, ultimately resulting in alterations in the activity of transcription factors and the expression of target genes [5,6,7,8]. Subsequently, these modifications in transcription factor activity and target gene expression further exert an impact on the functionality and physiological processes of the cells.

In recent times, remarkable headway has been achieved in deciphering cell–cell communication from single-cell transcriptome data, attributable to the ceaseless progression of single-cell transcriptomic sequencing (scRNA-seq) technology [9,10,11,12]. A number of methods have been devised to analyze and fathom cell–cell communication. These methods can be generally classified into two categories of algorithms: those involved in deducing intercellular communication between diverse cell types [13,14,15,16] and those considering intracellular signaling [17,18,19,20]. The algorithms utilized for inferring intercellular signaling depend on pre-existing databases and communication scores to single out significant ligand-receptor pairs. Moreover, the methods for considering intracellular signaling can be further partitioned into two classifications: ones based on established signal pathway annotations [17, 18, 20] and ones based on de novo inference [19]. Notwithstanding the relative maturity of these methods founded on scRNA-seq data, there remains a necessity to further probe into how spatial information can be harnessed in the exploration of cell–cell communication. This is because cell–cell communication customarily takes place within a specific range, and incorporating spatial distance information can assist in alleviating potential false-positive outcomes that might emerge from non-spatial communication methods [21].

Spatial transcriptome technology offers a powerful tool for investigating cell–cell communication within the tissue microenvironment [22, 23]. It enables us to obtain a more comprehensive understanding of cell–cell interactions, uncover the spatial dependence of such communication in tissues, and further deepen our insights into cell functions and regulatory mechanisms in the tissue microenvironment. Recently, multiple methods have been developed to study spatial cell–cell communication based on spatial transcriptome data. For instance, SpaTalk [24], Giotto [25], and Niches [26] are designed to identify neighborhood communication. SpaCI [27] and GCNG [28] leverage graph convolutional networks to detect adjacent potential communication relationships. HoloNet [29] employs a graph network model to analyze the impact of communication on target gene expression, while COMMOT [30] utilizes the collective optimal transmission method to analyze the competitive relationships in communication.

Although there are some methods that do take downstream signaling into consideration in the context of spatial communication, the current ones [24, 29,30,31,32] fail to make full use of gene networks and gene expression information when dealing with downstream signaling. Traditional non-spatial communication inference methods also have the same limitations since they do not adequately incorporate gene expression information and the pathways of signal transduction. In the study of cellular communication and regulation, taking the downstream signaling process within the cell into account can help us understand the regulatory mechanism of cell communication in a more comprehensive manner [33].

Consequently, in order to conduct a more precise analysis of the downstream functions influenced by cell–cell communication and identify the corresponding spatial communication patterns, we have put forward SpaCcLink, a spatial cellular communication analysis method that takes into account the downstream signals of individual cells and systematically explores both the spatial cellular communication patterns and downstream signal networks. Leveraging the prior gene network, the graph attention network is employed to identify distinct groups of target genes that are highly correlated with the receptors. By combining the specific expressions of these target genes, the activity fraction of receptors in individual cells can be quantified to analyze the communication patterns. Additionally, the Cross Moran index is incorporated to identify spatially dependent signal molecules.

To assess the performance of SpaCcLink, we compare it with other methods using real datasets. The results reveal that the downstream signals identified by our model are of greater significance. Moreover, the receptor pairs identified through the Cross Moran index also demonstrate strong reliability. In summary, SpaCcLink represents a novel approach for the downstream functional analysis of spatial cell–cell communication. It facilitates the identification of key signaling molecules and communication patterns, thereby contributing to a deeper understanding of cell localization and function within tissues.

Results and discussion

Overview of SpaCcLink

The schematics of SpaCcLink are shown in Fig. 1. Firstly, the downstream influence of each receptor within individual cells is taken into account. A graph attention network model is employed to pinpoint a group of target genes that exhibit a high degree of correlation with the receptor (Fig. 1a). Subsequently, the downstream influence scores of each receptor within every cell are derived by integrating the gene expression specificity scores (Fig. 1c).

[IMAGE OMITTED: SEE PDF]

Secondly, the Wasserstein distance is calculated for all ligand-receptor pairs. Given that the existing prior ligand-receptor (LR) databases fail to offer detailed accounts of the interaction modalities between receptors and ligands, SpaCcLink resorts to an optimal transport-based approach to categorize LR pairs into short-range and long-range interactions. By making use of d_ratio and p-values (as detailed in the “ Methods” section), the interaction pairs are sorted into significant short-range, significant long-range, and the default medium-range classifications. Grounded on this classification, cross-Moran’s I index is applied to single out ligand-receptor pairs with conspicuous spatial dependencies (as illustrated in Fig. 1b and Fig. 1c). These two sequential steps collaboratively facilitate the identification of ligand-receptor pairs that exhibit both remarkable spatial dependencies and powerful downstream effects.

Finally, a downstream analysis is carried out on the identified ligand-receptor pairs (Fig. 1c). This analysis encompasses several key components: decomposing the receptor impact scores through non-negative matrix factorization (NMF) to obtain corresponding communication patterns, dissecting the intercellular communication among different cell types, identifying the downstream target genes that are primarily affected by the ligand-receptor pairs with the aid of Sankey diagrams, conducting Fisher’s exact test to analyze the crucial transcription factors and visualizing the entire signal transduction pathways from receptors to target genes.

Identification of communication patterns and significant LR pairs in human melanoma data

We initially applied our methodology to human melanoma data derived from the ST platform. This dataset encompasses 293 spots and comprises 7 distinct cell types (as depicted in Fig. 2a.). We computed the downstream influence scores for each receptor and integrated them with spatial communication analysis. This enabled us to identify receptors that possess both significant downstream influence and conspicuous spatial dependence on ligands. Through the execution of standard NMF on the downstream influence scores of receptors, two patterns of downstream signaling communication within the melanoma data were obtained (shown in Fig. 2b). It is evident that Pattern 0 predominantly aligns with the distribution area of melanoma cells, whereas Pattern 1 mainly corresponds to B cell regions and certain boundary zones between cancer-associated fibroblasts (CAF) and melanoma (as illustrated in Fig. 2a and b). To further validate the efficacy of identifying communication patterns via the downstream signaling scores of receptors, we carried out supplementary dimensionality reduction and Leiden clustering on the outcomes obtained through the traditional scoring approach (see the “ Methods” section). By regulating the resolution, we acquired results that were highly congruent with the communication patterns identified by the downstream scores (as presented in Fig. 2c).

[IMAGE OMITTED: SEE PDF]

We subsequently acquired the top genes for each pattern by relying on the decomposed H-matrix (Fig. 2d). This matrix delineates the significance or influence of each feature within the original data, principally aiming to gauge the extent of contribution to the pattern. In pattern 1, the genes that rank at the top and exhibit high specificity are ALOX5, CD79A, and CXCR12. ALOX5 pertains to the human lipoxygenases family and is robustly correlated with immune cell populations [34, 35].

It orchestrates the signaling pathways or immune responses that impact these cell populations. Previous investigations have unearthed a correlation between ALOX5 and the prognosis of melanoma [36, 37]. The expression of CD79a is highly specific to B cells, making it one of the characteristic molecules of B cells. The CD79a/b heterodimer formed by the binding of CD79a and CD79b plays a crucial role in the normal functioning of B cells, including signal transduction, cell development, and survival [38, 39]. The expression distribution of CXC12 and CXCR4 manifested an evident spatial dependence (Fig. 2e), which intimated that CXCR4 predominantly received signals from CXCL12. The communication between CXC12 and CXCR4 is chiefly concentrated between B cells and B cells as well as between B cells and CAF (Fig. 2f). The interaction between CXCL12 and CXCR4 is highly expressed in B cells and is postulated to play a pivotal role in cell migration and proliferation [40]. These processes are intimately intertwined with cancer progression and metastasis, thereby designating CXCL12 and CXCR4 as potential targets for the treatment of skin malignancies, especially melanoma [41, 42]. In accordance with the GOBP (Gene Ontology Biological Process) enrichment analysis carried out on the top 20 genes in pattern 1, the outcomes signified a substantial association between pattern 1 and biological processes such as leukocyte cell adhesion regulation, cellular activation, and inflammatory response (Fig. 2g). These discoveries imply that pattern 1 might be implicated in the modulation of the immune system and inflammatory processes. This holds significant ramifications for further research into immune cell function and the disease mechanisms of melanoma.

Identification of prognostic receptor genes and downstream signals in human breast cancer data

Subsequently, we applied our methodology to human breast cancer data acquired from the 10X Visium platform. This dataset comprises 3798 spots and encompasses 9 cell types, including four malignant cell types, namely basal_like_1, basal_like_2, luminal_AR, and mesenchymal cells (as shown in Fig. 3a). In a similar vein, we initially conducted pattern recognition on the downstream signal scores of receptor genes that manifested significant effects and exhibited clear spatial dependency on ligands. Through this process, we identified four communication patterns, among which patterns 1–3 display distinct regional specificities. Pattern 1 is mainly associated with the area of immune cells; pattern 2 is predominantly concentrated between the malignant cell type basal_like_1 and stroma cells, potentially being related to the invasiveness and aggressiveness of the tumor; pattern 3 is mostly clustered in the area of malignant cells (illustrated in Fig. 3b). We then visualized the top 10 receptor genes with high load scores in pattern 1. Among these genes, CD247, CD3D, and CD48 demonstrated the regional specificity of pattern 1. Utilizing TCGA human breast cancer data, we carried out survival analysis on the top 10 receptor genes. The expression of specific receptor genes in pattern 1 significantly influenced the survival time of patients, whereas non-specific genes showed no significant correlation (Additional file 1: Fig. S1).

[IMAGE OMITTED: SEE PDF]

Next, we analyzed the communication relationship between CD3D and its spatially dependent ligand HLA-B across all cell types. The communication between HLA-B and CD3D primarily takes place within basal_like_1 cells and between basal_like_1 cells and macrophage cells. In our analysis, we discovered that the high communication score within basal_like_1 is mainly attributable to its relatively larger cell population. Consequently, we focused on the signaling process from basal_like_1 to macrophages. The signaling from basal_like_1 to macrophages is primarily mediated by ligand HLA-C and ligand HLA-B. These molecules act as signaling initiators and interact with receptor molecules CD3D and LILRB2 on the cell surface.

As signal receivers, CD3D and LILRB2 receive the signals from basal_like_1 cells and mainly transmit them to downstream target genes UBAP2L, VAMP3, and VCAN (shown in Fig. 3f). UBAP2L is associated with macrophage polarization [43]. VAMP3 plays a vital role in macrophages by regulating processes such as adhesion, spreading, and sustained migration, thereby having a significant impact on their function [44, 45]. Moreover, VAMP3 is involved in membrane transport pathways related to cytokine secretion and pathogen engulfment [46]. VCAN is widely expressed during immune inflammation processes and has been found to regulate macrophage activation and function [47, 48].

We performed Gene Ontology Biological Process enrichment analysis on the genes involved in the downstream signaling of CD3D and found that they were mainly enriched in the following functions: immune response-activating cell surface receptor signaling pathway, cellular response to cytokine stimulus, and cytokine signaling in the immune system (as presented in Fig. 3g). These results imply that CD3D may play a significant role in immune response and cytokine signaling. The study of these signaling pathways is essential for understanding the interactions between cancer cells and immune cells as well as their crucial roles in the development of breast cancer and the configuration of the tumor microenvironment.

Identification of key TFs and complete downstream signaling pathway in mouse brain data

We also applied our methodology to mouse brain data obtained from the 10X Visium platform. This dataset covered 2688 spots and encompassed 15 cell types (as shown in Fig. 4a). In a similar manner, we initially carried out pattern recognition on the downstream effect scores of receptors. Four patterns were identified, which were closely related to cell type distributions. Pattern 0 was predominantly clustered in the mouse cortex and pyramidal layer regions (as illustrated in Fig. 4a and b). Subsequently, we visualized the top 5 genes with high loading scores for each pattern. In pattern 0, the specific genes were identified as Robo2, Cckbr, and Cnr1 (as depicted in Fig. 4d). For each specific gene and its ligands that exhibited significant spatial dependence, we conducted a cell-type communication network analysis (as shown in Fig. 4c). Cckbr and its ligand CcK primarily communicate within Cortex_3 cells and between Cortex_1 and Cortex_3, which is in line with previous studies [49]. Cholecystokinin (CCK) is a peptide hormone and one of the most abundant neuropeptides in the vertebrate brain. CCKBR plays a preponderant role in the central nervous system, mainly in the neocortex and limbic structures. CCK has been recognized as a central modulator in neuronal circuits, and along with its receptors, it participates in the regulation of neural processes such as feeding, memory, pain perception, and exploratory behavior [50].

[IMAGE OMITTED: SEE PDF]

Next, we employed Fisher’s exact test to identify the downstream TFs activated by each receptor gene. Subsequently, we visualized the complete signal transduction network from the receptor to the transcription factors and further to the target genes (as shown in Fig. 4e).

FOXO1 and FOXO3 belong to the FOXO family, and they play crucial roles in neural development, cell survival, and metabolic regulation [51]. Trp53 plays a vital part in the nervous system, regulating the process of apoptosis [52]. NFKB1 regulates inflammatory and immune responses in neuronal cells, and its involvement is associated with the development of neurological disorders and cognitive impairments [53, 54]. By conducting GOBP enrichment analysis on downstream genes, we can further confirm the association of the identified functional transcription factors with cellular responses to organonitrogen compounds, responses to peptides, and behavior-related processes. Identifying key transcription factors and studying the complete signaling pathways can offer a profound understanding of the functionality and regulatory mechanisms of the mouse brain. This contributes to deciphering the interactions between cells, the organization and functionality of neural networks, and the mechanisms underlying complex physiological processes such as behavior.

Compared with other methods

We compared SpaCcLink with SpaTalk, NicheNet, and CellCall, assessing their performance in deducing the intracellular downstream signaling pathways within different cell types. Due to the complex nature of inputting prior LR data and pathway data from diverse methods, we decided to employ the specific prior data furnished by each method. It was discovered that a considerable portion of the LR data we amassed and arranged overlapped with the data utilized by other methods (Fig. 5a). Although NicheNet presented a more extensive range of LR pair information, a significant part still did not coincide with the other methods. We conducted pathway enrichment analysis on the downstream target genes of the receptors using the KEGG and Reactome databases. The significance of pathway enrichment was gauged based on the adjusted q-values acquired from the enrichment analysis. According to the outcomes presented in Fig. 5b, it was found that the downstream genes identified by SpaCcLink possess a stronger capacity to embrace receptor-related biological processes or pathways.

[IMAGE OMITTED: SEE PDF]

Furthermore, we also contrasted the performance of other cell communication methods in inferring LR interactions, such as CellChat, Giotto, COMMOT, and SpaTalk. Among these methods, only CellChat does not make use of spatial information. Giotto and SpaTalk depend on K-nearest neighbor techniques to infer cell communication with adjacent neighbors, whereas COMMOT deduces cell communication within a specific range. We applied a unified prior LR database for all these methods. Although CellChat inferred the majority of the communications, it had a restricted overlap with other methods. In contrast, while SpaCcLink inferred a relatively smaller number of LR pairs, it displayed a larger overlap with other methods, especially with COMMOT (as shown in Fig. 5c. and Additional file 1: Fig. S2). In the human breast cancer dataset, it was noted that SpaCcLink had a lower overlap rate in predicting LR interactions in comparison to SpaTalk. This is because SpaTalk’s predictions mostly overlapped with CellChat, while there was less coincidence with other methods that take spatial information into account. Hence, this implies that SpaCcLink demonstrates higher reliability in inferring LR interactions, as it exhibits a greater overlap with other methods. This indicates that the predictions from SpaCcLink might be more dependable, offering valuable perspectives for understanding cell–cell communication.

Conclusions

Cell–cell communication not only entails the binding of ligands and receptors but also a cascade of downstream signaling reactions initiated by the receptors. Moreover, cell–cell communication typically occurs within a specific range, with a higher probability of communication between adjacent cells. The progress in spatial transcriptomics has empowered us to directly decipher cell–cell communication within the cellular microenvironment. However, among the existing methods for analyzing communication based on spatial transcriptomic data, only a few take into account the downstream signaling reactions within cells. Additionally, the methods that consider intracellular signaling fail to integrate the topological information of prior gene networks and gene expression, thereby restricting the comprehensive analysis of cell communication. Here, we present SpaCcLink, a spatial cell–cell communication approach that combines the downstream signaling reactions within cells and the spatial expression relationships. Case studies on three real datasets (human melanoma, human breast cancer, and mouse brain) illustrate that SpaCcLink can effectively assist us in identifying ligands and receptors with significant downstream impact and spatial expression dependencies, thereby uncovering their downstream intracellular signaling. Moreover, comparisons with other methods suggest that SpaCcLink can identify more reliable ligand-receptor relationships and more efficiently identify pathways associated with biological processes.

SpaCcLink initially focuses on the receptor genes that are more pertinent to downstream reactions and utilizes Moran’s I index for spatial communication analysis to guarantee that the identified receptor genes possess stronger downstream influence and display spatial dependence with ligand expression. Experimental investigations on human melanoma have shown that the downstream signaling patterns identified by SpaCcLink are in accordance with the communication patterns. Thus, analyzing the downstream signaling patterns can aid us in identifying receptor genes that have communication relationships and strong downstream influence. Analysis results on human breast cancer imply that SpaCcLink can help us identify receptor genes that are significantly correlated with patient survival, providing crucial cues for disease research and treatment.

Furthermore, we have constructed a complete signaling pathway from receptors to transcription factors and then to target genes. In the analysis of mouse brain data, SpaCcLink can identify key transcription factors in the signaling cascade and visualize the entire signal transduction process between receptors and target genes. By revealing the complete pathways of intercellular signaling, we can gain insights into the patterns of interactions between cells, signal amplification and cooperation mechanisms, leading to a better comprehension of the functions and regulatory mechanisms of biological systems.

Owing to the consideration of gene networks, SpaCcLink is indeed more appropriate for scenarios with a larger number of genes. In cases with fewer genes, its accuracy and reliability might be constrained. Additionally, the available prior LR databases do not furnish detailed descriptions of the interaction modes between receptors and ligands. Consequently, SpaCcLink adopts an optimal transport-based method to classify LR pairs into short-range interactions and long-range interactions. This step demands a longer computation time. If future research can supply additional information regarding the interaction modes, this step could potentially be simplified or omitted.

Furthermore, with the continuous advancement of single-cell resolution spatial transcriptomic sequencing technologies, we can anticipate exploring the complete process of cell communication at a higher resolution level. This will enable us to achieve a deeper understanding of the interactions and signal transduction between cells, resulting in more accurate analysis outcomes.

Methods

Collection of ligand-receptor pairs and gene signal network data

We have successfully collected ligand-receptor pairs from CellTalkDB [55], CytoTalk, and Connectome [56], including 3324 known ligand-receptor pairs for human and 2484 known ligand-receptor pairs for mouse. TF-TG regulatory data are obtained from TRRUST [57], HTRIdb [58] and RegNetwork [59]. Furthermore, pathways from receptors to TFs are gathered from KEGG [60] and Reactome [61].

To construct the gene signal network, we integrate TF-TG regulatory data and pathways, and only retain the pathways that goes from receptors to TG, forming a strongly connected graph.

Construction of training data

To identify target genes highly correlated with receptors within the genes network, we employ triplet loss [62] for model training. The triplet loss optimizes the model by triplets, i.e., < anchor, positive, negative > . Specifically, it brings the anchor and the positive as close as possible within the manifold space, while pushing the anchor and the negative as far apart as possible. Additionally, it ensures that different domains, which are composed of similar samples, remain as separate as possible. This enables our model to effectively learn the gene interaction relationships and distinguish between interacting and non-interacting gene relationships.

Pearson correlation coefficients for all associations in the gene network are calculated to construct triplets. A default threshold of 5% (which is learnable) is utilized to identify highly correlated links as positive samples, denoted as < anchor, positive > . The gene set not associated with the anchor is taken as negative samples, denoted as < anchor, negative > . Highly related gene pairs and unrelated gene pairs together form triplets, such as \(\text{<}{\text{g}}_{a}\text{,}{\text{g}}_{p}\text{,}{\text{g}}_{n}>\), where \(\text{<}{\text{g}}_{a}\text{,}{\text{g}}_{p}>\in P\), \(\text{<}{\text{g}}_{a}\text{,}{\text{g}}_{n}>\in N\), \({g}_{a}\) refers to anchor, \({g}_{p}\) refers to positive, \({g}_{n}\) refers to negative. In this manner, the trained model is capable of identifying gene relationships with stronger correlations within the collected gene network.

Model architecture

To better identify target genes that are more relevant to the receptor in terms of the gene signaling network structure and expression features, we employ a graph attention network (GAT) model to extract the latent characteristics of gene expression.

We use log-normalized expression data \(X\in {R}^{m\times n}\) and gene interaction network \(A\in {R}^{m\times m}\) as inputs, where n and m denote the number of cells and genes, respectively. The model comprises a two-layer graph attention network. In each layer of GAT, each node is adaptively weighted based on the feature information of its neighboring nodes to obtain a comprehensive embedded representation \(Z\).The embedded representation of the gene \(i\) at layer k (\({Z}_{i}^{(k)}\)) can be defined as follows:

$$Z_i^{(k)}=\sigma\;(\sum\nolimits_{j\in N_i}\alpha_{ij}^kW_kZ_j^{(k-1)})$$

(1)

$${e}_{ij}=Leak\,{Re}\,LU\,(\overrightarrow{a} [W{Z}_{i} \parallel W{Z}_{j}])$$

(2)

$${\alpha }_{ij}=soft{\,{max}}_{j}({e}_{ij})=\frac{\mathit{exp}({e}_{ij})}{{\sum }_{k\in {N}_{i}}\mathit{exp}({e}_{ik})}$$

(3)

where \({W}_{k}\) is the weight parameter matrix of the k layer, \({N}_{i}\) represents the neighbor of gene i in the signal network, \({\alpha }_{ij}^{k}\) is the correlation coefficient between genes i and j, \(\parallel\) is concatenation operation, and \(\overrightarrow{a}\) is the trainable weight vector. Consequently, each gene is capable of integrating information from neighboring genes, thereby facilitating a more accurate capture of the interactions between genes.

Loss function

For a given triplet \(<{\text{g}}_{a}\text{, }{\text{g}}_{p}\text{,}{\text{g}}_{n}>\), the optimization objective of the model adopts triplet loss:

$$L={\sum }_{(a,p,n)\in {\rm K}}\mathit{max}\{d(a,p)-d(a,n)+\alpha ,0\}$$

(4)

where \({\rm K}\) is the triplet list; \(d(a,p)\) is the absolute value of cosine distance between \({z}_{a}\) and \({z}_{p}\); \(\alpha\) is the margin parameter. The margin \(\alpha\) controls the difference in manifold space between interacting pairs and non-interacting pairs. Cosine similarity is used to measure the correlation between gene pairs in the manifold space.

Computation of intracellular influence score

By taking into account the expression of target genes that are highly correlated with the receptor, we assess the downstream influence of the receptor \(r\) in the cell \(i\) using the following definitions:

$${\text{IntraScor}}{\text{e}}_{r}^{i}={\sum }_{Tg\in {T}_{r}}PE{M}_{r}^{i}\cdot PE{M}_{Tg}^{i}\cdot |cor{r}_{r-Tg}|/\left|{\rm T}_{r}\right|$$

(5)

where \({\rm T}_{r}\) is the set of target genes highly correlated with the receptor \(r\); and \(cor{r}_{r-Tg}\) is the absolute value of cosine similarity between the latent feature \({Z}_{r}\) and \({Z}_{Tg}\); \(\left|*\right|\) is the cardinality of the set; \(PE{M}_{j}^{i}\) measures the expression specificity of gene \(j\) in cell \(i\).The definition of \(PE{M}_{j}^{i}\) is as follows:

$${\text{PE}}{\text{M}}_{j}^{i}={\mathit{log}}_{2}\frac{Exp{r}_{j}^{i}}{{e}_{i}^{A}}$$

(6)

$${e}_{i}^{A}={\sum }_{m=1}^{M}Mea{n}_{j}^{m}\cdot \frac{{s}_{*A}}{{\sum }_{m=1}^{M}{s}_{*m}}$$

(7)

where \(Exp{r}_{j}^{i}\) is the expression of gene \(j\) in cell \(i\); \({e}_{i}^{A}\) is the expected expression of gene \(i\) in cell type \(A\) under the null distribution across all \(M\) cell types; \(Mea{n}_{j}^{m}\) is the average expression of gene \(j\) in cell type \(m\); \({s}_{*A}\)is the total expression of all genes in cell type A. Thus, the intracellular score takes into account the expression specificity and association between the receptor and target genes.

Division of the type of ligand-receptor pairs

Cell–cell interaction typically occurs only within a specific range. Based on the interaction distance, cell–cell interaction can be categorized into four types: autocrine, paracrine, juxtacrine, and endocrine. However, it is challenging to assess endocrinology in local tissue, and the interaction range between autocrine and juxtacrine is relatively short and difficult to distinguish. Moreover, the physical distance between cells impacts signal diffusion and intensity, and different interaction types need to be evaluated at appropriate spatial scales. The curated LR data that we have collected does not provide information about the types of interactions. Inspired by previous research [63], we employ the Wasserstein distance to partition the types of interactions for the ligand-receptor pairs. The Wasserstein distance is a metric for quantifying the difference between two probability distributions. It is founded on optimal transport theory and computes the minimum cost between two distributions by contemplating how to transform one distribution into another.

In the classification of ligand-receptor interaction types, the Wasserstein distance is utilized to measure the minimum transportation cost between the distributions of ligand expression and receptor expression. By supposing that the ligand genes are expressed in m cells and the receptor genes are expressed in n cells, and using the Euclidean distance matrix between cells as the cost matrix \(D\in {\mathbb{R}}^{m\times n}\), the Wasserstein distance from the ligand expression distribution to the receptor expression distribution can be calculated by resolving the optimal transportation cost problem as follows:

$${\gamma }^{\prime}=\underset{\gamma \in \Gamma (\text{L,}R)}{argmin}<\gamma ,D>$$

(8)

And the corresponding Wasserstein distance is obtained as follows.

$$W(\text{L,R})=\underset{\gamma \in \Gamma (L,R)}{min}<\gamma ,D>=<{\gamma }^{\prime},D>$$

(9)

We then conduct multiple random permutations (with a default of 200 permutations) on the expression distributions of the ligand and receptor. We recalculate the Wasserstein distance and utilize the d_ratio (the ratio of the actual distance to the mean of the multiple perturbations) and the two-sided significance p-value to classify the types of interactions.

After sorting lr pairs based on d_ratio and p_value, we select the lr pairs that satisfy the condition d_ratio < 1 and p_value < 0.05 from the top 10% as short-range interaction pairs. We then select the lr pairs that satisfy the condition d_ratio > 1 and p_value < 0.05 from the bottom 10% as long-range interaction pairs, and the remaining LR pairs are categorized as medium-range interaction pairs.

Spatial communication analysis

To identify ligand-receptor pairs with spatial dependence, we calculate the cross-Moran’s I index for all interaction pairs and obtain significant candidate pairs through a permutation test (with a default of 500 permutations). In the permutation test, we randomly shuffle the locations of all spots and assess statistical significance by comparing the observed cross-Moran’s I index to the likelihood of obtaining such an extremely positive value under the permutation-based random model. Furthermore, we consider that the communication relationship with a p-value less than 0.05 is significant.

$$\text{Cross-Moran}^{\prime}\,\text{I} =\frac{N}{{\sum }_{i}^{N}{\sum }_{j}^{N}{W}_{ij}}\frac{{\sum }_{i}^{N}{\sum }_{j}^{N}{W}_{ij}({x}_{i}-\overline{x})({y}_{j}-\overline{y})}{\sqrt{{\sum }_{i}^{N}({x}_{i}-\overline{x}{)}^{2}}\sqrt{{\sum }_{j}^{N}({y}_{j}-\overline{y}{)}^{2}}}$$

(10)

where \({x}_{i}\) represents the expression of ligand \(x\) in cell \(i\), while \({y}_{j}\) represents the expression of receptor gene \(y\) in cell \(j\); N is the number of cells; \({W}_{ij}\) is the distance weight of cells i and j, which is calculated based on the k-nearest neighbor method. The exact choice of k should reflect the data. For all data sets analyzed in this study, which are derived from the 10X Visium platform, according to Squidpy [64], we set k = 6 for short-range interactions and multiplied on the medium-range and long-range upper parameters to capture different ranges of interactions. For other data sets, users can visualize the weight matrix for different parameters and choose the most appropriate scale for their data.

We further use local cross-Moran’s I index to identify important local interaction regions for each LR pairs. The local cross-Moran’s I for spot i is defined as:

$$\text{Local Cross-Moran}^\prime\,{\text{I}}_{i}\text{ = }\frac{{\sum }_{j}^{N}{W}_{ij}[({x}_{i}-\overline{x})({y}_{j}-\overline{y})+({y}_{i}-\overline{y})({x}_{j}-\overline{x})]}{\sqrt{{\sum }_{i}^{N}({x}_{i}-\overline{x}{)}^{2}}\sqrt{{\sum }_{j}^{N}({y}_{j}-\overline{y}{)}^{2}}}$$

(11)

Similarly, we utilize permutation testing for conducting a significance test (with a default of 500 permutations). After identifying reliable communication relationships, we further quantify the cell–cell communication strength between cell \(i\) and cell \(j\) for LR pair \(l-r\) using the following definition:

$$\text{CCI }{\text{Scor}}{\text{e}}_{lr}^{ij}={x}_{l}^{i}\cdot {x}_{r}^{j}\cdot {\text{IntraScor}}{\text{e}}_{r}^{j}\cdot (1-\frac{{d}_{ij}}{{d}_{thre}})$$

(12)

where \({x}_{l}^{i}\) is the expression of gene \(l\) in cell \(i\); \({d}_{ij}\)is the distance between cell \(i\) and \(j\); \({d}_{thre}\) is the median distance among all local interaction hotspots of the interaction pair.

Analysis of communication between cell types

To identify ligand-receptor pairs with significant spatial communication between clusters, we calculate the cross-Moran’s I index for the interaction pair (ligand \(x\) and receptor \(y\)) between cell type A and cell type B as follows:

$${\text{Celltype}}{\text{s}} \, \text{Cross-Moran}^\prime\,\text{I}=\frac{N+M}{2{\sum }_{i\in A}{\sum }_{j\in B}{W}_{ij}}\frac{{\sum }_{i\in A}{\sum }_{j\in B}{W}_{ij}({x}_{i}-\overline{x})({y}_{j}-\overline{y})}{\sqrt{{\sum }_{i\in A}({x}_{i}-\overline{x}{)}^{2}}\sqrt{{\sum }_{j\in B}({y}_{j}-\overline{y}{)}^{2}}}$$

(13)

where N and M are the number of cells belonging to cell types A and B, respectively. And we also permute the locations of spots within each cluster for significance test (with a default of 500 permutations).

After filtering out LR pairs with significant communication relationships between different cell types through cross-Moran index, we can further calculate cell–cell communication scores between different cell types for the LR pair \(l-r\).The cell-type-level signaling strength is quantified as.

$${\text{Celltype}} \, \text{CCI }{\text{Scor}}{\text{e}}_{\text{lr}}^{\text{AB}}=\frac{{\sum }_{(k,l)\in {I}_{AB}}\text{CCI Scor}{\text{e}}_{\text{lr}}^{\text{kl}}}{\left|{I}_{AB}\right|}$$

(14)

where \({I}_{AB}=\{(k,l)|k\in {S}_{A},l\in {S}_{B}\}\), \({S}_{A}\) is the cluster A.

Identification of downstream key transcription factors and signaling pathways

After identifying the reliable ligand-receptor relationship, we further identify the signals affected by the ligand-receptor. First, we utilize the inter-gene co-expression coefficient to further screen the receptor-target gene relationship for each cell type (with a threshold of 0.1) and then employ the Fisher’s test to identify the activity fraction of each transcription factor. The activity score of TF k in the cell type \(t\) is calculated by Fisher’s exact test as follows:

$$P=\frac{\left(\begin{array}{c}\text{a+b}\\ a\end{array}\right)\left(\begin{array}{c}\text{c+d}\\ c\end{array}\right)}{\left(\begin{array}{c}\text{a+b+c+d}\\ a+c\end{array}\right)}$$

(15)

where \(a\) is the number of the set of target genes that are highly associated with the receptor \(r\) and TF \(k\); \(b\) is the number of the set of target genes that receptor r can act on through the prior gene network excluding \(a\); \(c\) is the number of the set of target genes that is highly correlated with receptor \(r\) but not associated with TF \(k\); \(d\) is the number of all genes excluding a, b, and c.

After identifying the key transcription factors, we search for the optimal signal pathway from the receptor to the target gene based on the co-expression coefficient between the genes.

Comparison with other methods

In our experiments, we conducted a comparison of our method against CellChat, COMMOT, SpaTalk, and Giotto using three real datasets. For each of the methods, our collected data was utilized as the prior ligand-receptor (LR) database, and default parameters were employed to predict LR pairs. We regarded the interaction pairs between clusters with a significant communication p-value of less than 0.05 as being of importance.

In addition to the aforementioned methods, we also made a comparison with methods (SpaTalk, CellCall, NicheNet) that take downstream signaling into account. Given the complexity of their prior databases, we opted to use their original databases. All methods were implemented with default parameters. Inspired by SpaTalk, we employed Fisher’s exact test to compare the significance of predicted downstream signaling pathways involved in KEGG and Reactome among different methods.

Data availability

The code of the model and example dataset can be downloaded from GitHub (https://github.com/LiangYu-Xidian/SpaCcLink, and Zenodo DOI: https://doi.org/https://doi.org/10.5281/zenodo.14504149 [65]). All data analyzed during this study are included in previously published datasets. All of these datasets have been annotated with cell types.

Human melanoma dataset

Li Z, Wang T, Liu P, Huang Y. SpatialDM for rapid identification of spatially co-expressed ligand–receptor and revealing cell–cell communication patterns. https://ndownloader.figshare.com/files/40178320. (2023).

Human breast cancer dataset

Li H, Ma T, Hao M, Guo W, Gu J, Zhang X, et al. Decoding functional cell–cell communication events by multi-view graph learning on spatial transcriptomics. https://cloud.tsinghua.edu.cn/f/dd941f0d12214e6583bb/?dl=1. (2023).

Mouse brain dataset

Palla G, Spitzer H, Klein M, Fischer D, Schaar AC, Kuemmerle LB, et al. Squidpy: a scalable framework for spatial omics analysis. https://ndownloader.figshare.com/files/26098397. (2022).

Abbreviations

scRNA-seq:

Single-cell transcriptomic sequencing

L:

Ligand

R:

Receptor

LR:

Ligand-receptor

NMF:

Non-negative matrix factorization

M:

Intermediate

TF:

Transcription factor

TG:

Target gene

CAF:

Cancer-associated fibroblasts

GOBP:

Gene Ontology Biological Process

GAT:

Graph attention network

CCK:

Cholecystokinin

References

Brücher BL, Jamall IS. Cell-cell communication in the tumor microenvironment, carcinogenesis, and anticancer treatment. Cell Physiol Biochem. 2014;34(2):213–43.

Kozar K, Ciemerych MA, Rebel VI, Shigematsu H, Zagozdzon A, Sicinska E, et al. Mouse development and cell proliferation in the absence of D-cyclins. Cell. 2004;118(4):477–91.

Armingol E, Officer A, Harismendy O, Lewis NE. Deciphering cell–cell interactions and communication from gene expression. Nat Rev Genet. 2021;22(2):71–88.

Yang BA, Westerhof TM, Sabin K, Merajver SD, Aguilar CA. Engineered tools to study intercellular communication. Adv Sci. 2021;8(3):2002825.

Almet AA, Cang Z, Jin S, Nie Q. The landscape of cell–cell communication through single-cell transcriptomics. Curr Opin Syst Biol. 2021;26:12–23.

Combarnous Y, Nguyen TMD. Cell communications among microorganisms, plants, and animals: origin, evolution, and interplays. Int J Mol Sci. 2020;21(21):8052.

Weidemüller P, Kholmatov M, Petsalaki E, Zaugg JB. Transcription factors: bridge between cell signaling and gene regulation. Proteomics. 2021;21(23–24):2000034.

Dai C, Jiang Y, Yin C, Su R, Zeng X, Zou Q, et al. scIMC: a platform for benchmarking comparison and visualization analysis of scRNA-seq data imputation methods. Nucleic Acids Res. 2022;50(9):4877–99.

Liu M, Chen W, Zhao J, Zheng C, Guo F. scTSSR-D: gene expression recovery by two-side self-representation and dropout information for scRNA-seq data. Curr Bioinform. 2023;18(4):285–95.

Wang Z, Ding H, Zou Q. Identifying cell types to interpret scRNA-seq data: how, why and more possibilities. Brief Funct Genomics. 2020;19(4):286–91.

Zhao M, He W, Tang J, Zou Q, Guo F. A hybrid deep learning framework for gene regulatory network inference from single-cell transcriptomic data. Brief Bioinform. 2022;23(2):bbab568.

Xu J, Xu J, Meng Y, Lu C, Cai L, Zeng X, et al. Graph embedding and Gaussian mixture variational autoencoder network for end-to-end analysis of single-cell RNA sequencing data. Cell Rep Methods. 2023;3(1):100382.

Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. Cell PhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes. Nat Protoc. 2020;15(4):1484–506.

Hou R, Denisenko E, Ong HT, Ramilowski JA, Forrest AR. Predicting cell-to-cell communication networks using NATMI. Nat Commun. 2020;11(1):5011.

Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan C-H, et al. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12(1):1088.

Noël F, Massenet-Regad L, Carmi-Levy I, Cappuccio A, Grandclaudon M, Trichot C, et al. Dissection of intercellular communication using the transcriptome-based framework ICELLNET. Nat Commun. 2021;12(1):1089.

Browaeys R, Saelens W, Saeys Y. NicheNet: modeling intercellular communication by linking ligands to target genes. Nat Methods. 2020;17(2):159–62.

Cheng J, Zhang J, Wu Z, Sun X. Inferring microenvironmental regulation of gene expression from single-cell RNA sequencing data using scMLnet with an application to COVID-19. Brief Bioinform. 2021;22(2):988–1005.

Hu Y, Peng T, Gao L, Tan K. CytoTalk: De novo construction of signal transduction networks using single-cell transcriptomic data. Sci Adv. 2021;7(16):eabf1356.

Zhang Y, Liu T, Hu X, Wang M, Wang J, Zou B, et al. Cell Call: integrating paired ligand–receptor and transcription factor activities for cell–cell communication. Nucleic Acids Res. 2021;49(15):8520–34.

Longo SK, Guo MG, Ji AL, Khavari PA. Integrating single-cell and spatial transcriptomics to elucidate intercellular tissue dynamics. Nat Rev Genet. 2021;22(10):627–44.

Cheng M, Jiang Y, Xu J, Mentis A-FA, Wang S, Zheng H, et al. Spatially resolved transcriptomics: a comprehensive review of their technological advances, applications, and challenges. J Genet Genomics. 2023;50(9):625–40.

Wang R, Jiang Y, Jin J, Yin C, Yu H, Wang F, et al. DeepBIO: an automated and interpretable deep-learning platform for high-throughput biological sequence prediction, functional annotation and visualization analysis. Nucleic Acids Res. 2023;51(7):3017–29.

Shao X, Li C, Yang H, Lu X, Liao J, Qian J, et al. Knowledge-graph-based cell-cell communication inference for spatially resolved transcriptomic data with SpaTalk. Nat Commun. 2022;13(1):4429.

Dries R, Zhu Q, Dong R, Eng C-HL, Li H, Liu K, et al. Giotto: a toolbox for integrative analysis and visualization of spatial expression data. Genome Biol. 2021;22:1–31.

Raredon MSB, Yang J, Kothapalli N, Lewis W, Kaminski N, Niklason LE, et al. Comprehensive visualization of cell–cell interactions in single-cell and spatial transcriptomics with NICHES. Bioinformatics. 2023;39(1):btac775.

Tang Z, Zhang T, Yang B, Su J, Song Q. spaCI: deciphering spatial cellular communications through adaptive graph model. Brief Bioinform. 2023;24(1):bbac563.

Yuan Y, Bar-Joseph Z. GCNG: graph convolutional networks for inferring gene interaction from spatial transcriptomics data. Genome Biol. 2020;21(1):1–16.

Li H, Ma T, Hao M, Guo W, Gu J, Zhang X, et al. Decoding functional cell–cell communication events by multi-view graph learning on spatial transcriptomics. Brief Bioinform. 2023;24(6):bbad359.

Cang Z, Zhao Y, Almet AA, Stabell A, Ramos R, Plikus MV, et al. Screening cell–cell communication in spatial transcriptomics via collective optimal transport. Nat Methods. 2023;20(2):218–28.

Li H, Pang Y, Liu B. BioSeq-BLM: a platform for analyzing DNA, RNA, and protein sequences based on biological language models. Nucleic Acids Res. 2021;49(22):e129.

Li H, Liu B. BioSeq-Diabolo: Biological sequence similarity analysis using Diabolo. PLoS Comput Biol. 2023;19(6):e1011214.

Wang J, Chen Y, Zou Q. Inferring gene regulatory network from single-cell transcriptomes with graph autoencoder model. PLoS Genet. 2023;19(9):1010942.

Hofheinz K, Kakularam KR, Adel S, Anton M, Polymarasetty A, Reddanna P, et al. Conversion of pro-inflammatory murine Alox5 into an anti-inflammatory 15S-lipoxygenating enzyme by multiple mutations of sequence determinants. Arch Biochem Biophys. 2013;530(1):40–7.

Zeng N, Ma L, Cheng Y, Xia Q, Li Y, Chen Y, et al. Construction of a ferroptosis-related gene signature for predicting survival and immune microenvironment in melanoma patients. Int J Gen Med. 2021;14:6423–38.

Rao Y, Zhu J, Zheng H, Dong W, Lin Q. A novel melanoma prognostic model based on the ferroptosis-related long non-coding RNA. Front Oncol. 2022;12:929960.

Xu C, Chen H. A ferroptosis-related gene model predicts prognosis and immune microenvironment for cutaneous melanoma. Front Genet. 2021;12:697043.

Luger D, Yang Y-A, Raviv A, Weinberg D, Banerjee S, Lee M-J, et al. Expression of the B-cell receptor component CD79a on immature myeloid cells contributes to their tumor promoting effects. PLoS ONE. 2013;8(10):e76115.

Liu Y, Chen Y, Hu X, Meng J, Li X. Development and validation of the B cell-associated Fc receptor-like molecule-based prognostic signature in skin cutaneous melanoma. Biomed Res Int. 2020;2020(1):8509805.

Nagy N, Busalt F, Halasy V, Kohn M, Schmieder S, Fejszak N, et al. In and out of the bursa—the role of CXCR4 in chicken B cell development. Front Immunol. 2020;11:1468.

Werner U, Künstner A, Drenckhan M, Pries R, Bruchhage K-L, Busch HS, et al. Linking complement C3 and B cells in nasal polyposis. J Immunol Res. 2020;2020(1):4832189.

Mitchell B, Mahalingam M. The CXCR4/CXCL12 axis in cutaneous malignancies with an emphasis on melanoma. 2014.

Pan Y, Jin K, Xie X, Wang K, Zhang H. MicroRNA-19a-3p inhibits the cellular proliferation and invasion of non-small cell lung cancer by downregulating UBAP2L Retraction in/10.3892/etm. 2022.11557. Exp Ther Med. 2020;20(3):2252–61.

Mishima S, Sakamoto M, Kioka H, Nagata Y, Suzuki R. Multifunctional regulation of VAMP3 in exocytic and endocytic pathways of RBL-2H3 cells. Front Immunol. 2022;13:885868.

Röhl J, West ZE, Rudolph M, Zaharia A, Van Lonkhuyzen D, Hickey DK, et al. Invasion by activated macrophages requires delivery of nascent membrane-type-1 matrix metalloproteinase through late endosomes/lysosomes to the cell surface. Traffic. 2019;20(9):661–73.

Veale KJ, Offenhäuser C, Lei N, Stanley AC, Stow JL, Murray RZ. VAMP3 regulates podosome organisation in macrophages and together with Stx4/SNAP23 mediates adhesion, cell spreading and persistent migration. Exp Cell Res. 2011;317(13):1817–29.

Huang XY, Liu JJ, Liu X, Wang YH, Xiang W. Bioinformatics analysis of the prognosis and biological significance of VCAN in gastric cancer. Immun Inflamm Dis. 2021;9(2):547–59.

Ren L, Huang D, Liu H, Ning L, Cai P, Yu X, et al. Applications of single-cell omics and spatial transcriptomics technologies in gastric cancer (Review). Oncol Lett. 2024;27(4):152.

Nishimura S, Bilgüvar K, Ishigame K, Sestan N, Günel M, Louvi A. Functional synergy between cholecystokinin receptors CCKAR and CCKBR in mammalian brain development. PLoS ONE. 2015;10(4):e0124295.

Rehfeld JF, Friis-Hansen L, Goetze JP, Hansen TV. The biology of cholecystokinin and gastrin peptides. Curr Top Med Chem. 2007;7(12):1154–65.

Keskin-Aktan A, Akbulut KG, Abdi S, Akbulut H. SIRT2 and FOXO3a expressions in the cerebral cortex and hippocampus of young and aged male rats: antioxidant and anti-apoptotic effects of melatonin. Biol Futur. 2022;73(1):71–85.

Baburamani AA, Sobotka KS, Vontell R, Mallard C, Supramaniam VG, Thornton C, et al. Effect of Trp53 gene deficiency on brain injury after neonatal hypoxia-ischemia. Oncotarget. 2017;8(7):12081.

Long J, Tian L, Baranova A, Cao H, Yao Y, Rao S, et al. Convergent lines of evidence supporting involvement of NFKB1 in schizophrenia. Psychiatry Res. 2022;312:114588.

Viola TW, Creutzberg KC, Zaparte A, Kestering-Ferreira É, Tractenberg SG, Centeno-Silva A, et al. Acute neuroinflammation elicited by TLR-3 systemic activation combined with early life stress induces working memory impairments in male adolescent mice. Behav Brain Res. 2019;376:112221.

Shao X, Liao J, Li C, Lu X, Cheng J, Fan X. Cell TalkDB: a manually curated database of ligand–receptor interactions in humans and mice. Brief Bioinform. 2021;22(4):bbaa269.

Ramilowski JA, Goldberg T, Harshbarger J, Kloppmann E, Lizio M, Satagopam VP, et al. A draft network of ligand–receptor-mediated multicellular signalling in human. Nat Commun. 2015;6(1):7866.

Han H, Cho J-W, Lee S, Yun A, Kim H, Bae D, et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res. 2018;46(D1):D380–6.

Bovolenta L, Acencio M, Lemke N. HTRIdb: an open-access database for experimentally verified human transcriptional regulation interactions. Nature Precedings. 2012;13:1–10.

Liu Z-P, Wu C, Miao H, Wu H. RegNetwork: an integrated database of transcriptional and post-transcriptional regulatory networks in human and mouse. Database. 2015;2015:bav095.

Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

Milacic M, Beavers D, Conley P, Gong C, Gillespie M, Griss J, et al. The Reactome Pathway Knowledgebase 2024. Nucleic Acids Res. 2024;52(D1):D672–8.

Schroff F, Kalenichenko D, Philbin J, editors. Facenet: a unified embedding for face recognition and clustering. Proceedings of the IEEE conference on computer vision and pattern recognition; 2015.

Liu Z, Sun D, Wang C. Evaluation of cell-cell interaction methods by integrating single-cell RNA sequencing data with spatial information. Genome Biol. 2022;23(1):1–38.

Palla G, Spitzer H, Klein M, Fischer D, Schaar AC, Kuemmerle LB, et al. Squidpy: a scalable framework for spatial omics analysis. Nat Methods. 2022;19(2):171–8.

Liu J, Yu L. SpaCcLink: Exploring downstream signaling regulations with graph attention network for systematic inference of spatial cell-cell communication. Zenodo (https://doi.org/10.5281/zenodo.14504149). 2024.