Computational exploration of cellular communication in skin from emerging single-cell and spatial transcriptomic data

Communication between cells is often mediated by various types of soluble and membrane-bound factors, such as ligands, receptors, extracellular matrix (ECM), integrins and junction proteins [1]. Cell–cell communication is critical for cell fate decisions and tissue organization in multicellular organisms [2–4]. Traditionally, cell–cell communication studies have relied on a large number of traditional experiments, including histological section analysis of tissues, in vitro co-cultures, and in vivo genetic manipulations [4]. However, these experiments are often limited to investigate communications among a very small number of cell types. Recently, single-cell RNA sequencing (scRNA-seq) has emerged as a revolutionary method to reveal cellular heterogeneity in tissues with unprecedented resolution and scale. It is increasingly clear that vast amounts of scRNA-seq data collected and published to date inherently contain gene expression information on signaling crosstalk between cells. This offers an unprecedented opportunity to analyze the cell–cell communication events driving diverse cell fate decisions across different tissues, and can leapfrog the field over the barrier typically associated with using traditional experimental tools.

A number of computational tools have been developed to infer cell–cell communication by integrating scRNA-seq data and prior knowledge of the interactions between signaling ligands, receptors and their cofactors [2,5]. However, scRNA-seq does not capture the spatial distribution and local environment within tissue [3]. Since signals often have a limited range between neighboring cells due to finite spatial diffusivity of the ligands, spatial information is essential for understanding intercellular communication and tissue function. Recently developed imaging techniques such as spatial transcriptomics, which can physically localize gene expression to specific cell subsets in tissues, have increased our understanding of the roles of cell–cell communication in tissue [6]. Due to the low resolution and coverage of spatial transcriptomics, several methods have been developed to characterize cell–cell communication by integrating scRNA-seq data with spatial information from imaging methods [3].

Because of its complex cellular composition and the interplay among its different cell types, the skin is a highly suitable and accessible model to study signaling mechanisms common to many other tissues [7–9]. In this review, after introducing the crucial role of cell–cell communication in cell fate decisions, we first summarize the computational approaches for cell–cell communication inference and analysis using scRNA-seq data and emerging approaches for using spatial transcriptomics data. Then, we describe the biological insights that can be gained by applying cell–cell communication analysis to the development, homeostasis, disease and aging of skin tissue. Finally, we discuss future directions in the field.

The skin consists of two distinct layers, including the upper epidermis and the dermis below it. Both layers host different types of cells, including keratinocytes, melanocytes and Langerhans cells in the epidermis, while the dermis contains fibroblasts, endothelial cells and various immune cells such as T cell, B cell, macrophages and dendritic cells. Interactions between epidermal and dermal cell populations influence cell fate decision in the skin [8]. For example, both epidermal Wnt signaling and dermal extracellular matrix cause fibroblast proliferation in adult mice [10]. During early hair follicle morphogenesis and development, WNT, EDA and FGF play vital roles in placode and dermal condensate fate commitment [11]. Dysregulation of these pathways leads to abnormal development. Previous studies showed that ablation of Fgf20 in mice results in the failure of dermal condensate formation [12]. Although many putative ligands and receptors involved in epidermis-dermis interactions have been identified using histology [8,9,13], our systematic understanding of the complex mechanisms driving cell fate decisions remains obscure. The recent advance of single-cell RNA sequencing technology can systematically assess gene expression, providing an unprecedented opportunity to understand complex cell–cell communications systematically and comprehensively with the help of emerging computational tools for cell–cell communication inference and analysis [2].

