Mapping and Validation of scRNA-Seq-Derived Cell-Cell Communication Networks in the Tumor Microenvironment.

Recent advances in single-cell technologies, particularly single-cell RNA-sequencing (scRNA-seq), have permitted high throughput transcriptional profiling of a wide variety of biological systems. As scRNA-seq supports inference of cell-cell communication, this technology has and continues to anchor groundbreaking studies into the efficacy and mechanism of novel immunotherapies for cancer treatment. In this review, we will highlight methods developed to infer inter- and intracellular signaling from scRNA-seq and discuss how they have contributed to studies of immunotherapeutic intervention in the tumor microenvironment (TME). However, a central challenge remains in validating the hypothesized cell-cell interactions. Therefore, this review will also cover strategies for integration of these scRNA-seq-derived interaction networks with existing experimental and computational approaches. Integration of these networks with imaging, protein secretion measurements, and network analysis and mathematical modeling tools addresses challenges that remain with scRNA-seq to enhance studies of immunosuppressive and immunotherapy-altered signaling in the TME.

Introduction

Single-cell technologies provide the resolution necessary to study intercellular heterogeneity. Specifically, single-cell RNA-sequencing (scRNA-seq) permits whole-transcriptome profiling of individual cells (1), which supports identification of rare cell subpopulations (2), inference of complex gene-gene regulatory networks (3), and tracing of developmental lineages (4). Further, scRNA-seq facilitates inference of cell-cell communication by measuring expression of genes encoding for corresponding ligands, receptors, intermediate signaling proteins, and intracellular targets across interacting cell types and under homeostatic and diseased conditions.

While scRNA-seq has been applied to a broad range of biological contexts, it has become an important tool to probe the tumor microenvironment (TME). The TME is a heterogeneous and constantly evolving milieu made up of populations of different cell types, cell-surface and secreted signals, and extracellular matrix (5, 6). These different cell types include cancer cells, stromal cells including fibroblasts and endothelial cells, and infiltrating immune cells, although this list is not exhaustive (7). Immune cells further partition into innate and adaptive immune cells, which are responsible for carrying out non-specific and tumor antigen-specific mechanisms of host defense, respectively (8). Innate immune cell types in the tumor predominantly include macrophages, dendritic cells (DCs), and neutrophils, while adaptive immune cell types include T cells and B cells. From the perspective of developmental lineage, it is worth noting that innate and adaptive immune cells do not derive from respective progenitors. Hematopoietic stem cells, which give rise to all immune cells, diversify into myeloid and lymphoid cell lineages (9). While myeloid cells can differentiate into the macrophages, DCs, and granulocytes of the innate immune system, lymphoid lineage cells contribute to both arms of immunity. For example, T cells and natural killer (NK) cells derive from a common lymphoid progenitor, although NK cells are considered innate immune cells due to their lack of tumor antigen-specific cell surface receptors (10).

Tumor cells promote their own survival by establishing an immunosuppressive and tolerogenic microenvironment. Through interactions with neighboring immune cells, tumor cells can manipulate immune functions to avoid immunosurveillance and coordinate delivery of key nutrients for growth (11). For example, tumor cells can take advantage of macrophage plasticity to reprogram them toward an anti-inflammatory phenotype (12). Anti-inflammatory tumor-associated macrophages (TAMs) have been shown to promote production of growth factors and neovascularization (13), which contribute to tumor growth and metastasis. These TAMs can also actively suppress adaptive T cell responses (14). Additionally, tumor cells can directly communicate with T cells via “immune checkpoints” (15), which act as gatekeepers to prevent an antitumor immune response from being mounted. Tumor cells can further act as architects of the spatial organization of the TME, orchestrating restriction of activated tumor-infiltrating lymphocytes to the periphery of the tumor bulk (16). This behavior yields an immune-excluded TME, which drastically limits the tumor-killing capacity of infiltrating immune cells. Broadly, immune cell infiltration into the tumor falls into three categories: immune-excluded, immune-inflamed (i.e., widespread presence in the tumor), and immune-desert (i.e., few to no cells in the tumor) (17). These patterns in infiltration are generally more applicable to the distribution of lymphoid lineage cells in the tumor, particularly T cells, while myeloid cells are commonly more ubiquitously distributed. Overall, the TME represents a spatially organized community of dysregulated cells, with communication between them giving rise to emergent behaviors including immune evasion and tumor growth.

Immunotherapy has emerged as a leading treatment option for advanced-stage tumors (18). By breaking immunosuppressive tumor-immune cell interactions and reestablishing anti-tumor immune signaling, this therapeutic option takes advantage of the body’s immune system to restore host defenses that were initially hijacked by the tumor. Immunotherapies have improved survival outcomes in numerous cancer types, and ongoing clinical trials serve to enhance the pool of effective interventions (19). Advancing our understanding of aberrant signaling in the TME directly supports the development of new immunotherapeutic strategies. Robust construction of hypothetical intercellular interaction networks from scRNA-seq is therefore critical to the design and testing of novel immunotherapeutic drugs that will ultimately continue to improve patient outcomes. This review will thus begin with an overview of existing computational methods for inference of communication between and within cells from scRNA-seq ( ). Particular attention will be paid to their contributions to the fields of tumor immunology and immunotherapy development.

Figure 1

Once scRNA-seq-derived cell-cell communication networks are mapped, spatial profiling, protein-level measurements, and computational pruning can aid validation and interpretation. Cell-cell communication can be inferred from scRNA-seq using a variety of computational methods, such as CellPhoneDB, CellChat, NicheNet, or CytoTalk. Experimentally, connections between cell populations can be validated with spatial profiling and protein-level measurements. Networks can also be further analyzed computationally with network analysis and mathematical modeling tools. Created with BioRender.com.

Despite these advances, a central challenge remains in the validation of the hypothesized cell-cell interactions. Computationally constructed communication networks often comprise hundreds to thousands of signaling interactions, which can make these networks impossible to probe experimentally in their entirety. Therefore, we will also review approaches to integrate these scRNA-seq-derived interaction networks with an array of existing experimental and computational methods to yield more meaningful biological conclusions ( ).

Initial Strategies to Infer Cell-Cell Communication From scRNA-Seq Focused on Matched Expression of Paired Receptor and Ligand

scRNA-seq affords the ability to measure expression of ligands and receptors to systematically decode how cells are communicating with one another. Evaluation of these interaction networks across homeostatic and diseased conditions facilitates identification of experimentally testable targets for tumor immunotherapy, among countless other applications. Many computational methods have been developed to examine cell-cell communication with scRNA-seq, and they exist along a spectrum with respect to their computational complexity and integration of additional signaling information. It should be noted that these algorithms require clustering of single cells into cell types or subsets of interest. There are many strategies for clustering and annotation of scRNA-seq data, which have been reviewed elsewhere (20–23).

