Amino Acid Physical Properties

Source data for the amino acid physicochemical properties were obtained from the proteomics resource at the Swiss Institute for Bioinformatics. Attempts were made to balance the types of physical data used and to obtain independent measurements of related physicochemical properties from two or more different studies. A total of 31 sets of physical parameters were tabulated (Additional File 2; Table S2) and the principal components (PC) computed from the correlation matrix. Results of this analysis are shown in Figure 1. The first three principal components account for approximately 90 percent of the variance in the 31 datasets. Table 1 is a listing of the 20 amino acids found in proteins, sorted independently by each of the first three principal components. It is apparent that the first PC is correlated with the polarity or hydrophobicity of the amino acid, and the second PC with molecular size. The underlying physical correlate of the third principal component is less obvious, but is generally considered to be electronic in nature [15]. The two sulfur-containing amino acids (cysteine and methionine) and histidine, with an imidazole ring, are ranked at one end of the third PC. Other amino acids with aromatic rings are ranked close behind. The actual z-scale numbers are slightly different from those originally reported by Wold et al [1, 16], likely because we used a significantly larger dataset to develop the PC.

PC are mutually orthogonal and are effectively uncorrelated proxies that embody the effect of the physical properties of amino acids found in proteins. An additional property characteristic of PC is that they are numerically appropriately weighted for their relative contributions to fitting a multivariate regression.

Human MHC-I and MHC-II Datasets

Extensive datasets of binding affinities of synthetic peptides have recently been made publicly available at the Immune Epitope Database (IEDB) [17, 18, 19]. Because we wanted to use the web-based programs NetMHCII and NetMHCIIPan at the Center for Biological Sequence Analysis (CBS) as comparators, we downloaded MHC-I and MHC-II databases from CBS [20]. Although the principles can readily be expanded to peptides of other sizes, the work described herein is restricted to 9-mers for MHC-I and 15-mers for MHC-II. The datasets downloaded contained binding affinity measurements as ic50 values for 29,336 9-mers for 35 alleles of MHC-I and 9,117 15-mers for 14 alleles of MHC-II (Additional File 3; Table S3a. and Additional File 4; Table S3b.). The binding data were natural logarithm-transformed ln(ic50) and the distributional characteristics of the binding data examined. The paper that is the primary source of the MHC-II dataset does not describe the experimental strategy for choosing the various peptides [19] but clearly the experimentalists designed the peptides in the datasets based on prior knowledge of binding characteristics. Several features of the datasets were noted. All of the datasets consist of several subsets. For example with the MHC-II allele datasets there are three major subsets: 232 peptides (232-subset) have been tested in all 14 HLA MHC-II supertypes, 167 (167-subset) with 11/14 supetypes, and 49 with 8/14 supertypes. In addition to these larger subsets there are subsets with fewer peptides that have been tested in various combinations of alleles and some of these have quite high average affinities. For example, one subset has an e^4.1 higher average affinity than the 232-subset. The amino acid combinatorials are different among the different subsets and they have statististically different means and variances. The experimentally measured average binding affinities of the 167 peptide subset is consistently greater than the 232 peptide subset in every one of the alleles where both have been tested. The difference between the means of the 232-subset and the 167-subset is e^2.2 (167-subset has statistically significant higher affinities) and among the different alleles the ratios of the means range from about 3-fold to 60-fold. A second feature apparent in the datasets is that the frequency distributions of the datasets have some anomalies; several are shown in Figure 2. The three highlighted bars in the histogram in Panel 2a. each contain a large fraction of the data points all ascribed exactly the same number 124/166 (1 nM), 701/800 (20,000 nM) and 168/195 (78,125 nM). Measurement of low affinity binding, where very small fractions of the total ligand are bound, are sometimes given identical scores. This is the assay experimental limit. Wang et al [19] do not indicate the source and history of the datasets; it is possible that they are composites of measurements made at different times and under different laboratory conditions. The reason for the cluster of high affinity binders is likewise unclear, but is also likely to represent an experimental limitation. Wang et al [19] state that they capped the binding affinities at 50,000 nM which is at odds with what is found in the datasets. Panels 2b. and 2c. are from the smaller MHC-II datasets; these have fewer anomalies than seen in Panel 2a. Clearly peptide synthesis was done by systematically varying different positions in the peptides, a process which likely gave rise to the similarities described by El-Manzalawy et al [14], who systematically evaluated the effects of similarity reduction on a number of commonly used datasets and concluded that regardless of the method employed, the performance of classifiers was substantially below that reported.

Use of Principal Components of Amino Acids in PLS Analysis of MHC Binding

