Computational Immune Monitoring Reveals Abnormal Double-Negative T Cells Present across Human Tumor Types Free

The immune system is a complex network of cells and tissue types, and it is increasingly important to simultaneously track cell subsets and understand the system as a whole. Longitudinal monitoring of changes in the immune system has provided insight into drug response and disease progression (1–3). Differences in response to perturbation can stratify clinical outcome (4, 5) and indicate mechanism of action (6). Challenges to the immune system, such as vaccination, infection, surgical intervention, or the emergence of a malignancy, can elicit detectable changes above the relatively stable basal state of each individual (7, 8). Cytomic approaches that can characterize all the cells in a system, such as high-dimensional flow and mass cytometry, are an area of active development in immunology (9). A common framework for data analysis will allow researchers using these cytomic approaches to compare and contrast immune systems from diverse research areas, including tumor immunology and treatment response (4, 10, 11), blood cancer (12, 13), bone marrow failure (14, 15), and human immune variation and autoimmunity (3, 7, 16).

Systematically monitoring the cells of the immune system and their features is especially powered when serial samples from an individual can be used as comparison points. However, this approach generates vast amounts of data. Simultaneously tracking entire systems and the parts that comprise them can become overwhelming in the context of clinical research, where cells of interest can be rare (4, 10), and the entire system can change quickly in unanticipated ways (15). Thus, in cytomic, system-level studies, the data analysis strategy is as important as the experiment design (4, 17, 18), and it is vital to track known reference populations and place new observations in the context of prior knowledge (19, 20). Tools exist for visualizing and gating (21–24), supervised population and biomarker analysis (25, 26), and describing cell population identity (19) within high-dimensional data sets. In contrast, a great need exists for automated analysis for the steps immediately after population identification (gating). Tools for cellular clustering are especially numerous, specific for each data type, and have been extensively addressed in prior work (27). This ensemble aims to quantify changes in the whole system, prior to cellular clustering, and to characterize how cellular populations change after the appropriate clustering tool has been applied to the system.

The growth and success of high-dimensional single-cell technologies relies on the ongoing development of data analysis tools needed to parse large amounts of data and place results into context (18–20). It remains relatively rare for researchers to choose data analysis approaches that explicitly incorporate time or other changes. Across longitudinal studies, four dynamic elements stand out as key features of cellular systems: (i) plasticity or stability of the system as a whole, (ii) changes in abundance of cell subsets, (iii) emergence of unexpected features on cell subsets, and (iv) emergence of novel or unexpected cell types. The central goal here was to simultaneously use distinct tools focused on these different data types as an ensemble and to organize this copious information into an automated "first pass" analysis that could be easily interpretable by an immunologist and that would highlight cell subsets for in-depth review. The components selected for the ensemble toolkit here included the Earth Mover's Distance (EMD) algorithm (28), t-distributed stochastic neighbor embedding (t-SNE; ref. 21), and Marker Enrichment Modeling (MEM; ref. 19).

To validate and challenge computational analysis tools, it is valuable to explore immunology problems representing differing levels of prior knowledge, amounts of change, and abundances of target cell types. Here, the first challenge (data set 1) comprised data from melanoma patients undergoing anti–PD-1 therapy with pembrolizumab. Data set 1 was chosen for the relevance of this therapeutic strategy and because the data set includes an unusual case from a patient with a distinct, unexpected immune system trajectory (15). This example was chosen to include an outlier case, which is unusual enough that training data sets would not normally include an example of it. Data set 2 comprised of a previously published data set of peripheral blood from acute myeloid leukemia (AML) patients undergoing chemotherapy (12). Data set 2 was chosen as an example of dynamic cellular populations shifting dramatically over therapy. This data set was also included to represent the challenge of tracking and characterizing treatment-refractory leukemic blasts, which did not converge on a single phenotype and, instead, shifted into different phenotypic compartments, none of which matched the phenotype of healthy cells (12). Data set 3 included serial melanoma tumor biopsies from patients treated with dabrafenib, a BRAF^V600 inhibitor (BRAFi), and trametinib, a MEK inhibitor (MEKi). A challenge of data set 3 was to apply multiple tools in a system that was relatively less well studied and included a diverse set of individuals and mass cytometry panels. By providing both high-level and detailed views of cellular systems changing over time in human patients, the ensemble approach revealed knowledge about immune system interactions in these three study types with contrasting changes and challenges.

In the case of peripheral blood collected from melanoma patients receiving pembrolizumab, the purpose of the study was to identify biological characteristics of melanoma occurring prior to treatment and at different time points following therapy for patients being treated with immune-based therapies. Each patient gave consent to an institutional review board (IRB)–approved research protocol. To be included in this study, patients met the following inclusion criteria: (i) pathologically proven diagnosis of melanoma, (ii) 18 years of age or older, (iii) treated with immune-based therapies, and (iv) willing to have several serial blood draws. Patients not receiving immune-based therapy or unwilling or unable to provide consent were not included. There was no blinding or randomization process in this study. Blood specimens were obtained from patients during the time of scheduled phlebotomy for routine clinical laboratory analysis. Peripheral blood draws were done on the day of therapy start and 21 days (±10 days), 84 days (±21 days), and 180 days (±21 days) following initiation of therapy.

