Optimising the β-parameter

The kernel is based on similarity scores between pairs of peptides. For each pair, a similarity score is composed of multiple sub-scores based on alignments between pairs of sub-sequences. The model parameter β regulates the relative influence that each sub-score will have on the cumulative score. In turn, it enables adjustment of the importance of sub-optimal alignments.

Experiments were undertaken to evaluate the effect on performance of varying the β-parameter using a simple test set. The MHCBench Set 4b was chosen for this purpose; it contains experimentally verified binders and non-binders of HLA-DRB1*0401. It consists of only natural peptides and an equal number of binders and non-binders (292 of each), which makes it well-suited for model testing.

The BLOSUM62 substitution matrix was used for S (·). This matrix is generally considered to be a good matrix for modelling evolutionary problems [34]. 10-fold CV was used to evaluate performance, and a rough search for a good β-parameter was undertaken. The effect on performance of varying β can be seen in Figure ​Figure1,1, which shows that the SKM is capable of distinguishing well between binders and non-binders. The best accuracy is 74.7% at β = 0.025. The best performance in terms of A_ROC is 0.827 at β = 0.035. Generally, the best results for most measures are found for β-values between 0.02 and 0.04. Higher β values of 0.2 to 5.0 were also evaluated, but as β becomes larger performance degrades for all measures.

Open in a separate window

Figure 1

Evaluating performance for varying the β-parameter using 10-fold CV. Graphs are plotted for accuracy (proportion of correct predictions), sensitivity (proportion of false negatives), specificity (proportion of true positives), and A_ROC (area under receiver operating characteristic curve).

Bootstrapping using case resampling [35] was performed to analyse the variance in results. 100 repetitions were undertaken, with data set sizes of 584. At a β-value of 0.025, which produced the best accuracy in the tests referred to above, the average bootstrapping accuracy was 73.1% with a standard deviation of 2.2%. This degree of variance was found throughout our experiments.

Interestingly, a low β indicates that the best solution is found when sub-optimal alignments have a large influence, as seen by the mathematical formulation below; when lowering the β-value, the relative contributions from scores of various alignments are evened out. This tendency was observed throughout the remaining experiments, with most "optimal" β-values being below 0.1. Interestingly, the same observation regarding the positive influence of sub-optimal alignments was also reported in Nielsen et al [27] using a very different method.

Selecting the substitution matrix

A substitution matrix was used in the calculation of the Smith-Waterman score to evaluate similarities between amino acids (see Eq. 0.3). The BLOSUM62 matrix initially used is regarded as a good descriptor for evolutionary problems [34]. However, using an alternative substitution matrix could prove more effective.

The AAIndex database [23] contains a large collection of substitution matrices produced during the last three decades. The matrices are based on numerous different measures, such as physicochemical properties and structural differences. An extensive search among the substitution matrices was conducted. Each substitution matrix was used instead of the BLOSUM62 matrix as above. Due to the scale of the experiment, only a crude search using 5 values of β was evaluated per substitution matrix. From the 83 matrices, the ten best performing substitution matrices with regard to A_ROC scores were retested with a refined search for the best value of β. The 3 best performing substitution matrices from this experiment are shown in Table ​Table22.

As can be seen from the table, the three matrices have very similar A_ROC values. The best performance was produced by a recently developed substitution matrix SM_THREADER_NORM, which is based on molecular mechanics force fields. [36] suggest that force fields can provide more reliable mutation matrices because of the incorporation of natural weighting of different physical contributions. Interestingly, the BLOSUM62 matrix is among the best three matrices out of 83. This suggests that the evolutionary rationale behind BLOSUM62 is also appropriate for MHC peptide similarity or that the chemical similarities underlying protein evolution also underlie peptide selectivity by the MHC. Such conjecture is supported in part by the fact that the SM_THREADER_NORM is also placed in the same family of substitution matrices when assessing the magnitude of distances between matrices [36]. In the following experiments, the SM_THREADER_NORM is used.

