Lessons learned in the analysis of high-dimensional data in vaccinomics

Abstract

The field of vaccinology is increasingly moving toward the generation, analysis, and modeling of extremely large and complex high-dimensional datasets. We have used data such as these in the development and advancement of the field of vaccinomics to enable prediction of vaccine responses and to develop new vaccine candidates. However, the application of systems biology to what has been termed "big data," or "high-dimensional data," is not without significant challenges-chief among them a paucity of gold standard analysis and modeling paradigms with which to interpret the data. In this article, we relate some of the lessons we have learned over the last decade of working with high-dimensional, high-throughput data as applied to the field of vaccinomics. The value of such efforts, however, is ultimately to better understand the immune mechanisms by which protective and non-protective responses to vaccines are generated, and to use this information to support a personalized vaccinology approach in creating better, and safer, vaccines for the public health.

Keywords: Data interpretation, statistical; Immunogenetics; Systems biology; Vaccination; Vaccines.

Figure 1

(a) Minus versus Average (MVA) plot demonstrating the effect of change in reagents and sequencing software. There is one data point for every feature measured on the assay. The x-axis is the average of each feature over all specimens in the study. Generally, the y-axis is the difference of each feature from the mean. Thus, if the observations are identical to the mean, all data points would lie on the y=0 line. Here, the y-axis is the difference between the before reagent change mean and the after reagent change mean. A reference line for y=0 as well as a loess smoother are included on the plot. If the smoother overlays the y-0 line, no normalization is needed. If the smoother is parallel to the y=0 line but shifted up or down, this indicates that between specimen biases are similar for all abundance levels and a linear normalization is needed. The nonlinear smoother demonstrates that nonlinear bias is present. (b). The same study shown after normalization and filtering out genes with median count <32. The fact that the smoother is now straight and lies on the y=0 line demonstrates that the bias has been removed.

Figure 2

Box-and-whisker plots showing global distribution of per-gene counts on the log scale (y-axis) by lane (x-axis) sorted by assay order. Top, mid-line and bottom of boxes indicate 75^th, 50^th and 25^th percentiles, respectively. (a) Pre-normalization. The total counts/lane increased from ∼150million to ∼200million after reagent and software upgrades. This is evident from the general shift up approximately two-thirds of the way across the plot. A failed specimen with median nearly half that of the neighboring specimens is evident about one-third of the way across the plot. The failed specimen was deleted in subsequent analyses. (b) Post-normalization. After normalization via Conditional Quantile Normalization (CQN)[82], the distributions of the specimens are aligned exactly at the maximum, 75^th and 50^th percentiles as expected. The lower counts are not exactly aligned since the smallest counts are not adjusted in CQN.

Figure 3

Over 450 PBMC specimens from healthy subjects aged 50-74 years old on were assayed on five bead array plates of the Illumina DNA methylation 450K assay. The assay utilizes two probe designs, each yielding an M and U intensity value (fluorescence intensity of methylated or un-methylated cells, respectively). These intensity values are mathematically combined to create an estimate of the percent methylation (β-value) in the specimen. (a) There is evidence of nonlinear between-specimen biases in the M and U expression intensities as demonstrated by these residual MVA plots. Each smoother represents one specimen. Nonlinearities are evident. (b) Between-specimen biases are near linear on the beta-value scale (left), are not large, and are essentially eliminated via this strategy (right).

Figure 3

Over 450 PBMC specimens from healthy subjects aged 50-74 years old on were assayed on five bead array plates of the Illumina DNA methylation 450K assay. The assay utilizes two probe designs, each yielding an M and U intensity value (fluorescence intensity of methylated or un-methylated cells, respectively). These intensity values are mathematically combined to create an estimate of the percent methylation (β-value) in the specimen. (a) There is evidence of nonlinear between-specimen biases in the M and U expression intensities as demonstrated by these residual MVA plots. Each smoother represents one specimen. Nonlinearities are evident. (b) Between-specimen biases are near linear on the beta-value scale (left), are not large, and are essentially eliminated via this strategy (right).

Figure 4

Power to detect genetic associations as a function of ordinal genotypic effect size for three different analyses, and with two different levels of significance. 1000 data sets were generated for each combination of parameters. Panel (a) shows statistical power for α=0.05 and panel (b) shows statistical power for a genome-wide significance threshold (α=5×10^-8).

Figure 5

Bias in estimating an ordinal genotypic effect, as a function of the simulated ordinal genotypic effect size for three different analytical approaches.