The re-analysis was performed on a dataset from the Cytobank project ‘Spitzer-et-al_Cell_2017’ (https://community.cytobank.org/cytobank/experiments#project-id=1025) obtained from the Cytobank repository (15). It was originally collected by Spitzer et al. (14), who used a common mouse carcinoma model to study four different treatments, of which two were effective and two were ineffective ( Figure 1 ). Briefly, the mice were implanted with tumor cells, which were left to grow until they reached a specific size. Then, the mice were shared between four treatment groups. Disease progression and treatment effect were monitored by blood and tissue sampling at different time points. Of these data, we re-analyzed blood samples taken from 12 mice on day three of treatment with 41 biological markers measured for each sample. Thus, a total of six animals were in the effective treatment and originally six animals were in the ineffective treatment. Each treatment group was assigned three animals: untreated (untr.), ineffective treatment (ineff.; anti-PD-1), effective treatment1 (effective 1; IFN-γ + anti-CD40 + CD1-allo-IgG) and effective treatment2 (effective 2; IFN-γ + anti-CD40 + B6-allo-IgG). However, one blood sample from the ineffective treatment (untr. blood3) was excluded from our analysis because of quality reasons (please see below).

In order to use PRI to visualize and analyze the dataset, some necessary steps to prepare the data must be performed. In Figure 2A we briefly summarize all steps that were part of the workflow to achieve the results as we present them here. The following paragraphs describe these steps in more detail.

Precleaned debarcoded data from blood samples were obtained from the public dataset (14). The blood samples were further processed in FlowJo™ v9.9.6 software as is exemplarily shown in the gating strategy in Figure 2B .

Here, we gated first on live single cells (191Ir193Ir⁺ and 195Pt- cells) and removed internal standard beads added before CyTOF measurement by using unique 140Ce⁺⁺ signal. Total leukocytes were gated based on CD45⁺ expression. T cells were further selected by CD3⁺CD19- expression. CD4 T helper cells used for the re-analysis in this study were determined by their CD4⁺CD8- expression.

The quality of the datasets was assessed with conventional metric multidimensional scaling (MDS, see Supplemental Figure 3A ) performed using the Bioconductor package ‘limma’ (16). Pairwise Euclidean distances of the median expression values for the 41 biological markers in all samples were calculated. It was evident that one blood sample from the untreated group differed in its SI medians from the other samples in this group. A similar behavior of this sample was seen using Ward’s hierarchical clustering (17) from the ‘stats’ package that is part of R ( Supplemental Figure 3B ). Both results and a very low cell count prompted us to exclude sample “untr. blood3” from further examination with PRI.

The extracted CD4⁺ T cell populations are imported as FCS files in PRI, where the SI are transformed with inverse hyperbolic sine (arcsinh), commonly used for flow cytometric analysis (18), with a cofactor of 0.1.

The signal intensities of two markers were divided into bins of width 0.2 x 0.2 asinh along the x and y axis for bin plot visualization. For all cells in each bin, different statistical values were calculated and displayed in a color-coded manner as heat maps. Low values are represented as shades of blue, median values as yellow, and high values as red color. The patterns of mean signal intensity of parameter z⁺ cells (MSI⁺) and the mean signal intensity of parameter z for all cells per bin (MSI) are analyzed. Additional statistics are displayed in black and red percentage numbers in each quadrant. Percentages in black indicate the frequency of cells relative to total cells. The frequency of z+ cells is given in red, and is relative to the cells inside the respective quadrants. All events (cells) are included in the analysis but only bins containing five or more cells are displayed, balancing the impact of identifying comparatively rare subpopulations while retaining statistical robustness.

In this publication we present the core functionality of PRI (for details visit https://github.com/InesHo/PRI-demonstration; also available as R-package within the repository). We used the R package ‘flowCore’ (19) to read the FCS files and R-core (20) functionalities to create the plots. The plots for MDS and Ward’s hierarchical clustering presented in the Supplementary Material have been created using the R-packages ‘ggplot2’ (21) and ‘dendextend’ (22).

Two-group comparisons (n = 5 vs. n = 6) were statistically assessed by non-parametric Mann-Whitney U test using GraphPad Prism8 Software. Median values were used. P < 0.05 was considered significant, with the numbers of asterisks indicating: *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001.

To evaluate the added value of the PRI approach (13) used in combination with other strategies for the analysis of high-dimensional mass cytometric data, we selected a publicly available dataset (14). Spitzer et al. studied the effectiveness of different immunotherapies in a mouse model of breast carcinoma. From the comprehensive data, we selected the three-day blood dataset, including six mice, from which the authors had described a mixture of CD4⁺ T cell clusters expanded by two different anti-tumor therapies, as revealed by a combination of Scaffold Maps (6) and Citrus analysis (7) ( Figure 1 ).

We attempted to use our bin-based approach PRI to further analyze, constrain and visualize the specificities of the described cluster mixture (CD44⁺ CD69^- CD62L^- CD27low CD90⁺ Tbet⁺). PRI uses an x-y plane divided into equally sized bins and color-codes the intensity or frequency of a z parameter. In this process, PRI semi-continuously analyzes and visualizes the combinatorics of mainly 3 parameters in so-called bin plots. Each axis is informative and intuitively understandable. PRI operates reproducibly and provides combinatorial features such as bin intensity ranges or maximum bin intensity values that are threshold-independent. In cases where frequencies are desired either in bins, quadrants, or cells, it is necessary to manually define thresholds for parameters (e.g., between negative and positive or between low and high values). All these features are useful to compare and classify samples.

Figure 2 summarizes and describes the work-flow ( Figure 2A ), including the gating strategy ( Figure 2B ). To ensure the quality of the analysis, we excluded one sample from the untreated group (see details in the Methods part).

Using the clustering information of Spitzer et al. together with background immunological knowledge and manual inspection of many 3-parametric overview plots, we decided to use CD90 and CD44 to define the x-y parameter plane. Evidence for CD90 and CD44 as possible x and y markers also comes from the non-redundancy score (NRS) ( Supplementary Table 1 ) (23). The selection criteria included a broad intensity range and sufficient grouping of Ki67^high cells, as well as Tbet cells. The first bin analysis statistic of the maximum bin MSI values of selected T cell parameters confirmed this choice and revealed statistically significant differences between the two classification groups (5 mice with untreated/ineffective treatment and 6 mice with effective treatments 1/2) for the five z markers CD86, CD27, Foxp3, KLRG1, and PDL1 ( Figures 3B , 4B ).

Furthermore, the manual inspection of the MSI bin plots with these z parameters revealed different patterns in the case of CD86 and CD27 ( Figure 3A ). For the visualization of these patterns, we used dynamic MSI bin plots, meaning that the z parameter of each is color-coded relatively to its own individual minimum-maximum range of MSI ( Supplemental Figures 1 , 2 ). Thresholds were defined for x and y to discriminate CD44^- and CD44⁺ cells, as well as CD90^-/low and CD90high cells.

It is obvious that there are higher CD86 MSI-bins in effective treatment (bin-range right in each plot). In addition, the biggest differences between untreated/ineffective treatment and effective treatments are mainly found in quadrant 3 (Q3). The highest CD27 bins are found in Q2 only in the samples from the untreated/ineffective treatment groups with one exception and the differences in the patterns are not as evident as in CD86.

However, lower cell frequencies in Q2 are also statistically significant in the effective treatment group ( Figure 3B left). Setting manual thresholds for the z parameters based on histograms similar to those in FlowJo allowed the determination of the frequency of CD86⁺ and CD27⁺ cells in the bin plots (red frequencies in the corners of each quadrant). High frequencies of CD86⁺ in Q3 and low frequencies of CD27⁺ in Q2 were found to be significant for the classification of the samples.

To simulate multiparametric visualization and characterize the therapy-expanded cell subsets grouped in Q3 and Q4, we plotted seven additional z markers side by side on the same x-y plane as in Figure 3 ( Figure 4A ). However, in contrast to Figure 3A , we took both exemplary samples (untr. blood1 and effective2 blood2 (IFN-γ + anti-CD40 + B6-allo IgG)) for each z marker and determined a common lower and upper bin-MSI value from both samples together. This provides a direct visualization of differences in intensity of both samples. The PRI feature statistics of all 11 samples ( Figure 4B ) reveals some classification value for Ki67 (frequencies in Q3 and Q4), KLRG1 (maximum bin MSI and frequencies in Q4), Foxp3 (maximum bin MSI and frequencies in Q3), and PDL1 (maximum bin MSI and frequencies in Q3). Interestingly, Tbet shows a tendency of different median values in the frequencies of Q3, although it was a major characteristic of the expanded cell cluster mix in the original work of Spitzer et al.

Next, we analyzed differences in expression level between bins taking only z⁺ cells into account ( Figure 5A ). Bins with less than five z⁺ cells are shown in gray. A distinctive pattern is apparent for the two transcription factors (TF) Tbet and Foxp3. The highest and lowest intensities are clearly divided between quadrants Q3 and Q4 and in opposite ways. Ki67 expression, on the contrary, shows no clear tendency to higher level in either of the two upper quadrants.

To investigate the co-expression of Ki67 with the two major TF, we compared the bins for two z parameters per cell in the two exemplary samples. The pie charts clearly show that effective therapy mainly increases the Ki67/TF co-expression sectors of Ki67⁺Tbet⁺ in Q3 and Ki67⁺Foxp3⁺ in Q4 ( Figure 5B ).