# SingleCellSignalR: inference of intercellular networks from single-cell transcriptomics

Simon Cabello-Aguilar, Mélissa Alame, Fabien Kon-Sun-Tack, Caroline Fau, Matthieu Lacroix, Jacques Colinge

► To cite this version:

Simon Cabello-Aguilar, Mélissa Alame, Fabien Kon-Sun-Tack, Caroline Fau, Matthieu Lacroix, et al. SingleCellSignalR: inference of intercellular networks from single-cell transcriptomics. Nucleic Acids Research, 2020, 48 (10), pp.e55-e55. 10.1093/nar/gkaa183. hal-04126980

HAL Id: hal-04126980

https://hal.science/hal-04126980

Submitted on 13 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# SingleCellSignalR: inference of intercellular networks from single-cell transcriptomics

Simon Cabello-Aguilar<sup>1,2,3</sup>, Mélissa Alame<sup>1,2,3,4</sup>, Fabien Kon-Sun-Tack<sup>1,2,3</sup>,
Caroline Fau<sup>1,2,3</sup>, Matthieu Lacroix<sup>1,2,3</sup> and Jacques Colinge<sup>1,2,3,*</sup>

Received December 12, 2019; Revised January 23, 2020; Editorial Decision March 06, 2020; Accepted March 10, 2020

## ABSTRACT

Single-cell transcriptomics offers unprecedented opportunities to infer the ligand-receptor (LR) interactions underlying cellular networks. We introduce a new, curated LR database and a novel regularized score to perform such inferences. For the first time, we try to assess the confidence in predicted LR interactions and show that our regularized score outperforms other scoring schemes while controlling false positives. SingleCellSignalR is implemented as an open-access R package accessible to entry-level users and available from https://github.com/SCA-IRCM. Analysis results come in a variety of tabular and graphical formats. For instance, we provide a unique network view integrating all the intercellular interactions, and a function relating receptors to expressed intracellular pathways. A detailed comparison of related tools is conducted. Among various examples, we demonstrate SingleCellSignalR on mouse epidermis data and discover an oriented communication structure from external to basal layers.

## INTRODUCTION

In multicellular organisms, cells engage in a large number of interactions with adjacent or remote partners. They do so to coordinate their fate and behavior from early developmental stages to mature tissues (1–4), in healthy and diseased (5) conditions. Although other mechanisms may play a role such as tunneling nanotubes, secreted vesicles or ion fluxes, a significant part of cellular communications is carried over by secreted ligand and cell surface receptor physical interactions (6). In the particular case of tumors, cancer cells can reprogram their micro-environment through secreted factors, turning neutral or anti-tumor cells into tumor supportive elements (7,8). The emergence of single-cell RNA sequencing (scRNA-seq) technologies (9–11) has pro-

vided researchers with powerful means of learning which cells compose specific tissues. The different cell populations present in a sample can be determined by applying unsupervised clustering (12). Further tools exist to infer intracellular pathway activity (13–15), i.e. internal cell states. While such compositional descriptions are essential, deciphering individual cell contributions in tissues requires the unraveling of cellular interactions. Recent scRNA-seq-based studies illustrated how ligand-receptor (LR) interaction mapping might provide a better insight in tissue development and homeostasis or tumor biology. For instance, Puram et al. (16) studied head and neck squamous cell carcinomas. They were able to identify an LR interaction, TGFB3-TGFBR2, involved in the communication between cancer cells undergoing (partial) epithelial to mesenchymal transition, at the leading edge of the tumor and cancer-associated fibroblasts. The latter interaction was shown to be necessary for invasiveness. Additional examples illustrating the discovery potential of LR maps can be found in the literature (17–21). Those results highlight the need for a systems biology tool that would assist investigators in portraying cellular networks by inferring confident putative LR interactions for follow-up validation.

SingleCellSignalR is the first such tool available in R. It relies on a comprehensive database of known LR interactions, which we called LRdb. It also introduces a new regularized product score aimed at adapting to variable levels of depth in single-cell datasets, i.e. the prevalence of censored read counts or dropouts. LRdb is the result of integrating and curating existing sources plus manual additions; to the best of our knowledge, it is the largest database of this kind. The new scoring approach has the advantage of facilitating the use of stable thresholds on LR interaction scores to control false positives (FPs) and not only rank LR interactions. SingleSignalR can start from raw read count matrices and use integrated data normalization, clustering and cell-type calling solutions before inferring LR interactions between cell populations, i.e. cell clusters. Alternatively, those preliminary steps can be substituted by any other tools or

<footnote n="1"><sup>1</sup>Institut de Recherche en Cancérologie de Montpellier, Inserm, F-34298 Montpellier, France, <sup>2</sup>Institut régional du Cancer Montpellier, F-34298 Montpellier, France, <sup>3</sup>Université de Montpellier, F-34090 Montpellier, France and</footnote>

<footnote n="4"><sup>4</sup>Département d'Hématologie Biologique, CHU Montpellier, Hôpital Saint Eloi, F-34090 Montpellier, France</footnote>

<footnote>*To whom correspondence should be addressed. Tel: +33 467 612 392; Fax: +33 467 613 787; Email: jacques.colinge@inserm.fr</footnote>
framework, and SingleSignalR used for LR interaction inference only, its primary purpose. In order to facilitate the interpretation of LR interactions and put them in context, a range of visualization and complementary analysis tools are provided, e.g. the inference of intracellular networks rooted at the receptors—or ligands—expressed by a particular cell type. LRdb contains human genes, but we can accommodate murine datasets by translating murine genes to their human orthologs. An application example illustrates this functionality on mouse skin data.

Several authors proposed tools for mapping LR interactions (18,22–26) from scRNA-seq data. They exploited different reference lists of potential LR interactions and LR interaction scoring schemes to rank those interactions. A systematic comparison of scoring sensitivity and selectivity, and tool features is presented.

## MATERIALS AND METHODS

### LRdb—a curated database of ligand–receptor interactions

We decided to favor interpretable inferences that are supported by the literature or experimental data. Accordingly, we compiled the content of existing databases that contain such LR pairs: FANTOM5 (6), HPRD (27), HPMR (28), the IUPHAR/BPS Guide to Pharmacology (29) and UniProtKB/Swissprot (30) annotations related to families of ligands or receptors covered in the previous databases. We also extracted LR pairs, with respective participants annotated as ligand or receptor in GO (31), from Reactome (32) pathways. We required GO Cellular Compartment (GOCC) annotation ‘receptor complex’ (GO:0043235) for receptors, and ‘extracellular space’ (GO:0005615) or ‘extracellular region’ (GO:0005576) for ligands. Reactome pathways were downloaded from Pathway Commons (33). That yielded 3191 pairs. Inspection of the latter revealed 106 dubious cases, i.e. ligand–ligand or receptor–receptor, or non-secreted gene products. After elimination, we obtained 3085 reliable LR pairs, which we extended with 166 pairs extracted from cellsignaling.com maps and related literature manually to reach 3251 LR pairs (Figure 1A and B). All the gene symbols were converted into the latest HUGO version as part of the curation process.

