key: cord-0806808-yv6sxkan authors: Lin, Yingxin; Loo, Lipin; Tran, Andy; Moreno, Cesar; Hesselson, Daniel; Neely, Greg; Yang, Jean Y.H. title: Characterization of cell-cell communication in COVID-19 patients date: 2020-12-30 journal: bioRxiv DOI: 10.1101/2020.12.30.424641 sha: f4d88dacdf8fc295a30bc387dc240944fa307633 doc_id: 806808 cord_uid: yv6sxkan COVID-19 patients display a wide range of disease severity, ranging from asymptomatic to critical symptoms with high mortality risk. Our ability to understand the interaction of SARS-CoV-2 infected cells within the lung, and of protective or dysfunctional immune responses to the virus, is critical to effectively treat these patients. Currently, our understanding of cell-cell interactions across different disease states, and how such interactions may drive pathogenic outcomes, is incomplete. Here, we developed a generalizable workflow for identifying cells that are differentially interacting across COVID-19 patients with distinct disease outcomes and use it to examine five public single-cell RNA-seq datasets with a total of 85 individual samples. By characterizing the cell-cell interaction patterns across epithelial and immune cells in lung tissues for patients with varying disease severity, we illustrate diverse communication patterns across individuals, and discover heterogeneous communication patterns among moderate and severe patients. We further illustrate patterns derived from cell-cell interactions are potential signatures for discriminating between moderate and severe patients. healthy human lung scRNA-seq datasets (4, 26, 27) including 189,967 cells and 44 cell types. For the Chua dataset with 5 healthy controls, 14 moderate and 13 severe samples ( Supplementary Fig. 1A-B) , scClassify is able to identify additional cell types and provide further refinement as illustrated in Fig. 1C . For example, the original "outliers epithelial" cluster was refined to "ciliated cells", and "secretory" cells to "goblet" cells. By accounting for such refinement, the new annotation recapitulates 78% of the original published analysis. The classified cell types are clearly identified by known markers (Supplementary Fig. 1C) , and further clustering of the Chua dataset generates 50 subclusters. Similar reannotation is applied to the Liao dataset (Supplementary Fig. 2A-B ) resulting in 15 cell types and 52 subclusters. To identify cell-cell interaction (CCI) differences with respect to disease severity in each dataset, we calculate the CCI scores that represent the communication probabilities among all pairs of subclusters across all ligand-receptor pairs (see Materials and Methods for details). Our group-specific CCI scores (CCI group ) aggregate the scores across all different pathways between each major cell type pair for different disease severity groups, represented as a network graph with thicker edges indicating stronger cell-cell interaction. Our results highlight that in healthy controls, most cell-cell interactions are between basal, ciliated, and goblet cells of the lung epithelium, with dendritic cells providing immune surveillance ( Fig. 2A) . As disease severity increases, cell-cell interactions become dominated by interactions between the lung epithelium and proinflammatory players within the immune compartment ( Fig. 2A-C, Supplementary Fig. 3A-E) . Overall, we observe significantly less communication (fewer edges in Fig. 2A ) in healthy individuals compared to moderate (Fig. 2B ) and severe patients (Fig. 2C ). We next study the cell-cell interaction in all three datasets (18, 28, 29) . Applying our developed workflow, we unified the annotation across all three datasets ( Supplementary Fig. 4A ). We used the Wilk dataset (with 44,721 cells and 20 cell types) as a reference (28) to reannotate the cell types for the Zhang dataset and the Arunachalam dataset. Cells with the same annotation were well integrated with scMerge (30) (Supplementary Fig. 4A) , indicating that similar cell types exist in the three studies. We observe distinct cell type compositions of these three studies, potentially due to different sampling or cell isolation procedures ( Supplementary Fig. 4B ). In all three PBMC studies, we observe that in monocyte-related interactions, COVID-19 patients generally have more cell-cell interactions than healthy controls (Supplementary Fig. 5A -C). The results across all five datasets suggest that cell type composition alone from single-cell experiments may not sufficiently discriminate between patients with different disease severity. This, along with the varying cell-cell interaction patterns across disease groups supports the further examination of the association between cell-cell interaction with disease outcomes and progression. Focusing on immune cell interactions in the upper airway, we observe different cell type interaction patterns across patients with different disease severity (Fig. 2D, Supplementary Fig. 3F ) in the Chua dataset. Compared to moderate patients, severe patients have less T cell to T cell interaction (illustrated by the thick blue edge on the T cell node in Fig. 2D ) and more macrophage to macrophage interaction (marked by thick red edge on the macrophage node in Fig. 2D ). Between different cell types, we see higher interactions between monocyte/macrophage or T cells towards neutrophils in severe patients, consistent with previous findings (31) . Similar patterns are observed in the Liao dataset where T cell to T cell interaction was higher in moderate patients, while interactions involving neutrophils and monocytes increased in severe patients ( Supplementary Fig. 2C ). The enrichment of monocyte dominated interactions in severe patients is also observed in all three PBMC datasets ( Supplementary Fig. 5D-F) . Interestingly, two of the PBMC datasets (Zhang and Arunachalam) illustrate that severe patients have a decrease in CD8 memory T cell interactions. The lack of T cell interactions in severe patients suggests potential T cell depletion, a finding supported by previous studies that lymphopenia is associated with severe disease (32) (33) (34) . Together, these data provide validation that our workflow can confirm known mechanisms and highlight new biology for further investigation. Focusing on individual pathways, Fig. 3A illustrates that all pathways can be broadly grouped into six large clusters. In particular, two of these pathway-clusters (pathway-cluster 2 marked by orange and pathway-cluster 4 marked by pink) are dominated by inflammatory pathways and these have significantly higher interaction between monocytes and neutrophils in severe patients compared to moderate ( Fig. 3A-B) . This is consistent with findings that in the productive immune response to SARS-CoV-2 infection, alveolar macrophages recognize and phagocytize apoptotic cells; however, under a dysfunctional immune response, excessive activation and accumulation of monocytes/macrophages and neutrophils leads to the overproduction of inflammatory cytokines which then damages the lung and other organs (35, 36) . To further delineate differences between moderate and severe patients observed in Fig. 2D (shown by thick red edges between monocytes and neutrophils), we investigated which subpopulations of monocytes actively interact with neutrophils ( Supplementary Fig. 6A-B) . The two inflammatory pathway-clusters mentioned above show that different cellular subtypes of monocytes in severe patients have significantly higher interaction scores in different pathways (Fig. 3B ). More specifically, we found that in severe patients, cellular subtype "monocyte 1" interacts with neutrophils through IL1, and CCL pathways, whereas interactions in moderate patients are dominated instead by TNF. Supplementary Fig. 6C -D shows that "monocyte 1" is marked by genes IL1B, IL1RN, IL8, TNFRSF1B and CCL4 and characterized by gene ontology terms "regulation of inflammatory response" as well as "regulation of apoptotic signaling pathways". The cellular subtype "monocyte 2" (marked by highly expressed IFI27), interacts primarily with neutrophils through pathways ANNEXIN and GALECTIN, which could suggest a role for this cluster in phagocytizing dying neutrophils. The cellular subtype "monocyte 3" expressing IFIT2, IFIT3, CCL8, CXCL10, and CXCL11 shows strong signatures of type 1 interferon cell-cell signaling ( Fig. 3B and Supplementary Fig. 6C ), suggesting equal support for antiviral immunity in moderate and severe patients. Alternatively, proinflammatory signaling via CXCL interactions is mainly through cellular subtype "monocyte 4", which highly expresses CCL2, CXCL1, CXCL2 and CXCL5 ( Fig. 3B and Supplementary Fig. 6C ). over time (Fig. 3C) . Interestingly, ANNEXIN downregulates across time in severe patients, but increases in moderate patients (Fig. 3C) . We observe heterogeneous interaction patterns from goblet cells to immune cells across patients and pathways (Fig. 4A ). Goblet cells are found to express high levels of genes associated with innate and antiviral immune functions indicating that the nasal epithelial cells interacting with immune cells may play an important role in reducing early viral load and this is also consistent with recent literature (37) . We observe one subgroup of severe patients (n = 3; including one deceased patient) showing clear differences in cell-cell interaction within the pathway-cluster 1 compared to moderate patients. In particular, they show a lack of interaction in the collection of pathways which includes immune signaling and costimulation pathways such as CD40, CD80, CD23 and CD86 inflammatory pathways IL6, IFN-II, and Th2 cytokines IL-4/IL10 (Fig. 4A ). Another subgroup of severe patients (n = 6) show clusters with a small subgroup of moderate patients that has low cell-cell interaction for antigen presentation (MHC-II), signaling pathways PTN and NPR2, and this subgroup is also lacking the Th2 cytokine IL-4 and the B cell activating factor BAFF (Fig. 4A ). Together, these results point to cohort heterogeneity within severe patients implicating immune costimulation or T cell polarizing pathways may contribute to disease severity with context. Focusing on a specific cellular subtype of epithelial cells (goblet 5), we observe a number of increased activities in moderate patients in the ANNEXIN pathway. This is most evident between cellular subtypes "goblet 5" and "monocyte 5". This epithelial to immune cell interaction within the ANNEXIN pathway also shows an increase in patients under moderate conditions (Fig. 4B ). Annexin plays a role in phagocytic uptake of dying cells, can drive neutrophil detachment and apoptosis, and plays a predominant role in immune resolution (38) . Remarkably, glucocorticoids, which are effective at treating COVID-19 patients, act via upregulating Annexin I (39), suggesting that natural moderate symptoms for COVID-19 may be linked to effective endogenous immune management, or that patients that respond to glucocorticoid drugs elevate Annexin cell communication pathways that then limit further inflammation, and this response is detectable via single-cell analysis. We have developed and provide an interactive resource (http://shiny.maths.usyd.edu.au/CovidCellInteraction/) to enable further investigation of cell-cell interaction at different resolutions from aggregate interaction between two major cell types to visualize expression values for specific ligand-receptor pairs. Finally, we found that information from cell-cell interactions provides a discriminating signal for patients with different disease progression. Table 1 ). The accuracy rate of leave-one-out cross-validation (LOOCV) based on the first three PCs using k nearest neighbor classification (k = 3) is 84.4%, highlighting the ability of cell-cell interaction features to predict the degree of severity of patients. By repeating our workflow on the three PBMC datasets, we further demonstrate that using CCI features can achieve higher LOOCV accuracy rate than using cell type composition as features. Despite the limited samples, repeating our workflow on a smaller dataset within BALF tissues in the Liao dataset demonstrates similar findings that cell-cell communication patterns between goblet cells to immune cells has potential discriminating power. A better understanding of virus and host cell interaction at the cellular level is an important component in understanding infectious disease progression and is critical for developing a treatment for the disease. In this paper, we provide a comprehensive workflow to integrate and examine multiple COVID-19 single-cell RNA-seq datasets to identify differential cell-cell interaction (CCI) pathways with respect to disease. Our initial results in upper airway tissues show strong intra-epithelial communication in the healthy lung, whereas the immune system then dominates communication pathways during COVID-19. We then discover that despite a higher cell-cell interaction (tCCI score) in severe patients compared to moderate patients between immune and neutrophil cells, the CCI scores between epithelial and immune cells are heterogeneous among severe patients, with a subpopulation illustrating lower CCI score when compared to moderate patients. Furthermore, features extracted from cell-cell interactions are potential signatures for discriminating between moderate and severe patients. In most multi-omics profiling in patients with COVID-19, strong acute inflammatory responses are commonly found in most of the cell types as expected. Since the airway epithelium is the primary site of infection for SARS-CoV-2 causing disease, investigating how epithelial cells interact with immune cells differentially leads to a better understanding of the initial host reaction to viral infection. Therefore, examining cell-cell communication offers an analytical approach to characterize specific cell type interaction and identify potential immune response drivers that results in different degrees of disease severity. The importance of using a workflow that accounts for cohort heterogeneity in examining severe and moderate patients is clearly illustrated when we examine the interaction pattern between nasal epithelium as a ligand and various receptors in immune cells. This is a different approach to the one taken by Chua and colleagues (25) , where they measure the overall/aggregate interaction between epithelial and immune cells and found a higher overall/aggregate interaction between epithelial and immune cells in severe patients. Here, when we examine the cell-cell interaction relationships at the individual sample level, we observe clear cohort heterogeneity among severe patients, and we are able to discover a subgroup of the moderate patients with higher interaction between epithelial and immune cells. In this study, we focus on the cell communication within COVID-19 patients via ligand-receptor signaling. Several methods have been developed recently to infer such cell-cell interaction from scRNA-seq data, such as CellPhoneDB, SingleCellSignalR, NicheNet and CellChat (20) (21) (22) 40) . Most of these methods aim to identify the significant ligand and receptor gene pair between two cell populations with the most recent method CellChat (22) that accounts for additional signaling factors. In addition, CellChat systematically categorizes the ligand-receptor pairs based on their signaling pathways, providing a comprehensive interpretation of cell-cell communication from single-cell RNAseq. There are also other types of cell communication like physical cell interaction that can be further investigated. Technology to sequence physically interacting cells like PIC-seq has been used to investigate epithelial-immune interaction and infectious disease in mice (41) . Application of such technology in COVID-19 research will potentially allow characterization of differential physical intercellular interaction at high resolution. Our analysis suggests the heterogeneity of cell-cell interaction patterns within patients, even if they have similar degrees of symptoms. One key variability is the sampling time since the onset of symptoms, as this may not fully capture the true underlying disease progression within each individual. Other potential factors that lead to the variability include age, gender, comorbidities and viral load. Currently, with the limited number of samples from patients with similar clinical characteristics, accounting for these uncertainties in modelling is challenging. Towards the future, as more large single-cell profiling resources in COVID-19 become publicly available, integrative analysis and metaanalysis of these studies by incorporating patient diversity to our workflow will provide a more comprehensive characterization of cell-cell interaction patterns in COVID-19 patients. Nevertheless, using the current databases our workflow supports that cell-cell interactions provide more meaningful predictions of disease progression ( Figure 4C ). In summary, our novel workflow enables integrative analysis of five different COVID-19 scRNA-seq data sets with a total of 415,856 cells and 85 samples. This generalizable workflow was built on the latest single-cell analytical methods and enables the identification of differential cell-cell interaction across disease progression. We discover clear cohort heterogeneity among the severe patients in the interaction between epithelial and immune cells, with signatures that can be linked with patient outcome. Together we provide a validated workflow for integration and analysis of diverse single-cell sequencing data to pinpoint communication networks that control disease outcome. [C] Wilk dataset -The raw count matrices of single-cell RNA-seq data from PBMC with metadata were downloaded from the COVID-19 Cell Atlas: https://www.covid19cellatlas.org/#wilk20 (28) . This data contains 6 healthy controls, 3 moderate patients and 4 severe patients. [D] Arunachalam dataset -The raw count matrices of single-cell RNA-seq data from PBMC and the clinical information were downloaded from GEO under accession number GSE155673. This data has 5 healthy controls, 3 moderate patients and 4 severe patients (29) . The cells with more than 20% mitochondrial proportion and UMI count greater than 50,000 are removed from the downstream analysis. [E] Zhang dataset -The raw sequence files of single-cell RNA-seq data from PBMC are downloaded from the Genome Sequence Archive of the Beijing Institute of Genomics (BIG) Data Center, BIG, Chinese Academy of Science using the accession code HRA000150 (18) . Cell Ranger (v3.0.2) with human reference version GRCh38 were used to generate the raw count matrices. The data includes 5 healthy controls, 7 moderate patients and 4 severe patients. Only the cells retained from the original study are used. Processing: For each dataset, we performed size factor standardization and log transformation on the raw count expression matrices using the logNormCount function in the R package scater (version 1.16.2) and generated log transformed gene expression matrices for analysis. Step 1 -Cell type annotation -For a given dataset, we perform a cell type identification using the scClassify framework (19) . Specifically, to identify the cell types from the Chua dataset and the Liao dataset, we performed a modified version of the joint classification from scClassify that incorporates the concept of iterative supervised learning. The initial model is built from four reference datasets including annotated cell information from healthy human lungs (4, 26, 27) . The final cell type labels were determined by the majority vote from individual classification labels using each single reference. An additional scClassify model based on the assigned cells was then built to predict the cells that are classified as "intermediate" or "unassigned" in the previous step. To identify cell types from the PBMC datasets, we used the Wilk dataset as a reference (28) to build the model and use it to predict the cell types for the Zhang dataset and the Arunachalam dataset. To prevent over clustering, we follow a similar workflow described in clusterExperiment to collapse the identified subclusters (23) . Hierarchical clustering is first performed on the aggregated average expression of each subcluster to construct a cluster hierarchy, and then from the bottom to top, the clusters of the same branches are merged if less than 10 genes are differentially expressed (log fold change > 1, FDR < 0.01). Note that we identified some cellular subtypes (ionocytes and squamous) that are inconsistently annotated between the original Chua dataset and scClassify (classified as goblet cells). In this instance, based on marker expression, we manually reannotated these two cell types using the original annotation for the downstream analysis. For a given individual sample and a pair of subclusters (i.e. cellular subtypes) obtained in Step 2, we calculate the aggregated ligand-receptor interaction score based on CellChat (22) The implementation is available as R code stored at the GitHub, https://github.com/SydneyBioX/COVID_CCI_analysis and as a web shiny application at http://shiny.maths.usyd.edu.au/CovidCellInteraction/. The output of the cell-cell interaction analysis can be considered as a three-dimensional array representing the cell-cell interaction (CCI) score. Let !"# denote the cell-cell interaction (CCI) score generated from the computational workflow for a pair of cellular subtypes , where ∈ (defined below as a set consisting of all pairs of cellular subtypes), pathway with = 1, . . . , , and individual sample with = 1, . . . , . In general, an individual sample represents the sample from one individual collected at a specific time point. For major cell types, we denote them by the sets $ , % , . . . , & and within a given major cell type ' consisting of ' cellular subtypes '$ , '% , . . . , '( ! we can write ' = { ') | = 1, . . . , ' }. We can represent a pair of cellular subtypes as = 6 ') , *+ 7, where , = 1, . . . , ; = 1 , . . , ' and = 1, . . . , * . Here, we consider ') as the sender cellular subtype within major cell type ' and *+ as the receiver cellular subtype from major cell type * . The collection of all pairs of cellular subtypes is written as = <6 '), *+ 7| , = 1, . . . , ; = 1, . . . , ' ; = 1, . . . , * =. We further denote -! ,-" as a subset of containing only pairs of cellular subtypes from major cell type ' to * which is For a given sample , the following measures of interest are explored: • Subtype cell-cell interaction (sCCI) between a pair of cellular subtypes for an individual sample is calculated as sCCI ( , )= @ !"# " . This measure totals the cell-cell interaction score across all pathways. Calculating this score for each pair of cellular subtypes and each individual sample is the same as totaling the array across the pathways resulting in a | | × two-dimensional matrix. • Pathway specific cell-cell interaction (pCCI) from the major cell type ' to the major cell type * for a pathway and an individual sample is pCCI 6 ' , * , , 7 = @ where -! ,-" is defined as above. This is a measure that sums the cell-cell interaction scores across all cellular subtypes between any two major cell types. For each pair of ( ' , * ), calculating this statistic for each pathway and individual sample results in a × matrix (see Figure 3A ). • Total CCI (tCCI) from major cell type ' to major cell type * for an individual sample is defined as tCCI6 ' , * , 7= F @ x !"# " = ∑ sCCI( , ) , where -! ,-" is defined as above. This is a measure that sums the cell-cell interaction scores across all cellular subtypes between two major cell types and across all pathways. For each individual sample , calculating the tCCI statistic for each pair of ( ' , * ) will result in a × matrix that can be visualized as a heatmap or network graph. • Suppose represents a set of pathways belonging to the same cluster termed as a pathwaycluster (see Clustering in Methods). The pathway-cluster cell-cell interaction for an individual sample between a pair of cellular subtypes is defined as psCCI( , , )= We calculate a group specific cell-cell interaction (CCI group ) between two cellular subtypes where groups represent any treatment of interest. Here it refers to control and disease progression such as moderate and severe patients. Let group denote a set of individual samples under the same condition of interest, where | group | indicates the size of the set. For example, the total number of samples having moderate response to COVID-19 in the dataset (see Fig. 2 A-C). The CCI group from the major cell types ' to the major cell types * can be calculated by #∈ group , where ( ) = 3456'((4)9 36:;(4)56'((4)9 is a scaling function to scale between individual samples. In practice, the differential CCI from ' to * between moderate (CCI moderate ) and severe (CCI severe ) patients can be calculated by CCI severe − CCI moderate measuring the differential patterns of the cell-cell interaction across different disease severity (see Figure 2D ). The pathway-cluster cell-cell interaction (used in Figure 3B ) for a group of individuals ℒ between a pair of cellular subtypes is simply the sums of psCCI across individual with a group ℒ and can be written as . For a pair of cellular subtypes , calculating this statistic results in a | | × | group | matrix. Suppose we have multiple samples collected from the same individual at different time points, say early and late then the cell-cell interaction across disease progression is the log-ratio of cell-cell interaction (illustrated in Fig. 3C ) between these two time points for a given pair of cell types (sender cell type ' , receiver cell types * within a pathway is log SpCCI6 ', * , late 7/pCCI6 ', * , early 7T. We group various pathways based on the similarity of intercellular communication patterns using hierarchical clustering with Euclidean distance and ward.D2 agglomerative method implemented in the function hclust in R. We integrate the three PBMC datasets using a modified version of scMerge (30) . Here, cell types annotated by scClassify are used as an input to scMerge to construct pseudo-bulk expression profiles. The resulting profiles are used to identify mutual nearest subgroups as pseudo-replicates and to estimate parameters of the scMerge model. To select the cell-cell interaction features that discriminate across samples under different conditions, we performed a Kruskal-Wallis rank sum test on pathway-specific cell-cell interaction (pCCI) to select the pathways that are significantly different across samples from healthy controls, moderate patients and severe patients. Feature selection is based on pCCI features with an adjusted p-value less than 0.1 for the Chua dataset, less than 0.2 for the Wilk and Zhang datasets and less than 0.4 for the Arunachalam dataset, we termed these selected features as "Top CCI". For the Chua dataset, we also selected the top pCCI from the cell-cell interaction between the two major epithelial cell types (Goblet and Ciliated) and the immune cell types (B cells, dendritic cells, macrophages, monocytes and T cells), termed as "Epi-Immune CCI". We further considered cell type proportion as another type of feature. The classification model to predict the samples' condition is built with linear discriminant analysis (LDA) and random forest (RF) on the selected features (Top CCI, Epi-Immune CCI, and cell type proportion) as well as k nearest neighbor classification (with k = 1, 3) using the first 3 principal components of the pCCI matrix. The classification performance was determined by leave-one-out cross-validation. Differential gene expressions were identified using moderated t-statistics implemented in the R package limma (version 3.44.3). The gene set over-representation analysis for the significant DE genes (top 100 genes selected) with biological process (BP) gene ontology is measured using the "enrichGO" function in the R package clusterProfiler (version 3.16.0) (43) . Significant GO term is defined by q-value < 0.1. To facilitate the interpretation of the complex data set, we have created an online interactive tool which allows researchers to explore different parts of the data. The first tab of the tool contains four columns. The first column allows the user to select two groups (or individual samples) to compare and it displays the associated cell-cell interaction network. The second column shows the difference between Airborne transmission of SARS-CoV-2: The world should face the reality Construction of a human cell landscape at single-cell level Single-Cell Transcriptomic Analysis of Human Lung Provides Insights into the Pathobiology of Pulmonary Fibrosis A cellular census of human lungs identifies novel cell states in health and in asthma Estimation of the asymptomatic ratio of novel coronavirus infections (COVID-19) Clinical course and outcomes of critically ill patients with SARS-CoV-2 pneumonia in Wuhan, China: a single-centered, retrospective, observational study Massachusetts Consortium for Pathogen Readiness, SARS-CoV-2 viral load is associated with increased disease severity and mortality Coronavirus Disease 2019 in patients with inborn errors of immunity: an international study Dysregulation of Immune Response in Patients With Coronavirus Dysregulated Type I Interferon and Inflammatory Monocyte-Macrophage Responses Cause Lethal Pneumonia in SARS-CoV-Infected Mice High expression of ACE2 receptor of 2019-nCoV on the epithelial cells of oral mucosa The protein expression profile of ACE2 in human tissues Mineralocorticoid receptor blocker increases angiotensin-converting enzyme 2 activity in congestive heart failure patients Alveolar macrophage dysfunction and cytokine storm in the pathogenesis of two severe COVID-19 patients The cytokine storm in COVID-19: An overview of the involvement of the chemokine/chemokine-receptor system Clinical features of patients infected with 2019 novel coronavirus in Wuhan Single-Cell Sequencing of Peripheral Mononuclear Cells Reveals Distinct Immune Response Landscapes of COVID-19 and Influenza Patients Single-cell landscape of immunological responses in patients with COVID-19 scClassify: sample size estimation and multiscale classification of cells using single and multiple reference CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes SingleCellSignalR: inference of intercellular networks from single-cell transcriptomics Inference and analysis of cell-cell communication using clusterExperiment and RSEC: A Bioconductor package and framework for clustering of single-cell and other large gene expression datasets Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19 COVID-19 severity correlates with airway epithelium-immune cell interactions identified by single-cell analysis A single-cell atlas of the human healthy airways A molecular cell atlas of the human lung from single-cell RNA sequencing A single-cell atlas of the peripheral immune response to severe Systems biological assessment of immunity to mild versus severe COVID-19 infection in humans scMerge leverages factor analysis, stable expression, and pseudoreplication to merge multiple single-cell RNA-seq datasets The trinity of COVID-19: immunity, inflammation and intervention Lymphopenia predicts disease severity of COVID-19: a descriptive and predictive study Clinical and immunological features of severe and moderate coronavirus disease 2019 Reduction and Functional Exhaustion of T Cells in Patients With Coronavirus Disease 2019 (COVID-19) Complement as a target in COVID-19? The potential danger of suboptimal antibody responses in COVID-19 HCA Lung Biological Network, SARS-CoV-2 entry factors are highly expressed in nasal epithelial cells together with innate immune genes Annexin A1 and glucocorticoids as effectors of the resolution of inflammation Antiinflammatory action of glucocorticoids--new mechanisms for old drugs NicheNet: modeling intercellular communication by linking ligands to target genes Dissecting cellular crosstalk by sequencing physically interacting cells Feature selection and dimension reduction for single-cell RNA-Seq based on a multinomial model He, clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters The authors declare no competing interests.Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.