FCS files and de-identified metadata from intracellular cytokine staining (ICS) analysis of CD4+ and CD8+ T cell responses was provided from the HVTN from HVTN 105 a Phase 1 preventative vaccine trial (ClinicalTrials.gov NCT02207920). Primary ICS analysis was previously reported (2). Details regarding the study design, participants, sample and data acquisition are included in the primary study manuscript (2). Briefly, participants were randomly assigned to 1 of 4 groups with an allocation ratio of 1:1:1:1 ( Figure 1 ). Participants received different combinations of AIDSVAX B/E, DNA-HIV-PT123, and placebo, administered intramuscularly. AIDSVAX B/E consisted of 300 μg of subtype B (MN) HIV gp120 glycoprotein and 300 μg of subtype A/E (A244) HIV gp120 glycoprotein adsorbed onto aluminum hydroxide gel adjuvant and administered into the right deltoid muscle. DNA-HIV-PT123 contained a mixture of 3 DNA plasmids: (a) clade C ZM96 gag, (b) clade C ZM96 gp140, and (c) clade C CN54 pol-nef, delivered at a total dose of 4 mg administered into the deltoid muscle via needle and syringe. Serum for humoral assays was obtained from serum-separating tubes (SSTs) and frozen at –80°C. Peripheral blood mononuclear cells (PBMCs) for cellular assays were isolated and cryopreserved from within 6 hours of venipuncture, as described previously (8). Flow cytometry was used to examine HIV-1–specific T cell responses using a validated intracellular cytokine staining (ICS) assay. The peptide pools evaluated were vaccine matched (ZM96 gp140-Env1, ZM96 gp140-Env2, 92TH023-Env, and ZM96 Gag), covering Env and Gag. Previously cryopreserved PBMCs were stimulated with the synthetic peptide pools. As a negative control, cells were not stimulated. Serum HIV-1–specific IgG, IgG3, IgG4, and IgA responses were measured with a custom HIV-1–binding antibody multiplex assay (BAMA) as previously described (9, 10) using gp120 proteins and V1V2 antigens detailed previously (11).

The set of fluorescent dimensions ℱ in z C were transformed using the “log-like” inverse hyperbolic sine, sinh^-1, in conjunction with a set of F-dimensional cofactors [α ₁, α ₂,…,α_F ] for each dimension j ∈ ℱ . Each vector z j C was divided by its corresponding cofactor α_j prior to transformation, which effectively removed the artifactual bimodality introduced by the raw sinh^-1 transformation.

To determine a suitable set of cofactors, each vector z j C was first transformed by sinh^-1 (Equation 1) and its intensity histogram was examined.

Each α_j was defined as the hyperbolic sine, sinh, of half the magnitude of the distance between the positive P ⁺ and negative P ^- peaks ( Equation 2 ) nearest zero in the intensity histogram of each sinh − 1 z j C ,

and because P ⁺ and P ^- were defined in the transformed space, sinh was required to convert values back to the raw data space. The cofactors were then applied as follows (Equation 3),

Note that scatter dimensions are typically not sinh^-1 transformed, but for convenience we refer to the full data (scatter included) after sinh^-1 transformation simply as z T .

To identify saturated events, all raw data vectors Zj    were transformed by ( Equation 3 ) with α_j = 100 to yield z j T . Then each z j T was allocated to 1024 uniformly-spaced bins, denoted s_bin. Each minimum bin was defined by (Equation 4),

and each maximum bin was defined by (Equation 5),

where R_j was the channel-specific keyword-value range parameter ($PnR) from the TEXT section of the FCS file. To determine the saturated event threshold h_j , we first examined a window w_j of the top-most 61 bins (Equation 6),

Then the median and robust standard deviation of the differences between consecutive bins were used to identify bins that contained extreme differences (Equation 7),

where w^D was the difference between consecutive bins w_b – w_{b -1} for b ∈ { 2 , 3 , … , 61 } , w^M was the median difference, w^σ was the robust standard deviation of differences, and w^X was a vector of 1’s and 0’s that indicated the presence or absence (respectively) of extreme differences. If no extreme differences were found, the examination window was shifted by -1 bin, w = s_bin _{[963,964,…,1023],} and re-examined. This process was performed iteratively until at least 1 extreme difference was found. Then the lowest s_bin_Xj that contained an extreme difference was identified by (Equation 8),