A proper comparison of SingleCellSignalR with existing tools is presented in ‘Results and Discussion’ section. At this stage, it is nonetheless interesting to note that some tools (23,24) relied on LR pairs deduced from STRING with principles analogous to what we did with Reactome. That yielded a large number of reference LR pairs, >30 000 typically, but the overlap with LRdb extensive collection of known LR interactions remained modest (50.2%) (Figure 1C). Although that approach might provide additional discovery potential, we decided to stay with literature-derived LR interactions.

### SingleCellSignalR implementation

The software package was implemented in R, following Bioconductor standards. Basic usage examples are provided in Supplementary Material and the package documentation. In principle, we encourage users to apply advanced data normalization, clustering and cell-type calling tools

(12,34), and to submit preprocessed count matrices to SingleCellSignalR to perform LR interaction inference and visualization. For convenience, we implemented simple data preparation steps enabling users to start from a raw read count matrix. We normalize individual cell transcriptomes according to their 99th read count percentile. In case a cell has its 99th percentile equal to zero, it is discarded. Normalized read counts (+1 to avoid zeros) are log-transformed. Clustering can be obtained either by chaining principal component analysis for dimension reduction and K-means (35), or by using the advanced SIMLR model (36). Cell-type calling is implemented based on a list of gene signatures following a format identical to PanglaoDB (37) exports such that users can easily add cell types from this rich source or provide their own. Our algorithm computes the average gene signature expression across all the cells and for all the signatures to obtain a signature × cell matrix. This matrix is normalized and a threshold is iteratively adjusted to maximize the number of cells assigned to a single cell type. Full details and an example are in Supplementary Material.

It is possible to infer paracrine or autocrine only interactions, or both types (Supplementary Figure S1); details in Supplementary Material. For annotation purposes, differentially expressed genes between each cluster and all the other clusters pooled, are successively searched with edgeR functions glmFit and glmRT (38). LR interactions with both the ligand and the receptor significantly enriched in their respective cell types are labeled ‘specific’ (Supplementary Figure S1).

In order to relate receptors to intracellular signaling, we make use of Reactome and KEGG (39) interactions downloaded from Pathway Commons. Interactions are assigned to several types that we simplified to facilitate the display of networks afterward. Interaction types ‘interactwith’ and ‘in-complex-with’ were assigned to the simplified type ‘complex.’ The interaction types ‘chemical-affects,’ ‘consumption-controlled-by,’ ‘controls-expression-of,’ ‘controls-phosphorylation-of,’ ‘controls-production-of,’ ‘controls-state-change-of,’ ‘controls-transport-of’ and ‘controls-transport-of-chemical’ were simplified as ‘control.’ The interaction types ‘catalysis-precedes,’ ‘reacts-with’ and ‘used-to-produce’ were simplified as ‘reaction.’ The simplified type ‘control’ was considered directional whereas ‘complex,’ and ‘reaction’ were considered undirected.

### Datasets

Different single-cell datasets were used to illustrate and benchmark SingleCellSignalR. MELANOMA is a metastatic melanoma dataset covering several patients (40), SMART-Seq2 protocol, GEO GSE72056. 10×PBMC is a peripheral blood monoclonal cell dataset from healthy donor (41), Chromium 2 protocol. 10×T is a pan T-cell dataset from healthy donor (42), Chromium 2 protocol. 10×PBMC and 10×T were downloaded from 10× Genetics web site. HNSCC is a head and neck squamous cell carcinoma (primary and metastatic) dataset (16), SMART-Seq2 protocol, GEO GSE103322. PBMC is a second, deeper peripheral blood monoclonal cell dataset (43), SCRB-seq protocol, GEO GSE103568.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** SingleCellSignalR databases and workflow. (A) LRdb sources. (B) Approximate overlap of sources. (C) Overlap (50.2%) of LRdb with an example of LRdb derived from STRING (here the database of scTensor). (D) SingleCellSignalR general workflow with input transcript expression matrix both normalized and clustered by either independent tools or by SingleCellSignalR basic built-in algorithms.</figcaption></figure>

## Mouse skin immunolabeling

Fresh adult mouse skin samples were embedded in OCT (Sakura) and cryosectioned. Immunolabeling was performed on unfixed 10 µm cryosections using the following antibodies: anti-PSEN1 (SAB4502423, Sigma) and CD44 (MA1-10225, ThermoFisher). For IF, we used anti-rat conjugated 488 (A-11006, ThermoFisher) and anti-rabbit conjugated 594 (A-21207, ThermoFisher) antibodies.

## RESULTS AND DISCUSSION

### Workflow overview

Independent of the chosen scRNA-seq platform, data come as a table of read or unique molecule identifier (UMI) counts, one column per individual cell and one row per gene. The prediction of LR interactions between cells requires that scRNA-seq data are normalized and clustered, with each cluster corresponding to a cell type (12). Ideally, the cell types should be called, e.g. using gene signatures, to facilitate the interpretation of the LR interactions. The design of SingleCellSignalR is such that those preliminary steps can be accomplished within the package, using built-in or integrated solutions or realized with other tools according to the user preference. The latter option is preferable to benefit from the latest or most advanced techniques. Once data have been prepared, the inference of LR interactions is performed by successively considering each possible

cell type—or cluster—couple, e.g. CD8+ T cells versus macrophages (Figure 1D). The reality of each potential LR pair according to LRdb is estimated by the computation of a score based on gene expression in the respective cell types. It is also possible to infer autocrine interactions. The output interaction lists are provided in various formats (tables, different plots and networks), and complementary functions were designed to help to interpret these lists, e.g. linking receptors to intracellular signaling pathways from Reactome and KEGG.

### Scoring ligand-receptor pairs

To infer LR interactions between cell types A and B, we interrogate LRdb and score each LR pair found with average ligand expression $l > 0$ in A, and average receptor expression $r > 0$ in B or vice versa. Existing tools often consider the product $lr$ and sometimes the average $\frac{l+r}{2}$. Alone, such scores can rank candidate LR pairs, but users are left without any clue where to cut off likely FPs. The computation of a score should hence be accompanied by a procedure to determine a threshold below which scores are deemed unreliable. A common choice is to shuffle cell-type assignments multiple times and to obtain a score null distribution to estimate score P-values. Although intuitive and statistically sound, this solution does not address the real question. It rather identifies LR pairs that are significantly specific to a
given couple of cell types. That should imply that the inter-action is likely to exist but it might result in poor P-values for real LR interactions that are shared by many couples of cell types. For instance, in tumors, it is common that many immune cell populations express immune checkpoints and their ligands, e.g. PD-1/PD-L1, forming an immunosup-pressive micro-environment. In such a case, many combi-nations of cell types would involve the PD-1/PD-L1 inter-action resulting in bad P-values if an insufficient number of other cell populations, not expressing these molecules, would be present in the dataset. To address the points above, we introduced a regularized product score:

$$LRscore = \frac{\sqrt{lr}}{\mu + \sqrt{lr}},$$

where $\mu$ = mean(C) and C is the normalized read count matrix. In this empirical score, the mean $\mu$ acts as a scaling factor, and the square roots are meant to keep the $lr$ products and $\mu$ on the same scale. The LRscore is bounded by 0 and 1, independent of the dataset depth.