To facilitate the cell–cell communication exploration and analysis, an increasing number of computational tools have been developed to systematically infer cell–cell communication [5,14–22] (Figure 1). Computational inference relies on prior knowledge of ligand–receptor interactions based on relevant literature and public databases. Currently, over 15 ligand–receptor databases have been built [2] (Figure 1). However, curation of such information is laborious and requires careful examination of the role of different signaling molecules in each context, as well as consideration of interactions from multiple other sources in the system. These databases are often built from different resources, leading to incomplete and inherently biased information of ligand–receptor pairs. Indeed, although these databases share the same original resources, including KEGG, STRING, Reactome and Guide to Pharmacology, different databases often include distinct proportions of the secreted and cell–cell contact interactions, and may under-represent certain signaling pathways and categories [14]. For example, innate immune pathways and Hedgehog are under-represented in ICELLNET [21] and CellPhoneDB [14]. Compared with other databases that use only one ligand/one receptor gene pairs, CellChatDB [15], CellPhoneDB and ICELLNET take into account complexes with multimeric ligands and receptors, leading to an accurate representation of known heteromeric molecular complexes. For example, signaling via BMP, IL, Interferon, and TGFβ pathways requires the presence of more than one membrane-bound receptor subunits [23]. A lack of expression of any subunit blocks ligand–receptor interactions and the resulting downstream communication cascade. Recently, the integrated database OmniPath was built by combining all interactions from the above databases as well as additional resources [24], leading to a more comprehensive resource for cell–cell communication analysis. Although such a unified resource is particular useful for benchmarking various methods, it also raises the possibility of false-positive data due to experimentally unvalidated ligand–receptor interactions as well as false-negative data due to cross-talk between unvalidated ligands not being accounted for. In addition, most of the databases are curated for human and mouse systems to the detriment of other established and emerging models. Future efforts should be made to address this, in particular for well-established vertebrate models like zebrafish, a paradigm in developmental biology. Since inference methods are sensitive to the quality of the databases, a comprehensive, high-quality resource is urgently required.

In addition to the varying databases of ligand–receptor interactions, different methods employed different computational strategies for cell–cell communication inference [2,5,14–16,25,26]. The core assumption of these methods is that the chance of a sender cell (i.e. a cell as signal source) communicating with a receiver cell (i.e. a cell as signal target) is positively correlated with the expression levels of ligand and cognate receptor in the sender and receiver cell, respectively. Most of the methods infer cell–cell communication between interacting cell clusters (i.e. a group of cells) instead of individual cells [5]. Individual cell-based methods such as SoptSC [27] are useful for studying cell–cell communication at single-cell resolution, but suffer from noise and high data dimensionality and computational costs due to the large number of cells and ligand–receptor pairs involved in computation. For methods that predict cluster–cluster communication, both CellChat [15] and CellPhoneDB [16] infer statistically significant interactions by randomly permuting the cluster labels of cells. More importantly, compared with other methods that only consider the expression of one ligand/one receptor per gene pair, CellChat, CellPhoneDB and ICELLNET [28] consider the expression of the members of the heteromeric complex, which highlights the importance of subunits in cell signaling. For instance, soluble ligands from the TGFβ pathway signal via heteromeric complexes of type I and type II receptors [29]. Tgfbr1 or Tgfbr2 knockout mice exhibited impaired phenotypes in diverse biological processes, such as head and neck carcinogenesis [30] and female reproductive tract [31]. The ligand–receptor interaction also depends on soluble and membrane-bound cofactors [15]. For example, WNT-related cofactors positively and negatively modulate WNT signaling. More recently, a new modeling framework was presented in CellChat to integrate all these known interactions.

For a given signaling pathway, a receiver cell responds to signals from a sender cell by triggering downstream gene responses, including altered transcription factor activity and target gene expression. While most of the existing methods infer intercellular communication based only on expression of ligand–receptor pairs, several methods have recently taken into account the intracellular downstream response in receiver cells [27,32–37] (Figure 1). These methods could reduce false-positive communications by considering the expression of genes for each receptor, transcription factor (TF) and target genes involved in a particular pathway. However, an assumption of these methods is that the TF activity can be approximated by the gene expression level. Because TF activation has multiple molecular mechanisms, it is essential that TF activity can be measured directly or estimated based on the expression levels of its target genes. Thus, these methods potentially lead to additional false positives. Including additional information from emerging technologies such as INs-seq [38] which records scRNA-seq and intracellular protein activity, can improve inference.