and its corresponding histogram value v j T was inverse-transformed back to a raw intensity by (Equation 9),

The saturated event threshold h_j was set to the raw intensity v_j (or 80% of the maximum data range, whichever was higher) as follows (Equation 10),

Finally, all events with raw intensities above the saturated event threshold were removed.

To identify time defect events, corrected fluorescence data were sorted by time. Then each z j C was allocated to B non-uniformly-spaced bins, denoted t_bin, and each contained the same bin_size number of events as follows (Equation 11),

Then the median event value m was determined for each bin. The vector of bin median event values within each dimension j were Z-score standardized by (Equation 12),

Then any bins containing time defects were defined by (Equation 13),

and all events within each t _ b i n D j were removed.

The censoring process (described above) identified and removed:

1. raw fluorescent events that saturated above the limits of detection (saturated events).

2. corrected fluorescent events that contributed to inconsistent signals over time (time defects).

Following the removal of saturated events and time defects, new FCS files were generated from the remaining data. Supplementary Figure 1 shows the number of cells per sample before and after censoring for all samples.

The quality of compensation matrices was assessed in FlowJo, and any sub-optimal compensation values were manually corrected. The optimized compensation matrices were inserted into the FCS files.

All marker-fluor combinations were consistent across the entire dataset. However, some FCS files contained channel names that did not match other files. Any mismatched channel names were corrected, and new FCS files were generated.

To remove variation due to experimental batches, while maintaining as much biological variation as possible, swiftReg (7) was used to register batches. This approach first registered each batch separately to the same reference batch, then applied the resulting batch-specific shifts to all individual samples in that batch.

To do this, first a SWIFT cluster template was produced from a concatenate of antigen-stimulated samples from a single reference batch, and similar concatenates from each of the other batches were registered by NDCR to the reference template. The resulting batch registration template contained batch-specific maps of cluster movement vectors that specify the value-adjustments necessary to bring that batch’s clusters into alignment with reference clusters. All individual samples in each batch were then registered using these batch-specific cluster movement vectors. This process generated new batch-registered FCS files.

To enhance detection of rare, biologically-significant populations and reduce computational burden, all batch-registered samples were randomly sub-sampled and combined into a single concatenated FCS file that was then clustered by SWIFT. The resulting SWIFT cluster template was used to identify debris clusters in FSC-A and SSC-A, as well as non-CD4 T cells. New FCS files were generated from non-debris CD4+ events.

Detection of positive markers was selectively enhanced by smoothly increasing intensity values about a user-specified inflection point. The smooth increase was achieved by multiplying intensity values within a channel by a sigmoid function ( Supplementary Figure 2 ) as follows (Equation 14),

where r was a 1×601 vector of values between -3.00 and 3.00 with intervals of size 0.01, P was the degree of overlap in the expanded region (default P = 0.5), L was the user-specified inflection point where expansion occurred in that channel, W was the width of the expanded region in decades (default W = 1 ), s_j was a vector of scaling values that were multiplied with z j C element-wise to produce expanded data z j E , normcdf is a MATLAB function that returned a cumulative standard normal distribution, and interp1 is a MATLAB function that returned interpolated values of the function y = f(x) at specific query points z j C by spline interpolation (option = 'spline').

Because the Env-1-ZM96 and Env-2-ZM96 peptide pools constituted the non-overlapping peptides covering the Env-ZM96 Env sequence, the total Env-ZM96 response for each sample was calculated by combining Env-1-ZM96 + Env-2-ZM96 and subtracting negctrl1. Note that Negctrl1 was subtracted here once to account for the additional background contribution of combining raw Env-1-ZM96 + Env-2-ZM96 counts. The 92TH023-ENV samples were stimulated with peptides covering the whole 92TH023-ENV sequence (2).

The cell counts for AnyEnvNeg1 were then defined as the maximum of the cell counts for 92TH023-ENV or Env-ZM96 + Env-2-ZM96, minus the background from negctrl1. Because the 92TH023-ENV and ZM96 ENV sequences have some homology, it is very likely that some peptides, presented by the MHC alleles of some participants, will be cross-reactive between the two ENV peptide sets. However, the extent of cross-reaction cannot be estimated from this dataset, and so we used a conservative definition of the “Any-Env” response as the larger of the response against either ENV sequence. This uses the conservative assumption that all T cells cross-reacted, and therefore the total response is revealed by the higher of the two anti-Env responses.