Initial approaches generated hypotheses about cell-cell communication by quantifying matched expression of corresponding ligand and receptor. Zhou et al. (24) constructed the first cell-cell communication network of its kind from scRNA-seq collected from 19 patients with primary and metastatic melanoma (25). Interactions between malignant, immune, stromal, and endothelial cells were hypothesized by identifying paired receptor- and ligand-encoding genes from an author-curated list for which respective expression in two cell type populations was three standard deviations above the average across all populations. While the authors were then able to characterize tumor-immune, tumor-stromal, and tumor-endothelial cell crosstalk through specific signaling axes, they recognized the need to link this information with biological outcomes of interest, like tumor growth.

To address this limitation, Kumar et al. (26) developed a computational technique to quantitatively relate scRNA-seq-derived interaction networks from mouse and human metastatic melanomas to pathophysiological characteristics of the TME. They quantified cell-cell communication by calculating the product of average ligand and receptor expression across interacting cell type populations (referred to as interaction scores), with literature-supported ligand-receptor pairs sourced from the FANTOM5 database (27). To filter out interactions with inappropriately high scores due to noisy gene expression, the authors used statistical comparisons against null networks. Correlating interaction scores against tumor growth rate and anti-tumor immune response hypothesized ECM- and adhesion-related interactions, in addition to certain cytokine and chemokine interactions, as potential contributors to tumor progression in mice. This approach has thus helped to lay the foundation for identifying experimentally testable targets for immunotherapy from scRNA-seq-derived cell-cell communication networks.

Emergence of Open-Source Tools for Interaction Network Construction From scRNA-Seq Has Advanced the Development of Cancer Immunotherapies

Despite the progress posed by the aforementioned methods, scalable application to existing scRNA-seq datasets remained a challenge as user-friendly platforms for construction of interaction networks did not exist. CellPhoneDB (28, 29) presented one of the first open-source tools for inference of cell-cell communication from single-cell transcriptomic data, which was originally applied to scRNA-seq from the human placenta. It quickly became one of the most utilized algorithms for the task, with their manuscripts amassing over 1,200 citations collectively since publication, and their online resource (www.cellphonedb.org) being used by over 500 users per month as of July 2020 . As with previously discussed approaches, CellPhoneDB probes cell-cell interactions by investigating matched expression of ligand-receptor pairs sourced from the UniProt (30, 31), Ensembl (32), Protein Data Bank (PDB) (33, 34), International Molecular Exchange (IMEx) consortium (35), and IUPHAR (36–38) databases. Critically, CellPhoneDB considers subunit architecture for both ligands and receptors, as opposed to the binary representation adopted by most methods.

CellPhoneDB has been widely used to understand cell-cell communication in the TME, particularly to characterize pro-tumor crosstalk between tumor, immune, and stromal cells. Macrophages are a type of immune cell that are especially vulnerable to exploitation by tumor cells because their activation is heavily context-dependent (39). Therefore, uncovering the signaling mechanisms by which tumor cells reprogram macrophages toward an anti-inflammatory pro-tumor phenotype, along with the mechanisms by which macrophages then carry out their pro-tumor behavior, will provide new targets for immunotherapeutic intervention. In both hepatocellular carcinoma (40) and esophageal squamous cell carcinoma (ESCC) (41), CellPhoneDB implicated the SPP1-CD44 signaling axis, among other receptor-ligand pairs, as a potential reprogramming interaction from tumor cells to macrophages. This axis has been previously implicated as an immune checkpoint in human cancers (42), and tumor cell signaling to macrophages and monocytes through the CD44 receptor has been shown to inhibit their anti-tumor response in in vitro co-culture models (43). In colorectal cancer (CRC), Zhang et al. (44) used CellPhoneDB to characterize the anti-inflammatory signaling from macrophages that had already shown evidence of pro-tumor reprogramming. Their analysis linked interactions involving SDC2, SPP1, and FN1 ligands between macrophages and cancer-associated fibroblasts, myofibroblasts, and endothelial cells to tumorigenesis and metastasis in CRC.

Beyond identifying potential therapeutic targets by characterizing pro-tumor signaling in the TME, CellPhoneDB has aided studies in understanding immunotherapy-altered signaling from single-cell transcriptomics. Immune checkpoint blockade (ICB) is the current standard of care for immunotherapeutic intervention for many cancers (45). Over time, anti-tumor effector immune cells in the tumor become “unresponsive” or “exhausted,” and become negatively regulated by ligands expressed in the TME. ICB aims to prevent or revert that immunosuppression by blocking these inhibitory interactions. One of the common targets of ICB is the programmed cell death protein 1 (PD-1) checkpoint. Signaling through PD-1 has been shown to negatively regulate T cell effector activity (46), and thus its blockade should limit the ability of tumor cells to evade immune surveillance and promote a successful anti-tumor response. However, mixed response rates to ICB in patients highlight the need for a deeper understanding of the mechanism of this therapy in the microenvironment. Studies in patients with breast cancer (47), advanced renal cell carcinoma (48), and basal cell carcinoma (49) have leveraged CellPhoneDB to characterize cell-cell communication associated with response and resistance to ICB. Across tumor types, this analysis identified signaling to CD8+ T cells via HAVCR2-LGALS9 (TIM3-Galectin9) as enhanced in non-responding and resistant patients. As TIM3 has been known to inhibit CD8+ cytotoxic activity and therefore can contribute to tumor progression (50), these studies emphasize the promise of CellPhoneDB to uncover conserved biomarkers for therapy resistance from scRNA-seq.

It is worth noting additional open-source packages for inference of cell-cell communication from scRNA-seq that use algorithms similar to CellPhoneDB. Network Analysis Toolkit for Multicellular Interactions (NATMI) (51) quantifies potential communication by calculating the product of average expression of matched ligand and receptor across two interacting cell type populations. The novelty of this method primarily derives from the curation of its reference database of ligand-receptor pairs, which borrows from competing methods and only includes pairs with primary literature support. NATMI does not consider subunit architecture for ligands or receptors. This method has been used to identify prognostic communication between tumor, stromal, and immune cells in patients with pancreatic ductal adenocarcinoma (PDAC) (52). CellTalker (53) probes matched expression of paired ligands and receptors from the FANTOM5 (27) database. This algorithm also does not take heteromeric complexes into account. While CellTalker does limit its analyses to genes that are differentially expressed across cell populations, which supports identification of population-specific signaling, it does not assign quantitative scores to potential interactions which restricts data-driven comparisons across conditions of interest. Finally, ICELLNET (54) also assigns interaction scores by multiplying the average expression of paired ligand and receptor, with an adjustment to this formula intended to account for receptor subunits. The ligand-receptor databases referenced by ICELLNET largely overlap with its competitors, including CellPhoneDB. This method has yet to be applied to investigate the TME in a published study.

