Single-cell systems immune monitoring approaches offer new ways to compare how an individual patient’s cells respond to treatment or change during infection (Chattopadhyay et al., 2014; Davis et al., 2017; Greenplate et al., 2016a; Schultze, 2015). However, the computational analysis tools used widely in this type of systems immunology study were primarily designed to track major cell populations representing >1% of a sample. Viral immune and cancer immunotherapy responses can include mechanistically important and extremely rare T cells that proliferate rapidly over the course of days but as an aggregate exist as <0.1% of blood CD3⁺ T cells at their peak. These cells can be tracked genetically through clonal expansion, but may be lost in computational analyses focused on describing the global landscape of phenotypes. The specific expansion or contraction of phenotypically distinct cells may be a hallmark feature of key immune effectors and could reveal these cells without the need for prior knowledge of their identity or specialized tracking reagents like MHC tetramers.

The datasets tested here were all suspension flow cytometry, a data type where it is typical to have multiple snapshot samples of cells over time; however, an ongoing challenge in the field is to match or register cells to their phenotypic cognates between samples (Irish, 2014; Pyne et al., 2014; Weber and Robinson, 2016). Analysis algorithms typically rely on aggregate statistics for clustered groups of cells, but the process of grouping the cells works best with larger, established populations (Diggins et al., 2015; Irish et al., 2006; Saeys et al., 2016) and can depend on pre-filtering of cells by human experts (Greenplate et al., 2016a; Greenplate et al., 2019). Cytometry tools like SPADE (Bendall et al., 2011; Qiu et al., 2011), FlowSOM (Van Gassen et al., 2015), Phenograph (Levine et al., 2015), Citrus (Bruggner et al., 2014), and RAPID (Leelatian et al., 2020) generally work best to characterize cell subsets representing >1% of the sample and are less capable of capturing extremely rare cells or subsets distinguished by only a fraction of measured features. Tools like t-SNE (Davis et al., 2013; Krijthe et al., 2015), opt-SNE (Belkina et al., 2019), and Uniform Manifold Approximation (UMAP) (Becht et al., 2018; McInnes et al., 2018) embed cells or learn a manifold and represent these transformations as algorithmically generated axes. For a biologist, these tools provide a way to organize cells according to phenotypic relationships that span multiple measured features, such as the proteins quantified on each of millions of cells in the datasets here. In addition to assisting with data visualization, these tools frequently reveal unexpected cells and facilitate their identification through manual or automated clustering (Davis et al., 2013; Becher et al., 2014; Diggins et al., 2015; Diggins et al., 2017; Gandelman et al., 2019; Leelatian et al., 2020). Sconify (Burns et al., 2018) is one such tool that applies k-nearest neighbors (KNN) to calculate aggregate statistics for the immediate phenotypic neighborhood around a given cell on a t-SNE plot representing data from multiple cytometry samples. This approach to creating a population around every cell was a key inspiration for the Tracking Responders EXpanding (T-REX) tool presented here, which applies KNN to every cell to pinpoint rare cells in phenotypic regions of significant change. In addition to combining UMAP, KNN, and Marker Enrichment Modeling (MEM) in a rapid, unsupervised analysis workflow for paired samples from one individual, T-REX contrasts with prior approaches in its specific focus on the regions of great difference between samples. This T-REX design is based on the observation that, in the absence of a perturbation such as disease or infection, adults tend to have a stable signature of blood cell abundances over weeks to months (Greenplate et al., 2019; Lakshmikanth et al., 2020; Mathew et al., 2020), and the hypothesis that short-term, dramatic changes in rare immune cell subsets will identify cells associated with exposure to an immunogenic agent, such as a virus.

Data types used to challenge the T-REX algorithm here included a new spectral flow cytometry study (Dataset 1) and three existing mass cytometry datasets (Dataset 2, Dataset 3, and Dataset 4). Mass cytometry is an established technique for human immune monitoring where commercial reagents presently allow 44 antibodies to be measured simultaneously per cell (Greenplate et al., 2016a; Mistry et al., 2019; Spitzer and Nolan, 2016). Spectral flow cytometry is gaining attention in human immune monitoring as it generates data that compares well to mass cytometry (Ferrer-Font et al., 2020; Mistry et al., 2019). Spectral flow cytometers collect cells at around 10-fold the number of cells per second as mass cytometers. While the availability of spectrally distinct antibody-fluorochrome conjugates imposes some practical limits on spectral flow cytometry at present, established panels like the one in Dataset 1 measure ~30 features per cell with excellent resolution, and that capacity is expected to roughly double in the next few years as recent work has demonstrated 40 features (Park et al., 2020). Spectral flow cytometry is thus well-matched to studies of very low-frequency cells, as was the case in Dataset 1, where a goal was to computationally pinpoint hundreds of virus-specific T cells in datasets of over 5 million collected cells.

Datasets 1 and 2 were from individuals infected with two different respiratory viruses, rhinovirus (RV) or SARS-Cov-2, respectively. Respiratory viruses are ubiquitous, and while some, like rhinovirus, are generally benign, they nonetheless pose risks to patients with underlying chronic health conditions. The common colds associated with rhinovirus are characterized by shifts in very rare virus-specific cells in the blood (Muehling et al., 2020; Muehling et al., 2018). In contrast, novel respiratory viruses, such as SARS-CoV-2, the coronavirus causing COVID-19, continue to emerge that enact high morbidity and mortality, even among healthy subjects. Understanding the immune response to such viruses is vital to treatment and vaccine design, and there has been rapid progress applying human immune monitoring to COVID-19 patients (Mathew et al., 2020; Rodriguez et al., 2020). An ongoing challenge in the field is to quantitatively compare novel diseases, like COVID-19, to other disease states and immune responses. T cells are pivotal to such responses. Severe COVID-19 has been linked to a pathogenic ‘cytokine storm’ in which cellular immune responses likely play a crucial role (Ragab et al., 2020). Nonetheless, in the case of both rhinovirus and COVID-19, it is clear that host factors are a key determinant of the degree of the T cell response (Mathew et al., 2020; Muehling et al., 2020). Datasets 3 and 4 were from cancer patients that included melanoma (MB) patients being treated with α-PD-1 checkpoint inhibitor therapies or acute myeloid leukemia (AML) patients undergoing induction chemotherapy. By tracking the CD4⁺ T cells that expand rapidly during infections and respond to immunotherapy, it may be possible to pinpoint or therapeutically guide cells into helpful vs. harmful roles or niches. Overall, a goal of this study was to develop an automated, quantitative toolkit for immune monitoring that would span a wide range of possible immune changes, identify and phenotype statistically significant cell subsets, and provide an overall vector of change indicating both the direction and magnitude of shifts, either in the immune system as a whole or in a key cell subpopulation.

We report here T-REX, a novel unsupervised machine learning algorithm for characterizing cells in phenotypic regions of significant change in a pair of samples (Figure 1). The primary use case for developing the T-REX algorithm was a new dataset from individuals infected with rhinovirus, where changes in the peripheral immune system are expected in very rare memory cells responding directly to the virus (Dataset 1). Infection with rhinovirus is known to induce expansion of circulating virus-specific CD4⁺ T cells in the blood, and a key feature of the new rhinovirus dataset here is that rare and mechanistically important virus-specific CD4⁺ T cells were marked with MHC II tetramers in the context of multiple other T cell markers. The T-REX algorithm was blinded to tetramers during analysis so that they could subsequently be used to test algorithm performance. In addition, T-REX was tested with paired samples from patients with moderate or severe COVID-19 (Dataset 2), melanoma patients being treated with α-PD-1 checkpoint inhibitor therapies (Dataset 3), and acute myeloid leukemia patients undergoing induction chemotherapy (Dataset 4). These datasets were used to determine whether the T-REX algorithm functions effectively across a spectrum of human immune monitoring challenges and to see how the algorithm performs when changes are restricted to rare cell subsets, as in Dataset 1 and Dataset 3, or when many cells may be expanding or contracting, as in Dataset 2 and Dataset 4.

Figure 1

Download asset Open asset

T-REX identifies cells in phenotypically distinct regions of significant change