The number of CD4+ T cells producing IL-2+ and/or IFNγ+ was expressed as a percentage of the total live cells in the corresponding sample. Percentages below 0.005 were thresholded to a minimum of 0.005.

To identify responders, the variance of cell counts was first stabilized across clusters. Cluster-specific scaling factors were defined as half the median of cell counts across all negctrl samples for each cluster, with low scaling factors thresholded to a minimum of 10. All counts were then transformed by inverse hyperbolic sine (asinh) after division by the cluster-specific scaling factor. Transformed stimulated counts (TSC) were obtained by subtracting each cluster’s transformed background count from its pairwise transformed stimulated count (for values above the threshold, this is analogous to a log ratio).

For each sample group defined by Treatment, Stimulation, Visit, and Cluster, the standardized pairwise background variances (SPBV) between Negctrl1 and Negctrl2 were defined as the square root of the sample-mean of the squares of their pairwise differences.

Then for each sample, the p-values (with a Null hypothesis of “no difference”) were determined by applying the normal distribution survival function (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html) to the ratio of TSC over SPBV. This was performed separately for the stimulated sample with each background (Neg1 or Neg2). The final reported probability of a sample being a responder was then 1 minus the mean of its two p-values. All samples with a final probability of ≥98% were considered to be responders.

Antibody levels measured by binding antibody multiplex assay (BAMA) for HVTN 105 were obtained from HVTN. To compare the CD4 responses to related antibody levels, the Spearman correlation coefficients were calculated for related antigens using log transformed antibody abundances. Specifically, CD4 responses upon ZM96 gp140-Env1 and ZM96 gp140-Env2 stimulations were compared with antibody responses to 96ZM651.D11gp120.avi, gp41, gp70–96ZM651.02 V1v2 antigens. CD4 responses to 92TH023-Env were correlated with antibody response to A244 gp120 gDneg/293F/mon, AE.A244 V1V2 Tags/293F and gp41 antigens. Finally, CD4 responses to ZM96 gag were compared with antibody responses to p24.

The HVTN 105 dataset comprised 3,200 .FCS files representing 24 batches, with accompanying compensation matrices for each batch. In general, the Visit 5 (V5), V7 and V9 PBMC samples for one participant were all analyzed in the same batch, whereas the V11 PBMC samples were analyzed in separate sets of batches. As described previously (2), if PBMC samples did not meet quality control criteria, those samples were re-analyzed in a subsequent batch, resulting in duplicate analyses. After curation according to these rules, the complete dataset potentially comprised four vaccine groups each containing 26 participants, eight in vitro antigen stimulations, and four time points, for a total of 3,328 flow cytometry samples. The study design did not include a placebo group receiving no HIV antigens, and the T cell data did not include a baseline sample, i.e. before vaccination. Therefore the important negative controls are the pairs of “negctrl” samples that did not receive in vitro stimulation with any antigens. Our re-analysis focused on six of the eight antigen stimulations: two negative control samples (negctrl2 and negctrl2); E92TH023_ENV; Env_1_ZM96; Env_2_ZM96; and Gag_ZM96. This resulted in a total of 2,496 potential samples. Due to some dropouts and missing negctrl replicates, the final number of flow cytometry samples in our re-analysis was 2,393. The samples, batches, repeated samples and final analyzed samples are shown in detail in Supplementary Figure S3 .

A consensus .FCS file was produced by concatenating sub-samples of all HIV Ag-stimulated samples from all batches at the V9 time point. V9 was chosen because Visit 9 was the pre-determined immunogenicity time point for the HVTN105 trial, and for most groups and antigen stimulations in vitro, this was also the strongest response (see below). V9 samples were therefore enriched for the rare, activated T cells, facilitating capture of these cell populations in the cluster template. Because this concatenate included samples representing all HIV antigen stimulations, this is an objective way to include potential cell phenotypes induced by any of the HIV antigens in any treatment group.