Peptides of 9-mers and 15-mers were converted into 27 (9 × 3) and 45 (15 × 3) vectors of principal components. This type of descriptor is commonly used in QSAR analysis where it is known as the "z"-scale [3, 21]. These z-scales were used in a partial least squares (PLS) analysis of the binding data. Each amino acid in a 9-mer is replaced by three z-scale descriptors. {z₁(aa₁), z₂(aa₁), z₃(aa₁)}, {z₁(aa₂), z₂(aa₂), z₃(aa₂)}... {z₁(aa₉), z₂(aa₉), z₃(aa₉)}. A 15-mer for MHC-II analysis has a correspondingly larger set of descriptors.

PLS Binding Affinity Predictions

The z-scale descriptors were used in PLS to develop prediction equations of peptide binding by MHC-I and MHC-II. This is essentially the approach pioneered by Hellberg et al [22], which was the genesis of the use of QSAR techniques for peptide binding. For evaluation of multivariate datasets such as these, the variable importance projection (VIP) is a very useful method to reduce the dimensionality of the data and to produce a single metric that ranks the relative importance of a particular predictor in the overall response (2,3). The VIP is an output generated using all the experimental measurements in the MHC datasets (29,336 MHC-I and 9,117 MHC-II) and the relationships among them. Any VIP metric with a value greater than 1 is considered to be of importance in the overall prediction. A by-product of the symmetry between the statistical models for the MHC binding interaction and the underlying binding reaction itself is that it possible to visualize and conceptualize the physicochemical properties of the interaction. The matrices of VIP values for all the alleles of both sets of HLA molecules are found in Additional File 5; Tables S4a and Additional File 6; Table S4b. The PLS models for all of the MHC-II alleles show evidence of significant latent factors, with twelve of the alleles having evidence of two latent factors and two of the alleles having three [2, 3, 23]. Through experimentation (data not shown) it appeared that the latent factors arise as a consequence of the clustering of different subsets having different means as discussed above. It is expected that the least squares processes of PLS should be less sensitive to the anomalies in the datasets than PSSM based methods.

Neural Network Regression using PCAA

As discussed by Bishop [24] a NN regression of data with continuous responses is effectively a non-linear PLS. In fact, the predecessors of the JMP^® platform we used, the statistical analysis programs of SAS http://www.SAS.com have had extensive neural network capabilities for several decades and treat neural networks as simply a special type of regression analysis. Bishop also discusses the use of principal components as a useful, if not essential, adjunct to building appropriately weighted reliable neural networks.

Thus, we submitted the MHC ln(ic50) binding data used in PLS above to the NN platform within JMP^® using the z-scale vectors to predict the ln(ic50)) in a simple three-layer perceptron (Figure 3). The non-linear equation that was used by the NN platform for fitting was the logistic function 1/(1+e^x) that is a well behaved smooth function that asymptotically approaches both a minimum and maximum and at the limit is practically a straight line. One of the key decisions in the design of the NN is to determine the appropriate number of hidden nodes. We chose to use a number of hidden nodes equal to the number of residues in the peptide that interacts with the MHC binding domain (9 hidden nodes for MHC-I and 15 hidden nodes for MHC-II). The rationale was that this symmetry would facilitate comparison of the fit obtained by this method with that described above for the linear PLS fit. It also provides symmetry between the statistical model and the underlying physicochemical model. It should be noted that this number of hidden nodes is actually a rather small number compared to that used by Nielsen and Lund [20] for their substitution matrix NN NetMHCII, used to benchmark our predictors. The JMP^® platform has a variety of mechanisms and statistical output for "training" of the NN, in order to control the underlying non-linear regression convergence, to assess the statistical reliability of the output, and to monitor and control overfitting through the use of an overfitting penalty coefficient. We systematically experimented with these control elements to evaluate the quality of the predictions through several cross validation strategies. We found that the presence of peptide subsets with different numbers of peptides, some having radically different mean affinities in the predictors (detected as latent factors in the PLS), are also somewhat problematic for random selection of training subsets during cross validation. The results of two different strategies are reported here. The two different models are referred to as Method 1 and Method 2.

In Method 1 multiple "tours" (different random seeds) of a random holdback strategy were used. Examination of the residuals in the various hyperplanes was used to examine the residuals of these fits. In as much as the three principal components we used for the model account for approximately 90% of variance in the underlying physical properties, we set the overfitting penalties to target an r² of 0.9. For benchmarking, the prediction models the IEDB datasets downloaded from CBS were contemporaneously submitted to the webservers for NetMHCII (version 2.0) and NetMHCIIPan (version 1.0) at CBS [25, 26, 27, 28]. The performance of Method 1 is compared to the PLS model and the output of the servers at CBS in Table 2. As described above for the PLS, both an r² comparing the fit and a categorical transformation were used to make the comparisons.