For the rhinovirus challenge study in Dataset 1, sample pairs available for T-REX included cells taken immediately prior to intranasal inoculation with virus (i.e., pre-infection, day 0), as well as those during (day 7) or following inoculation (day 28). Cells were subsampled equally from each timepoint and then concatenated for a single UMAP specific to each analysis pair. UMAP axes were labeled to indicate they were specific to a comparison for a given individual (Figure 2). Thus, each UMAP comparison was a new run of the algorithm. Although it is also possible to map all sample times or all individuals into a single UMAP for analysis, a key goal here was to imagine a minimal T-REX use case with only a pair of samples from one individual. The features selected for UMAP analysis were intentionally limited to surface proteins in order to test whether suitable features for live cell fluorescence-activated cell sorting (FACS) could be identified. Following UMAP, each cell was used as the seed for a KNN search of the local neighborhood within the UMAP axes (i.e., the KNN search was within the learned manifold, as with the analysis in Sconify [Burns et al., 2018] or RAPID [Leelatian et al., 2020]). The k-value for KNN was set to 60 as a starting point based on prior studies and later optimized. For each cell, the KNN region could include cells from either time chosen for analysis, and the percentage of each was calculated to determine the representation of each sampled time in a cellular neighborhood. When cells in regions of expansion (≥95% of cells in the KNN region from one sampling time) were clustered together in one phenotypic region of the UMAP, they were considered a ‘hotspot’ of significant change. Cells in change hotspots were aggregated and the phenotype automatically characterized using MEM (Diggins et al., 2017). MEM labels here indicated features that were enriched relative to a statistical null control on a scale from 0 (no expression or enrichment) to +10 (greatest enrichment). Ultimately, T-REX and MEM were used to reveal hotspots of ≥95% change and assign a label that could be used by experts to infer cell identity.

Figure 2 with 4 supplements see all

Download asset Open asset

In the human rhinovirus challenge study yielding Dataset 1, MHC class II tetramers were used to identify rhinovirus-specific CD4⁺ T cells with the goal of tracking phenotypic changes over the course of infection. Increases in tetramer+ cells on day 7 (Figure 2A) corresponded to the acute infection phase (Muehling et al., 2020). This tetramer tracking system for virus-specific T cells provided an opportunity to test whether the cells identified by T-REX were biologically significant by leaving the tetramer stain features out of the computational analysis (i.e., not using tetramers to make the UMAP or in other parts of T-REX) and then testing to see whether hotspots of cellular change identified by T-REX were statistically enriched for tetramer+, virus-specific cells. A hotspot region was considered to be enriched for tetramer+ cells if >5% of the cells in the region were tetramer+. In the example subject shown, the pairwise comparisons used in T-REX analysis included CD4⁺ T cells from day 0, immediately prior to rhinovirus infection, and day 7, a well-studied timepoint at which rare, virus-specific CD4⁺ T cells are observed at higher frequencies (Muehling et al., 2020). This trajectory of virus-specific cell expansion was confirmed by a peak in the log₂ fold change in the frequency of tetramer+ CD4⁺ T cells (Figure 2A). Applying T-REX to the rhinovirus data revealed that KNN regions with expansion from day 0 to day 7 were greatly enriched for tetramer+ cells, as compared to regions with less expansion (Figure 2B). UMAP axes were labeled as UMAP_RV001_7_0 to denote this UMAP analyzed day 0 and day 7 for individual RV001 (Figure 2C). Regions of contraction were observed but were not enriched for tetramer+ cells, except in the case of one individual, RV007, studied here (Figure 3). Notably, two of the eight study subjects challenged with rhinovirus were not infected (RV002 and RV003); all other individuals were infected (Supplementary file 1). As the focus of this study was virus-specific T cells, CD4⁺ T cells were the main populations analyzed in all subjects (Supplementary file 2). However, T-REX finds the largest tetramer+ cell population with the same phenotype even when using all 4.9 million live cells instead of only the 1.3 million CD4⁺ T cells (Figure 2—figure supplement 1).

Figure 3

Download asset Open asset

A key question for the T-REX algorithm is where to set a statistical cutoff for what is considered to be a biologically significant amount of expansion. Two change cutoffs were tested with subject RV001, ≥90% and ≥95% (Figure 2C). Using a cutoff of ≥95% identified 2/2 (100%) tetramer+ hotspots of change for RV001 and did not identify any additional regions that were not tetramer hotspots, whereas the ≥90% cutoff identified both tetramer+ hotspots and an additional tetramer– hotspot (Figure 2C, top). Thus, ≥95% represented a stringent cutoff that still captured biologically significant cells. An analysis of tetramer enrichment as a function of percentile of expansion from day 0 to day 7 (Figure 2B) showed that tetramer+ cells were not commonly observed to be in local neighborhoods around cells with change below 80% in their KNN region. In contrast, above 90% change, the median CD4⁺ T cell had 10% or more tetramer+ neighbors around it in the KNN region (Figure 2B). Thus, only regions of 80% or more expansion from day 0 to day 7 were enriched for tetramer+ CD4⁺ T cells in study individual RV001.

A k-value of 60 effectively identified immune hotspots in T-REX

A critical question for KNN analysis is the value of k, the number of neighbors to assess. While it is useful to have a lower k-value as the analysis will complete more quickly, increasing the k-value might better represent the phenotypic neighborhood or be more statistically robust. To assess how k-value impacted detection of cells in regions of change and the degree to which these cells were virus-specific in rhinovirus challenge Dataset 1, the k-value was systematically changed. In example case RV001, an optimal k was determined to be an inflection point in a graph of the average tetramer enrichment (y-axis, Figure 4) vs. increasing values of k (x-axis, Figure 4). To calculate this curve, a KNN search was repeated while increasing k in steps from 0 to 300 for every cell in each sampling. This analysis was performed for all tetramer+ cells from day 7 (dark purple, Figure 4), all tetramer+ cells from day 0 (light purple, Figure 4), and, as a negative control, random tetramer– cells from day 7 (black, Figure 4). Within each of these neighborhoods, tetramer enrichment was calculated. This approach identified the inflection point of the percent tetramer+ curve as k = 70 for RV001 (Figure 4). In further analysis of the remaining infected rhinovirus subjects, optimal k-values ranged from 30 to 80. A k-value of 60 was chosen and used in all other analyses of rhinovirus subjects (Figure 3), as well as Datasets 2, 3, and 4 described below.

Figure 4

Download asset Open asset

Regions of significant change contained rhinovirus-specific CD4⁺ T cells in Dataset 1

The association between regions of change and enrichment for virus-specific cells observed in the example subject shown (Figure 2B) was observed in five infected rhinovirus subjects; tetramer+ CD4⁺ T cells were not enriched in KNN regions around cells that had not expanded from day 0 to day 7 (one infected, two uninfected; Figure 2—figure supplement 2). This observation suggested that cutoffs at the 5th and 95th percentile would accurately capture cells representing phenotypic regions with significant change over time. In addition, 15th and 85th percentiles were chosen as cutoffs to capture a more moderate degree of change and track cells that might still be of interest but not from regions experiencing significant change. The remaining cells in phenotypic regions between the 15th and 85th percentiles were not considered to have not changed significantly in the context of these studies. Going forward, it was of interest to determine how often regions of significant change (i.e., the 95th and 5th percentile cutoffs) would contain tetramer+ CD4⁺ T cells in different individuals participating in the rhinovirus challenge study. Cells in regions of significant expansion (≥95%) were also from regions that were enriched for virus-specific cells in nearly all rhinovirus-infected individuals (4/6 at 95% cutoff, 5/6 at 85% cutoff) (Figure 2, Figure 2—figure supplement 2, Figure 3). Thus, by focusing specifically on cells in regions representing the most change over time, T-REX analysis revealed subpopulations containing virus-specific cells. This highlights the ability of T-REX to pinpoint such cells without the use of antigen-specific reagents. Following T-REX, MEM analysis was performed using all available features, including intracellular features not used to define the UMAP space (TCF1, TBET, and Ki-67). The phenotype of the regions of significant change enriched for virus-specific cells was quantitatively described with MEM scores (hotspot 1: ▲CD45R0⁺¹⁰ CD38⁺⁸ ICOS⁺⁶ CCR5⁺⁵ TCF1⁺⁵ CD27⁺⁴ PD-1⁺⁴ CXCR3⁺³ CD95⁺³ TBET⁺² CD25⁺²; hotspot 2: ▲CD45R0⁺¹⁰ CD38⁺⁸ ICOS⁺⁶ CD27⁺⁵ TCF1⁺⁵ CCR5⁺⁴ CXCR3⁺³ CD95⁺³ CCR7⁺³ PD-1⁺² CD25⁺² CXCR5⁺²). The change hotspots thus contained activated memory cells (CD45RO + CD38+) that were notable for their early memory/stem-like T cell signature (TCF1+ CD27+), as well as their expression of CCR5 and CXCR3, both of which are chemokine receptors found on rhinovirus-specific CD4+ T cells that respond during infection (Muehling et al., 2020; Muehling et al., 2018).