Moreover, identification of signaling changes across different conditions is important for understanding how distinct cell states respond to evolution, perturbations and diseases. Most of these methods only focus on the intercellular communications in one biological condition, while several methods [15,22,39,40] have also been developed recently to perform comparison analysis of cell–cell communication across conditions (Figure 1). Different from iTALK [22] and Connectome [39], which identify altered signaling based on differential expression analysis, CellChat [15] performs comparison analysis of both cell–cell communication structure and strength by performing joint manifold learning and quantitative contrasts. Differential expression-based methods have advantages in detecting context-specific signaling, but likely fail to identify shared interactions across distinct contexts. More recently, different from scTensor for a single condition [41], a sophisticated approach called Tensor-cell2cell has been presented to decipher complex cell–cell communication patterns across diverse conditions [42]. This method is attractive because it deciphers context-driven intercellular communication by simultaneously accounting for multiple conditions by utilizing a tensor decomposition framework. However, one limitation of this method is the requirement of the same cell types appearing across all experimental conditions.

Rather than the interaction score-based methods, more advanced approaches are highly needed to infer and compare cell–cell communication networks. Beyond inference, methods such as CellChat also include powerful visualization features, which greatly facilitate the analysis and interpretation of complex intercellular communication networks (Figure 1). Currently, few efforts have been made to identify general rules for cell communication and how multicellular organisms exploit these rules to coordinate tissue morphogenesis and function. While CellChat addresses this critical point by using a matrix factorization method, further improvements, like including niche neighbor information, will greatly improve accuracy.

Due to the finite spatial diffusivity of the soluble ligands in paracrine signaling and the physical contact between adjacent cells in juxtacrine signaling [3], spatial information is vital to study cell–cell communication. This spatial information of biological tissue is lost in scRNA-seq, but preserved in imaging-based technologies. The recent advances of spatial transcriptomics such as Visium, STARmap, MERFISH and seqFISH+ [6], greatly help to understand complex cell–cell communication by offering both spatial information and molecular profiles. Given that spatial transcriptomics techniques often do not yet provide gene expression profiles at single-cell resolution, integration of scRNA-seq and spatial transcriptomics data will be helpful in understanding cellular composition and communication within and across tissues. Incorporating spatial information will likely reduce false-positive inferred signaling links, because cells only communicate directly over a limited spatial distance. Computational methods that do not use spatial information may fail to detect certain expected ligand–receptor interactions [43]. Here, we discuss some emerging approaches that have been developed to infer cell–cell communication using spatial information [44].

The key for spatial-informed approaches lies in the identification of spatially proximal cell clusters [45], which are used to filter out spatially distant cell–cell communication events. Determination of spatial proximity starts by representing the spatial information through a spatial network. Giotto randomly permutes cell cluster labels [46], while Squidpy computes a co-occurrence score for clusters [47]. However, both methods solely utilize spatial coordinates while neglecting gene expression information. To address this limitation, Stlearn learns a joint representation by integrating both gene expression and image information. In contrast with other methods, Stlearn also identifies spatial regions where interactions between cell types are most likely to occur by grouping the spots with the most similar ligand–receptor co-expression values and calculating cell type diversity across the tissue [48]. The identification of spatial regions with strongly active cell–cell communication can potentially shed light into the signaling mechanisms responsible for certain biological processes. Recently, several machine learning-based approaches have been proposed [49–52], which provide additional insights of microenvironmental mechanisms from different spatial views. For example, NCEM reconciles variance attribution and communication modeling in a single model of tissue niches based on graph neural networks [52]. Several exciting new directions still not covered by existing approaches include the spatially dependent interplay of signaling networks and gene regulatory networks, the correlation of cell–cell communication inferred from different molecular modalities in spatial multiomics, and the incorporation of 3D neighborhood information of cells. Further development is also needed to infer cell–cell communication to fully characterize cell function and spatial cellular organization across a whole tissue [43]. While the rapid development of spatial transcriptomics has provided new opportunities, it also brings additional exciting challenges. One outstanding challenge of spatial transcriptomics data analysis is the inference of cellular compositions within each spot, which does not have single-cell resolution and often captures more than one cell type [6]. Accurate integration of spatial transcriptomics and scRNA-seq data is important for such deconvolution and to further improve the discovery of cell communication events between cells. Moreover, due to the complicated spatial distribution of cell types in tissues, the determination of the spatial proximity cell types remains challenges. Future work in this direction will help to improve the inference of spatial-informed cell–cell communication in tissues.