This cohort has been previously published and includes peripheral blood from AML patients undergoing chemotherapy (12). Briefly, patients consented to a blood draw through an IRB-approved research protocol. Peripheral blood was drawn pretreatment at the time of diagnosis, every 2 to 3 days within the first 2 weeks after the start of chemotherapy, at the 2-week time point, and at the time of recovery, if applicable.

In the case of tumors sampled sequentially from melanoma patients treated with targeted therapy, the objective was to identify biomarkers of response and resistance to B-RAF– and MEK-targeted therapy in melanoma. Patients with advanced, operable BRAF mutation–positive melanoma received GSK-2118436 (BRAF inhibitor) for 2 weeks, followed by the combination of GSK-2118436 and GSK-1120212 (MEK inhibitor) for 2 weeks, followed by surgical resection of the disease. Tumor biopsies were obtained prior to start of therapy and 2 weeks after combined GSK-2118436 and GSK-1120212 (29). To be included in this study, patients met the following inclusion criteria: (i) signed written, informed consent, (ii) between the ages of 18 and 90, (iii) patients with locally or regionally advanced melanoma being considered for resection of the lesion(s) for local–regional control and potential cure, (iv) BRAF V600 mutation–positive by snapshot molecular analysis, (v) measurable disease, (vi) all prior treatment-related toxicities CTCAE ≤ grade 1 at the time of enrollment, (vii) adequate baseline organ function, (viii) women of childbearing potential with a negative serum pregnancy test within 14 days of the first dose of study or men with female partner of childbearing potential must have had either had a prior vasectomy or agree to use effective contraception, and (ix) able to swallow and retain oral medication. There was no blinding or randomization process in this study.

All human samples were obtained in accordance with the Declaration of Helsinki following protocols approved by Vanderbilt University Medical Center IRB. The patient information for unpublished samples can be found in Supplementary Table S1. Healthy donor tonsil, adenoid, and blood were collected as “non-human subjects,” without gender or age information. Upon single-cell isolation, all cells were cryopreserved in 88% fetal calf serum plus 12% DMSO. Cells from human samples were collected and isolated as follows.

Peripheral blood mononuclear cells (PBMC) were collected, isolated, and cryopreserved from approximately 20 mL of freshly drawn blood as previously described (15). Briefly, peripheral blood was drawn into sodium heparin anticoagulant, and PBMCs were isolated by centrifugation after layering on top of a Ficoll-Paque PLUS (GE Healthcare Bio-Sciences) gradient.

Melanoma tumors, glioblastoma tumors, and nonmalignant human adenoid and tonsil tissue were resected from patients with consent and in accordance with the Declaration of Helsinki. All solid tissue samples were dissociated into live, single-cell suspensions, and samples were cryopreserved using a previously documented protocol (30, 31). Solid tissue samples were first manually dissociated using a scalpel. The minced tissue was then incubated in RPMI 1640 (Corning/Mediatech) plus 10% FBS, collagenase II (1 mg/mL; Sigma-Aldrich), and DNase (0.25 mg/mL; Sigma-Aldrich) for 1 hour in a 37°C incubator with 5% CO₂. Cells were then strained through a 70-μm and 40-μm filter prior to cryopreservation.

Renal cell carcinoma samples were processed and stored as described by Siska and colleagues (32). Briefly, malignant tissue was removed from patients and processed by mechanical dissociation and enzymatic digestion. Single-cell suspensions were achieved by passing dissociated tissue through a 70-μm strainer. Finally, red blood cells were lysed, and the single-cell suspension was cryopreserved.

Reprogramming of human neonatal foreskin fibroblast cells (strain BJ; ATCC no. CRL2522) was induced by transduction with CytoTune Sendai virus (Life Technologies). All experiments were performed under the supervision of the Vanderbilt Institutional Human Pluripotent Cell Research Oversight Committee. Induced pluripotent stem cells were grown in feeder-free conditions in plates coated with Matrigel (BD Biosciences) and maintained in mTESR1 media (Stem Cell Technologies) at 37°C with 5% CO₂. iPSCs were generated in the lab of Vivian Gama, PhD and were passaged 30 to 35 times prior to staining by mass cytometry. The iPSCs were first characterized by live staining with Tra1-60 or Tra1-81 antibodies that recognize undifferentiated iPSCs. Every 30 passages, iPSCs were characterized by gene-expression analysis using the TaqMan hPSC score card panel (Life Technologies) and karyotyping. Cells were checked daily for differentiation and were passaged every 3 to 4 days using Gentle dissociation solution (Stem Cell Technologies). iPSCs were treated with 0.5% EDTA prior to staining with mass cytometry antibody panel described below. iPSCs used in this study were mycoplasma negative.