To determine the sensitivity of this method, all tetramer+ regions were next reviewed, including those that did not meet the criteria for hotspots of significant change (Figure 2C, bottom). In analysis of RV001, 66.6% (2/3) of tetramer+ regions were captured, meaning there was one region with lower change that contained tetramer+ cells. However, there were only 87 cells in these missed regions compared to 896 cells and 826 cells in the regions with ≥95% expansion, confirming that T-REX captured the majority of virus-specific cells in the dataset.

In the case of an emerging infectious disease, it may not be possible to have a pre-infection sample and it would be useful to know whether T-REX analysis of change between a peak of infection and a later time might also reveal virus-specific T cells. To test this idea, pairwise comparisons were performed with cells from day 7 following rhinovirus inoculation and at day 28 after inoculation (Figure 5). Strikingly, cells in phenotypic regions of significant change again were enriched for virus-specific, tetramer+ CD4⁺ T cells. The MEM values for these cells further identified them as CD45R0⁺ memory cells enriched for CD38, ICOS, CD27, TCF1, CXCR5, PD-1, and CD95 expression, a phenotype matching that of the cells identified in the day 0 to day 7 analysis for this individual (RV001, Figure 5).

Figure 5

Download asset Open asset

T-REX tracking of direction and degree of change contextualizes diverse immune responses

The next goal was to test T-REX with additional data types and contextualize the results from rhinovirus infection (Dataset 1) with changes that might be observed in other immune responses, such as another respiratory infection (Dataset 2), cancer immunotherapy (Datasets 3), or cancer chemotherapy (Dataset 4). To accomplish this, metrics for degree of change as well as direction of change in each sample were devised (Figure 6A). Degree of change was calculated as the sum of the percent of cells in the 5th and 95th hotspots of change. Therefore, degree of change is the percentage of the sample that has changed significantly (≥95%) between the two samples (i.e., the percentage of the sample that was dark red or dark blue in the T-REX plot). A degree of change of 100% would mean that the KNN region around every cell was populated only by cells from the same sample, which is interpreted as meaning the pair of samples had changed completely (half the sample contracting and half expanding). Direction of change is the ratio of ≥95% expansion to contraction in the sample and was calculated as the difference between the number of cells in the 95th and 5th hotspots of change divided by the sum of the number of cells in the 95th and 5th hotspots of change. A negative direction of change value indicates that there was significant contraction that was more common than significant expansion. The values for directionality range from 0.00 (evenly contracting and expanding) to ±1.00 (all contracting or expanding). This way of looking at the data provided a method for comparing changes across different systems immune monitoring settings (Figure 6). Rhinovirus subjects had small changes in samples over time with a median of 0.019% and an interquartile range (IQR) of 0.0006%. Rhinovirus also had large directionality across all subjects either up or down, with a median of 0.029 and an IQR of 2.00. Thus, rhinovirus infections resulted in an extremely low magnitude of change, as very rare cell subsets were responding, and the direction of this change was typically fairly high or low for a given individual (i.e., the changes were not balanced and tended to represent marked expansion or contraction in the rare subsets that changed; Figure 6).

Figure 6 with 2 supplements see all

Download asset Open asset

T-REX outperforms other algorithms based on accuracy and speed

The T-REX workflow combines multiple algorithms, and each step represents a choice of tool and settings. T-REX was next tested against (1) T-REX with different choices of dimensionality reduction or no dimensionality reduction, (2) T-REX with a different clustering method from KNN, and (3) other algorithms. The results were evaluated based on speed (overall runtime, including any dimensionality reduction), whether regions of significant change were detected, and whether regions of change were also tetramer+ hotspots (Figure 7).

Figure 7 with 1 supplement see all

Download asset Open asset

Throughout this work, we have presented a version of T-REX that utilizes UMAP for dimensionality reduction and KNN for clustering. This standard version of T-REX ran in 0.35 hr on the 1.3 million cell dataset from RV001 (Figure 7, left column) and identified a cluster that was 22% tetramer+. An altered version of T-REX, in which t-SNE replaced UMAP as the method of dimensionality reduction, resulted in a much longer runtime of 23 hr and comparable capture of areas of change, containing 22% tetramer+ cells (Figure 7). Finally, we tested T-REX without dimensionality reduction, running KNN directly on the 14 original surface protein features previously utilized when performing either UMAP or t-SNE dimensionality reduction. This KNN on the original feature space was both slower (1.7 hr) and less accurate, in that it identified too many putative regions of change, most of which were not tetramer+ hotspots (Figure 7).

We also compared T-REX with other commonly used dimensionality reduction and clustering algorithms (Weber and Robinson, 2016), including SPADE, Phenograph, FlowSOM, and Citrus. These algorithms were designed to run on the original measured features and not on t-SNE or UMAP axes, although we have previously used FlowSOM on t-SNE as part of the RAPID algorithm (Leelatian et al., 2020). When SPADE was tested using the 14 original features, no regions of change ≥95% were identified (Figure 7). Of the 160 cell clusters requested and produced in SPADE, only two clusters had a percent change between 85% and 90% and contained 7% and 10% tetramer+ cells. The other clusters changed little (average of 53% ± 10.5%). Next, the KNN-based Phenograph algorithm was tested using the 14 original features and a k-value of 60 (as in T-REX, Figure 4). This test of Phenograph identified 13 clusters, none of which displayed significant change (average of 49% ± 6.4%) and none of which contained >5% tetramer+ cells. Phenograph was also tested using the UMAP coordinates and a k-value of 60 (as with standard T-REX in Figure 2). This approach identified 113 clusters, one of which was between 85% and 95% change, 8% tetramer+, and contained cells that were phenotypically similar to those identified by T-REX using MEM (Figure 7). The remaining clusters identified by Phenograph did not significantly change (average of 50% ± 6.6% change in the KNN region). Next, FlowSOM was tested both on UMAP axes and original features using an optimized cluster number, as determined by RAPID (Leelatian et al., 2020). RAPID cluster optimization (testing 0–90 clusters) ran in 0.93 hr and identified an optimal value of 66 clusters. FlowSOM did not identify regions of change using either UMAP or original features (average of 50% ± 5%) and did not identify clusters with >5% tetramer+ cells. Finally, we used 10 subsamples of each day to test Citrus on RV001 as Citrus requires two groups and several samples per group. This test of Citrus required 3 hr and found 156 total clusters, with >30 being significantly associated with either day 0 or day 7 abundance. Of these significant clusters, only one had >5% tetramer+ cells. This cluster also closely matched the T-REX tetramer+ cell phenotype (MEM label: ▲ CD45RO⁺¹⁰ CD38⁺⁸ ICOS⁺⁷ CD25⁺⁷ Ki67⁺⁷ CD27⁺⁶ CCR5⁺⁶ PD1⁺⁶ Tcf-1⁺⁶ CCR6⁺⁵ CD95⁺⁵ CCR7⁺⁵ CXCR3⁺⁴ CXCR5⁺⁴ CD127⁺³ Tbet⁺²). Importantly, phenotypes were comparable for all methods that identified at least one hotspot of >85% change, as quantified by RMSD analysis of MEM labels (Greenplate et al., 2019) (ranged from 97.1 to 98.7, with 100 being the highest degree of similarity).