Cell–cell communication tools have been successfully applied to a diverse range of biological systems to dissect mechanisms of cell fate decisions during tissue homeostasis, development and disease [2,3]. Here, we discuss how cell–cell communication contributes to skin development, wound healing, disease and aging. Examples of studies were summarized in Table 1.

Skin development requires coordination between different cell types in epidermis and dermis [8]. Cell–cell communication analysis using scRNA-seq data has increased our understanding of the roles of novel signaling pathways and cell subpopulations in skin development. For example, CellChat predicts a novel role of Edn3 signaling in stimulating directed migration of melanocytes into placodes during mouse hair follicle formation [15]. Application of SoptSC to scRNA-seq data of human neonatal epidermis showed distinct cell–cell communication patterns for heterogeneous basal stem cell subpopulations and basal cell populations as crucial signaling hubs to maintain epidermal communication [53]. CellPhoneDB analysis of scRNA-seq data from fetal human skin development predicted an active role of double positive αβγδ T cells during fetal skin immune responses via interactions with other immune cells, including myeloid cells and natural killer cells, and interactions with fibroblasts and endothelial cells [54].

Extensive cell–cell communication occurs between the diverse compartments of the skin in response to wound healing after skin injury [8,65,66]. Understanding of signaling mechanisms underpinning tissue regeneration and repair is critical for wound-healing therapeutics and translation[66]. CellChat analysis of scRNA-seq data from day 12 mouse skin wound healing showed complicated TGFβ signaling with redundant ligand sources targeting fibroblasts. On the contrast, Wnt5a was predicted to mediate signaling from fibroblasts to fibroblasts, endothelial and myeloid, which exhibits a non-redundant network architecture. CellPhoneDB analysis of mouse scRNA-seq data predicted ephrin signaling-mediated epithelial-mesenchymal crosstalk that enables mesenchymal competence for regeneration in skin wound healing[55]. Another application of CellPhoneDB to skin injury revealed that αvβ8 expression in skin Tregs helps to activate TGF-β in neighboring keratinocytes, which act directly on epithelial cells to promote CXCL5 production and neutrophil recruitment, suggesting that αvβ8-expressing Tregs contributes to innate inflammation and delayed epidermal repair after inducing skin injury [56]. In addition, NicheNet analysis was applied to identify the signals driving the response of OTULIN-deficient keratinocytes to inflammation [57]. NicheNet [32] predicted that cytokine IL-1β and chemokine MCP-1 production by infiltrating immune cells contributes to the inflammatory skin phenotype in OTULIN-deficient mice [57].

Acne, atopic dermatitis (AD), psoriasis and rosacea are the four most common skin diseases. Combining scRNA-seq with cell–cell communication analysis have provided critical insights into the disease pathogenesis and uncover potential opportunities for therapeutic interventions. Comparison analysis of non-lesional and lesional human skin from patients with AD using CellChat discovered major signaling changes in response to disease [15]. CellChat identified ligand–receptor pair CCL19–CCR7 as the most significant signaling that was specifically active in lesional skin, contributing to the communication from Inflammatory fibroblasts to Inflammatory dendritic cells. To investigate if F13A1+ macrophage subset and SNCG+ and ACKR1+ vascular endothelial cells were interacting with each other or other immune cells to coordinate leukocyte migration, CellPhoneDB was applied to assess cell–cell interactions in healthy, AD and psoriasis skin. CellPhoneDB predicted a significant enrichment for ACKR1 on VE3 to interact with CXCL8 (IL-8) on Mac2 and an enhanced interaction between VE3 and Mac2 with lymphocytes in AD and psoriasis compared with healthy skin, supporting a role for these cells in lymphocyte recruitment into inflamed skin [58]. Customized ligand–receptor analyses of scRNA-seq data collected from psoriatic and vitiligo skin have also identified signaling changes specific to diseases. For example, ligand–receptor analysis revealed the regulatory potential from resident epidermal/mesenchymal cells to dendritic cells during psoriasis [59]. More recently, cell–cell communication analysis of scRNA-seq of human vitiligo revealed that cell type–specific signaling programs and CCR5–CCL5 signaling was critical to effector CD8+ T cell and Treg function in vitiligo, implying the potential role of chemokine circuits in driving lymphocyte localization [60].