The predictions produced by Method 1 and its ability to generalize in the training sets compared favorably to NetMHCII (Table 2) evaluated either as a continuous fit or as a categorical classifier. The statistical metrics associated with the model suggested that some overfitting was likely occuring with this model and therefore a second method (Method 2) was developed.

Having five prediction methods based on different underlying predictors, substitution matrices for NetMHCII and NetMHCIIPan and physical properties of amino acids for PLS, Method 1 and Method 2 described above provided an opportunity to examine the comparative performance of the different prediction methods with both the IEDB training sets as well as with other peptides. This was done by creating a test set of 1000 15-mer peptides selected at random from the proteome of Staphylococcus aureus COL (Genbank NC_002951). This random test set was submitted to each of prediction tools and the results tabulated for comparison. Figure 4 shows the results of comparisons of the different methods with Method 2 as the base method, using the Pearson correlation coefficient of the predictions as the metric for comparison for the training sets (see additional file 7 for detail). Method 1, NetMHCII and NetMHCIIPan all produce highly correlated predictions, the highest correlations being between Method 2 and NetMHCII. The results of evaluation using categorical predictors gave comparable results (not shown).

As with the training set, the correlated response of between Method 2 and Method 1 is also seen for the random peptide set. Table 3 also shows the comparison of Method 2 with both the training set and the random set. Interestingly, with the random set the correlation with PLS is substantially better than for the training set, however the correlation between Method 2 and both NetMHCII and NetMHCIIPan is diminished. Also, the correlation coefficients of the later two prediction methods show a higher degree of variability.

MHC-I predictions

Similar comparisons to those described above for MHC-II were conducted with the MHC-I datasets (results not shown). Distributional issues are very common in the MHC-I datasets, like those shown in Figure 2A., in which a clustering of identical ic50 values occur probably as a result of experimental limitations. For the MHC-I, final prediction r² values are somewhat lower, in the range of 0.7 - 0.85. Visualization of the VIP from the PLS analysis is shown below. As seen with the MHC-II alleles, the PLS analysis suggested the presence of several latent factors in the datasets.

Visualization of the Variable Importance Projection

Computation of the VIP provides a means of estimating the relative importance of predictors to the overall model and as such is a very useful aspect of a PLS analysis. In the case of the peptide datasets used here, where the predictors are conceptual physicochemical proxies, the VIP provides a statistical view of the underlying physical interactions between the peptide and the MHC binding domains. As a visualization paradigm for the VIP we use a "heat plot", commonly used in other areas of systems biology. In these relatively simple plots, each cell represent tens of thousands of statistical comparisons reduced to a single metric and displayed in a uniformly scaled graphic system that facilitates visual comprehension [29]. All of the plots in Figures 5, 6, 7, and 8 are produced from two VIP matrices (the matrices can be found in the Additional File 3; Tables S3a and Additional File 4; Table S3b), one for MHC-I binding and the other for MHC-II, that have been colorized in different ways to accentuate different potential physicochemical properties and to summarize all binding reactions of all peptides in the datasets used to derive the predictions described above. Any value in the VIP matrices >1 are interpreted as being the most relevant in explaining the binding affinity. In Figures 5, 6, 7, and 8 the matrices are colored by three different methods which provide different views of the underlying physicochemical properties and interactions between the peptide and binding site. Each of the underlying physical properties are on individual plots; Additional File 8; Figures S5, S6, S7, and S8 provide copies of Figures 5, 6, 7, and 8 that include the "thermometers" for the heat plots that detail how the numerical information in the VIP matrix was color scaled. We designate the MHC binding pocket to consist of 9 zones of interactions given the label P1 to P9 based on MHC-I. Peripheral to this, the MHC-II molecules comprises 3 additional residues on each side (N-3, N-2, N-1 and C+1, C+2 and C+3) to make up the 15-mer. In MHC-II molecules the ends of the binding pocket (groove) are open and longer peptides can interact with additional residues.