Now, to define a score threshold is not a simple task since, by essence, all LR pairs listed in LRdb are correct in a particular context, i.e. defining false LR interactions must rely on external knowledge specific to the couple of cell types considered each time. Ideally, single-cell transcrip-tomes would be available for a tissue where all possible LR pairs would have been tested experimentally, e.g. by immunofluorescence, eventually including functional evi-dence. Such data are obviously not available.

To estimate an appropriate LRscore threshold and to compare LRscore to other scoring schemes nonetheless, we decided to use an *ad hoc* benchmark that would be unbi-ased regarding the scoring schemes and accurate enough bi-ologically. A first opportunity was provided by Ramilowski *et al.* (6) data, which are part of FANTOM5 and included a table reporting the expression (in TPM) of many ligands and receptors over 144 primary cell types. This table, which we denote $T_{ref}$, was obtained by sequencing techniques that were much deeper than scRNA-seq, i.e. virtually devoid of dropout for our purpose. The authors considered that ex-pression above 10 TPM could be taken as a conservative gene expression basal limit. We thus designed an evaluation where candidate LRdb pairs, restricted to ligands, recep-tors and cell types covered in $T_{ref}$, could be considered cor-rect if both the ligand and the receptor TPM were above 10. If both were below 10 TPM, the pair was deemed false, and mixed cases (one >10 and the other $\le 10$) were ignored. Ad-mittedly, ligand and receptor concomitant expression does not guarantee a functional interaction, but since those lig-ands and receptors are known to interact for at least one combination of cell types, we considered the benchmark above sufficiently accurate for our purpose. A second op-portunity was provided by a proteomics study of 28 primary human hematopoietic cell populations in steady and acti-vated states (44), which enabled a similar benchmark relat-ing censored scRNA-seq data to their proteomics counterpart. Namely, each cell population was measured in quadruplicate and we stored the average peptide counts in a table $P_{ref}$ (similar to $T_{ref}$). This second evaluation was conducted

as above with a threshold of an average spectral counts $\ge 2$ for an LR pair to be true.

We used ROC curves to compare methods and determine score thresholds at 5% FPs. We already mentioned LR pair scores equal to the product $lr$ or the mean $\frac{l+r}{2}$ as well as their P-values. A further proposed scheme was to select in each cell type, those genes whose expression was specific. Subsequently, whether the reference database would con-tain pairs of such genes was checked (22,24). In the case of Zhou *et al.* (22), gene selection relied on a z-score calcula-tion. To obtain a ROC curve, we varied the coefficient ap-plied to the standard deviation at that gene selection stage. Another selection-based tool (24) did not report any score. It is thus not included here but covered in the tool compar-ison below only. scTensor (23) applied non-negative Tucker decomposition to model ligand and receptor read counts as sums over the many LR pairs they might contribute to. Despite considerable efforts, we could not use it beyond its pre-packaged example. It is therefore absent from the pre-sented comparison.

We applied the two benchmarks above to five datasets covering several cell types in $T_{ref}$ and in $P_{ref}$ ('Materials and Methods' section). Every couple of cell types, both in $T_{ref}$ or $P_{ref}$, yielded ROC curves such as the example fea-tured in Figure 2A for 10×PBMC data against $T_{ref}$. Areas under the curves (AUCs) over all the datasets and all the cell-type couples represented in $T_{ref}$ are presented in Figure 2B; individual curves and dataset-specific AUCs are in Sup-plementary Figures S2–6. We see that the best AUCs were achieved by LRscore, the product, and the average. As ex-pected, P-value calculations yielded inferior performance. Zhou *et al.* selection mechanism was also inferior to the best scores. We next asked how variable a threshold aimed at cutting at a certain FP rate would be; we imposed 5% FPs (Figure 2C and Supplementary Figures S2–6 for dataset-specific plots). Naturally, due to its design, LRscore thresh-olds turned out to be much more stable. In Figure 2D, we show how a threshold common to all the single-cell datasets could be determined, imposing <5% FPs in 75% of all the ROC curves of all the five datasets. We found LRscore > 0.4. Repeating the analysis for the proteomics reference $P_{ref}$ gave the results in Figure 2E and F (and Supplementary Figures S7–11) that are very similar. The computation of a common threshold to achieve <5% FPs in 75% of the ROC curves resulted in LRscore > 0.6 (Supplementary Fig-ure S12). This more stringent threshold can be explained by a potential lower sensitivity (~10 000 proteins) compared to FANTOM5 transcriptomics and by the obvious differ-ences between transcriptome and proteome. These consid-erations indicate that a universal threshold that would guar-antee, e.g. 5% FPs cannot be determined, but score regular-ization is already a major step toward more rigorous cut-offs. To illustrate this, we imposed LRscore > 0.5 to all the datasets against both $T_{ref}$ and $P_{ref}$. From Figure 2G, we see that we managed to maintain the FPs in a reasonable range, and, most importantly, the variability of the FP rate within each dataset is modest. That is, a threshold can be set and FP rates on LR interactions between any combinations of two cell types remain comparable. Repeating this analysis with the product score and common threshold $10^{-2}$ found as above, we note large discrepancies between datasets.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<table><caption>**Figure 2.** Statistical analysis. (A) Representative ROC curve. (B) AUCs over all the ROC curves of the five datasets; transcriptomics reference. (C) Relative variability (with respect to the median value) of the thresholds required to achieve 5% FPs in each ROC curve (all couples of cell types, all the datasets). (D) FP rate upon various LRscore thresholds. A small box plot is featured for each threshold value and the figure inset shows all the ROC curves. LRscore threshold (blue dashed line) such that 75% of the ROC curves would yield FPs below 5%. (E and F) Same as B and C, but for the proteomics reference. (G) FP rates (±sd) on each dataset against the transcriptomic and the proteomic references when imposing a consensus LRscore threshold of 0.5. Same results for an equivalent product score threshold consensus in the two right-most columns.</caption><thead><tr><th rowspan="2">Dataset</th><th colspan="2">LRscore</th><th colspan="2">Product score</th></tr><tr><th>FPR for threshold = 0,5</th><th>+/-</th><th>FPR for threshold = 10<sup>-2</sup></th><th>+/-</th></tr></thead><tbody><tr><td>Transcriptomics</td><td></td><td></td><td></td><td></td></tr><tr><td>10xT</td><td>0,0%</td><td>0,0%</td><td>0,0%</td><td>0,0%</td></tr><tr><td>10xPMMC</td><td>0,5%</td><td>1,0%</td><td>0,0%</td><td>0,0%</td></tr><tr><td>PMMC</td><td>0,5%</td><td>1,0%</td><td>30,8%</td><td>9,0%</td></tr><tr><td>HNSCC</td><td>2,5%</td><td>1,6%</td><td>14,6%</td><td>4,4%</td></tr><tr><td>MELANOMA</td><td>3,7%</td><td>3,0%</td><td>20,8%</td><td>7,7%</td></tr><tr><td>Proteomics</td><td></td><td></td><td></td><td></td></tr><tr><td>10xPMMC</td><td>4,3%</td><td>0,7%</td><td>1,8%</td><td>0,9%</td></tr><tr><td>10xT</td><td>5,7%</td><td>1,2%</td><td>0,5%</td><td>0,4%</td></tr><tr><td>PMMC</td><td>10,3%</td><td>2,5%</td><td>65,9%</td><td>8,9%</td></tr><tr><td>HNSCC</td><td>9,4%</td><td>1,6%</td><td>43,8%</td><td>6,9%</td></tr><tr><td>MELANOMA</td><td>8,4%</td><td>1,7%</td><td>36,7%</td><td>4,0%</td></tr></tbody></table>
Representing ligand–receptor interactions