Finally, as an additional test of T-REX stability, sample RV001 was subsampled into 10 groups, each containing 10% of the cells. T-REX was then run using both the same original (common) UMAP (Figure 7—figure supplement 1) and subsample-specific (unique) UMAPs. In all cases of subsampling, T-REX consistently identified tetramer+ hotspots based on ≥90 % or ≥ 95 % change across the 10 subsamples. To summarize, T-REX using either UMAP or t-SNE performed the best at identifying regions of change that were greatly enriched for rare tetramer+ cells, and T-REX with UMAP was faster than T-REX with t-SNE (Figure 7). Thus, T-REX with UMAP and MEM provides an automated approach to quickly reveal and quantitatively characterize cells in regions of great expansion or contraction.

A signature feature of the immune system is the ability of rare cells to respond to a stimulus by activating and proliferating, leading to rapid expansion of highly specialized cells that may share both a distinct phenotype and a clonal origin. Here, we report the design of a new algorithm that can identify and reliably phenotype biologically relevant cell types, including very rare cells, which respond in human disease. The T-REX algorithm was designed to capture phenotypic regions where significant change was occurring between a pair of samples from one individual. The fact that T-REX was able to identify the phenotype of cells whose regions were greatly enriched for virus-specific T cells in rhinovirus Dataset 1 highlighted its ability to pinpoint rare cells responding to infection and closely matches what would be expected based on a current understanding of rhinovirus immunology (Muehling et al., 2020; Muehling et al., 2018). Specifically, after rhinovirus challenge, expanded regions displayed molecular signatures consistent with activated memory (CD45RO + CD38+) and tissue trafficking (e.g., CCR5 enrichment in the MEM labels) that aligned with our previous findings for rhinovirus-specific CD4+ T cells using manual gating methods and a limited marker panel (Muehling et al., 2020; Wisniewski et al., 2018). The algorithm also reliably identified memory phenotypes of cells responding to rhinovirus infection, thereby revealing its potential to track the evolution of memory responses in vivo, in addition to defining candidate signatures that might be probed in functional assays. Although comparative phenotyping across time was beyond the scope of this study, it will be of high interest in the future to determine whether the vector of change in specific subsets correlates with additional aspects of disease or complicating host factors, such as allergy and asthma.

The T-REX algorithm also revealed potential new research directions, as there were cells that one might predict would be virus-specific, based on the T-REX enrichment analysis, but which were not enriched for the specific tetramers available here (e.g., Figure 3, RV007). Genetic analysis for the clonal origin of cells in such regions might help to determine whether these cells correspond to a clonal response for which a tetramer was not available, or else another type of CD4⁺ T cell response that may or may not be related to rhinovirus infection, such as a ‘bystander’ response. Additionally, it will be important to test whether this type of finding holds true for other well-studied viruses for which tetramers are available, such as influenza (Turner et al., 2020), and whether these findings extend to MHC class I tetramers and CD8 T cells. It was also striking that in the comparisons of day 7 to either day 0 (Figures 2 and 3) or day 28 (Figure 5), only the expanding cells (red) were in regions that were also tetramer hotspots. However, despite the focus on expansion in the T-REX acronym, contracting cells will likely also be of biological significance in different disease settings (as with AML) or potentially at different timepoints during the course of an infection, for example, as a result of egress from the circulation in the acute phase, or else transitions in memory and tissue-homing subsets that occur later. This aspect would also be expected to translate to different disease settings such as AML.

Extending the use of T-REX algorithm beyond rhinovirus further highlighted its ability to identify responding cells in a consistent manner across different subjects and different disease settings. Indeed, it is notable that regardless of the disease context the patient served as an effective baseline for comparison and allowed T-REX to find phenotypically similar cells in individuals with different starting immune profiles. A central question in systems immune monitoring is to place newly emerging diseases into the context of other well-studied diseases and immune responses. In working to compare COVID-19 and rhinovirus, it became clear that a summary of change indicating both the direction and magnitude of shifts was needed (Figure 6). This framework represents a way to summarize both broad populations of immune cells, like all CD45⁺ leukocytes, and key cell subpopulations, like CD4⁺ T cells. The striking changes observed in patients with moderate and severe COVID-19 were far beyond the subtle changes observed in individuals with rhinovirus and more closely matched the immune reprogramming observed in melanoma patients receiving checkpoint inhibitor therapy (Figure 6). A primary finding of T-REX analysis of Dataset 2 from the blood of COVID-19 patients was that some patients experienced very large changes in the blood, and that these changes were typically associated with more decreases than increases (Figure 6). This finding closely matches reported findings from others who observed a systematic reprogramming in many immune cell populations in severe COVID-19 patients (Mathew et al., 2020). Also observed were T cell subsets with enrichment of CD38, PD-1, and CD95, as has also been reported. While disease severity is not available for individual patients from Dataset 2, it is known that all these cases were at least moderate or severe (Rodriguez et al., 2020). It will be of interest to test the hypothesis that the more severe cases will be one of the two groups, either the patients with very little change and just expansion of cells, or those with more marked change and a general decrease of cells (Figure 6).

Notably, CD147/Basigin was highly expressed on many cells that changed during infection and was observed to change greatly on some populations over time. CD147 has been proposed in pre-prints as both a binding partner for SARS-CoV-2 spike protein and a potential mechanism of cellular entry, although evidence is needed to support this controversial hypothesis (Shilts and Wright, 2020). In the study of Dataset 2, the authors noted that immune responses were dominated by cells expressing CD38 and CD147 (Rodriguez et al., 2020). In the T-REX analysis of the same Dataset 2, for the cells that were changing greatly, CD147 was sometimes present on cells from day 0 that decreased greatly and was lower or absent on cells that emerged only at later times (Figure 6). An example of this was seen in cells from patient COV40, for which the authors noted CD147 expression on effector subsets at 1 week and onwards. The cells pinpointed by T-REX as emerging at day 6 included B cells that expressed CD147 (e.g., CX3CR1⁺⁸ CD9⁺⁸ CD29⁺⁸ CD147⁺⁵ IgD⁺³ CD99⁺³ CD33⁺¹ CD11c⁺¹ HLA-DR⁺¹ CD24⁺¹, Supplementary file 3), but the level of enrichment was lower than on myeloid cells from day 0 that decreased over time (e.g., CD147⁺⁸ CD29⁺⁶ CD55⁺⁵ CD38⁺⁵ CD99⁺⁴ CD64⁺³ CD62L⁺² CD45⁺¹ CD33⁺¹ CD14⁺¹, Supplementary file 3). This pattern of decreased enrichment of CD147 on cells emerging after day 0 was seen on other patients (Supplementary file 3) and is consistent with multiple explanations. Overall, there was a strong downward trend in many of the markers and cell subsets in COVID-19 patients, suggesting either selection against cells expressing a high level of proteins, downregulation of expression of key surface markers like CD147, expansion of immature or abnormal cells, or extreme trafficking of cells into tissues. These potential outcomes cannot be distinguished from each other with the analysis here. The utility of the T-REX algorithm is primarily in generating these hypotheses automatically and in pinpointing cells with extreme behavior within the context of the patient as their own baseline. Given the large amounts of change (Figure 6) and the generally lower numbers of T cell subsets observed in COVID-19 than in healthy individuals (Supplementary file 3), it may be the case that therapeutic stabilization of the immune system will be needed before virus-specific T cells will be identifiable with the T-REX method. It will be especially interesting to explore more mild cases of COVID-19 with this approach and determine whether the hotspots of change are truly virus-specific, analogous to the scenario with rhinovirus.

For the melanoma and AML cases presented here, the cohort sizes were not large enough to allow robust statistical comparison of patient response to degree or direction of change, although this information is available in the original studies (Ferrell et al., 2016). Of the AML patients, those with a high magnitude of change (Figure 6) were also those that had a high blast count and were complete responders to induction therapy, suggesting that the change represents the overall ‘reset’ of the immune system following chemotherapy. It will be of high interest to ask whether the identification of virus-specific T cells extends to populations of cells on checkpoint inhibitor therapy. The dynamics of regulatory cells may also be of interest, especially for autoimmunity, and it is possible, but not known, whether these cells will follow the same pattern as the CD4⁺ T cells in rhinovirus infection.