This concatenate was clustered using SWIFT to establish a high-resolution cluster template of all cell sub-populations. All samples were assigned to the resulting cluster template, establishing the number of cells in each cluster, in each sample. All cluster membership information was then condensed to two dimensions using UMAP (Uniform Manifold Approximation and Projection) (12). The results in Figure 2A show batches encoded by colors, stimulations by symbols, and visit number by symbol size. The strongest contribution to diversity was clearly the batch - most members of each batch are clustered together, and the batches are substantially resolved. This is particularly true for the V11 batches on the right, that were analyzed in a different set of batches from V5, V7 and V9. The presence of batch effects is not surprising in samples analyzed over a period of months - we have seen batch effects in all such datasets that we have examined. The HVTN105 batch effects were relatively minor, and so registration could be used to reduce batch effects and improve the comparison of the vaccine groups.

We have previously developed swiftReg (7), an automated registration tool that builds on the SWIFT clustering algorithm to perform high-resolution alignment of samples at the single-cluster level. The HVTN 105 batches were registered by producing a SWIFT cluster template from Batch 2204, producing consensus samples from each batch, and then registering each batch consensus sample to the Batch 2204 consensus cluster template. This generated, for each batch, a map of registration shifts that were then applied to each individual sample in the respective batch. This procedure registers the overall batch trends, without altering the differences between individual samples within each batch that might carry biological information.

A new SWIFT cluster template was generated from a consensus of all registered V9 HIV antigen-stimulated samples. After assignment of all registered samples to the resulting cluster template, the cell numbers per cluster were reduced to two dimensions by UMAP, and Figure 2B shows that the registered batches were intermingled. ‘Micro-aggregates’ of samples from the same batch were still visible - focusing on just 15 participants for clarity, each micro-aggregate comprised samples from a single participant (including different stimulations and time points). These tended to group in close proximity on the UMAP projection ( Figure 2C ), even though the Visit 11 samples were analyzed in different batches from the Visit 5, 7, 9 samples. Overall, the samples included time points spanning 18 months. Thus, most individuals are sufficiently diverse for SWIFT analysis of flow cytometry data to identify a unique ‘fingerprint’ of cell populations in different participants. We have observed this pattern in other studies (unpublished). The proximity of the registered V11 data points to the V5, V7 and V9 points from the same participant reinforces the interpretation that the HVTN 105 batch effects have been substantially reduced by the registration process. Examination of individual parameters by the same approach identified parameters, e.g., CD4 and CD8, that contributed to these batch differences ( Figures 2D–G ). Interestingly, the groupings of similar batches were variable between different parameters ( Figure 2 ; Supplementary Figure 4 ).

We then used the model-based SWIFT multidimensional clustering algorithm (5, 6) to generate an unbiased cluster map from a sample constructed by concatenating a random subset of events from samples across all batches. The SWIFT algorithm is particularly useful for detecting rare populations (13), possibly because these were the type of samples used during SWIFT development (5, 6). Preliminary analysis of the cluster map indicated that the antigen-specific responses in many samples were small, consistent with the previous analysis (2). To maximize the sensitivity of detecting all cytokine producing sub-populations, we produced a new SWIFT cluster template from a concatenate of random subset of events from all antigen-stimulated samples, using only the scatter, live/dead, and CD4 parameters. All samples were assigned to the resulting template. As described previously (2), the non-replicating vaccines induced almost no CD8 T cell responses, and so further analysis focused on CD4 T cells. Clusters containing CD4 T cells were selected, and all the events in this set of clusters were saved, for each flow cytometry sample, as reduced-size .FCS files for further analysis. This “cluster gating” (6) allowed subsequent analysis to focus more clearly on the cells of interest, because clustering could then be performed on a full concatenate of the entire dataset. The resulting .FCS files are more amenable to analysis by SWIFT, other automated algorithms, and manual analysis.

A large concatenate was then produced from all cluster-gated events in all samples stimulated with HIV antigens (Env, Gag), from all groups at Visit 9, which was the time point that showed the highest responses overall. A second concatenate was produced from the corresponding negative control samples. SWIFT cluster templates were created from each of these two large concatenates, using all parameters for high-resolution clustering. The two resulting cluster templates were combined, and all individual samples from all groups, all visits, all stimulations were assigned to the resulting combined template (total clusters 2,246). This cluster competition approach (7) sharpens the differences between the two groups represented by the two templates, in this case stimulated and unstimulated cell populations. Note that each concatenate included samples from all vaccine groups, so the competition process should not affect the resolution or statistical analysis of any study group differences.