Thawed samples were first incubated with cisplatin (25 μmol/L, Enzo Life Sciences). After incubation with cisplatin, cells were washed in PBS containing 1% BSA. Staining occurred in 50 μL PBS/1% BSA for 30 minutes at room temperature using the antibodies listed in Supplementary Table S2. Cells were then washed twice with PBS/1% BSA and fixed with a final concentration of 1.6% paraformaldehyde (PFA, Electron Microscopy Sciences). Cells were washed again, using PBS, and then resuspended in ice-cold methanol to permeabilize. Cells were incubated at −20°C overnight before being washed twice in PBS and stained with iridium DNA intercalator (Fluidigm Sciences). Purified, carrier-free antibodies were purchased from the listed provider and labeled with the listed metal using the protocol provided by Fluidigm. Stained samples were collected at Vanderbilt University Flow Cytometry Shared Resource on a CyTOF 1.0 mass cytometer (Fluidigm Sciences). All events were normalized prior to analysis using Fluidigm normalization beads.

Data (FCS files) were collected and stored in the online analysis platform Cytobank. Data analysis was performed in Cytobank and statistical programming environment R (version 3.4.0) via R Studio.

EMD was calculated between each pair of populations using the "transport" library for R (refs. 28, 33; https://cran.r-project.org/web/packages/transport/citation.html). The parent population (e.g., live CD45⁺ events) was gated in Cytobank, followed by the creation of a viSNE map in Cytobank. A viSNE analysis with two output dimensions was performed, equally sampling 5,000 events per file, with 1,000 iterations, perplexity equal to 30, and theta equal to 0.5. The events with their viSNE axes were then downloaded from Cytobank, and the EMD was calculated between each pair of files using the "transport" library for R. The "wpp" object was used to represent each set of points in the two viSNE axes, and the "wasserstein" function was called on each pair of point sets to produce a distance matrix. Each point was assigned unit weight.

Because calculating a matrix with the EMD between each set of 5,000 events from the viSNE analysis is computationally expensive, four optimizations were performed. (i) Each file was further down-sampled to 1,000 out of the original 5,000 events per file in the viSNE analysis. Each event was still assigned unit weight, and each point set, therefore, still had an equal total mass of 1,000. (ii) The "shortsimplex" method was used for the "wasserstein" function in the "transport" library, which accepted no other parameters besides the pair of weighted point sets (34). (iii) Each population was automatically assigned a zero EMD compared with itself, and EMD scores already computed across the diagonal were simply copied because EMD is a metric. (iv) The "parallel" library was used to parallelize the computation of each row of the matrix, in addition to the above, using the number of cores detected from the "detectCores" function in the "parallel" library.

EMD values computed by “emdist” were compiled in a CSV file and used to create a heat map, in R, for visualization. Statistical comparisons of EMD values between groups were done in Excel using a Student t test. CSV file and heat map are each produced as an output.

The frequency of immune populations was determined in Cytobank and exported into CSV files prior to reorganization. For data sets 1 and 3, populations were identified by traditional biaxial gating. For data set 2, populations were identified by first running a viSNE on nucleic acid expressing events from all patients at all time points and then running a SPADE on the t-SNE axes. Fifteen nodes (15) were identified with 5% downsampling. The following equation was used to determine the change in frequency for all data sets where FREQ_t is equal to the frequency of a population at a given time point and FREQ_pre is the frequency of that same population prior to the start of therapy. The addition of 0.01 to both the numerator and the denominator is to account for the appearance of new populations over the course of therapy.

Change in frequency = ln((FREQ_t + 0.01)/(FREQ_pre + 0.01))

R was used to conduct a paired Student t test to compare samples from the same patient at different time points of treatment. R script provided by Carr and colleagues was used to create boxplots in R (7). In the case of data set 1, a Bonferroni correction was used for multiple hypothesis testing.

MEM creates a quantitative label of cell identity for given populations (19), and the MEM equation is implemented in R. MEM labels were either created for the indicated populations using the bulk, nonpopulation as the reference, except, where indicated, when iPSCs or hematopoietic stem cells were stained and run on mass cytometry as a respective common reference (19). Median MEM labels were created by taking the median MEM score of each marker for each population. Standard deviation is shown. ΔMEM scores are calculated by subtracting the MEM score of the pretherapy sample from the MEM score of the indicated time point.

Root-mean-square deviation (RMSD) and hierarchical clustering were used to compare MEM labels, as previously described (19). The MEM vectors for each nonreference population were calculated over phenotype channels which were shared across all nonreference populations and the single reference population. Each MEM vector contained the population's MEM score, calculated for each of the common phenotype channels, in reference to the single reference population. The MEM RMSD between pairs of nonreference populations was then calculated using the Euclidean distance between these MEM vectors.

Heat maps representing population similarity were generated from each distance matrix using the "heatmap.2” function of the "gplots" library for R. The distance matrix was normalized by the maximum nonnormalized distance d_max between any pair of populations, then multiplied by 100, then subtracted from 100. The result was that zero entries in the original distance matrix would receive a similarity score of 100, whereas the pair of populations with greatest distance in the original distance matrix would receive a similarity score of 0. Thus, two populations with the exact same enrichment score would have 100% similarity (19). To compare populations, median RMSD scores were compared using a two-tailed, Student t test.

Original data sets were provided as FCS files in Flow Repository. The software for calculating EMD and displaying it as a heat map is available as Supplementary Software. The software for generating MEM scores is available in the Supplementary Software by Diggins and colleagues (ref. 19; http://mem.vueinnovations.com/).

Peripheral blood from melanoma patients treated with anti–PD-1 and healthy peripheral blood controls is available as FCS files in Flow Repository. SNaPshot genotyping was done in the clinic on tumors resected from each patient. Forty-eight mutations in NRAS, BRAF, KIT, CTNNB1, and GNAQ were monitored (35).

Peripheral blood from AML patients treated with chemotherapy was generated by CyTOF analysis as described by Ferrell and colleagues (12) and is available as FCS files in Flow Repository (http://flowrepository.org/id/FR-FCM-ZZMC). Patient characteristics and treatment details are available in Ferrell and colleagues.

Serially biopsied melanoma tumors from patients treated with BRAFi and MEKi were generated in separate mass cytometry experiments. Patients MP-034, MP-029, MP-031, MP-032, MP-055, and MP-059 were stained with the mass cytometry panel described in Supplementary Table S2. Patients MP-019, MP-023, MP-054, MP-052, and MP-062 were stained with the panel described by Doxie and colleagues (36). FCS files are available in Flow Repository. SNaPshot genotyping was done as described above.

Data for comparison across cancer types were generated by us in separate mass cytometry studies characterizing untreated melanoma tumors, glioblastoma tumors, and nonmalignant tonsil and adenoid. FCS files are available in Flow Repository. Renal cell carcinoma tumors RC-29, RC-37, and RC-52 were published by Siska and colleagues (32), and Staphylococcal enterotoxin B (SEB)–stimulated PBMCs were published by Nicholas and colleagues (37). Briefly, PBMCs were stimulated with SEB (final concentration 1 μg/mL in in 200 μL) at a cell concentration of 10 × 10⁶ cells/mL. Cells were incubated at 37°C for 16 hours before being washed twice with PBS and stained for mass cytometry.

Data sets used in the RMSD heat map are described by Diggins and colleagues (19).

We sought to investigate how the overall composition of circulating immune populations changes over the course of anti–PD-1 therapy. To do this, we generated viSNE maps, which reduce the dimensionality of the data by plotting cells in two dimensions on the basis of their high-dimensional similarity. Changes in the distribution, or topography, of these maps are indicative of changes in the abundance or profile of individual cell populations. Quantifying changes that exist between viSNE maps relies on either deconstructing the system into populations or qualitative assessment of the overall “shape” of the map. Approaches for quantifying differences in multidimensional cytometry distributions include Quadratic Form (38) and EMD (28). Here, we used EMD to quantify the similarity of two viSNE maps, where differences were quantified as the amount of work required to make the images like another. In this instance, “work” was defined as the number of units, or cells, times the distance moved. To quantify differences in viSNE maps, all samples must be run together on the same axes. It is imperative that measures be taken to reduce batch effects, such as the use of normalization beads here.

Blood was drawn from melanoma patients immediately prior to therapy, and 3 weeks, 12 weeks, and 6 months after the start of anti–PD-1 therapy (data set 1). Mass cytometry was used to characterize PBMCs from each patient at each time point (Supplementary Tables S1 and S2). To monitor and quantify the phenotypic plasticity of the peripheral immune system as a whole, the EMD algorithm was used to quantify differences between viSNE maps as described above (21, 28, 33). One viSNE map was created for 8 patients analyzed at 4 clinical time points and 8 healthy controls (Fig. 1A; Supplementary Table S1). EMD was then used to quantify the differences between each viSNE map, and the numerical results were displayed in a heat map (Fig. 1B). Low EMD scores indicated that the maps were similar, whereas larger EMD scores indicated divergent maps. To determine whether the peripheral immune systems of each patient remained stable or had increased phenotypic plasticity over the course of anti–PD-1 therapy, intrapatient EMD values, those generated by comparing viSNE maps within a single patient over therapy, were compared with interpatient EMD values, those generated by comparing viSNE maps across patients. In 7 of 8 melanoma patients receiving anti–PD-1 therapy, the intrapatient EMD score was lower than the interpatient EMD score, indicating that each patient's peripheral blood immune system was more similar to itself than to that of any of any other patient, regardless of any ongoing therapy response (Fig. 1C). The exceptional patient, MB-009, did not conform to this pattern (median intrapatient EMD value ± standard deviation: 4.22 ± 2.69, versus median interpatient EMD value: 3.99 ± 2.48). This patient was diagnosed with myeloid dysplastic syndrome 8 months after starting anti–PD-1 therapy and was known to have an expansion of mature and blasting myeloid cells and a decrease in all other major cell types in the periphery (15). Thus, combining EMD and viSNE allowed for an automated approach to quantify stability and plasticity of a system over the course of therapy.

Although viSNE and EMD provided a quantitative, first glance at system stability, it is possible that shifts in abundance or phenotype of small, biologically relevant populations may be overlooked when using a single summary statistic. To avoid missing small but crucial cell subsets, cells should be broken into smaller populations manually or computationally and compared. The cells were next manually split into classic immune populations using traditional, biaxial gating (Supplementary Fig. S1). Although biaxial gates were used for cell population identification, it is possible to integrate other automated clustering tools, such as Phenograph (23), VorteX (39), and FlowSOM (40), into the ensemble at this step. The change in frequency of these populations was calculated as a fold change over pretherapy, allowing for visualization of all patients and populations on one graph.

Twenty-eight cell populations were identified, and the change in frequency of those populations compared with the pretherapy time point was calculated (Fig. 2A and B; Supplementary Fig. S2). Two populations showed significant changes in frequency. CD4⁺PD-1⁺ and CD8⁺PD-1⁺ T cells had a significant change in frequency [N = 10; P < 0.0001 (pre- vs. 3 weeks) and P = 0.0011 (pre- vs. 12 weeks) for CD4⁺PD-1⁺ T cells, and P = 0.00057 (pre- vs. 12 weeks) for CD8⁺PD-1⁺ T cells] over the course of therapy, in which both decreased (Fig. 2A–C). PD-1 gates were drawn with nonmalignant tonsil as a reference and considerations discussed and explored by Nicholas and colleagues (37). This decrease in relative frequency is most likely the result of receptor occupancy by pembrolizumab, as opposed to a biological decrease in these cell populations (41, 42).

The next step in the ensemble systems immune analysis pipeline was to automatically quantify enrichment of measured parameters and determine how those enriched parameters changed during treatment. MEM was used to identify signature features of each population at each time point following established methods (19). MEM scores can range from +10 (maximum enrichment) through 0 (no enrichment) to −10 (maximum lack) and here are reported as the median ± standard deviation in MEM value for the cell population. Enrichment quantifications for each population are a vector to be used for additional computational analysis. In this ensemble workflow, vectors describing the high-dimensional phenotype of populations are compared using RMSD (19).

PD-1⁺ CD8⁺ and CD4⁺ T cells from pretherapy samples were enriched for canonical identity makers CD3, and CD8 and CD4, respectively (Fig. 3A and B, top). To indicate changes in enrichment patterns, changes in MEM score (ΔMEM) were calculated by subtracting the median MEM score for each parameter at the pretherapy time point from the indicated time point after the start of therapy. After 3 weeks and 6 months of anti–PD-1 therapy, PD-1⁺CD8⁺ T cells lost enrichment, but not expression, for CD3 (ΔMEM of –5; median ± standard deviation, CD3 MMI 40.9 ± 11 (3 weeks) and 48.3 ± 14.1 (6 months; Fig. 3A, bottom; Fig. 3C). This was not the case for PD-1⁺CD4⁺ T cells (Fig. 3B, bottom). Given a loss of enrichment of the TCR CD3 subunit on these peripheral blood PD-1⁺CD8⁺ T cells, it was informative to determine their novelty by comparing them to other subsets, including PD-1⁺CD8⁺ T cells from tumors or healthy donors. To assess this population, MEM labels were created for PD-1⁺CD4⁺ and PD-1⁺CD8⁺ T cells gated from (i) blood from melanoma patients during anti–PD-1 therapy, (ii) human melanoma tumors, (iii) blood from healthy donors, and (iv) tonsil or adenoid tissue from healthy donors. IPSCs, analyzed by mass cytometry, were used as a control cell reference for MEM calculations. Similarity in MEM labels was then compared using RMSD (Fig. 4A and B). B cells gated from healthy donor blood and tonsils were also included for contrast because of their distinctly different enrichment profiles. As expected, B cells clustered separately from the T-cell populations. PD-1⁺CD4⁺ T cells from the blood of healthy donors and melanoma patients receiving anti–PD-1 therapy clustered together, whereas PD-1⁺CD4⁺ T cells from melanoma tumors and healthy donor tonsils formed a different cluster (Fig. 4A). Similarly, PD-1⁺CD8⁺ T cells from melanoma tumors clustered with those from healthy donor tonsils, whereas PD-1⁺CD8 T cells from healthy donor blood and melanoma blood formed two, intermixed clusters. These results indicated PD-1⁺ T cells found in the blood were distinct from those found in the tumor or healthy tonsil.

Median MEM labels, which display enrichment scores for each measured feature, are shown for each tissue's PD-1⁺ T-cell populations (Fig. 4B). PD-1⁺CD8⁺ T cells from the blood of melanoma patients were enriched for CD43 protein expression and trafficking markers, such as CCR4 and CXCR3, and specifically lacked activation markers CD38 and CD69 compared with PD-1⁺CD8⁺ T cells in melanoma tumors (Fig. 4C, left). PD-1⁺CD4⁺ T cells were enriched for CD4 (+3) compared with their melanoma tumor counterparts. No difference in enrichment between PD-1⁺CD4⁺ T cells in the peripheral blood of melanoma patients and healthy donors (Fig. 4C, right). PD-1⁺ T cells from the blood of melanoma patients and healthy donors were, therefore, phenotypically similar. Thus, MEM, as part of this ensemble, automated quantitative comparisons of high-dimensional data from different tumor and donor types.

The ensemble systems immune monitoring pipeline was applied to a previously published data set of peripheral blood from AML patients undergoing chemotherapy (data set 2; ref. 12) in order to describe and dissect system-wide changes. Five patients with AML were consented for peripheral blood draws over the course of chemotherapy, which were then characterized by mass cytometry. All PBMCs, including blasts and nonblasts, were identified using the gating scheme published by Ferrell and colleagues (12). EMD on t-SNE revealed lack of intrapatient stability in 3 of 5 patients, indicated by no significant difference between intra- and interpatient EMD values (Fig. 5A; Supplementary Fig. S3). The remaining 2 patients showed intrapatient stability, most likely due to the presence of leukemic blasts throughout treatment. Unlike in data set 1, where cellular populations were identified by traditional biaxial gates, 15 populations were defined automatically using the SPADE algorithm to cluster on t-SNE axes (Fig. 5B), as previously described (17, 24). These populations were then fed into the remaining 3 steps of the ensemble analysis pipeline.

Of the 15 populations identified, 7 showed significant changes at some point during chemotherapy (Fig. 5C; Supplementary Fig. S4). MEM was used to label the automatically characterized signature features on each population (Fig. 5D; Supplementary Fig. S4). Two populations of blasts were defined by HLA-DR enrichment and were observed to decrease over the course of chemotherapy. In contrast, CD4⁻ T cells and two subsets of CD4⁺ T cells were observed to expand, relative to other populations, following chemotherapy. Population 13 was enriched for CD4 (+5), CD7 (+4), and CD45 (+2), while specifically lacking expression of HLA-DR (−3) and CD123 (−2). Combined with the expression of CD3 (Supplementary Fig. S5), population 13 is likely a population of T cells. Next, the similarity of population 13 to CD4⁺ T cells from healthy donors was assessed. To do this, a common reference population of hematopoietic stem cells was used to create MEM labels for population 13 and for previously published healthy donor CD4⁺ T cells and B cells (19). RMSD was used to compare MEM labels and determine whether population 13 shared a phenotype similar to that of previously characterized CD4⁺ T cells. Population 13 clustered with CD4⁺ T cells from healthy donors, suggesting that they are phenotypically similar (Fig. 5E). Taken together, these results represent an automated way to dissect immune responses to therapy and suggest that chemotherapy resulted in a relative increase of a nonmalignant population of CD4⁺ T cells.

Samples from a cohort of melanoma patients with BRAF^V600E mutations treated with targeted therapies dabrafenib and trametinib (n = 11; data set 3) were characterized by mass cytometry, and the data were analyzed by the ensemble systems immune monitoring pipeline. Combining EMD and viSNE to quantify stability of the immune compartment showed that each patient remained more similar to itself than to other melanoma tumors or healthy tonsils over the course of therapy. However, EMD run on t-SNE axes created from analysis of T cells revealed that the T-cell compartment did not have lower intrapatient EMD values compared with interpatient EMD values (Fig. 6A). Thus, significant immune plasticity followed therapy. Immune populations were defined using traditional biaxial gates, and the change in frequency for those populations was calculated. After 4 weeks of treatment, a statistically significant change in frequency of 5 immune populations was observed. One population of interest included double-negative (DN) T cells lacking expression of CD4 or CD8 (DN T cells) that comprised 7.23% ± 17.18% in pretherapy tumors, 26.27% ± 16.36% in tumors 4 weeks after therapy, and 3.57% ± 1.52% of a healthy lymph node (median ± standard deviation; Fig. 6B; Supplementary Figs. S6–S7). Pretherapy MEM labels showed that DN T cells were enriched for CD3 (+2), CD45 (+2), CD45RO (+1), CD4 (+1), and CD28 (+1) but specifically lacked CD8 (–3) and CD45RA (–2) when compared with all cells found within the melanoma tumors. ΔMEM scores indicated that, over the course of therapy, DN T cells became more enriched for CD45RO (+1) and CD44 (+1) but lost enrichment of CD69 (–1), CD43 (–1), CD27 (–1), and HLA-DR (–1; Fig. 6C).

Given that T cells can downregulate expression of CD4 or CD8 if activated (43), DN T cells from melanoma tumors were compared with acutely activated peripheral blood T cells stimulated through the TCR by Staphylococcal enterotoxin B (SEB) (37). Additionally, we sought to understand whether similar T cells were found in other cancers, such as glioblastoma (GBM) and renal cell carcinoma (RCC; Fig. 6D). SEB-stimulated T cells formed their own cluster, apart from all other DN T cells, and shared 87.2% (±2.88%, n = 78 comparisons) similarity with DN T cells from melanoma tumors. In contrast, DN T cells from melanoma, GBM, and RCC clustered together, with a similarity score of 93.6% (±1.06%, n = 78 comparisons). Of the T cells activated by SEB, close to half were CD69⁺ (41.26% ± 7.41%), and most retained CD4 or CD8 coreceptor expression (Fig. 6E). In contrast, DN T cells from posttreatment melanoma tumors contained fewer CD69⁺ cells (22.45% ± 8.76%; Fig. 6F). Thus, DN T cells from melanoma, GBM, and RCC consisted of a phenotypically distinct population of DN T cells not observed in resting nondiseased blood and tonsils or acutely stimulated blood. Taken together, these data suggest that an unusual population of T cells emerged in tumors from melanoma patients treated with BRAFi and MEKi and that this DN T-cell population occurs in multiple types of tumors.

The article describes an analysis suite designed to be a common starting point for immunologists tracking cells over time, provides three reference data sets for testing tools designed to discover and characterize cell subsets, and reveals unexpected cells in the context of cancer therapies. This ensemble data analysis strategy was designed specifically for systems immune monitoring in longitudinal, clinical studies, but it could be applied for any system with changes. We envision adapting it to study experimental perturbations to map signaling networks and drug responses (44). Although we expect the algorithms in the ensemble will change and improve over time, the four cellular properties identified should be considered essential features for immune monitoring with any single-cell platform. Cancer immune monitoring strategies must expect the unexpected and be prepared for novel phenotypes (9, 45, 46). The first property, systems plasticity, automatically quantifies the state of a system compared with baseline. In the cohort described here, the peripheral blood of melanoma patients remained stable over the course of therapy. By monitoring systems plasticity with EMD and t-SNE, a patient who experienced bone marrow failure (15) was identified without gating or clustering. In this case, quantifying systems plasticity allowed for the unbiased identification of a patient experiencing large, biologically relevant changes during therapy. Caveats exist for using EMD to quantify differences between viSNE maps. All samples must be embedded on the same viSNE map in order to be appropriately compared, and precautions, such as barcoding and bead normalization, must be taken to ensure that differences in viSNE maps are attributed to biological differences and not batch effects. In addition to EMD, the quadratic form (38) could be used to quantify differences in t-SNE axes. Looking forward, we envision automating a process in which t-SNE and a plasticity test such as EMD are first run on the system as a whole, as described here, and then run iteratively on increasingly refined populations in order to pinpoint the populations undergoing the greatest change or stability.

The remaining components of the ensemble work together to describe changes or stability revealed in the first step, systems plasticity. To accomplish this, the system is divided into subpopulations, and the frequency, signature features, and novelty of each subpopulation are quantified. Identification of cells into groups can be accomplished in whichever way is most appropriate, or in multiple ways, prior to ensemble analysis (22–24, 27). For example, FlowSOM (40) could be used to identify cellular populations, and those populations could be fed into the remaining parts of the ensemble. Because this ensemble toolkit does not rely on known populations, it remains independent from methods of automated population identification (17, 27) and can be used to compare and communicate analysis results from teams relying on computational approaches, immunologists, and bioinformatics experts. The ensemble approach was especially adept at capturing shifts over time and was able to identify a shift in relative abundance of a subset of T lymphocytes, their phenotype, and intrapatient recovery of these phenotypes, all of which could have substantial implications for maintenance of remission and clinical outcomes (47, 48). A human immune monitoring strategy using this ensemble toolkit might provide new insight into the immune system's interaction with leukemia remission and relapse.

This ensemble approach was robust across multiple, contrasting studies in quantifying changes and stability in immune system cells. This approach also provided a detailed analysis of the abundance and tractable quantitative phenotype of populations that comprised each system. We found a population of CD4⁻CD8⁻ DN T cells observed to have a common phenotype in three human tumor types, melanoma, RCC, and glioblastoma that had a phenotype distinct from that of resting CD4⁻CD8⁻ DN and SEB-activated T cells from healthy individuals. This result is consistent with reports characterizing the phenotypic similarity of T cells infiltrating melanoma and colon cancer mouse tumor models (49) and provides evidence that common changes to immune cell mechanisms are shared across human tumor types and play a role in response to targeted therapy.

The ensemble toolkit detected known biological occurrences, as well as identified potential mechanisms for further study. The ensemble toolkit detected an overall loss of PD-1⁺ T cells in the peripheral blood over the course of therapy. This has previously been described and attributed to receptor occupancy by the drug itself (41, 50). Previous work revealed an expanded population of Ki67⁺ T cells in the peripheral blood during immunotherapy (10, 51). Although our comparable study did not measure Ki67 in the blood and, thus, could not be compared directly, we did observe a trend toward increased activated CD4⁺ and CD8⁺ T cells, defined by expression of HLA-DR, at 3 weeks after the start of therapy, and CD4⁺ effector memory T (T_EM) cells at all time points after the start of therapy. In addition to capturing known clinical events, the ensemble tool kit identified a novel population of DN T cells present in unexpectedly high frequencies in melanoma tumors before and after treatment with BRAFi and MEKi compared with healthy lymphoid tissue. Expanded DN T-cell populations have previously been reported in metastatic lymph nodes of melanoma patients (52). However, their deep phenotype or frequency in response to inhibitor therapy has not yet been described. Previous work has described the expansion of a regulatory CD3⁺ T-cell population lacking both CD4 and CD8 after TCR and cytokine stimulation (53–55). MEK inhibitors can support antitumor T-cell function by blunting TCR-induced apoptosis (56). Therefore, it is possible that high frequencies of DN T cells after treatment with MEKi and BRAFi indicate an accumulation of T cells derived from tumor-resident T cells, although additional mechanistic studies are required.

Overall, the melanoma tumor data set (data set 3) highlights the complexity of tumor-associated T cells in human malignancies and provides further evidence that phenotypically diverse populations of tumor-resident T cells can be found across multiple, distinct tumor types (49). The data presented here will join a common immunology reference set (19, 57, 58) that can be mined further to characterize and understand changes in immune and cancer cell populations in diverse disease settings and build a reference of cellular identity.

C.E. Roe reports receiving commercial research funding from Incyte Corporation, Pharmacyclics LLC, and Janssen Pharmaceutica. M.C. Kelley reports providing Medical Professional Liability Reviews unrelated to the content of this research manuscript. D.B. Johnson reports receiving commercial research funding from Bristol-Myers Squibb and Incyte and is a consultant/advisory board member for Array, Bristol-Myers Squibb, Genoptix, Merck, and Novartis. J.M. Irish is co-founder of Cytobank Inc.; reports receiving commercial research funding from Janssen, Incyte, and Pharmacyclics; and has ownership interest in Cytobank Inc. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A.R. Greenplate, D.B. Johnson, J.M. Irish

Development of methodology: A.R. Greenplate, D.D. McClanahan, B.K. Oberholtzer, D.B. Doxie, J.M. Irish

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A.R. Greenplate, C.E. Roe, N. Leelatian, M.C. Kelley, P.J. Siska, J.C. Rathmell, P.B. Ferrell, D.B. Johnson

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): A.R. Greenplate, D.D. McClanahan, K.E. Diggins, N. Leelatian, P.B. Ferrell, J.M. Irish

Writing, review, and/or revision of the manuscript: A.R. Greenplate, D.D. McClanahan, B.K. Oberholtzer, D.B. Doxie, C.E. Roe, N. Leelatian, M.C. Kelley, V. Gama, P.B. Ferrell, D.B. Johnson, J.M. Irish

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A.R. Greenplate, D.B. Doxie, C.E. Roe, M.L. Rasmussen, V. Gama

Study supervision: J.C. Rathmell, J.M. Irish

Other (generation and characterization of induced pluripotent stem cells from human skin fibroblasts to be used as controls for analysis): M.L. Rasmussen

This study was supported by NIH/NCI R00 CA143231 (J.M. Irish), F31 CA199993 (A.R. Greenplate), U01 CA196405 (D.D. McClanahan and J.M. Irish), R25 GM062459 (D.B. Doxie), T32 CA009592 (D.B. Doxie), K23 HL138291 (P.B. Ferrell), and K23 CA204726 (D.B. Johnson); the Vanderbilt-Ingram Cancer Center (VICC, P30 CA68485); a Vanderbilt University Discovery Grant (J.M. Irish and N. Leelatian); a VICC Provocative Question award (M.C. Kelley and J.M. Irish); and VICC Ambassadors (J.M. Irish and A.R. Greenplate).

The authors thank Holly Crandall and Kate Carson for their help with consenting patients and collecting patient samples.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.