A major strength of the algorithm is that once cell regions of change are identified, the key features highlighted by T-REX and MEM can be used in lower parameter flow cytometry or imaging panels to provide further information, confirm findings, and physically isolate cells by FACS (Figure 2—figure supplement 3). Thus, low parameter cytometry approaches may rely more on manual analysis methods and cell signatures that are determined a priori, and T-REX may provide a useful tool for narrowing in on such features using exploratory high-dimensional data. The computational approach here emphasizes unsupervised UMAP and KNN clustering and uses statistical cutoffs to guide the analysis. Further optimization of the algorithm could include a stability testing analysis where the stochastic components of the algorithm are repeated to determine whether clusters or phenotypes are stable (Leelatian et al., 2020; Melchiotti et al., 2017). Overall, the unsupervised approach aims to diminish investigator bias and reveal novel or unexpected cell types. While unsupervised analysis tools have impacted high-dimensional cytometry for at least a decade (Davis et al., 2013; Becher et al., 2014; Bendall et al., 2011; Diggins et al., 2015; Saeys et al., 2016), T-REX is designed to capture both very rare and very common cell types and place them into a common context of immune change. The extremely rare T cells identified here were overlooked by other tools, likely due to these tools typically seeking clusters of cells representing at least 1% and generally more than 5% of the sample. This observation was expected and is consistent with past tests of algorithms on detection of rare subsets, which found that over half the algorithms tested returned poor F1 accuracy scores for detecting known populations that were 0.8% or 0.03% of the cells in a sample (Weber and Robinson, 2016). Notably, FlowSOM was a better performing algorithm in those examples. Here, FlowSOM did not effectively detect the rare 0.02% of virus specific cells in analyses of original features or dimensionality reduction results (Figure 7). An aspect of data structure in this study is the phenotypic similarity within the CD4+ T cell population, which means that the 0.02% of target cells is connected to a larger population and thus less apparent in the raw data. Thus, it is the threshold of change in T-REX removing the vast majority of cells (i.e., the cells that did not expand or contract ≥95%) that makes the virus-specific cells appear as separate clusters and greatly simplifies the task of distinguishing cluster groups.

T-REX was designed for use on a pair of data files, such as the pairs of timepoints tested here, but it can also be used for other paired comparisons. In the rhinovirus examples here, T-REX revealed virus-specific cells based on change over time only. However, the dataset timepoints in the rhinovirus study were selected carefully using prior knowledge to best reveal antigen-specific T cells (i.e., it was known that 7 days would be a good timepoint). In other disease settings, one may not know ideal timepoints and so it is not necessarily the case that T-REX will always reveal antigen-specific cells. It may also be the case that the level of change is so great that it obscures mechanistically significant rare cell subsets. This is an interpretation of the COVID-19 infection examples presented here, wherein a high degree of overall change was apparent in the immune system and many cells are changing that are not antigen-specific T cell subsets. This could potentially be due to timepoint selection or the natural history of the disease. However, we would hypothesize that COVID-19 vaccine responses might be more specific and localized to key cell subsets and, if this is the case, that T-REX would pinpoint key cell populations responding to the vaccine. For studies involving multiple timepoint tests, when it is possible to run a large common UMAP on a pool of the samples (i.e., there are not batch effects), T-REX can be used to do systematic paired comparisons of every individual against the pool. Data from different panels should not be used as the paired samples must have all channels in common. T-REX is also very sensitive to batch effects, and issues like these need to be addressed before analysis, or samples need to be analyzed in a batch-specific manner. We expect that one advantage of using MEM labels to summarize the phenotype of T-REX populations is that MEM labels have been shown to allow quantitative comparison across batches, instruments, and data types (Diggins et al., 2017). T-REX performance may depend on the number of cells and the algorithm was designed for studies with hundreds of thousands to millions of cells. Analyses with fewer than 1000 cells are unlikely to be productive, and analyses with fewer cells may require testing and adjusting k more than analyses with >100,000 cells. In testing subsamples, tetramer+ hotspots were still detected with 100,000 cells per sample (Figure 7—figure supplement 1). While equal subsampling is part of the algorithm of T-REX, subsampling tested here did not significantly impact the results. So while sampling and rarity of cells may impact the ability of T-REX to find subsets, the method was stable with far fewer cells (10%) than were originally planned as needed in the original study (Figure 7—figure supplement 1). The parameters that are chosen by the user are the k-value, the markers used to create the UMAP, the markers used for MEM, and the DBSCAN settings. Additional optimization of the k-value and cutoffs based on a formal statistical test is something that should be explored in further study of T-REX and related algorithms.

It is a central goal of systems immunology to map people with vastly contrasting immune system changes onto a common plot of change (as in Figure 6). The approach here goes beyond prior single measures of systematic change, such as Earth Mover’s Distance (Greenplate et al., 2019; Orlova et al., 2016), by including both direction and magnitude of change in one view of an individual’s immune response. This improvement proved useful for comparing settings with great change in many cell types (COVID-19 infection, AML chemotherapy responders) to settings with rare cells that specifically expanded or contracted (rhinovirus infection, melanoma checkpoint inhibitor therapy). This sensitivity of T-REX for extremely rare cells allowed the algorithm to reveal virus-specific CD4⁺ T cells without prior knowledge of their phenotype. T-REX should now be tested further to determine whether cells identified in SARS-CoV-2 also share a clonal origin. In addition, it is likely that T-REX will be useful beyond immunology settings in paired comparisons of quantitative single cell data, such as discovery screening or paired analysis of tumor cells.