Cluster gating (6) was then used to narrow down the cell populations of interest. During cluster gating, all cells are assigned their cluster medians in all dimensions, so that the two-dimensional gating shown in Figure 3A takes advantage of all the information in all dimensions. Activated CD4 T cell clusters were identified as live, singlet, CD3+ CD4+ CD154+ TNF+ clusters ( Figure 3A ). Additional marker intensities for all parameters are shown in Supplementary Figure S5 . These activated CD4 T cell clusters were then examined by testing the significance of differences between antigen-stimulated and negative control clusters in all participants at visit 9. A Wilcoxon test was followed by the Benjamini-Hochberg correction for multiple measurements, because of the number of clusters examined. Figure 3B shows the ratios and magnitudes of differences between antigen-stimulated and negative control cultures in a volcano plot. All clusters that were significantly increased in the antigen-stimulated samples (green shaded area) were chosen for further analysis. To facilitate comparisons with previous analysis (2), the SWIFT clusters were aggregated into four groups: IL-2+IFN-γ+, +/-, -/+ and -/- cells (15, 5, 3 and 4 clusters, respectively). The heatmaps ( Figure 3C ) show the marker characteristics of each cluster.

The samples showing significant responses to each antigen, at each time point, were then evaluated as described in Methods, using the aggregated cluster data for all clusters producing IL-2 and/or IFN-γ. Figure 4 shows the results for each time point, each vaccine treatment, and five antigen stimulations, or combinations of stimulations: AnyEnv (Env92 or Env1/2), Env92 (Env92TH023 only), Env1/2 (Env1 plus Env2), GAG-ZM96, and the negctrl2. Background values (negctrl1) were subtracted from all antigen-stimulated values (similar conclusions were obtained if the negative controls were reversed). All samples with >98% probability of being genuine responders are shown in red. As expected, very few negctrl samples were evaluated as responders (at a confidence level of 98%, a small number of false positives are expected). As Gag antigen was only included in the DNA vaccine, Treatment 1 uniquely lacks immunization with Gag for the first blood sample evaluated, at Visit 5. Consistent with this, only Treatment group 1 lacks a response to Gag at Visit 5. At a very high confidence level of 99.9%, there were still high rates of responders (up to 88%) but no responders in any negative controls ( Supplementary Figure S6 ). Supplementary Figure S7 shows an alternative layout of the responder data to emphasize the time course within each group.

Several combinations of vaccine treatments and times induced responses in the great majority of participants, particularly in Treatment group 3 at Visits 7 and 9. The numbers of responders were generally higher than evaluated previously (2), possibly because the extensive pre-processing and the competitive cluster templates used in our analysis provided sharper distinction between antigen-stimulated versus background cells producing cytokines. The magnitude of the net anti-HIV T cell responses was well-correlated between the original analysis and the re-analysis ( Supplementary Figure S8 ). There is a general trend towards higher magnitudes detected by SWIFT (compared to the 1:1 reference line), possibly due to the effectiveness of high-dimensional definition of populations, as well as the sharper signal:noise discrimination by focusing on the clusters that were significantly increased by antigen stimulation.

The quality of the cytokine response to protein or DNA-derived immunogens was assessed between the different vaccine treatments by comparing the ratio of T cells producing IFN-γ vs. T cells producing IL-2 but not IFN-γ. The anti-Gag response is easiest to interpret, as this is induced only by the DNA vaccine. Figure 5A shows that this response is biased towards IFN-γ production, consistent with a previous report (14). The response to the ZM96 clade C peptides, primed by DNA, also showed a tendency towards an IFN-γ-biased response. In contrast, the response to clade E 92TH023 protein immunization was biased more towards IL-2-only responses, consistent with our previous demonstration (15) that viral infections tend to induce more Th1/IFNγ responses, whereas protein vaccines tend to produce responses biased towards IL-2-producing central memory (16) cells. Figure 5B summarizes these results, including the results for the minority IL-2- IFN-γ- and IL-2- IFN-γ+ responses.