After scoring and application of a threshold, Single-
CellSignalR outputs LR interactions in various formats.
10×PBMC raw UMI counts were submitted to our de-
fault pipeline, which clustered cells in six populations: B
cells, T cells, regulatory T cells (Tregs), neutrophils, cyto-
toxic cells and macrophages (Figure 3A). A summary chord
diagram can be generated that indicate the number of LR
interactions between each cell population couple (Figure
3B). Chord diagrams for individual interactions between
two populations are possible as well (Figure 3C). Inspection
of the expression of the ligand and the receptor involved in
a LR interaction is implemented in mixed or dual 2D pro-
jections (Figure 3D–E).

Chord diagrams are convenient for looking at specific
couples of cell types but not to account for the whole
dataset. We hence implemented two integrated views as ei-
ther a tabular plot (Figure 3F) or a network exported in
graphML format. The latter can be imported in tools such
as Cytoscape (45). We exemplify the network view with
MELANOMA patient 89 data. In Figure 3G, we picture
an overview of the intercellular network (full details in Sup-
plementary Figure S13). The complexity and partial re-
dundancy of this network are particularly noteworthy, with
multiple cell types expressing the same ligands or receptors.
That is typical of immunosuppressive TMEs.

Putting ligand–receptor interactions in context

To help to understand the functional consequences of LR
interactions, we relate the receptors in each cell type with
downstream biological pathways (Figure 3H). This requires
a global reference of all possible pathways, including links
from the receptors and a method to generate receptor-
related, intracellular regulatory networks that are restricted
to each cell type. We decided to use KEGG (39) and Re-
actome (32) pathways taken from Pathway Commons (33)
as a global reference. This choice guaranteed interpretable,
literature-supported network interactions. Surprisingly, we
found 100 receptors in LRdb that were devoid of down-
stream interactions in this pathway ensemble. By manually
checking their respective UniProtKB/Swissprot entries, we
found 1–12 known interactions for 61 of them, which were
added along with literature references (176 interactions in
total).

Construction of a ligand-related internal network for
each cell type was obtained by the following algorithm.
Given a list of receptors, e.g. all the receptors involved in
LR interactions for the cell type at hand or a subset, we
first identify all the pathways including these receptors. The
union of all such pathways is intersected with the set of
genes expressed by the cell type. That gene list is obtained
by selecting the genes expressed by more than a propor-
tion *p* of the individual cells after normalization of the raw
read count matrix, that is matrix *C* (*p* default is 20%; it can
be adjusted by the user). To the intersected pathways, we
add direct receptor/expressed gene interaction that would
not be included in Reactome or KEGG, which are those
we manually added. The functionality can be applied to lig-
ands as well to obtain information about their pathways of