Figure 5 shows the VIP of the PLS regression prediction of ln(ic50) of peptide binding by using the first three principal components of the amino acids in each of the amino acids in the 9-mer as predictors for MHC-I (5a) and 15-mer for MHC-II (5b). In these Figures the coloration of the plot is uniform over all cells in the matrix (complete scaling information is found in the Additional File 8; Figure S5). The colors compare the relative importance of the particular numbered residue of the binding domain among all of the MHC alleles indicated. As the selection of peptide sequences were highly non-random it is likely that the coloration seen is a combination of two factors, the prior knowledge of the experimentalists use in designing the amino acid combinatorials in the peptides and the actual physicochemical properties. The coloration is consistent with the amino acid composition thought to provide the best fit into the binding groove. The VIP provides an insight into the potential physicochemical interactions in the binding groove of MHC molecules which have not been studied as intensely. The relative importance of the polarity of position 9 for the MHC-I locus A binding is seen. PC3, the third ranked principal component electronic physical correlate, emerges as important in position 2 of MHC-I locus B alleles. In Figures 6 and 7 the coloration emphasizes a different aspect of the binding interaction of MHC-I (Figure 6) and MHC-II (Figure 7). In these two Figures the color scaling compares the relative importance of the particular numbered residue of the binding domain among the all of the alleles indicated. Again, cells in the VIP matrix with a value >1 are the most relevant in explaining the binding affinity. Coloration is column-relative for each position in the binding domain. One position of note in Figure 7 is the strong signal of PC1 (the polarity correlate) for C+1 position for DRB3*0101. Giving credence to the properties visualized, a recent paper of Parry [30] using an entirely different methodology, suggested a major role of P10 (the binding zone located one amino acid position toward the C-terminus and equivalent to C+1 in our notation) in the binding of peptides to DRB3*0101. For different MHC-II alleles the possible role of each of the three physical properties in the overall binding reaction of those N-terminal and C-terminal residues outside the core binding pocket is suggested by these Figures.

Figure 8 provides a visualization only for MHC-I; output for MHC-II is included only in the Additional File 8; Figure S8. Coloration scaling is column-relative for each position in the binding domain. The comparison is within MHC-I allele and is repeated for each of the 3 principal components. The colors compare the relative importance of the physicochemical property of each particular numbered residue of the binding domain within each of the MHC-I alleles indicated.

Using the JMP^® "contour profiler" it is possible to segregate and visualize of the interaction between any two of the input variables and the output prediction. An example is shown in Figure 9, however a two-dimensional representation does not adequately illustrate the fully interactive demonstration of the physicochemical interactions involved in the binding reactions. Figure 9 is a plot showing a hyperplane of the first principal component predictions for amino acid P9 and C+1 (also called P10 by Parry [30]). It shows not only how the polarity of these two residues might play a role in binding, but also how the experimentalist might systematically vary the amino acids in these positions and what would be the expected impact on binding affinity. In the perspective depicted in Panel A the non-randomness of the amino acids used in the training sets is readily apparent. Panel B provides a means of examining of the relative scatter of the actual datapoints around the fitted hyperplane in multidimensional space.

Data sources

Amino acid physical properties were retrieved from the repository at the Swiss Institute of Bioinformatics: http://expasy.org/tools/protscale.html

MHC-I and MHC-II peptide binding datasets were retrieved from the repository at Center for Biological Sequence Analysis http://www.cbs.dtu.dk and IEDB http://www.IEDB.org.

A random dataset of 1000 15-mer peptides was selected from the surfome and secretome of Staphylococcus aureus COL (Genbank ID NC_002951) (Additional File 9; Table S6).

Internet data analysis and manipulation

Web servers for benchmark comparisons NetMHCII version 2.1, NetMHCIIPan 1.0 http://www.cbs.dtu.dk/services/NetMHCII/. Some testing was done with the MHCserver but there were no extensive or systematic comparisons to NetMHC 3.0.

Workstation applications

The Partial Least Squares and Neural Network platforms in JMP^® 8.02 and JMP^® Genomics 4.2 were used for all local calculations and data manipulation (SAS Institute Inc. Cary, NC).

The statistical background for the JMP^® platforms is found in: JMP^® 8 Statistics and Graphics Guide, SAS Institute Inc., Cary, NC, USA (ISBN 978-1-59994-923-9).

Variable influence

For a given PLS dimension a, (VIN)_ak ² is equal to the squared PLS weight (w_ak)² of that term, multiplied by the explained sums of squares (SS) of that PLS dimension [2, 3]. The accumulated (overall PLS dimensions) value,

is then divided by the total explained SS by the PLS model and multiplied by the number of terms in the model. The final VIP is the square root of that number. The formula can also be expressed as follows:

The SS of all VIP is equal to the number of terms in the model and hence the average VIP is equal to 1. The VIP of one term can be compared to others and terms with the larger VIP, larger than 1, are the most relevant for explaining the response. The VIP were computed within the PLS platform of JMP^®.