The flow cytometry panel included several cytokines, including IL-21 (produced by Tfh and some other cells) and IL-4 (produced by Th2 cells). Manual examination of the concatenated results suggested that antigen stimulation appeared to induce a small IL-21 response in a relatively low number of TNFa+ IL-2+ CD4 T cells. However, the IL-21 staining was weak, and did not result in a clearly separated sub-population of positive cells. As the SWIFT clustering algorithm uses a criterion of multidimensional unimodality to define individual sub-populations (6), the putative IL-21+ cells were initially difficult to identify by clustering. We therefore used a ‘stretching’ modification that slightly broadened the cell distribution across the expected junction between IL-21- and IL-21+ cells. Clustering the resulting data in SWIFT allowed the reproducible detection of IL-21+ clusters ( Figure 6 ). In contrast, applying the same stretching modification to the IL-4 channel did not result in the detection of IL-4+ clusters, consistent with the manual examination of the IL-4 data ( Figure 6 ). The IL-21+ clusters were activated memory CD4 T cells (CD154+ CD45RAlo CD4+), but interestingly, did not express the CXCR5 chemokine receptor that is characteristic of circulating T follicular helper (Tfh) cells ( Figure 6 ), perhaps due to down-regulation of CXCR5 on the in vitro activated cells as we have observed previously (S. De Rosa, unpublished). In contrast to the IL-2+ IFN-γ+ versus IL-2+ IFN-γ- skewing described above, the IL-21+ cells were observed in all treatment groups, and did not show obvious biases towards particular antigens or immunization strategies ( Figure 7 ).

Binding antibody (Ab) responses were the major correlates of risk (CoR) identified in the RV144 Trial (3). Subsequently, at V9, 2 weeks after the final immunization we assessed the relationship of IL-2/IFN-γ and IL-21 defined clusters with the contemporary HIV-specific plasma Abs ( Figure 8 ). Overall T1 ( Figures 8E, J ) had low antibody responses at this timepoint compared to the other groups as expected due to the boosting immunizations with DNA alone. T2 overall exhibited the greatest number of significant correlations with the IgG response, primarily associated with responses to AE.A244 ( Figure 8B ) which most closely matches the protein component of the vaccine regimen. T3 and T4 overall had significant associations relatively balanced between AE.A244 ( Figures 8C, D ) and 96ZM651 ( Figure 8H ) Ab responses which most closely matches the DNA component of the vaccine regimen, and is consistent with T3 and T4 receiving 4 doses of DNA. T3 had the greatest number of significant correlations between IL-21+ and Ab responses ( Figures 8C, H ), consistent with T3 having the overall greatest IL-21+ response. Supplementary Figure S9 shows the magnitude of IgG responses for the T cell responders identified by the original analysis or the new SWIFT analysis.

Both total IgG specific for the V1V2 region of gp120 and IgG3 specific for V1V2 were inverse CoR in RV144 (3). IL-2-IFN-γ+ cells were significantly correlated with IgG AE.A244 V1V2 in T4 ( Figure 8D ), with IL-2+IFN-γ+ and IL-2-IFN-γ also significantly correlating with IgG AE.A244 V1V2 in T2 ( Figure 8B ). IL-2+IFN-γ+ also, and IL-21+ also significantly correlated with IgG gp70–96ZM51 V1V2 in T3 ( Figure 8H ). V1V2 responses of the specific IgG subclass, IgG3 which is known to be a potent mediator of Fc-effector functions such as antibody dependent cellular cytotoxicity were inverse CoR in RV144, and only IL-2+IFN-γ+ in T1 was significantly correlated with IgG3 gp70–96ZM51 V1V2 ( Figure 8F ). IgA specific for gp120 overall as well as the V1V2 region was a CoR in RV144, suggested to compete with the binding of protective IgG3 (3, 17). Only T2 and T3 had measurable IgA responses to AE.A244 gp120 or AE.A244 V1V2 ( Figure 8E ), with IL-2-IFN-γ- in T2 significantly correlating with IgA AE.A244 gp120 and IgA AE.A244V1V2 ( Figure 8B ). IL-21+ in T3 was significantly correlated with both IgA AE.A244 gp120 and IgA AE.A244V1V2 ( Figure 8C ). Overall these results indicate subtleties in the association of CD4+ T cell responses and plasma Ab responses, that are impacted by vaccine regimen and may provide insight into efficacy outcomes.