Performance on HLA-DRB1*0101 and HLA-DRB1*0301

Two MHC Class II alleles from the MHCBN database [2] were evaluated. The MHCBN database contains 475 binders and 105 non-binders for HLA-DRB1*0101 and for HLA-DRB1*0301 contains 219 binders and 150 non-binders. Duplicates and peptides with 75% or more Alanines were removed. 5-fold cross-validation was undertaken (5 fold CV used instead of 10-fold CV for comparison with [21]), and a crude optimisation of the β-parameter was performed as described in the methods section. Results are shown in Table ​Table3.3. The table reports comparison results from two other methods, a linear programming model "LP_top2" [21] and TEPITOPE, a quantitative matrix method [37].

The results in Table ​Table33 show that the SKM method performs significantly better than LP_top2 and TEPITOPE. The relative improvements in A_ROC scores are 8% and 41%, and the improvements in A_OVER-ROC are 37% to 60%. Other methods have also been evaluated on the data sets. In [25], results on two methods using Hidden Markov models combined with successive state splitting are reported. The best 10-fold CV results were 0.85 (S-HMM1) and 0.89 (S-HMM2), which are close to but still lower than the 0.91 A_ROC of SKM using 5-fold CV; using 10-fold CV the SKM performance increases to 0.93 A_ROC.

Performance on allele HLA-DRB1*0401

MHCBench [3] contains 8 data sets of binders and non-binders for HLA-DRB1*0401. Again, the SKM was evaluated using 10-fold cross validation, and a crude search for optimal values of β was performed for each set. A_ROC performance from all 8 sets are reported in Table ​Table4,4, which also includes results of PERUN, a method based on TEPITOPE [38], Neural Networks [7], Gibbs Sampler, a method based on Metropolis sampling [27], and LP_top2, a linear programming method [21].

The results on Table ​Table44 show that the SKM is significantly better on 6 out of 8 benchmark sets (Set 1 to Set 4b). On sets 5a and 5b, another method, LP_top2, scored the best results (0.815 A_ROC and 0.859 A_ROC, respectively). The relative lower accuracies of the SKM method on these two data sets may be due to the small training set sizes (117 and 85); in training, the SKM method selected nearly all training samples as Support Vectors, indicating there may not be sufficient samples to properly describe the model space.

Local alignment kernels

The principle of local alignments has been shown to provide a powerful approach to detecting relationships between sequences, using the optimal local alignment via the Smith-Waterman algorithm [46] and it's efficient (PSI-)BLAST approximations [47]. Therefore, we utilise a complex kernel comprising several sub-kernels based on the Smith-Waterman algorithm [32]. On a protein homology detection problem, this approach was found to significantly outperform scores based solely on optimal alignments [31].

Local Alignment Kernels are convolution kernels [48] consisting of a number of simple sub-kernels:

$\begin{array}{l}{K}_1\cdot K_2\cdot \mathrm{...}\cdot K_p(x,y)=\hfill \\ {\displaystyle \sum _{x=x_1\mathrm{...}{x}_p,y=y_1\mathrm{...}{y}_p}{K}_1(x_1,y_1)\cdot \mathrm{...}\cdot K_p(x_p,y_p)}\hfill \end{array}$

Eq. 0.1

Where the components K₁...K_p consists of three different kernels: (a) A constant kernel K_const, (b) a kernel for measuring the difference between aligned letters K_align, and (c) a kernel for penalizing gaps K_gap:

$\begin{array}{l}{K}_{const}(x,y)=1\hfill \\ K_{align}(x,y)=\{\begin{array}{cc}0& ,if\left|x\right|\ne 1or\left|y\right|\ne 1\\ e^{\beta \ast S(x,y)}& otherwise\end{array}\hfill \\ K_{gap}(x,y)=e^{\beta (g(\left|x\right|)+g(\left|y\right|))}\hfill \end{array}$

Eq. 0.2

where x and y are the amino acid sequences, S (x, y) the Smith-Waterman score, g (·) the gap penalty function, and β a scaling parameter to adjust importance of gaps and sub-optimal alignments.

The Smith-Waterman score SW_S,g(π) is calculated as:

$\begin{array}{l}{S}_{S,g(\pi )}={\displaystyle \sum _{i=1}^{\left|\pi \right|}S\left(x_{\pi 1(i)},y_{\pi 2(i)}\right)}-{\displaystyle \sum _{i=1}^{\left|\pi \right|-1}[g({\pi }_1(i+1)-{\pi }_1(i))+g({\pi }_2(i+1)-{\pi }_2(i))]}\hfill \\ S{W}_{S,g(\pi )}(x,y)=\underset{\pi \in \Pi (x,y)}{\max }{S}_{S,g(\pi )}\hfill \end{array}$

Eq. 0.3

Where π is the alignment between two sequences x and y, S (·) is a substitution matrix and g (·) a gap penalty function.

The component sub-kernels are combined by convolution to represent a kernel for an alignment of length n. The Local Alignment score is the sum over all possible alignments in the sequence:

$\begin{array}{l}{K}_{(n)}(x,y)=K_{const}\cdot {(K_{align}\cdot K_{gap})}^{(n-1)}\cdot K_{align}\cdot K_{const}\hfill \\ K_{LA}(x,y)={\displaystyle \sum _{i=0}^N{K}_i(x,y)}\hfill \end{array}$

Eq. 0.4

where N is the number of possible alignments.

The above formulation results in a computational complexity exponential with |x| and |y|, and is thus not a feasible solution for this problem. Hence, a dynamic programming algorithm by [32] is used, which is a slight modification of the Smith-Waterman algorithm [46]. The kernel computation is done in O (n²·|x|·|y|), where n is the number of samples, and |x| and |y| are the lengths of the peptide strings. In the context of MHC-peptide binding, the gap penalty term must be maximised since gaps are not possible within bound peptides.

Calculating the sub-kernels requires the following two parameters: the substitution matrix S (·) in Eq. 0.3, and the β parameter in Eq. 0.2. S (·) quantifies a similarity between pairs of amino acids, and is a well-known term in bioinformatics with numerous substitution matrices designed from evolutionary, physicochemical or structural properties (a list can be found at [23]). The β parameter regulates the effect of individual contributions from alignments, and allows adjustment of the relative importance between low- and high scoring alignments. When β is low, the model will increase the importance of low scoring (sub-optimal) alignments to quantify the similarity between sequences. Similarly, as β → ∞ the contributions from sub-optimal values is reduced.

All kernels must be symmetric and positive semi-definite. Some values of β and substitution matrices S in Eq. 0.2 resulted in invalid kernels, and caused convergence problems. A trick used in [32] subtracts the smallest negative Eigenvalue from the diagonal of the kernel to ensure kernels are positive semi-definite. Symmetry of the resulting kernel, K_LA, is guaranteed as long as substitution matrix S (·) is symmetric.

Data set

Multiple data sets were used in the experiments. Eight benchmark data sets with samples of known binders and non-binders of the HLA-DRB1*0401 allele were taken from MHCBench [3]. Within the data sets, peptide strings are assigned binding strengths of level 0 (non-binders) to 4 (strong binders), and are collected from multiple sources, mainly MHCPEP [1]. The sets are derived from the same base set of peptides, created with varying levels of curation. Set 1 includes all peptides whereas Set5b is a homology reduced set containing only natural peptides. Data set sizes range from 85 to 1017 samples with peptide lengths of 9 to 33. As these sets have been used in many published experiments [21,26,27], we use them in preference to alternatives.

In addition to the 8 sets, two data sets from specific alleles HLA-DRB1*0101 and HLA-DRB1*0301 were taken from MHCBN [2]. The data sets separate peptides into binders (low, moderate, and high), and non-binders (peptides having IC₅₀ values of more than 50,000 nm). Finally, binders from multiple alleles from the MHCPEP database [1] were used. In the database, peptides are labelled as having low, moderate or high binding affinity.