Skin cancer, including melanoma, basal cell carcinoma, and squamous cell carcinoma, is a disease in which malignant (cancer) cells form in the tissues of the skin. The advance understanding of known and new ligand–receptor pairs in the context of tumor-immune cell interaction within a tumor is extremely important for the further development of immunotherapies [67]. The tumor microenvironment (TME) comprises non-immune cells such as fibroblasts, blood and lymphatic endothelial cells, and numerous immune populations, aiding the growth and development of malignant cells. Application of CellChat to scRNA-seq data of esophageal squamous cell carcinoma (ESCC) highlighted the role of TNS1high fibroblasts in TME, in particular the interaction with cytotoxic T cells. Of note, TNS1high fibroblasts was associated with immune exclusion phenotype and poor prognosis of ESCC patients [61]. CellPhoneDB was applied to systematically study interactions within the TME in murine melanoma and identified stromal-immune interactions, such as C3–C3AR1, CXCL12–CXCR4 and CSF1–CSFR1 with macrophages as primary target [62]. Recently, cell–cell communication analysis has also been applied to spatial imaging data in skin cancer. A recent study revealed an immunosuppressive tumour-specific keratinocyte (TSK) subpopulation as a hub for intercellular communication by integrating scRNA-seq and spatial transcriptomics data in human squamous cell carcinoma [63]. NicheNet analysis predicted this subpopulation modulates TME-specific cell-type signatures at the leading edge and revealed TSK signaling to cancer-associated fibroblasts (CAFs) via MMP9–LRP1 and TNC–SDC1 and to endothelial cells via PGF–FLT1, PGF–NRP2, and EFNB1–EPHB4. Applying the computational pipeline STRISH to two types of skin cancer — basal cell carcinoma and squamous cell carcinoma — revealed considerable interaction of IL34–CSF1R around the areas where the cancer nests were located in both cancers, particularly in the epidermal compartments [43].

Altered intercellular communication is one of the nine candidate hallmarks of aging [68], which can directly affect tissue homeostasis and function. Single-cell sequencing allows us to comprehensive investigate altered cell–cell communication underlying aging processes, and offers novel therapeutic concepts to combat aging-associated skin diseases [69]. By comparing the inferred cell–cell communication in the young and old human skin using CellPhoneDB, one study showed that aging causes a substantial reduction in the interactions between dermal fibroblasts and other skin cells, including undifferentiated keratinocytes at the dermal-epidermal junction [64]. We anticipate there is growing insight into the roles of cell–cell communication in the aging process with the advance of various single-cell sequencing and imaging technologies.

The importance of cell–cell communication in tissue development and homeostasis has sparked a growing number of computational methods for inferring cell–cell communication from scRNA-seq data or spatial transcriptomics data. Recently systematical comparisons of ligand–receptor interaction resources and computational tools observed uneven coverage in terms of pathways and biological categories among the resources and varying predictions among the tools [14], making it difficult for researchers to choose an appropriate resource and method. Since computational inference is based on prior knowledge of ligand–receptor interactions, a unified resource similar to ominPath [24], with comprehensive literature-supported ligand–receptor interactions, is highly needed. Moreover, these issues also pose an urgent need to comprehensively benchmark the accuracy and robustness of each method using datasets with different biological variabilities, tissues, protocols, and platforms. While few studies have attempted to generate simulated datasets with ground truth against which various methods could be benchmarked [41,42,51], simulation of cell–cell communication with spatial information remains challenging due to the inherent complexity and redundancy of signaling mechanisms across multicellular systems. Existing case studies usually validate and compare inference results based on findings from biological experiments [15], bulk tissue data [20] or spatial transcriptomics data [34]. Therefore, we first need to establish a methodological framework that enables the simulation of cell–cell communication in different scenarios. Second, a curation of ground truth datasets would faithfully reflect differences in the performance of methods while avoiding the shortcomings of simulated data. Third, a set of quantitative metrics needs to be established that can better reflect the accuracy and robustness of methods for inferring cell–cell communication. Finally, a framework like ligrec_decoupler [14] should be established to facilitate a comparative assessment of methods and to provide a unified interface for users to applying diverse methods to their own data. We anticipate that the growing number of datasets from spatial transcriptomics will greatly boost the benchmarking of these computational tools.