Datasets analyzed in this manuscript are available online, including at FlowRepository. COVID-19 Dataset 2 (https://ki.app.box.com/s/sby0jesyu23a65cbgv51vpbzqjdmipr1), melanoma Dataset 3 (https://flowrepository.org/id/FR-FCM-ZYDG), and AML Dataset 4 (https://flowrepository.org/id/FR-FCM-ZZMC) were described and shared online in the associated manuscripts. Rhinovirus Dataset 1 is a newly generated dataset created at the University of Virginia available on FlowRepository (FR-FCM-Z2VX available at: https://flowrepository.org/id/FR-FCM-Z2VX). Transparent analysis scripts for all four datasets and all presented results are publicly available on the CytoLab Github page for T-REX [(https://github.com/cytolab/T-REX) copy archived at https://archive.softwareheritage.org/swh:1:rev:5e7ae8d512b36bef58b9e0df2b36a6775e82f734] and include open source code and commented Rmarkdown analysis walkthroughs.

1. 1. Barone SM

   (2020) FlowRepository

   ID FR-FCM-Z2VX. Rhinovirus data from: Unsupervised machine learning reveals key immune cell subsets in COVID-19, rhinovirus infection, and cancer therapy.

   https://flowrepository.org/id/FR-FCM-Z2VX

1. 1. Rodriguez L

   (2020) COVID-19_project Brodin lab accession

   ID sby0jesyu23a65cbgv51vpbzqjdmipr1. COVID-19 data from: Systems-Level Immunomonitoring from Acute to Recovery Phase of Severe COVID-19.

   https://ki.app.box.com/s/sby0jesyu23a65cbgv51vpbzqjdmipr1
2. 1. Greenplate AR

   (2016) FlowRepository

   ID FR-FCM-ZYDG. Melanoma data from: Mass Cytometry of Peripheral Blood from Melanoma Patients Receiving anti-PD-1.

   http://flowrepository.org/id/FR-FCM-ZYDG
3. 1. Ferrell PB

   (2016) FlowRepository

   ID FR-FCM-ZZMC. AML data from: High-Dimensional Analysis of Acute Myeloid Leukemia Reveals Phenotypic 684 Changes in Persistent Cells during Induction Therapy.

   http://flowrepository.org/id/FR-FCM-ZZMC

1. Software

   1. Barone SM
   2. Irish JM

   (2021a) T-REX (Tracking Responders EXpanding), version swh:1:rev:5e7ae8d512b36bef58b9e0df2b36a6775e82f734

   Software Heritage.

   https://archive.softwareheritage.org/swh:1:dir:7a6565a4ef9c1c888734e8f2665a8ad0d0aafbd0;origin=https://github.com/cytolab/T-REX;visit=swh:1:snp:153f999274aa8a1b4cd07b3030cdf91a0b9bf7a6;anchor=swh:1:rev:5e7ae8d512b36bef58b9e0df2b36a6775e82f734
2. Software

   1. Barone SM
   2. Diggins K
   3. Iris JM

   (2021b) MEM v3, version swh:1:rev:fc72a290c706c1268678b6300007eb59183af2f9

   Software Heritage.

   https://archive.softwareheritage.org/swh:1:rev:fc72a290c706c1268678b6300007eb59183af2f9;origin=https://github.com/cytolab/mem;visit=swh:1:snp:097b2ae3302b4f73f2cbc5457a9f62fbb4b94909
3. 1. Becher B
   2. Schlitzer A
   3. Chen J
   4. Mair F
   5. Sumatoh HR
   6. Teng KWW
   7. Low D
   8. Ruedl C
   9. Riccardi-Castagnoli P
   10. Poidinger M
   11. Greter M
   12. Ginhoux F
   13. Newell EW

   (2014) High-dimensional analysis of the murine myeloid cell system

   Nature Immunology 15:1181–1189.

   https://doi.org/10.1038/ni.3006

   • PubMed
   • Google Scholar
4. 1. Becht E
   2. McInnes L
   3. Healy J
   4. Dutertre CA
   5. Kwok IWH
   6. Ng LG
   7. Ginhoux F
   8. Newell EW

   (2018) Dimensionality reduction for visualizing single-cell data using UMAP

   Nature Biotechnology 1:nbt.4314.

   https://doi.org/10.1038/nbt.4314

   • PubMed
   • Google Scholar
5. 1. Belkina AC
   2. Ciccolella CO
   3. Anno R
   4. Halpert R
   5. Spidlen J
   6. Snyder-Cappione JE

   (2019) Automated optimized parameters for t-distributed stochastic neighbor embedding improve visualization and analysis of large datasets

   Nature Communications 10:5415.

   https://doi.org/10.1038/s41467-019-13055-y

   • PubMed
   • Google Scholar
6. 1. Bendall SC
   2. Simonds EF
   3. Qiu P
   4. Amir ED
   5. Krutzik PO
   6. Finck R
   7. Bruggner RV
   8. Melamed R
   9. Trejo A
   10. Ornatsky OI
   11. Balderas RS
   12. Plevritis SK
   13. Sachs K
   14. Pe’er D
   15. Tanner SD
   16. Nolan GP

   (2011) Single-cell mass cytometry of differential immune and drug responses across a human hematopoietic continuum

   Science 332:687–696.

   https://doi.org/10.1126/science.1198704

   • PubMed
   • Google Scholar
7. 1. Bruggner RV
   2. Bodenmiller B
   3. Dill DL
   4. Tibshirani RJ
   5. Nolan GP

   (2014) Automated identification of stratifying signatures in cellular subpopulations

   PNAS 111:E2770.

   https://doi.org/10.1073/pnas.1408792111

   • PubMed
   • Google Scholar
8. Preprint

   1. Burns TJ
   2. Nolan GP
   3. Samusik N

   (2018) Continuous Visualization of Differences between Biological Conditions in Single-Cell Data

   bioRxiv.

   https://doi.org/10.1101/337485

   • Google Scholar
9. 1. Chattopadhyay PK
   2. Gierahn TM
   3. Roederer M
   4. Love JC

   (2014) Single-cell technologies for monitoring immune systems

   Nature Immunology 15:128–135.

   https://doi.org/10.1038/ni.2796

   • PubMed
   • Google Scholar
10. 1. Davis KL
    2. Tadmor MD
    3. Simonds EF
    4. Levine JH
    5. Bendall SC
    6. Shenfeld DK
    7. Krishnaswamy S
    8. Nolan GP

    (2013) VISNE enables visualization of high dimensional single-cell data and reveals phenotypic heterogeneity of leukemia

    Nature Biotechnology 31:545–552.

    https://doi.org/10.1038/nbt.2594

    • PubMed
    • Google Scholar
11. 1. Davis MM
    2. Tato CM
    3. Furman D

    (2017) Systems immunology: Just getting started

    Nature Immunology 18:725–732.

    https://doi.org/10.1038/ni.3768

    • PubMed
    • Google Scholar
12. 1. Diggins KE
    2. Ferrell PB
    3. Irish JM

    (2015) Methods for discovery and characterization of cell subsets in high dimensional mass cytometry data

    Methods 82:55–63.

    https://doi.org/10.1016/j.ymeth.2015.05.008

    • PubMed
    • Google Scholar
13. 1. Diggins KE
    2. Greenplate AR
    3. Leelatian N
    4. Wogsland CE
    5. Irish JM

    (2017) Characterizing cell subsets using marker enrichment modeling

    Nature Methods 14:275–278.

    https://doi.org/10.1038/nmeth.4149

    • PubMed
    • Google Scholar
14. 1. Diggins KE
    2. Gandelman JS
    3. Roe CE
    4. Irish JM

    (2018) Generating quantitative cell identity labels with marker enrichment modeling (MEM). Current protocols in cytometry / editorial board

    J Paul Robinson 83:10.

    https://doi.org/10.1002/cpcy.34

    • Google Scholar
15. 1. Ferrell PB
    2. Diggins KE
    3. Polikowsky HG
    4. Mohan SR
    5. Seegmiller AC
    6. Irish JM

    (2016) High-dimensional analysis of acute myeloid leukemia reveals phenotypic changes in persistent cells during induction therapy

    PLOS ONE 11:e0153207.

    https://doi.org/10.1371/journal.pone.0153207

    • PubMed
    • Google Scholar
16. 1. Ferrer-Font L
    2. Mayer JU
    3. Old S
    4. Hermans IF
    5. Irish J
    6. Price KM

    (2020) High-dimensional data analysis algorithms yield comparable results for mass cytometry and spectral flow cytometry data

    Cytometry. Part A 97:824–831.

    https://doi.org/10.1002/cyto.a.24016

    • PubMed
    • Google Scholar
17. 1. Gandelman JS
    2. Byrne MT
    3. Mistry AM
    4. Polikowsky HG
    5. Diggins KE
    6. Chen H
    7. Lee SJ
    8. Arora M
    9. Cutler C
    10. Flowers M
    11. Pidala J
    12. Irish JM
    13. Jagasia MH

    (2019) Machine learning reveals chronic graft-versus-host disease phenotypes and stratifies survival after stem cell transplant for hematologic malignancies

    Haematologica 104:189–196.

    https://doi.org/10.3324/haematol.2018.193441

    • PubMed
    • Google Scholar
18. 1. Greenplate AR
    2. Johnson DB
    3. Ferrell PB
    4. Irish JM

    (2016a) Systems immune monitoring in cancer therapy

    European Journal of Cancer 61:77–84.

    https://doi.org/10.1016/j.ejca.2016.03.085

    • PubMed
    • Google Scholar
19. 1. Greenplate AR
    2. Johnson DB
    3. Roussel M
    4. Savona MR
    5. Sosman JA
    6. Puzanov I
    7. Ferrell PB
    8. Irish JM

    (2016b) Myelodysplastic syndrome revealed by systems immunology in a melanoma patient undergoing anti-pd-1 therapy

    Cancer Immunology Research 4:474–480.

    https://doi.org/10.1158/2326-6066.CIR-15-0213

    • PubMed
    • Google Scholar
20. 1. Greenplate AR
    2. McClanahan DD
    3. Oberholtzer BK
    4. Doxie DB
    5. Roe CE
    6. Diggins KE
    7. Leelatian N
    8. Rasmussen ML
    9. Kelley MC
    10. Gama V
    11. Siska PJ
    12. Rathmell JC
    13. Ferrell PB
    14. Johnson DB
    15. Irish JM

    (2019) Computational Immune Monitoring Reveals Abnormal Double-Negative T Cells Present across Human Tumor Types

    Cancer Immunology Research 7:86–99.

    https://doi.org/10.1158/2326-6066.CIR-17-0692

    • PubMed
    • Google Scholar
21. 1. Irish JM
    2. Kotecha N
    3. Nolan GP

    (2006) Mapping normal and cancer cell signalling networks: Towards single-cell proteomics

    Nature Reviews. Cancer 6:146–155.

    https://doi.org/10.1038/nrc1804

    • PubMed
    • Google Scholar
22. 1. Irish JM

    (2014) Beyond the age of cellular discovery

    Nature Immunology 15:1095–1097.

    https://doi.org/10.1038/ni.3034

    • PubMed
    • Google Scholar
23. Website

    1. Krijthe J
    2. van der Maaten R
    3. Rtsne L

    (2015) L. T-distributed stochastic neighbor embedding using a barnes-hutimplementation

    Accessed March 20, 2020.

    https://cran.r-project.org/web/packages/Rtsne/citation.html
24. 1. Lakshmikanth T
    2. Muhammad SA
    3. Olin A
    4. Chen Y
    5. Mikes J
    6. Fagerberg L
    7. Gummesson A
    8. Bergström G
    9. Uhlen M
    10. Brodin P

    (2020) Human immune system variation during 1 year

    Cell Reports 32:107923.

    https://doi.org/10.1016/j.celrep.2020.107923

    • PubMed
    • Google Scholar
25. 1. Leelatian N
    2. Sinnaeve J
    3. Mistry AM
    4. Barone SM
    5. Brockman AA
    6. Diggins KE
    7. Greenplate AR
    8. Weaver KD
    9. Thompson RC
    10. Chambless LB
    11. Mobley BC
    12. Ihrie RA
    13. Irish JM

    (2020) Unsupervised machine learning reveals risk stratifying glioblastoma tumor cells

    eLife 9:e56879.

    https://doi.org/10.7554/eLife.56879

    • PubMed
    • Google Scholar
26. 1. Levine JH
    2. Simonds EF
    3. Bendall SC
    4. Davis KL
    5. Amir ED
    6. Tadmor MD
    7. Litvin O
    8. Fienberg HG
    9. Jager A
    10. Zunder ER
    11. Finck R
    12. Gedman AL
    13. Radtke I
    14. Downing JR
    15. Pe’er D
    16. Nolan GP

    (2015) Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis

    Cell 162:184–197.

    https://doi.org/10.1016/j.cell.2015.05.047

    • PubMed
    • Google Scholar
27. 1. Mathew D
    2. Giles JR
    3. Baxter AE
    4. Oldridge DA
    5. Greenplate AR
    6. Wu JE
    7. Alanio C
    8. Kuri-Cervantes L
    9. Pampena MB
    10. D’Andrea K
    11. Manne S
    12. Chen Z
    13. Huang YJ
    14. Reilly JP
    15. Weisman AR
    16. Ittner CAG
    17. Kuthuru O
    18. Dougherty J
    19. Nzingha K
    20. Han N
    21. Kim J
    22. Pattekar A
    23. Goodwin EC
    24. Anderson EM
    25. Weirick ME
    26. Gouma S
    27. Arevalo CP
    28. Bolton MJ
    29. Chen F
    30. Lacey SF
    31. Ramage H
    32. Cherry S
    33. Hensley SE
    34. Apostolidis SA
    35. Huang AC
    36. Vella LA
    37. UPenn COVID Processing Unit
    38. Betts MR
    39. Meyer NJ
    40. Wherry EJ

    (2020) Deep immune profiling of COVID-19 patients reveals distinct immunotypes with therapeutic implications

    Science 369:abc8511.

    https://doi.org/10.1126/science.abc8511

    • PubMed
    • Google Scholar
28. Preprint

    1. McInnes L
    2. Healy J
    3. Melville J

    (2018) UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction

    arXiv.

    https://arxiv.org/abs/1802.03426

    • Google Scholar
29. 1. Melchiotti R
    2. Gracio F
    3. Kordasti S
    4. Todd AK
    5. de Rinaldis E

    (2017) Cluster stability in the analysis of mass cytometry data

    Cytometry. Part A 91:73–84.

    https://doi.org/10.1002/cyto.a.23001

    • PubMed
    • Google Scholar
30. 1. Mistry AM
    2. Greenplate AR
    3. Ihrie RA
    4. Irish JM

    (2019) Beyond the message: Advantages of snapshot proteomics with single-cell mass cytometry in solid tumors

    The FEBS Journal 286:1523–1539.

    https://doi.org/10.1111/febs.14730

    • PubMed
    • Google Scholar
31. 1. Muehling LM
    2. Mai DT
    3. Kwok WW
    4. Heymann PW
    5. Pomés A
    6. Woodfolk JA

    (2016) Circulating memory CD4+ T cells target conserved epitopes of rhinovirus capsid proteins and respond rapidly to experimental infection in humans

    Journal of Immunology 197:3214–3224.

    https://doi.org/10.4049/jimmunol.1600663

    • PubMed
    • Google Scholar
32. 1. Muehling LM
    2. Turner RB
    3. Brown KB
    4. Wright PW
    5. Patrie JT
    6. Lahtinen SJ
    7. Lehtinen MJ
    8. Kwok WW
    9. Woodfolk JA

    (2018) Single-Cell Tracking Reveals a Role for Pre-Existing CCR5+ Memory Th1 Cells in the Control of Rhinovirus-A39 After Experimental Challenge in Humans

    The Journal of Infectious Diseases 217:381–392.

    https://doi.org/10.1093/infdis/jix514

    • PubMed
    • Google Scholar
33. 1. Muehling LM
    2. Heymann PW
    3. Wright PW
    4. Eccles JD
    5. Agrawal R
    6. Carper HT
    7. Murphy DD
    8. Workman LJ
    9. Word CR
    10. Ratcliffe SJ
    11. Capaldo BJ
    12. Platts-Mills TAE
    13. Turner RB
    14. Kwok WW
    15. Woodfolk JA

    (2020) Human Th1 and Th2 cells targeting rhinovirus and allergen coordinately promote allergic asthma

    The Journal of Allergy and Clinical Immunology 146:555–570.

    https://doi.org/10.1016/j.jaci.2020.03.037

    • PubMed
    • Google Scholar
34. 1. Orlova DY
    2. Zimmerman N
    3. Meehan S
    4. Meehan C
    5. Waters J
    6. Ghosn EEB
    7. Filatenkov A
    8. Kolyagin GA
    9. Gernez Y
    10. Tsuda S
    11. Moore W
    12. Moss RB
    13. Herzenberg LA
    14. Walther G

    (2016) Earth Mover’s Distance (EMD): A True Metric for Comparing Biomarker Expression Levels in Cell Populations

    PLOS ONE 11:e0151859.

    https://doi.org/10.1371/journal.pone.0151859

    • PubMed
    • Google Scholar
35. 1. Park LM
    2. Lannigan J
    3. Jaimes MC

    (2020) OMIP-069: Forty-color full spectrum flow cytometry panel for deep immunophenotyping of major cell subsets in human peripheral blood

    Cytometry. Part A 97:1044–1051.

    https://doi.org/10.1002/cyto.a.24213

    • PubMed
    • Google Scholar
36. 1. Pyne S
    2. Lee SX
    3. Wang K
    4. Irish J
    5. Tamayo P
    6. Nazaire M-D
    7. Duong T
    8. Ng S-K
    9. Hafler D
    10. Levy R
    11. Nolan GP
    12. Mesirov J
    13. McLachlan GJ

    (2014) Joint modeling and registration of cell populations in cohorts of high-dimensional flow cytometric data

    PLOS ONE 9:e100334.

    https://doi.org/10.1371/journal.pone.0100334

    • PubMed
    • Google Scholar
37. 1. Qiu P
    2. Simonds EF
    3. Bendall SC
    4. Gibbs KD Jr
    5. Bruggner RV
    6. Linderman MD
    7. Sachs K
    8. Nolan GP
    9. Plevritis SK

    (2011) Extracting a cellular hierarchy from high-dimensional cytometry data with spade

    Nature Biotechnology 29:886–891.

    https://doi.org/10.1038/nbt.1991

    • PubMed
    • Google Scholar
38. 1. Ragab D
    2. Salah Eldin H
    3. Taeimah M
    4. Khattab R
    5. Salem R

    (2020) The COVID-19 Cytokine Storm; What We Know So Far

    Frontiers in Immunology 11:1446.

    https://doi.org/10.3389/fimmu.2020.01446

    • PubMed
    • Google Scholar
39. Preprint

    1. Rodriguez L
    2. Pekkarinen P
    3. Tadepally LK
    4. Tan Z
    5. Rosat Consiglio C
    6. Pou C
    7. Chen Y
    8. Habimana Mugabo C
    9. Nguyen Quoc A
    10. Nowlan K

    (2020) Systems-Level Immunomonitoring from Acute to Recovery Phase of Severe Covid-19

    medRxiv.

    https://doi.org/10.1101/2020.06.03.20121582

    • Google Scholar
40. 1. Saeys Y
    2. Van Gassen S
    3. Lambrecht BN

    (2016) Computational flow cytometry: Helping to make sense of high-dimensional immunology data

    Nature Reviews. Immunology 16:449–462.

    https://doi.org/10.1038/nri.2016.56

    • PubMed
    • Google Scholar
41. 1. Schultze JL

    (2015) Teaching “big data” analysis to young immunologists

    Nature Immunology 16:902–905.

    https://doi.org/10.1038/ni.3250

    • PubMed
    • Google Scholar
42. Preprint

    1. Shilts J
    2. Wright GJ

    (2020) No Evidence for Basigin/CD147 as a Direct Sars-Cov-2 Spike Binding Receptor

    bioRxiv.

    https://doi.org/10.1101/2020.07.25.221036

    • Google Scholar
43. 1. Spidlen J
    2. Breuer K
    3. Rosenberg C
    4. Kotecha N
    5. Brinkman RR

    (2012) Flowrepository: A resource of annotated flow cytometry datasets associated with peer-reviewed publications

    Cytometry. Part A 81:727–731.

    https://doi.org/10.1002/cyto.a.22106

    • PubMed
    • Google Scholar
44. 1. Spitzer MH
    2. Nolan GP

    (2016) Mass cytometry: Single cells, many features

    Cell 165:780–791.

    https://doi.org/10.1016/j.cell.2016.04.019

    • PubMed
    • Google Scholar
45. 1. Turner JS
    2. Lei T
    3. Schmitz AJ
    4. Day A
    5. Choreño-Parra JA
    6. Jiménez-Alvarez L
    7. Cruz-Lagunas A
    8. House SL
    9. Zúñiga J
    10. Ellebedy AH
    11. Mudd PA

    (2020) Impaired cellular immune responses during the first week of severe acute influenza infection

    The Journal of Infectious Diseases 222:1235–1244.

    https://doi.org/10.1093/infdis/jiaa226

    • PubMed
    • Google Scholar
46. 1. Van Gassen S
    2. Callebaut B
    3. Van Helden MJ
    4. Lambrecht BN
    5. Demeester P
    6. Dhaene T
    7. Saeys Y

    (2015) FlowSOM: Using self-organizing maps for visualization and interpretation of cytometry data

    Cytometry Part A: The Journal of the International Society for Analytical Cytology 87:636–645.

    https://doi.org/10.1002/cyto.a.22625

    • PubMed
    • Google Scholar
47. 1. Weber LM
    2. Robinson MD

    (2016) Comparison of clustering methods for high-dimensional single-cell flow and mass cytometry data

    Cytometry. Part A 89:1084–1096.

    https://doi.org/10.1002/cyto.a.23030

    • PubMed
    • Google Scholar
48. 1. Wisniewski JA
    2. Muehling LM
    3. Eccles JD
    4. Capaldo BJ
    5. Agrawal R
    6. Shirley DA
    7. Patrie JT
    8. Workman LJ
    9. Schuyler AJ
    10. Lawrence MG
    11. Teague WG
    12. Woodfolk JA

    (2018) Th1 signatures are present in the lower airways of children with severe asthma, regardless of allergic status

    The Journal of Allergy and Clinical Immunology 141:2048–2060.

    https://doi.org/10.1016/j.jaci.2017.08.020

    • PubMed
    • Google Scholar

Author details

1. Sierra M Barone

   1. Department of Cell and Developmental Biology, Vanderbilt University, Nashville, United States
   2. Vanderbilt-Ingram Cancer Center, Vanderbilt University Medical Center, Nashville, United States

   Contribution

   Conceptualization, Formal analysis, Software, Validation, Visualization, Writing – original draft, Writing – review and editing

   Contributed equally with

   Alberta GA Paul and Lyndsey M Muehling

   Competing interests

   none

   "This ORCID iD identifies the author of this article:" 0000-0001-5944-750X

2. Alberta GA Paul

   Allergy Division, Department of Medicine, University of Virginia School of Medicine, Charlottesville, United States

   Contribution

   Data curation, Methodology, Writing – review and editing

   Contributed equally with

   Sierra M Barone and Lyndsey M Muehling

   Competing interests

   became an employee of Cytek Biosciences, Inc after performing this research at University of Virginia.

   "This ORCID iD identifies the author of this article:" 0000-0002-9318-3760

3. Lyndsey M Muehling

   1. Allergy Division, Department of Medicine, University of Virginia School of Medicine, Charlottesville, United States
   2. Department of Microbiology, Immunology, and Cancer Biology, University of Virginia School of Medicine, Charlottesville, United States

   Contribution

   Data curation, Formal analysis, Visualization, Writing – review and editing

   Contributed equally with

   Sierra M Barone and Alberta GA Paul

   Competing interests

   None

   "This ORCID iD identifies the author of this article:" 0000-0003-3203-3264

4. Joanne A Lannigan

   Department of Microbiology, Immunology, and Cancer Biology, University of Virginia School of Medicine, Charlottesville, United States

   Contribution

   Data curation, Methodology, Writing – review and editing

   Competing interests

   became a paid consultant of Cytek Biosciences, Inc after performing this research at University of Virginia.

   "This ORCID iD identifies the author of this article:" 0000-0002-3981-8681

5. William W Kwok

   Benaroya Research Institute at Virginia Mason, Seattle, United States

   Contribution

   Methodology, Writing – review and editing

   Competing interests

   none

   "This ORCID iD identifies the author of this article:" 0000-0003-4843-4599

6. Ronald B Turner

   Department of Pediatrics, University of Virginia School of Medicine, Charlottesville, United States

   Contribution

   Methodology, Writing – review and editing

   Competing interests

   none

7. Judith A Woodfolk

   1. Allergy Division, Department of Medicine, University of Virginia School of Medicine, Charlottesville, United States
   2. Department of Microbiology, Immunology, and Cancer Biology, University of Virginia School of Medicine, Charlottesville, United States

   Contribution

   Data curation, Investigation, Methodology, Project administration, Resources, Supervision, Writing – review and editing

   For correspondence

   jaw4m@virginia.edu

   Competing interests

   none

   "This ORCID iD identifies the author of this article:" 0000-0002-8915-4334

8. Jonathan M Irish

   1. Department of Cell and Developmental Biology, Vanderbilt University, Nashville, United States
   2. Vanderbilt-Ingram Cancer Center, Vanderbilt University Medical Center, Nashville, United States
   3. Department of Pathology, Microbiology and Immunology, Vanderbilt University Medical Center, Nashville, United States

   Contribution

   Conceptualization, Formal analysis, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review and editing

   For correspondence

   jonathan.irish@vanderbilt.edu

   Competing interests

   was a co-founder and a board member of Cytobank Inc and received unrelated research support from Incyte Corp, Janssen, and Pharmacyclics.

   "This ORCID iD identifies the author of this article:" 0000-0001-9428-8866

• 3,073

  views

• 345

  downloads

• 16

  citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.