Functional Understanding of Cell-Cell Communication may Require More Than Paired Expression of Receptor and Ligand

CellPhoneDB in particular has advanced the study of intercellular signaling in the tumor microenvironment from scRNA-seq with its scalable and accessible platform. However, functional cell-cell communication extends beyond the interaction of extracellular ligand and receptor, as signaling cofactors can upregulate or antagonize downstream signal propagation, and ligands from “sender” cells can influence intracellular signaling cascades and target gene expression in “receiver” cells.

CellChat (55) is an algorithm for cell-cell communication inference from scRNA-seq that was developed to address the lack of modulation by signaling cofactors, such as soluble agonists, antagonists, and membrane-bound co-receptors, considered in previous methods. A mass action-based model underlies the interaction probability assigned to each ligand-receptor pair for each communicating set of cell populations. As opposed to inferring cell-cell communication directly from gene expression, CellChat first projects the transcriptomic profiles onto experimentally validated protein-protein interaction networks (56) using a network propagation approach (57). These projected profiles are then used to calculate interaction probability scores, with the scoring equation taking advantage of Hill functions to model positive and negative modulation of an interaction by signaling agonists and antagonists, respectively. This intermediate projection step is intended to account for the role of protein-protein interactions in ligand-receptor binding, which cannot be reflected in sparse gene expression profiles from scRNA-seq data without prior knowledge of how the proteins encoded by these genes can influence each other’s expression. Using the law of mass action addresses the fact that these protein-protein interactions can saturate.

First applied to the wound healing microenvironment, CellChat has shown promise in the characterization of malignant and immunotherapy-altered signaling in glioblastoma (58), esophageal squamous cell carcinoma (41), gastric cancer (59), and breast cancer (60). When compared directly to CellPhoneDB (55) across four mouse skin scRNA-seq datasets, CellChat achieved a slightly better true positive rate, a lower false positive rate, and higher accuracy. Accuracy in this context was defined as having more predictions that overlapped with the predictions of one or more methods, which in this study included the methods iTALK (61) (which has been published as a preprint) and SingleCellSignalR (62) in addition to CellPhoneDB. CellChat’s advantage over CellPhoneDB may also relate to the fact that this method limits its inquiries to genes that are differentially expressed (DEGs) across all cell groups in the data, as opposed to maintaining the transcriptome-wide gene set, which simplifies comparison across conditions. Including this same DEG-limiting step in the CellPhoneDB workflow, which was implemented in the newest version of the CellPhoneDB toolkit (CellPhoneDB v.3.0) (63), has improved its predictions as well. Further, CellChat and CellPhoneDB were found to similarly distinguish spatially adjacent from spatially distant interactions (55), unlike competing algorithms that only consider expression of paired receptor and ligand. This result supported the assertion that including subunit architecture and signaling cofactors serves as a proxy to help distinguish spatially relevant interactions from scRNA-seq. CellChat can also be distinguished from many of its competitors by its built-in tools for visualization and downstream analyses, which use concepts from graph theory, pattern recognition, and manifold learning to identify conserved and context-specific patterns across inferred signaling networks, among other important interpretations.

Other recently published algorithms NicheNet (64) and CytoTalk (65) have recognized the need to integrate intracellular signaling into scRNA-seq-derived communication networks. NicheNet was the first computational method to consider prior knowledge of intracellular signaling networks in its analysis of scRNA-seq by leveraging databases describing interactions between ligands, receptors, intracellular intermediaries, transcriptional regulators, and target genes. Strong correlation between expected patterns in gene expression, based on links between ligands and target genes in NicheNet’s prior model, and observed patterns in a user’s own scRNA-seq data identifies potential communication between interacting cell types. Similar to CellChat, NicheNet also takes advantage of DEGs to construct differential signaling networks across conditions of interest.

The use of a prior model with ligand-to-intracellular target gene axes to infer cell-cell communication from scRNA-seq data represented a departure from the interaction scoring approach matching ligand and receptor expression adopted by most preceding methods. Algorithmically, preceding methods follow a common procedure ( , column 1). Taking annotated scRNA-seq data as input, they calculate interaction scores by combining average expression of ligand and receptor across two interacting cell populations with an elementary mathematical operation, generally multiplication. For example, Kumar et al.(26) and NATMI (51) rely on a simple product ( , column 1, ii. A), while ICELLNET (54) includes receptor subunits in its formula ( , column 1, ii. B). SingleCellSignalR (62) uses a regularized product, which accounts for total expression across cells ( , column 1, ii. C). CellChat (55) technically also uses a regularized product, although its interaction scoring equation takes inferred protein expression as input as opposed to mRNA counts as competing algorithms do. Downstream analyses then generally focus on comparing these interaction scores statistically to null scores (to distinguish signal from noise) and across experimental conditions. As previously described, NicheNet instead correlates patterns in observed expression from annotated scRNA-seq data against expected expression in its prior model ( , column 2). Expected expression of each gene in the prior model, given a particular stimulating ligand, is mathematically derived based on predicted propagation of signal from ligand to target gene through possible intermediate signaling genes. Downstream analyses with NicheNet generally further explore links between predicted ligands, receptors, and target genes in the prior model.

Table 1

Comparison of matched ligand-receptor (L-R) expression algorithms versus correlation against a prior model (e.g., NicheNet) for inference of cell-cell communication from scRNA-seq.