origin. In every case, the resulting network edges are anno-
tated with Reactome/KEGG data plus a simplified interac-
tion typology that is convenient for visualization ('Materials
and Methods' section).

The neutral and straightforward procedure above for call-
ing expressed genes was retained to be independent of the
cell population size, on the one hand, and to avoid interfer-
ence with raw read count preprocessing, on the other hand.
In particular, users could apply imputation strategies (46) to
reduce dropouts and/or advanced data normalization and
regularization (12,47).

Mouse interfollicular epidermis application

Mouse interfollicular epidermis (IFE) is a multilayered ep-
ithelium in which proliferating cells reside in the basal layer
(IFE B), where they undergo regulated cell division. Their
daughter cells move upward into the suprabasal layers while
further differentiating until they reach the outermost layer
(IFE K2) (Figure 4A). In this application, we demonstrate
SingleCellSignalR mouse data functionality, which is im-
plemented by internally mapping mouse genes to their hu-
man orthologs according to Ensembl (48) to exploit LRdb.
Murine gene names are preserved in the different outputs.

From Joost et al. data (49), we selected IFE cells exclu-
sively (658 cells, Figure 4B and C). To confront these data
with our FP analysis above, we considered LR interactions
between the keratinized layer (IFE K2 and K1 cells pooled)
and the suprabasal layer (IFE D1 and D2 cells pooled).
The pooling was motivated by minor differences between
cells and a facilitated comparison with human epidermis
(see below). We found 248 LR interactions with LRscore >
0.5 (Supplementary Table S1). Among the top-ranked inter-
actions, several cases involved Presenilin-1 (Psen1). Based
on novelty and specific antibody availability, we decided
to validate in mouse epidermis the interactions Presenilin-
1/CD44, which is likely involved in cell differentiation and
tissue organization. From Figure 4D, we observe correct lo-
calization of those two proteins in the IFE layers (left and
middle pictures), Presenilin-1 being both cytoplasmic and
secreted in agreement with the literature. In the merged pic-
ture (right), we see an overlap between the diffusive, extra-
cellular Presenilin-1 signal and suprabasal cell plasma mem-
branes harboring CD44 in some regions, e.g. at the two loca-
tions indicated by white asterisks. Using the Human Protein
Atlas (HPA) (50) as a systematic resource to search human
orthologs expression in the corresponding epidermis layers,
literature (51–53) and the above experimental validation, we
confirmed 158 out of the 176 inferred LR interactions for
which we had data (Supplementary Table S1). That is 10.2%
FPs, in line with our ROC curve analyses above compared
to the proteomics reference.

Top layer IFE K2 cells supposedly maintain elementary
activities only, which was supported by the limited num-
ber of expressed genes found in those cells (Figure 4E).
The number of expressed genes coding for ligands or recep-
tors was directly correlated with the total number of genes
in Figure 4E, reflecting the proportion they occupy in the
mouse genome. This configuration changed considering lig-
ands and receptors that were involved in reliable LR inter-
actions only (LRScore > 0.5). The largest numbers of lig-
<figure><img src="image_9.png" /><figcaption>**Figure 3.** Graphical representations. (A) 10×PBMC data with SIMLR clusters. (B) Summary chord diagram of the paracrine interactions; largest number of interactions from regulatory T cells toward macrophages and neutrophils. (C) Paracrine interactions from neutrophils toward cytotoxic cells. (D and E) Joined, and separated expression plots over the t-SNE map to assess LR interaction specificity and prevalence. (F) Integrated tabular view of the N most variable LR pairs with LRscore > 0.5 in one cell-type couple at least. (G) Integrated network of MELANOMA data patient 89 intercellular interactions. Overview and chosen interactions. The full network is in Supplementary Figure S13. (H) Example of intracellular signaling downstream CTLA-4 in T cells. Node sizes represent the gene expression level (arbitrary scale).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 4.** Mouse IFE. (A) Schematics of the IFE (arrow = axis of differentiation). (B) t-SNE plot of the IFE cells with the underlying black arrow representing the axis of differentiation. Some incompletely differentiated non IFE B cells lie among the IFE B subpopulation. (C) Cell type calling with actual subpopulations in the top color band. (D) Presenilin-1 (left panel) and CD44 (middle panel) immunostainings of mouse skin sections. Sections were counterstained with DAPI (epi = epidermis, \* = co-localization). (E) Number of genes with average expression different from 0 in each IFE layer. (F) Number of ligands and receptors involved in LRscore > 0.5 interactions. (G) Flow diagram representing the number of interactions between the cell subpopulations. (H and I) Chord diagrams of adjacent and remote IFE K2 cell–cell interactions. (J) Immune marker gene expression. Il18 is expressed by macrophages, Ifi30 by antigen-presenting cells and Prdm1 by T cells (genecards.org).</figcaption></figure>
ands were found in upper layers and the receptors followed a reversed pattern (Figure 4F). The overall LR interactions are quantified in Figure 4G, which shows an IFE K2 layer sending most signals and receiving the least. The number of LR interactions with a ligand contribution diminish from the outer layer to the most basal one. That is, the application of SingleCellSignalR unraveled upper IFE layer cells that produce much information but receive little. The basal layers behaved anti-symmetrically. LR interactions between IFE K2 and K1 and more remotely between IFE K2 and D1 are reported in Figure 4H and I. Pathway analysis showed that the interactions issued from IFE K2 cells were mostly involved in growth (e.g. GAB1 signalosome, or signaling by EGFR) and differentiation (EPH-ephrin mediated repulsion of cells, or Ephrin signaling). That illustrates the notion of a tight long-distance connection between outer and basal layers, which might be of particular interest in understanding the regulation of the epidermis constant self-renewal.

Interestingly, we noticed that LR pairs involved in wound healing regulation (54) were expressed by IFE K2 and D1 (Figure 4I). That suggests a mechanism for rapid reaction upon wounding, bringing K2 and D1 cells in contact. Although we excluded the immune cells present in Joost et al. data, we found 146 genes associated with immune cells. Selected genes with a clear pattern, as well as the average of all immune-related genes, are reported in Figure 4J. IL18 and PRDM1 were found with epidermis keratinocyte expression in human skin according to HPA (Supplementary Figure S14). This striking observation strengthens the potential for an immune function of the epidermis (55,56) and may have essential implications in inflammatory skin diseases such as psoriasis or atopic dermatitis.

## Comparison with other tools

The general characteristics of SingleCellSignalR and related tools appear in Table 1. iTALK (25) infers LR pairs from a compilation of public LR databases using a product score limited to the most abundantly expressed genes. It requires a read count matrix, where cell type calling was performed beforehand. iTALK includes a feature to deal with multiple datasets (time course, different conditions) such that variability and trends in ligand and receptor expression across datasets can be added to the LR plots. CellPhoneDB (26) is both an online tool and a Python package that can be downloaded. It scores LR pairs after P-values of the mean score. A modification is operated in case of a multimeric receptor and/or ligand, requiring that all the subunits are expressed. This unique and biologically sound feature might cause false negatives due to dropouts in scRNA-seq data. Kumar et al. (18) employed a product score. Instead of simple data shuffling, they applied a statistical test (Wilcoxon) to assess that over several tumors, the score was significantly different from zero. That procedure cannot be reproduced when analyzing a single dataset. Zhou et al. (22) proposed to select LR pairs based on specific expression of the ligand and the receptor in two cell subpopulations, referring to a database of known LR pairs. Specific gene expression was tested, requiring that the average expression of a gene, in a subpopulation, was above its mean expression over the whole count matrix plus three standard deviations. scTen-

sor (23) aims at inferring LR pairs through non-negative Tucker decomposition, explaining the observed ligand and receptor read counts as sums of the contributions of all the interactions they would engage in. It exploits a database of potential LR pairs generated from STRING interactions and Swissprot annotations (secreted/membrane) automatically, which yields a considerable number of putative LR pairs and does not cover many known interactions (Figure 1C). One attractive scTensor feature is the LR reference available for multiple organisms. PyMINEr (24) is a Python program that can perform k-means plus clustering based on log-transformed and normalized data. Cell-type calling is based on subpopulation specific gene pathway enrichment; application of characteristic gene signatures is not available. By contrast, SingleCellSignalR relies on gene signatures to perform cell-type calling (see Supplementary Methods for details and Supplementary Figure S15 and Table S7 for an illustration). It is possible to prepare and cluster data externally with PyMINEr. LR pairs between two cell subpopulations are inferred from pairs of genes found to be specific to each subpopulation, in interaction in STRING (57), and respectively classified as receptor and ligand in GO. There is no scoring *per se*, only this selection process. PyMINEr further infers a so-called co-expression graph, linking genes correlated over the whole read count matrix. This graph can be exploited to project different data such as the expression of genes in one specific subpopulation. Like scTensor, PyMINEr relies on automatically generated LR references from STRING.

As explained above, we could not use scTensor beyond its prepackaged example dataset (58) despite considerable efforts. We nonetheless tried to obtain some elements of comparison by processing the same data with SingleCellSignalR. We found two reliable LR pairs whereas scTensor returned 14 pairs, none in common, see Supplementary Table S2. Not being able to test scTensor on more examples did not allow us to perform a real comparison. We analyzed 10×PBMC data with PyMINEr and observed a substantial discrepancy with our results. Without imposing any threshold on PyMINEr output (no guidance provided by the authors), PyMINEr typically retrieves 10–20 times more paracrine interactions, and twice the number of pooled paracrine and autocrine interactions, with little interactions with our LRscore > 0.5 results (Supplementary Tables S3 and 4). The ROC curve analysis we conducted and such enormous lists of interactions inferred by PyMINEr suggest contamination with FPs. Obviously, PyMINEr output should be imposed a threshold in practice, which would reduce the FPs but would not improve the overlap with the known interactions in LRdb. We next submitted the same 10×PBMC data to CellPhoneDB website. The smaller LR reference database, strict rules regarding multimeric ligands or receptors, and the impact of a re-shuffling-based P-value calculation caused reduced sensitivity compared to our solution by a factor 8.5 on average (Supplementary Table S5). Inspection of CellPhoneDB output showed that the interactions we missed were absent from LRdb and imported in CellPhoneDB from sources that we did not consider, such as InnateDB, MINT, IntAct, MatrixDB or manual curation. At last, we installed iTALK R package and applied it to 10×PBMC data as downloaded from iTALK GitHub
<table><caption>**Table 1.** Software tool comparison</caption><thead><tr><th></th><th>SingleCellSignalR</th><th>iTALK</th><th>PyMINER</th><th>CellPhoneDB</th><th>scTensor</th><th>Zhou et al.<sup>a</sup></th><th>Kumar et al.</th></tr></thead><tbody><tr><td>LR database size</td><td>3251</td><td>2648</td><td>65 910<sup>b</sup></td><td>1396</td><td>34 449</td><td>2558</td><td>1800</td></tr><tr><td>Complete pipeline</td><td>Y</td><td>N</td><td>N</td><td>N</td><td>N</td><td>N</td><td>N</td></tr><tr><td>Accept pre-processed data</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td></tr><tr><td>Multiple samples</td><td>N</td><td>Y</td><td>N</td><td>N</td><td>N</td><td>N</td><td>N</td></tr><tr><td>Intracellular signaling</td><td>Y</td><td>N</td><td>Y</td><td>N</td><td>N</td><td>N</td><td>N</td></tr><tr><td>Scoring approach</td><td>regularized product</td><td>product</td><td>selection</td><td>modified<sup>c</sup> mean P-value</td><td>linear decomposition</td><td>selection</td><td>product<sup>d</sup></td></tr><tr><td>Export types</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Tables</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>Y</td><td>NA</td><td>Y</td></tr><tr><td>Circular plots</td><td>Y</td><td>Y</td><td>N</td><td>N</td><td>N</td><td>NA</td><td>N</td></tr><tr><td>graphML or equiv.</td><td>Y</td><td>N</td><td>Y</td><td>N</td><td>N</td><td>NA</td><td>N</td></tr><tr><td>Num. species</td><td>2</td><td>1</td><td>238</td><td>1</td><td>12</td><td>1</td><td>2</td></tr><tr><td>Platform/language</td><td>R</td><td>R</td><td>Python</td><td>website or Python</td><td>R</td><td>NA</td><td>MATLAB</td></tr></tbody></table>

repository. As for PyMINER, iTALK scores and orders the LR interactions but does not propose any cutoff. Contrary to PyMINER, iTALK relies on a selective database of LR pairs. We hence decided to consider the 118 top-scoring pairs of its output since LRscore > 0.5 gave us an average of 118 selections. Results are in Supplementary Table S6 and we observe a good overlap with our larger selections. Dissecting the data showed that all of iTALK unique LR interactions were caused by LR pairs with deprecated HUGO gene symbols; they were present in LRdb with the current symbols. For CellPhoneDB and iTALK, it was not possible to clearly separate paracrine from autocrine, and we hence compared results on the pooled paracrine and autocrine predictions. Overall, we found significant discrepancies with PyMINER and scTensor, while CellPhoneDB and iTALK were more comparable with our tool. In the latter two cases, differences originated from reference databases and scoring.

## CONCLUSION

In the nascent topic of intercellular network mapping, SingleCellSignalR is an R software package that facilitates the transformation of complex data into higher-order information. The package comes with a large variety of graphical representations and export formats to accommodate users and to support the downstream analysis of data. Of particular note are the abilities to represent a complete intercellular network and to import the latter in systems biology tools such as Cytoscape, and to explore receptor downstream signaling by integrating Reactome and KEGG pathways.

For the first time, we discussed the question of the significance of inferred LR interactions. We showed that SingleCellSignalR regularized score achieves better control of the FP compared to other solutions, independent of the retained single-cell platform. That is an important feature to explore intercellular communication networks, where each cell type might be involved in a different number of interactions. It allowed us to evidence cells that emit more sig-

nals than they receive in mouse IFE (Figure 4). The estimation of false/true positive rates entailed the construction of a benchmark relying on cell line deep RNA-seq and proteomic data to compensate for the absence of an ideal and systematic functional validation of LR interactions in specific tissues. Results obtained in mouse IFE supported those estimations.

A detailed comparison with existing tools revealed the specific features contributed by each solution. The particular design choices we made, e.g. to rely on well-documented sources of interaction data and to control false inferences, clearly produced more or more reliable information. SingleCellSignalR is open to other software packages; UMI or read count matrices prepared with other tools can be imported, or preparatory steps can be accomplished with SingleCellSignalR built-in procedures. Although SingleCellSignalR was designed for scRNA-seq data, it could be used with emerging single-cell proteomics technologies such as CyTOF (59) or SCoPE-MS (60). At last, our R library does not require advanced R skills from users.

## DATA AVAILABILITY

The R package and LRdb are available from https://github.com/SCA-IRCM under the GPL v3 license; submitted to Bioconductor.

## SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.

## ACKNOWLEDGEMENTS

We thank our colleagues Emmanuel Cornillot, Gabriel Jimenez-Dominguez and Andreï Turtoi for their insightful comments.

## FUNDING

Labex EpiGenMed Postdoctoral Fellowship [ANR-10-LABX-12-01 to J.C., S.C.A.]; Fondation ARC pour la

<footnote><sup>a</sup>Not available as a software package.</footnote>

<footnote><sup>b</sup>This number may vary since the reference interaction list is generated from STRING automatically.</footnote>

<footnote><sup>c</sup>In case a complex is considered, all subunits must be selected in a significant interaction with the partner.</footnote>

<footnote><sup>d</sup>The authors added a Wilcoxon test to assess that the LR pair score was different from zero in multiple tumors.</footnote>
Recherche sur le Cancer [PJA 20141201975 to J.C.]; Région Occitanie, programme Recherche et Société(s), and European Union (FEDER) to the project biomarqueursMETCP (to J.C.). Funding for open access charge: Inserm. *Conflict of interest statement.* None declared.

## REFERENCES

1. Potente, M., Gerhardt, H. and Carmeliet, P. (2011) Basic and therapeutic aspects of angiogenesis. *Cell*, **146**, 873–887.

2. Kopan, R. and Ilagan, M.X.G. (2009) The canonical notch signaling pathway: unfolding the activation mechanism. *Cell*, **137**, 216–233.

3. Orkin, S.H. and Zon, L.I. (2008) Hematopoiesis: an evolving paradigm for stem cell biology. *Cell*, **132**, 631–644.

4. Batlle, E., Henderson, J.T., Beghtel, H., van den Born, M.M., Sancho, E., Huls, G., Meeldijk, J., Robertson, J., van de Wetering, M., Pawson, T. (2002) Beta-catenin and TCF mediate cell positioning in the intestinal epithelium by controlling the expression of EphB/ephrinB. *Cell*, **111**, 251–263.

5. Gorelik, L. and Flavell, R.A. (2000) Abrogation of TGFbeta signaling in T cells leads to spontaneous T cell differentiation and autoimmune disease. *Immunity*, **12**, 171–181.

6. Ramilowski, J.A., Goldberg, T., Harshbarger, J., Kloppmann, E., Lizio, M., Satagopam, V.P., Itoh, M., Kawaji, H., Carninci, P., Rost, B. et al. (2015) A draft network of ligand–receptor-mediated multicellular signalling in human. *Nat. Commun.*, **6**, 7866.

7. Hanahan, D. and Weinberg, R.A. (2011) Hallmarks of cancer: the next generation. *Cell*, **144**, 646–674.

8. Junttila, M.R. and de Sauvage, F.J. (2013) Influence of tumour micro-environment heterogeneity on therapeutic response. *Nature*, **501**, 346–354.

9. Tang, F., Barbacioru, C., Wang, Y., Nordman, E., Lee, C., Xu, N., Wang, X., Bodeau, J., Tuch, B.B., Siddiqui, A. et al. (2009) mRNA-Seq whole-transcriptome analysis of a single cell. *Nat. Methods*, **6**, 377–382.

10. Pan, X., Durrett, R.E., Zhu, H., Tanaka, Y., Li, Y., Zi, X., Marjani, S.L., Euskirchen, G., Ma, C., Lamotte, R.H. et al. (2013) Two methods for full-length RNA sequencing for low quantities of cells and single cells. *Proc. Natl. Acad. Sci. U.S.A.*, **110**, 594–599.

11. Klein, A.M., Mazutis, L., Akartuna, I., Tallapragada, N., Veres, A., Li, V., Peshkin, L., Weitz, D.A. and Fournier, M.W. (2015) Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. *Cell*, **161**, 1187–1201.

12. Luecken, M.D. and Theis, F.J. (2019) Current best practices in single-cell RNA-seq analysis: a tutorial. *Mol. Syst. Biol.*, **15**, e8746.

13. Aibar, S., González-Blas, C.B., Moerman, T., Huynh-Thu, V.A., Imrichova, H., Hulselmans, G., Rambow, F., Marine, J.-C., Geurts, P., Aerts, J. et al. (2017) SCENIC: single-cell regulatory network inference and clustering. *Nat. Methods*, **14**, 1083–1086.

14. Fan, J., Salathia, N., Liu, R., Kaeser, G.E., Yung, Y.C., Herman, J.L., Kaper, F., Fan, J.-B., Zhang, K., Chun, J. et al. (2016) Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. *Nat. Methods*, **13**, 241–244.

15. Pont, F., Tosolini, M. and Fournié, J.J. (2019) Single-cell signature explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets. *Nucleic Acids Res.*, **47**, e133.

16. Puram, S.V., Tirosh, I., Parikh, A.S., Patel, A.P., Yizhak, K., Gillespie, S., Rodman, C., Luo, C.L., Mroz, E.A., Emerick, K.S. et al. (2017) Single-cell transcriptomic analysis of primary and metastatic tumor ecosystems in head and neck cancer. *Cell*, **171**, 1611–1624.

17. Skelly, D.A., Squiers, G.T., McLellan, M.A., Bolisetty, M.T., Robson, P., Rosenthal, N.A. and Pinto, A.R. (2018) Single-Cell transcriptional profiling reveals cellular diversity and intercommunication in the mouse heart. *Cell Rep.*, **22**, 600–610.

18. Kumar, M.P., Du, J., Lagoudas, G., Jiao, Y., Sawyer, A., Drummond, D.C., Lauffenburger, D.A. and Raue, A. (2018) Analysis of single-cell RNA-seq identifies cell-cell communication associated with tumor characteristics. *Cell Rep.*, **25**, 1458–1468.

19. Camp, J.G., Sekine, K., Gerber, T., Loeffler-Wirth, H., Binder, H., Gac, M., Kanton, S., Kageyama, J., Damm, G., Seehofer, D. et al. (2017) Multilineage communication regulates human liver bud development from pluripotency. *Nature*, **546**, 533–538.

20. Costa, A., Kieffer, Y., Scholer-Dahirel, A., Pelon, F., Bourachot, B., Cardon, M., Sirven, P., Magagna, I., Fuhrmann, L., Bernard, C. et al. (2018) Fibroblast heterogeneity and immunosuppressive environment in human breast cancer. *Cancer Cell*, **33**, 463–479.

21. Vento-Tormo, R., Efremova, M., Botting, R.A., Turco, M.Y., Vento-Tormo, M., Meyer, K.B., Park, J.-E., Stephenson, E., Polański, K., Goncalves, A. et al. (2018) Single-cell reconstruction of the early maternal-fetal interface in humans. *Nature*, **563**, 347–353.

22. Zhou, J.X., Taramelli, R., Pedrini, E., Knijnenburg, T. and Huang, S. (2017) Extracting intercellular signaling network of cancer tissues using ligand-receptor expression patterns from whole-tumor and single-cell transcriptomes. *Sci. Rep.*, **7**, 8815.

23. Tsuyuzaki, K., Ishii, M. and Nikaido, I. (2019) Uncovering hypergraphs of cell-cell interaction from single cell RNA-sequencing data. *bioRxiv* doi: https://doi.org/10.1101/566182, 04 March 2019, preprint: not peer reviewed.

24. Tyler, S.R., Rotti, P.G., Sun, X., Yi, Y., Xie, W., Winter, M.C., Flamme-Wiese, M.J., Tucker, B.A., Mullins, R.F., Norris, A.W. et al. (2019) PyMINER finds gene and autocrine-paracrine networks from human islet scRNA-seq. *Cell Rep.*, **26**, 1951–1964.

25. Wang, Y., Wang, R., Zhang, S., Song, S., Jiang, C., Han, G., Wang, M., Ajani, J., Futreal, A. and Wang, L. (2019) iTALK: an R package to characterize and illustrate intercellular communication. *bioRxiv* doi: https://doi.org/10.1101/507871, 04 January 2019, preprint: not peer reviewed.

26. Efremova, M., Vento-Tormo, M., Teichmann, S.A. and Vento-Tormo, R. (2020) CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. *Nat Protoc.*, doi:10.1038/s41596-020-0292-x.

27. Prasad, T.S., Kandasamy, K. and Pandey, A. (2009) Human protein reference database and human proteinipedia as discovery tools for systems biology. *Methods Mol. Biol.*, **577**, 67–79.

28. Ben-Shlomo, I., Yu Hsu, S., Rauch, R., Kowalski, H.W. and Hsueh, A.J.W. (2003) Signaling receptome: a genomic and evolutionary perspective of plasma membrane receptors involved in signal transduction. *Sci. STKE*, **2003**, R9.

29. Harding, S.D., Sharman, J.L., Faccenda, E., Southan, C., Pawson, A.J., Ireland, S., Gray, A.J.G., Bruce, L., Alexander, S.P.H., Anderton, S. et al. (2018) The IUPHAR/BPS Guide to PHARMACOLOGY in 2018: updates and expansion to encompass the new guide to IMMUNOPHARMACOLOGY. *Nucleic Acids Res.*, **46**, D1091–D1106.

30. Wu, C.H., Apweiler, R., Bairoch, A., Natale, D.A., Barker, W.C., Boeckmann, B., Ferro, S., Gasteiger, E., Huang, H., Lopez, R. et al. (2006) The universal protein resource (UniProt): an expanding universe of protein information. *Nucleic Acids Res.*, **34**, D187–D191.

31. Dimmer, E.C., Huntley, R.P., Alam-Faruque, Y., Sawford, T., O'Donovan, C., Martin, M.J., Bely, B., Browne, P., Mun Chan, W., Eberhardt, R. et al. (2011) The UniProt-GO annotation database in 2011. *Nucleic Acids Res.*, **40**, D565–D570.

32. Fabregat, A., Jupe, S., Matthews, L., Sidiropoulos, K., Gillespie, M., Garapati, P., Haw, R., Jassal, B., Korninger, F., May, B. et al. (2018) The reactome pathway knowledgebase. *Nucleic Acids Res.*, **46**, D649–D655.

33. Cerami, E.G., Gross, B.E., Demir, E., Rodchenkov, I., Babur, O., Anwar, N., Schultz, N., Bader, G.D. and Sander, C. (2011) Pathway Commons, a web resource for biological pathway data. *Nucleic Acids Res.*, **39**, D685–D690.

34. Zhang, A.W., O'Flanagan, C., Chavez, E.A., Lim, J.L.P., Ceglia, N., McPherson, A., Wiens, M., Walters, P., Chan, T., Hewitson, B. et al. (2019) Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling. *Nat. Methods*, **16**, 1007–1015.

35. Zheng, G.X.Y., Terry, J.M., Belgrader, P., Rykvin, P., Bent, Z.W., Wilson, R., Ziraldo, S.B., Wheeler, T.D., McDermott, G.P., Zhu, J. et al. (2017) Massively parallel digital transcriptional profiling of single cells. *Nat. Commun.*, **8**, 14049.

36. Wang, B., Zhu, J., Pierson, E., Ramazzotti, D. and Batzoglou, S. (2017) Visualization and analysis of single-cell RNA-seq data by kernel-based similarity learning. *Nat. Methods*, **14**, 414–416.

37. Franzén, O., Gan, L.-M. and Björkegren, J.L.M. (2019) PanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data. *Database (Oxford)*, **2019**, doi:10.1093/database/baz046.
38. McCarthy,D.J., Chen,Y. and Smyth,G.K. (2012) Differential expression analysis of multifactor RNA-seq experiments with respect to biological variation. *Nucleic Acids Res.*, **40**, 4288–4297.

39. Kanehisa,M., Araki,M., Goto,S., Hattori,M., Hirakawa,M., Itoh,M., Katayama,T., Kawashima,S., Okuda,S., Tokimatsu,T. et al. (2008) KEGG for linking genomes to life and the environment. *Nucleic Acids Res.*, **36**, D480–D484.

40. Tirosh,I., Izar,B., Prakadan,S.M., Wadsworth,M.H., Treacy,D., Trombetta,J.J., Rotem,A., Rodman,C., Lian,C., Murphy,G. et al. (2016) Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. *Science*, **352**, 189–196.

41. 8k PBMCs from a Healthy Donor (2017) https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/pbmc8k.

42. 4k Pan T Cells from a Healthy Donor (2017) https://support.10xgenomics.com/single-cell-gene-expression/datasets/1.2.0/t.4k.

43. Bagnoli,J.W., Ziegenhain,C., Janjic,A., Wange,L.E., Vieth,B., Parekh,S., Geuder,J., Hellmann,I. and Enard,W. (2018) Sensitive and powerful single-cell RNA sequencing using mcSCRB-seq. *Nat. Commun.*, **9**, 2937.

44. Rieckmann,J.C., Geiger,R., Hornburg,D., Wolf,T., Kveler,K., Jarrossay,D., Sallusto,F., Shen-Orr,S.S., Lanzavecchia,A., Mann,M. et al. (2017) Social network architecture of human immune cells unveiled by quantitative proteomics. *Nat. Immunol.*, **18**, 583–593.

45. Smoot,M.E., Ono,K., Ruscheinski,J., Wang,P.L. and Ideker,T. (2011) Cytoscape 2.8: new features for data integration and network visualization. *Bioinformatics*, **27**, 431–432.

46. Arisdakessian,C., Poirion,O., Yunits,B., Zhu,X. and Garmire,L.X. (2019) DeepImpute: an accurate, fast, and scalable deep neural network method to impute single-cell RNA-seq data. *Genome Biol.*, **20**, 211.

47. Hafemeister,C. and Satija,R. (2019) Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. *Genome Biol.*, **20**, 296.

48. Flicek,P., Aken,B.L., Beal,K., Ballester,B., Caccamo,M., Chen,Y., Clarke,L., Coates,G., Cunningham,F., Cutts,T. et al. (2008) Ensembl 2008. *Nucleic Acids Res.*, **36**, D707–D714.

49. Joost,S., Zeisel,A., Jacob,T., Sun,X., La Manno,G., Lönnerberg,P., Linnarsson,S. and Kasper,M. (2016) Single-cell transcriptomics reveals that differentiation and spatial signatures shape epidermal and hair follicle heterogeneity. *Cell Syst.*, **3**, 221–237.

50. Persson,A., Hober,S. and Uhlen,M. (2006) A human protein atlas based on antibody proteomics. *Curr. Opin. Mol. Ther.*, **8**, 185–190.

51. Castrogiovanni,P., Mazzone,V. and Imbesi,R. (2011) Immunolocalization of HB-EGF in human skin by streptavidin-peroxidase (HRP) conjugate method. *Int. J. Morphol.*, **29**, 1162–1167.

52. Bauer,T., Zagórska,A., Jurkin,J., Yasmin,N., Köffel,R., Richter,S., Gesslbauer,B., Lemke,G. and Strobl,H. (2012) Identification of Axl as a downstream effector of TGF-β1 during Langerhans cell differentiation and epidermal homeostasis. *J. Exp. Med.*, **209**, 2033–2047.

53. Langan,E.A., Vidali,S., Pigat,N., Funk,W., Lisztes,E., Bíró,T., Goffin,V., Griffiths,C.E.M. and Paus,R. (2013) Tumour necrosis factor alpha, interferon gamma and substance P are novel modulators of extrapituitary prolactin expression in human skin. *PLoS One*, **8**, e60819.

54. Bünemann,E., Hoff,N.-P., Buhren,B.A., Wiesner,U., Meller,S., Bölke,E., Müller-Homey,A., Kubitza,R., Ruzicka,T., Zlotnik,A. et al. (2018) Chemokine ligand-receptor interactions critically regulate cutaneous wound healing. *Eur. J. Med. Res.*, **23**, 4.

55. Di Meglio,P., Perera,G.K. and Nestle,F.O. (2011) The multitasking organ: recent insights into skin immune function. *Immunity*, **35**, 857–869.

56. Guttman-Yassky,E., Zhou,L. and Krueger,J.G. (2019) The skin as an immune organ: tolerance versus effector responses and applications to food allergy and hypersensitivity reactions. *J. Allergy Clin. Immunol.*, **144**, 362–374.

57. Szklarczyk,D., Franceschini,A., Kuhn,M., Simonovic,M., Roth,A., Minguez,P., Doerks,T., Stark,M., Muller,J., Bork,P. et al. (2011) The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. *Nucleic Acids Res.*, **39**, D561–D568.

58. Li,L., Dong,J., Yan,L., Yong,J., Liu,X., Hu,Y., Fan,X., Wu,X., Guo,H., Wang,X. et al. (2017) Single-Cell RNA-Seq analysis maps development of human germline cells and gonadal niche interactions. *Cell Stem Cell*, **20**, 858–873.

59. Giesen,C., Wang,H.A.O., Schapiro,D., Zivanovic,N., Jacobs,A., Hattendorf,B., Schüffler,P.J., Grolimund,D., Buhmann,J.M., Brandt,S. et al. (2014) Highly multiplexed imaging of tumor tissues with subcellular resolution by mass cytometry. *Nat. Methods*, **11**, 417–422.

60. Budnik,B., Levy,E., Harmange,G. and Slavov,N. (2018) Mass-spectrometry of single mammalian cells quantifies proteome heterogeneity during cell differentiation. *Genome Biol.*, **19**, 161.