Investigation of each method's advantages and limitations will provide practical guidelines for users. Since cell–cell communication often relies on multi-subunit protein complexes [70], methods like CellChat, CellPhoneDB and ICELLNET, which take into account multimeric ligand–receptor complexes, could potentially better prioritize biologically relevant interactions. However, modeling multimeric structure of ligand–receptor complexes is challenging due to possible ‘zero expression’ of subunits caused by dropout events in scRNA-seq data. Fortunately, dropouts are unlikely to affect strong signals predicted by these methods because dropouts commonly happen for genes with low expression [71–73]. Integrating other type of data such as single-cell proteomics [74] and CITE-seq data [75] could further improve the modeling accuracy. In addition, users can also choose methods based on the unique characteristics and capabilities of these methods, as shown in Figure 1. We believe that comprehensive benchmarking will provide better guidance for users and accelerate newly developed methods for the community.

Single-cell sequencing and imaging technologies are rapidly developing, in particular for examining the spatial transcriptome [3,6]. While cell–cell communication analysis has advanced considerably in recent years, enhanced computational approaches are clearly warranted. First, because of the spatial restrictions inherent to ligand diffusivity in juxtracrine and paracrine communication, consideration of spatial information is critical for cell–cell communication inference. Gene regulatory networks often drive cell fate and decision-making; however, most of the existing cell–cell communication inference methods do not consider the downstream response within the target cells. Integrating scRNA-seq data with spatial transcriptomics as well as gene regulatory networks offers a unique opportunity to study the spatial patterns of cell–cell communication and gene regulatory networks as well as to investigate how cell–cell communication affects cell-specific signaling networks within the context of a tissue. Second, sophisticated methods like Tensor-cell2cell [42] are highly needed to discern context-shared and -specific signaling patterns across conditions. Particularly, how to extract the biologically relevant signaling while removing the batch effects remains challenging [76]. Third, newly emerging experimental modalities from single-cell multi-omics such as single-cell proteomics and epigenomics [77–79] can further improve inference and our understanding of cell–cell communication, given that cell signaling occurs at the protein level instead of the gene level. Fourth, developing new computational methods that can incorporate phenotype and clinical information may be important to identify intracellular communication driving disease progression and enhance the predictive prognostic power of cell–cell communication. Last but not least, combining cell–cell communication analysis and mechanism-based systems biology modeling will likely deepen our understanding of the role of environmental signals in cell lineage fate decisions.

• Cell–cell communication orchestrates tissue development, homeostasis and disease. Single-cell sequencing and imaging technologies provide a unique opportunity to study cell–cell communication systematically and comprehensively.

• Integrating single-cell genomics and spatial transcriptomics data will greatly increase our understanding of the roles of cell–cell communication and cell-based interceptive medicine.

• Linking phenotype and clinical information with cell–cell communication is important for identifying signaling mechanisms driving cell fate decisions and disease pathogenesis and enhancing the predictive prognostic power of cell–cell communication.

The authors declare that there are no competing interests associated with the manuscript.

This work is partially supported by the National Natural Science Foundation of China under Grant number 11801207.

S.J. and R.R. wrote the manuscript. S.J. prepared the figures and tables.

AD

atopic dermatitis

ESCC

esophageal squamous cell carcinoma

TF

transcription factor

TME

tumor microenvironment

TSK

tumour-specific keratinocyte