Matched L-R expression	Correlation against a prior model
i. Input: preprocessed scRNA-seq data annotated by cell type	i. Input: preprocessed scRNA-seq data annotated by cell type
ii. Calculate interaction score (IS) between cell types j and k, based on known ligand-receptor pairs (reference databases listed in Table 2), which takes on one of three general forms based on the method:	ii. Identify “target genes” in receiving cell populations (e.g., genes differentially expressed across conditions); remaining cells considered “background”
IS(L,R,j,k)={L_j∗R_k,L_j+R_k(A)(∏zL_j,z^z)∗(∏mR_k,m^m)(A) or (B)regularization term(C)(B)	iii. For each possible ligand i, calculate Pearson correlation (ρi) between observed expression of target genes j and expected expression of j given i (Ei,j) from prior model; ρi is penalized for correlations with background genes. Ei,j is defined as:
where L_j=(1n_j∑i=1n_jl_i),R_k=(1n_k∑i=1n_kr_i),	E_i,j=∑k=1n(PPR_i,k)*(GRN_k,j)
Or Lj , Rk are otherwise scaled representations of expression, e.g., z-scores, and z and m account for ligand and receptor subunits, respectively	where, based on the prior model, PPRi,k denotes the probability of ligand i signaling through intermediate gene k, and GRNk,j quantifies evidence that k regulates expression of target gene j. n is the number of possible intermediate genes k
iii. Downstream analyses: compare IS to null networks and across experimental conditions	iv. Downstream analyses: further explore ligand-target gene and ligand-receptor axes in the prior model

NicheNet has been cited over 200 times since its publication in 2020, with many of the citing papers investigating the TME. Cheng et al. (66) presented a prime example of NicheNet’s utility in studying myeloid signaling across tumor microenvironments. The authors investigated conventional dendritic cells (cDCs), mast cells, and tumor-associated macrophages (TAMs) using scRNA-seq data from 210 patients across 15 human cancer types. Following previous work from Maier et al. (67), cDCs were clustered across tumor types into two classical DC subsets (cDC1s and cDC2s) and a novel LAMP3+ mature cDC subset enriched for immunoregulatory genes (mregDCs). Cheng et al. were keen to understand mregDC development, as this novel DC subset has been shown to originate from both cDC1s and cDC2s, resulting in different roles in anti-tumor immunity. Using NicheNet, IL-4 and IL-12 were identified as drivers of the cDC1-derived mregDC signature, while IL-15 and IFNs were linked to the cDC2-derived mregDC signature (66). Modulating these NicheNet-identified signals in the TME has the potential to adjust the balance of cDC1- and cDC2-derived mregDCs in a way that supports an anti-tumor immune response. The authors also probed mast cells in the TME, as TNF+ and pro-angiogenic VEGFA+ mast cells have been thought to play dueling roles in the tumor (68). Both NicheNet and CellPhoneDB were used to uncover potential signaling underlying the TNF+ mast cell signature. Interestingly, while both methods implicated IL1β+ macrophages as the drivers of anti-tumor TNF+ mast cell behavior, the full suite of their top results offered little overlap. Whether the inclusion of intracellular signaling information in NicheNet’s method generates more accurate or biologically relevant hypotheses than algorithms like CellPhoneDB requires further experimental validation.

Although NicheNet advanced cell-cell communication inference from scRNA-seq by presenting a more mechanistic picture of interaction though ligand-to-target gene axes, its prior model does not include the additional prior knowledge that many signaling pathways are specific to a given cell type or tissue. Hu et al. (65) illustrated this heterogeneity in signaling downstream of fibroblast growth factor receptor 2 (FGFR2) between three pairs of interacting cell populations across mammary gland and skin. They therefore developed CytoTalk, which similarly constructs scRNA-seq-derived cell-cell communication networks through ligand-to-target gene axes, but with data-driven inference of signal transduction networks instead of a prior model. While NicheNet’s intracellular signaling networks were derived from published databases including OmniPath (69), CytoTalk constructs these networks de novo by quantifying relationships between genes from the user’s own data. CytoTalk has yet to be applied to the TME in a published study, although the utility of this method is clear considering the roles cell type and tissue specificity have been shown to play in signaling in cancer (70). Altogether, much progress has been made to extract increasingly accurate and biologically relevant information about cell-cell communication from scRNA-seq.

Underlying Reference Databases Play a Key Role in the Inference of Cell-Cell Communication

It is worth mentioning that, across methods, predictions about cell-cell communication rely heavily on the underlying reference databases summarizing known ligand-receptor pairs and intracellular signaling pathways. Many of the aforementioned methods source interaction information from a similar pool of databases ( , “Referenced databases”). The FANTOM5 database (27) specifically, which includes ~2,500 human ligand-receptor pairs, is referenced by nearly all of the selected algorithms in . Further comparison of the general features, advantages, and limitations of FANTOM5 against other commonly used reference databases, including those underlying CellPhoneDB, NATMI, ICELLNET, CellChat, iTALK, SingleCellSignalR, and NicheNet, has been reviewed elsewhere (89, 90).

Table 2

Selected summary of computational methods for inference of cell-cell communication.

Method	Application to TME?	Signaling components considered	Underlying mathematical model	Referenced databases	Primary advantages	Primary limitations*
Zhou et al. (24)	Melanoma	Paired ligand and receptor (L-R)	Matched significant upregulation of L and R	FANTOM5, DLRP (71), IUPHAR, HPMR (72), author additions	Laid foundation for interaction inference	Lacks other signaling components, statistical framework, links to emergent behaviors
Puram et al. (73)	Head and neck cancer	Paired L-R	Matched significant upregulation of L and R	FANTOM5	Laid foundation for interaction inference	Lacks other signaling components and statistical framework
Kumar et al. (26)	Melanoma	Paired L-R	Product of average L and R	FANTOM5, author additions	Analysis framework links interactions to emergent behaviors	Lacks other signaling components
Raredon et al. (74)	No	Paired L-R	Sum of average L and R	FANTOM5	Uses graph theory for cross-network comparison	Lacks other signaling components
CellPhoneDB (28, 29)	ESCC, CRC, breast cancer, among others	Paired L-R, including subunits	Matched significant upregulation of L and R	UniProt, Ensembl, PDB, IMEx, IUPHAR	Scalable, user-friendly platform; includes subunits	Lacks other signaling components
NATMI (51)	PDAC	Paired L-R	Product of average L and R	CellPhoneDB, SingleCellSignalR, ICELLNET, STRINGDB (56), RNA-magnet (75), author additions	Considers L-R specificity in interacting cell types; underlying database	Similar to predecessors, lacks other signaling components
CellTalker (53)	Head and neck cancer	Paired L-R	Matched expression of L and R	FANTOM5	Uses differentially expressed genes for focused exploration	Similar to predecessors, lacks other signaling components and statistical framework
ICELLNET (54)	No	Paired L-R, including receptor subunits	Product of average L and R	STRINGDB, Ingenuity, BioGRID (76), Reactome (77), CellPhoneDB	Includes receptor subunits	Similar to predecessors, lacks other signaling components
CellChat (55)	GBM, ESCC, breast cancer, among others	Paired L-R and signaling cofactors	Mass action-based model	KEGG (78), author additions	Includes signaling cofactors, uses graph theory for cross-network comparison	Lacks intracellular signaling information
iTALK (61)	Lung adeno-carcinoma (79)	Paired L-R	Matched significant upregulation of L and R	FANTOM5, DLRP, IUPHAR, HPMR, author curation of cytokine and chemokine interactions (80–82)	Author-curated reference interactions, built-in tools for visualization	Lacks other signaling components
SingleCellSignalR (62)	CRC (83), breast cancer (84), metastatic melanoma (85), among others	Paired L-R	Regularized product of average L and R	FANTOM5, HPMR, IUPHAR, UniProt, HPRD (86), GO (87), Reactome	Author-curated reference interactions	Lacks other signaling components
NicheNet (64)	Melanoma, CRC, breast cancer, among others	L-R and downstream intracellular signaling	Pearson correlation of user data with prior model	FANTOM5, IUPHAR, KEGG, OmniPath (69), MSigDB (88), among other PPI, gene regulatory interaction databases	Includes intracellular signaling components	Prior model doesn’t consider cell type or tissue specificity
CytoTalk (65)	No	L-R and intracellular signaling in both “sender” and “receiver”	Regularized sum of average L and R, mutual information between genes for inference of intracellular signaling	FANTOM5, author additions	Considers cell type and tissue specificity	Lacks extracellular signaling components (subunits, cofactors)

This column highlights limitations of communication inference methods when compared to each other. Limitations across methods (e.g., lack of spatial information) are discussed in the text.

Among many attributes of reference databases, there are two to consider when evaluating predictions from cell-cell communication inference methods. The first is the proportion of interactions in each database with primary literature evidence. For example, Hou et al. (51), in the manual curation of their own interaction database, excluded 94 ligand-receptor pairs from CellPhoneDB and 143 from the Human Protein Reference Database (HPRD) (86) for a lack of literature support. Ensuring primary experimental evidence for predicted interactions increases the accuracy and biological relevance of hypotheses generated by these cell-cell communication inference tools. The second consideration is the frequency with which these databases are updated with newly discovered or otherwise revised interactions. Whether a database of interactions reflects the current state of the literature further impacts a user’s interpretation of the predictions made by these methods on their own datasets. Altogether, there are many caveats to keep in mind when deriving cell-cell communication networks from scRNA-seq, with the underlying reference database of ligand-receptor and ligand-target gene axes being one of them.

Shortcomings Across Existing Computational Methods Derive From Challenges That Remain With scRNA-Seq

scRNA-seq experiments require tissue dissociation (91, 92), which prevents the preservation of critical spatial information. Although the aforementioned computational methods (summarized in ) leverage ligand-receptor pair databases and integrate critical information about intracellular and other relevant signaling, it remains difficult to overcome this key shortcoming in the method of data collection. Although some ligands are soluble, distances traveled by diffusion are limited; in other cases, both receptor and ligand are membrane-bound. Thus, inference of cell-cell communication would be greatly improved by incorporating knowledge of the physical proximity between two cell populations in the microenvironment. For example, tumors can disrupt successful immune responses by limiting T cell infiltration into their bulk (93), thus restricting the physical interaction between T cells and other cell types. In this case, integration of spatial information is crucial to filtering edges in scRNA-seq-derived communication networks describing these pathologies.

Additionally, heterogeneous signaling among broad cell types and cell subpopulations within a single cell type often differentially contribute to cancer progression. All cell-cell communication inference methods depend on how cell populations are defined. While clustering single cells into broad cell types has been extensively explored (20–23), subpopulation definition remains a contested topic. From scRNA-seq, unsupervised clustering of broader cell type lineages is often used to delineate cell subsets; this approach has been extensively applied to understand heterogeneity in the myeloid compartment of sarcoma (94) and non-small-cell lung cancer (95) tumors, among other malignancies (58, 66, 96, 97).

Because mRNA is not a functional measurement of the cell, defining cell subpopulations based on gene expression does not always separate cells into subsets that are relevant with respect to their functional signaling behavior. Directly measuring extracellular signaling events at the protein level partly overcomes this limitation. Single-cell secretomic technologies in particular have shown promise in the characterization of functional heterogeneity in the TME. For example, Perry et al. (98) used a single-cell secretion device to investigate myeloid heterogeneity in the melanoma TME, which identified specific functional subgroups of tumor-associated macrophages and myeloid cells (TAMMs) that aid orchestration of a pro-tumorigenic immune response. Inference of informative interaction networks from scRNA-seq thus also requires integration of protein-level measurements to confirm secretion of ligand and separate large cell-type populations into functionally relevant, context-dependent subpopulations.

Moreover, once the cell-cell communication networks are constructed, they need to be systematically compared across experimental conditions and interpreted with respect to emergent biological behaviors. As these networks are directed, weighted, and densely connected, data mining becomes an increasingly complex task. Therefore, this review will next discuss strategies for network validation and interpretation by leveraging imaging, protein secretion measurements, and network analysis and mathematical modeling tools. We anticipate that these integrations will aid generation of more biologically accurate and actionable hypotheses about aberrant and altered signaling mechanisms underlying tumor progression and as a result of immunotherapeutic intervention, respectively, from scRNA-seq.

Spatial Information Is Critical to Confirming the Possibility of Interaction Between Cells in the TME

Because scRNA-seq requires tissue dissociation, the spatial information necessary to fully infer intercellular communication, particularly between membrane-bound receptors and ligands, is lost. Knowing the spatial localization of cells is especially critical to studies of the TME, as tumors can become “immune-excluded” (99), restricting the dialog between cancer and immune cells and rendering typical immunotherapeutic strategies useless. For example, T cell exclusion from the tumor bulk places crucial limits on the efficacy of T cell-targeted treatments, like ICB or adoptive cell transfer therapies. The field of spatially resolved transcriptomics (SRT) has begun to address the issue of a lack of spatial information by enabling confirmation of the physical proximity of cells expressing transcripts for a corresponding ligand and receptor pair. The advent of SRT has been particularly exciting for cancer research, yielding preliminary spatial maps of the prostate (100) and melanoma (101) TMEs, among others. With SRT, characterizing the TME in situ has suddenly evolved from low throughput “snapshots” given by immunohistochemistry and in situ hybridization to unbiased molecular profiling without compromising spatial context. For many of these reasons, SRT was named Nature Method’s “Method of the Year” for 2020 (102).

Due to the complexity of preserving spatial information, SRT technologies do lag behind scRNA-seq in depth of sequencing, resolution, and accessibility. Some SRT methods require a preselected panel of genes for profiling (103), while others have yet to achieve single-cell resolution. As opposed to replacing scRNA-seq, SRT methods thus present an opportunity to extract more biologically relevant hypotheses about cell-cell communication in the TME through integration with scRNA-seq. This is particularly important given the opportunity to further mine the plethora of scRNA-seq datasets available in the literature as a result of both independent studies and collaborative efforts like the Human Cell Atlas (104). Recently, several excellent in-depth reviews of SRT have been published (103, 105–107). Here our goal will be to highlight how these SRT methods have been integrated with scRNA-seq to enhance the investigation of interactions between cell populations in the TME. Moncada et al. (108) presented one of the first examples of the integration of the recently developed spatial transcriptomics (ST) method (109), capable of transcriptome-wide sampling but at low resolution, with scRNA-seq to understand the tissue architecture of pancreatic ductal adenocarcinoma (PDAC). With ST, the authors were able to identify patterns in spatially restricted enrichments and co-enrichments of tumor, myeloid, and healthy ductal cells with other cell types in the TME. Although the ST array cannot cover the entire tissue, nor does it achieve single-cell resolution, integration with scRNA-seq supported a focused exploration of the relationship between the co-localization of specific cell types and their potential to communicate in the PDAC TME.

Two other examples of SRT methods, multiplexed error-robust fluorescent in situ hybridization (MERFISH) (110, 111) and GeoMX Digital Spatial Profiling (DSP) (112), have aided studies of spatial organization in the TME. MERFISH is an imaging-based approach capable of generating single molecule resolution spatial maps of thousands of preselected genes (i.e., not transcriptome wide, a fundamental trade-off) at high detection efficacy. In glioblastoma (GBM), MERFISH revealed an enrichment of macrophages in the neighborhood of a specific subpopulation of malignant cells in both mouse models and patient samples (113). Integration of this information with scRNA-seq lent support for five ligand-receptor axes through which macrophages could be inducing a mesenchymal-like state in this GBM cell subtype. As mesenchymal-like GBM has been associated with more aggressive, invasive, and inflammatory features of the TME (114), the integration of MERFISH and scRNA-seq by Hara et al. directly resulted in the identification of targets for immunotherapeutic intervention in GBM. MERFISH has also been used to spatially map all cell types and states identified by scRNA-seq in CRC patient samples (115), linking spatial organization in the tumor to previously predicted cell-cell communication networks.

GeoMX DSP is an in situ barcoding method for multiplexed profiling of a predefined panels of RNAs or proteins from fixed tumor tissues (112). This method has been used to identify spatially informed biomarkers of response and resistance to immunotherapies such as ICB across multiple cancer types, including melanoma (116), non-small cell lung cancer (NSCLC) (117) and head and neck squamous cell carcinoma (HNSCC) (118). In melanoma and NSCLC, PD-L1 expression in the CD45+CD68+ spatial compartment (i.e., macrophages) was found to correlate with prolonged progression-free survival and overall survival in patients. Similar to MERFISH, GeoMX DSP thus provides high resolution imaging of focused regions of the TME, supporting local rather than global evaluations of ligand, receptor, and cell type co-localization.

It is worth noting that additional spatial profiling methods continue to be developed with potential to advance the study of the TME. Deterministic barcoding in tissue for spatial omics profiling (DBiT-seq) is a microfluidics-based method for 10μm-resolution multi-omic profiling of fixed tissues (119). By delivering barcodes to the surface of tissue slides, the authors were able to generate two-dimensional mosaics of protein and transcriptome-wide mRNA expression in mouse embryos. DBiT-seq specifically capitalized on the shortcomings associated with the technically demanding nature of in situ hybridization methods like MERFISH and sequential fluorescence in situ hybridization (seq-FISH+) (120) by taking advantage of next generation sequencing (NGS) technology and resources. Unlike GeoMX DSP and similar barcoding strategies [e.g., Slide-seq (121)], DBiT-seq uses microfluidics for guided delivery of barcodes to tissue, which greatly improves spatial accuracy. As of yet, DBiT-seq has not been used to map the spatial organization of a tumor in a published study, although it is likely only a matter of time due to the versatility and ease of use of this technology. The transcriptomic measurements captured by DBiT-seq would directly support inference of intercellular interactions at nearly (but not quite) single-cell resolution. DBiT-seq, like many SRT technologies, is still limited by pixel size and the total tissue area able to be mapped (119). Taken together, SRT methods, while not yet replacing scRNA-seq due to limitations in multiplexing and resolution, have and can continue to advance investigations of cell-cell communication in the TME when analyzed in conjunction with scRNA-seq.

Beyond SRT methods, an array of complementary imaging techniques for both dynamic and static receptor-ligand and sender-receiver cell co-localization has been developed and used in the TME. Intravital microscopy (IVM) (122) has emerged as a leading tool to visualize in vivo dynamics and migration in tumor models at single-cell resolution. Following surgical implantation of an imaging window, live fluorescent imaging of the TME can be collected, generally over a period of hours to days. Luthria et al. (123) showcased the utility of IVM in the identification of dynamic physical interactions between immune and tumor cells in vivo across multiple cancer models. TAMs were found to often neighbor and intertwine themselves around tumor cells, confirming physical interaction between the two cell populations. Moreover, the authors linked TAM polarization state to microtubule dynamics in tumor cells, as quantified by IVM, elucidating an association between the anti-inflammatory or pro-tumorigenic macrophage state and increased cancer cell migration and metastasis. Integration with scRNA-seq-derived communication networks generated hypotheses about the role that autocrine signaling (i.e., TAM-TAM interaction) played in anti-inflammatory TAM polarization. Using the cell-cell communication inference strategy from Kumar et al. (26) and scRNA-seq from patient biopsies across cancers, the authors implicated the interleukin 10 (IL10)-IL10R(A,B) axis as a critical feedback loop in TAMs that restricts their reprogramming toward a pro-inflammatory or anti-tumor phenotype (123). Treatment with an IL10R-blocking antibody was further found to stunt tumor cell motility and migration both in vitro (in a co-culture model system) and in vivo. The cell-cell interaction networks inferred from scRNA-seq therefore directly provided an experimentally testable target for immunotherapeutic intervention in this study, although the authors did primarily focus on crosstalk between immune cell populations in these hypothetical communication networks. Investigation of TAM-tumor cell crosstalk as inferred from scRNA-seq was not presented in this work. The focus of their IVM data was to demonstrate physical interactions between TAMs and tumor cells in the TME and pinpoint these interactions as a catalyst for downstream alterations in microtubule dynamics. Reanalysis of the scRNA-seq data to instead concentrate on inferring extracellular signaling events between TAMs and tumor cells has the potential to identify additional mechanisms of cancer cell migration and anti-inflammatory TAM polarization for therapeutic targeting with direct support from the IVM work.

Further, IVM has helped to elucidate the mechanism of anti-PD-1 checkpoint inhibition underlying effective anti-tumor responses (124). IVM revealed co-localization of CD8+ cytotoxic T cells and dendritic cells (DCs), producing IFNγ and IL-12 signals, respectively, upon treatment with anti-PD-1. This T cell-DC crosstalk was thus found to form a positive feedback loop resulting in stimulation of anti-tumor T cell immunity that was indirectly downstream of anti-PD-1 binding. Garris et al. (124) collected paired scRNA-seq from the MC38 colon carcinoma TME to support their claim that IL-12 production was exclusive to DCs, although inference of cell-cell communication from this data was not reported. Reanalyzing this scRNA-seq data to predict differential signaling between immune and tumor cells upon checkpoint inhibition could further the conclusions of this work by contextualizing the IFNγ/IL-12-mediated signaling in the greater network of cell-cell communication at play in the tumor, with support for communication between certain cell type populations coming from their IVM co-localization data. Other studies, which characterized crosstalk between immune cell populations (125) and between tumor cells and healthy tissue (126) in the TME with IVM, similarly collected paired scRNA-seq but did not present scRNA-seq-derived interaction network analysis. These works highlight existing opportunities in the field for reanalysis of publicly available transcriptomic data to infer and validate specific mechanisms and global signaling underlying emergent behaviors. IVM has also shown promise in studying the targeting (127, 128), dynamics, and effects on angiogenesis (129) of nanoparticle-based therapies for immunomodulation in the TME. This information would aid both validation and interpretation of therapy-altered communication networks inferred from scRNA-seq.

While IVM affords in vivo measurement of interacting cell co-localization and dynamics in the TME, this method is limited by the number of fluorophores that can be measured simultaneously. The need for an imaging window also makes IVM more suitable for collecting spatial information from animal models than humans, restricting the validation of cell-cell communication networks inferred from patient-derived scRNA-seq, although human intravital microscopy (HIVM) is under development (130). Albeit from fixed tumor slices, imaging methods such as tissue-based cyclic immunofluorescence (t-CyCIF) (131) can visualize a greater multiplicity of fluorophores simultaneously, permitting spatial mapping of cell type lineage markers, signaling proteins, drug targets, and immune cell antigens concurrently in the TME. Jerby-Arnon et al. (132) leveraged t-CyCIF in addition to scRNA-seq to define a relationship between T cell exclusion and therapy resistance in ICB-treated melanomas. Malignant cells expressing a resistance-associated gene program were often found to localize to T cell-depleted niches in the TME. As T cell exclusion has been recognized as a main player in the mixed patient responses to ICB, t-CyCIF can shed light on whether the spatial organization of a TME might support or hinder certain immunotherapeutic strategies. t-CyCIF has also aided studies of myeloid-targeted immunotherapies in breast cancer (133). Characterization of BRCA1-associated triple negative breast cancer (TNBC) tumors as highly infiltrated by immunosuppressive macrophages with t-CyCIF directly led the authors to explore therapies intended to reprogram these cells toward a more pro-inflammatory state in this TME. Integration with scRNA-seq-derived interaction networks would help to identify mechanisms by which macrophages orchestrate immunosuppression and provide auxiliary targets for immunotherapy.

Altogether, many strategies for spatial mapping of the TME exist, offering dynamic and static measurements of receptor-ligand and sender-receiver cell co-localization, in vivo and ex vivo. All have merit in the confirmation of physically possible interactions in scRNA-seq-derived communication networks, depending on the context-specific parameters of the TME being studied. Further, updating the aforementioned cell-cell communication inference methods to directly take spatial profiling information as input will enhance the integration of imaging with transcriptomics technologies. For example, as shown in their documentation , CellPhoneDB v.3.0 (63) now allows the user to input additional information describing the localization of each cell type to different spatial niches in the microenvironment, as derived from prior knowledge, SRT, or imaging. Their algorithm then places restrictions on predicted interactions between cell types in different niches.

Functional Signaling in the Tumor Microenvironment Can be Validated by Direct Measurement of Protein Secretion

mRNA expression does not provide a functional cellular measurement, nor does it generally correlate with protein levels (134), which fundamentally limits inference of protein signaling activity from scRNA-seq. Direct measurement of protein secretion is particularly critical when characterizing signaling in the tumor microenvironment, as many cytokine and chemokine networks modulated by immune cells are responsible for orchestrating pro- or anti-tumor effects. There are many commonly used platforms for protein measurement in the TME. Two of the most popular, fluorescence-activated cell sorting (FACS) and enzyme-linked immunosorbent assays (ELISA), maintain inherent limitations. Both are limited in their multiplexing capabilities, and FACS in particular measures surface markers and intracellular proteins as opposed to the secreted proteins of interest to signaling studies. Measuring protein secretion in the tumor microenvironment to complement scRNA-seq-derived communication networks requires immunophenotyping technology capable of higher throughput and multiplexing and that is tuned to capture secreted signals only.

The advent of microfabricated single-cell secretomic analysis devices has begun to clear some of the aforementioned hurdles (135–137). The assay platform presented in Lu et al. (138, 139) in particular has laid the foundation for investigating macrophage secretion in vitro and ex vivo. Their device, referred to as the single-cell barcode chip (SCBC), integrates subnanoliter microchambers with high-density antibody barcode microarrays for simultaneous detection of up to 15 cytokines from over a thousand single cells in parallel. Beyond fundamental studies of macrophage proinflammatory dynamics in vitro (140, 141), the SCBC has been leveraged to measure secretion directly ex vivo in tumor-associated myeloid cells and macrophages (i.e., TAMMs, which includes dendritic cells) isolated from 8-week-old murine melanoma tumors (98). Upon combinatorial treatment with two myeloid-targeted immunotherapies, CD40 agonist (CD40ag) and CSF1R blockade, a subset of TAMMs were found to upregulate their co-secretion of proinflammatory cytokines TNF, IL-6, IL-12, and chemokine CCL3, along with anti-inflammatory factor Chi3l3. Other TAMM subsets decreased their secretion of immunosuppressive proteins MMP9, CCL17, and CCL22 after treatment relative to control. Many of these cytokines and chemokines are known to modulate the behavior of other cell types in the TME through cell-cell communication. For example, IL-12 represents a key link between innate and adaptive immunity, as it has a potent effect on T cell activation and Th1-type responses through induction of IFNγ production (142). Further, defining cell subpopulations based on differences in functional signaling is a key component of disentangling biologically relevant interaction networks in the TME. Regarding scRNA-seq-derived cell-cell communication networks, integration with multiplexed measurement of cytokine and chemokine secretion, particularly of key inflammatory mediators like IL-12, from paired experimental conditions would help to validate hypothesized sender-population signaling at single-cell resolution.

While the SCBC has proven its utility in measuring secretion from tumor-associated myeloid cells (namely macrophages and DCs), an up-to-42-plex version of this device has been commercialized with IsoPlexis (www.isoplexis.com), and it has been primarily tuned to investigate T cell secretion. Commercialization has facilitated access to this technology by over 100 pharmaceutical companies and medical centers worldwide . In particular, this T cell-tuned SCBC has been used to examine the relationship between T cell polyfunctionality and resistance to anti-PD-1 therapy in human patients across multiple cancer types (143, 144). Polyfunctionality in this context was defined as co-secretion of two or more cytokines. Co-secretion of effector proteins Granzyme B and IFNγ by CD4+ and CD8+ T cell subsets, as measured by the SCBC, was identified as a distinguishing factor between response and resistance to anti-PD-1. The commercialized device has also been used to understand the functional response to other T cell-targeted treatments, including chimeric antigen receptor (CAR) T cell therapy (145, 146). Altogether, microwell assay platforms like the SCBC have advanced functional measurement of immune cell signaling in the TME, which can greatly aid the validation of functional cell subsets and cell-cell communication inferred from scRNA-seq.

While these microwell devices provided marked advances regarding multiplexed measurement of protein secretion by single cells from the TME, they maintain two shortcomings. As with scRNA-seq, the aforementioned assays offer static, “snapshot” measurements of protein secretion, although multiple experiments could be conducted over a series of timepoints. Further, platforms like the SCBC require cells to be sorted from the tumor prior to measurement of secretion, which could alter their functional behavior. Real time in vivo measurement of protein secretion would address both of these limitations, as well as provide a more dynamic picture of functional signaling for validation of scRNA-seq-derived cell-cell communication networks. As of yet, no real time in vivo methods have been demonstrated in the TME.

“Pruning” Cell-Cell Communication Networks Computationally Narrows the Search for Experimentally Testable Targets

Another significant challenge when trying to reconstruct cell-cell communication networks from scRNA-seq data is to manage the inherent scale and complexity that these networks can attain. One approach is to use the aforementioned experimental approaches to collect spatial and protein-level information in order to restrict scRNA-seq-derived connections to those that can be experimentally confirmed. However, this can limit the potential of these computational inference methods to discover novel interactions. Because interaction networks are directed (i.e., unidirectional along each individual axis of communication), weighted (i.e., each axis has an associated score indicating whether the data supports the occurrence of that interaction), and densely connected, data mining and biological interpretation are complicated. However, they also present an opportunity to use concepts from graph theory and mathematical modeling to identify the most central, influential, or biologically relevant signaling axes for experimental validation.

Graph theory is a field of mathematics that views networks as a collection of nodes, or vertices, connected by edges. In the context of scRNA-seq-derived interaction networks, cell types or subtypes would be considered nodes, with ligand-receptor or ligand-target gene axes comprising the edges between them. Centrality is a graph theory concept concerned with identifying “important” nodes in a graph, where the definition of importance changes depending on the “angle” at which the graph is analyzed. Measures of centrality include degree, betweenness, closeness, pagerank, and eigenvalues. Degree is the simplest measure of centrality, counting the absolute number of edges pointing toward (i.e., in-degree) and away from (i.e., out-degree) a node. Betweenness quantifies a node’s influence on the flow of information in the network, which can identify “bottlenecks” or nodes that serve as “ambassadors” between modules of signaling. Similarly, closeness indicates how directly or indirectly each node in a graph can communicate with one another. Pagerank is a measure of “popularity,” as it answers the question of how often a particular node is communicated through during a random walk on the graph. Finally, eigenvalue centrality is similar to degree centrality in that it quantifies a node’s influence through its connections with others. However, this metric is better suited to a weighted network, as it takes strength of communication and influence of other nodes into account.

While centrality metrics are often reserved for analyzing scRNA-seq-derived gene regulatory networks (3, 147), they are becoming a useful tool to evaluate cell-cell communication. Raredon et al. (74) leveraged degree and eigenvector centralities, in addition to Kleinberg hub and authority scores, to quantitatively compare scRNA-seq-derived interaction networks from adult mammalian lungs across species. The Kleinberg hub and authority scores (148) were developed as extensions of eigenvector centrality for directed networks, where a “hub” is a node that communicates to other important nodes, and an “authority” is a node that is communicated to by other important nodes. With these measures of node centrality, Raredon et al. were able to confirm conserved patterns in signaling topologies in lung tissue across species. For example, alveolar type I (ATI) cells were unexpectedly characterized as a major signaling hub in mouse, rat, pig, and human lung, particularly for growth factor-related signaling families including vascular endothelial growth factor (VEGF) and semaphorins (SEMA). This work will lay the foundation for recapitulating homeostatic cell-cell communication in engineered lung tissue. A similar approach has yet to be applied to scRNA-seq from the TME, although this analysis is nearly directly transferable. Instead of revealing conserved mechanisms of signaling, comparing network centralities across homeostatic and cancerous conditions could identify changes in signaling topologies that give rise to aberrant communication underlying tumor growth, metastasis, and immune evasion. User-friendly tools for the interpretation of scRNA-seq-derived interaction networks from the perspective of graph theory are starting to become more popular in the literature, including the aforementioned CellChat (55).

Outside of graph theory, regression modeling presents an opportunity to interpret scRNA-seq-derived cell-cell communication networks with respect to relevant quantitative and qualitative biological behaviors, including tumor outcome, tumor growth rate, and immune cell infiltration. Elastic net regularization (149) is a regression modeling method that estimates the relationship between one or more features, or predictors, and a dependent outcome of interest, while accounting for overfitting and the instability imposed by high-dimensional data with few examples. Features are assigned coefficients representative of their relevance with respect to the outcome of interest, with coefficients for less predictive features shrinking to zero. Elastic net regularization has been previously applied to scRNA-seq to generate predictive gene signatures of immune cell identity (150), with the features in the regression problem being each gene in the transcriptome. Similarly, with features as ligand-receptor or ligand-target gene axis edges instead of genes, elastic net regularization could be used to select signatures of signaling predictive of aberrant or immunotherapy-altered tumor behaviors. This analysis would further filter network edges considered irrelevant in the context of the biological outcome of interest. Similar regression modeling approaches like Partial Least Squares Discriminant Analysis (PLS-DA) (151, 152) could further disentangle the directionality of the relationship between predictive signaling axes and emergent behaviors in the tumor. Overall, computational “pruning” strategies from graph theory and mathematical modeling are underappreciated and under-utilized resources for the distillation and interpretation of hypothesized interaction networks of the TME.

Conclusion

scRNA-seq has proven to be a powerful tool for high throughput prediction of cell-cell communication. Transcriptome-wide profiling of genes encoding for ligands, receptors, intermediary signaling proteins, and signaling targets, in conjunction with novel algorithms like CellPhoneDB, CellChat, NicheNet, and CytoTalk, has characterized signaling mechanisms underlying homeostasis and disease across biological systems. In particular, scRNA-seq-derived interaction networks have facilitated identification of immunosuppressive and immunotherapy-altered cell-cell signaling in the TME, which contribute to the design and development of more effective immunotherapeutic interventions in the clinic. Validation of these networks will only serve to strengthen the hypotheses that they generate. Integration with SRT and other imaging technologies will help to define which interactions are physically possible in the microenvironment through co-localization of both communicating cell types and paired ligand, receptor, and signaling intermediaries. Functional measurement of protein secretion will confirm outgoing communication from “sender” cells and organize broader cell populations into relevant subsets based on patterns in secreted signaling. Further computational analysis of interaction networks by borrowing from the fields of mathematics and statistics will focus the search for experimentally targetable signaling mechanisms underlying aberrant behaviors. Together, integration of scRNA-seq-derived communication networks with powerful experimental and computational tools will expand our understanding of the role that cell-cell communication plays in emergent biological behaviors in the TME and beyond, improving outcomes for patients across malignancies.

Author Contributions

KB and KM-J contributed equally to all aspects of this article. All authors read and approved the final manuscript.

Funding

This work was supported by the U.S. National Institutes of Health (U01CA238728 to KM-J).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.