Predictor for the Effect of Amino Acid Composition on CD4+ T Cell Epitopes Preprocessing

2.1. statistical dispersion

The comparison between the dispersion of diAA distributions in the peptide body and the cleavage sites was carried out by two steps:

• To get the diAA distributions, the observed/expected ratio of diAA probability was measured for each AA ( $\frac{P\left(\text{di}A{A}_{ij}\right)}{P\left(A{A}_i\right)·P\left(A{A}_j\right)}$, where i is to the left of the site and j to the right).

• The statistical dispersion of these distributions was then compared by using Average Absolute Deviation (AAD) from the median.

The average absolute deviation (AAD) of a data set is the average of the absolute deviations from the median and is a summary statistic of statistical dispersion or variability. Lower AAD means low dispersion and convergence to specific value.

2.2. Training and validation sets

Epitopes from the IEDB database were used for the positive datasets. We extracted 1027 different epitopes classified as ligand elution (Vita et al.) in the IEDB database (Vita et al.). For each peptide, we found its origin in the appropriate genome and extended it to include the flanking regions at its N and C termini. Peptides that had no explicit origin or peptides where an appropriate origin could not be found were excluded from the data set. As a negative data set, we extracted all the peptides in the source proteins which are not included in the epitope dataset. Out of this pool, we pulled a subset of 20,000 peptides with a length distribution similar to the length distribution of the positive data set (Supp Mat. Figure S1). No peptide was found in more than one dataset. The positive data set was divided into two parts: 856 sequences for the learning set and distinct 171 sequences for the validation set. We also used a database of 42 sequences of cleavage sites for Cathepsins L and S extracted from the MEROPS database (Rawlings and Barrett, 1999) as an in vitro cleavage positive learning set. We have also used an external test set for further analysis. This additional test set was not used for the learning process or for the choice of the optimal model. This dataset contains a set of 3862 epitopes inducing a T cell response from the IEDB (Vita et al.). The goal of this set is to see whether peptides measured to induce a T cell response are indeed naturally processed.

2.3. Simulated annealing

Optimal weights for each diAA were learnt using a Simulated Annealing (SA) process. The initial configuration of the learnt weights (20 AAs × 20 AAs = 400 parameters) was initialized at zero. The initial temperature was set to T₀ = 5, and was decreased exponentially once during every cycle (T_n =λⁿ·T₀, λ= 0.8). Each step was a random change of a weight by a random factor uniformly distributed between 2.5 and −2.5; the value of the weights was limited to the domain: 6 to −6. We performed 60 steps of decreasing the temperature; each step composed of 300 cycles per parameter.

2.4. SVM

For the classification of preprocessed peptides, we used a support vector machine (SVM) algorithm with linear and quadratic kernels (Boser et al., 1992; Müller et al., 2001). The SVM finds a maximal margin separator between two data sets with either a linear or quadratic distance kernel. In both cases, the error cost (C) values were set between 10⁻⁶ – 10⁰ and the positive to negative learning set size ratio was varied between 1:1 and 1:3 (through the increase of the negative learning set size). The observations used in the SVM were AA or diAA frequency vectors. The size of each vector was 1×20 and 1×400 respectively for each position. Each domain of the peptide has its unique frequency vector. The AA composition of the peptide was computed using an occupancy vector, where the appropriate AA (or diAA) value was the number of times it appeared, and all others had values of 0. For example, in the diAA formalism, an Ala-Cys combination was represented as a 400 dimensional vector with a value of 1 at position 21 (second row first column) and all other values as 0.

2.5. Performance evaluation

In order to assess the performance of the classifiers, we used the common statistical measures of sensitivity, specificity and accuracy. The sensitivity was computed as the true positive rate, i.e., the preprocessed peptides that were classified correctly (true positive) among all preprocessed peptides either classified correctly (true positive) or classified incorrectly (false negative). The specificity was computed as the true negative rate, i.e., the non preprocessed peptides that were classified correctly (true negative) among all non preprocessed peptides either classified correctly (true negative) or classified incorrectly (false positive). As the real prevalence is not known due to lack of experimental negative data, we used a prevalence of 0.5 in the accuracy calculations.

2.6. Propensity scales analysis

In order to identify properties that correlate with the weights of the cleavage score, we used the Pearson correlation coefficient to selected propensity scales extracted from the AA index database (Kawashima and Kanehisa, 2000). Multiple SVM based classifiers were tested. In all classifiers, we treated the C and N termini as one region by unifying the appropriate weight vectors. Matrices of 20×20 elements were reduced to 20 elements vectors by averaging columns and rows. The resulting averages were compared to the propensity scales to detect the biochemical elements affecting the classification.

2.7. Feature selection

We performed a feature selection using the following score: $\frac{|{\mu }_{\mathit{pos}}-{\mu }_{\mathit{neg}}|}{\sqrt{{\sigma }_{\mathit{pos}}+{\sigma }_{\mathit{neg}}}}$, where μ_pos is the average value of the positive set for a given position and σ_pos is the standard deviation in the same position, μ_neg and σ_neg are the parallel for the negative set. This criterion normalizes the absolute value of the distance between the mean frequencies of specific AA (or diAA) in the positive and negative datasets by dividing it with the square root of the sum of standard deviations of specific AA (or diAA) in the positive and negative datasets. Higher values are more meaningful for the classifier, since they represent clear differences between the positive and negative sets. Note that this score was obviously only computed on the learning set and not on the test set, to avoid any information leakage.

To further reduce the dimensionality, we used both Principle Component analysis (PCA) and Partial Least Squares (PLS) algorithms. Both methods project the data on an alternative space composed of components that have either the largest possible variance (PCA) or the maximal covariance with the labels (i.e., cleaved/not cleaved) (PLS).

3.1. Amino acids distribution

We have first tested if the amino acid usage of cleaved peptides contains information about the peptide cleavage probability. If the peptide and flanker sequences do not contain such information, sequence based classifiers will probably fail. The information content was estimated by comparing the differences between the AA composition in peptide bodies, at the flanking regions and outside peptides (Figure 1). The distributions of single AAs in each region were calculated as: A) their frequency inside peptides (except for the first and last positions); B) their frequency around the C and N termini of the peptide (the first position surrounding the cleavage site inside the peptide and K positions outside the peptide); C) the single AA frequency in the same proteins outside the cleaved peptides. The same proteins were used in order to avoid biases that may result from the choice of proteins. Similarly, the distribution of diAA was calculated as the frequency of diAA in the above mentioned regions. A clear difference is observed in the single AA distribution around the cleavage site (Figure 1 middle panel) (Chi square test p < 5·10⁻⁵ for all comparisons), within peptides and outside peptides (Figure 1 upper panel) (chi square test p <2·10⁻¹³¹). On the other hand, only minor differences are observed between the single AA composition at the second position after the cleavage site and in random regions farther away from the cleavage site (Figure 1 lower panel) (chi square, p< 6·10⁻⁷). Given the limited effect of these regions, we ignored all positions beyond the first flanking residue (K=1) in the learning score.

Beyond the differences in the single AA composition, a correlation between the diAAs may affect the cleavage probability. Such an effect can be determined through the dependence between the diAA composition distributions at neighboring sites. We measured this dependence using the mutual information (I) and Shannon Entropy (H) of each position. The mutual information between the AAs at the two sides of the cleavage site was higher than within the peptide (I_N= 0.2726 out of H_N = 5.3889 and I_C = 0.2904 out of H_C = 5.3230 compared with I_PB = 0.0600 out of H_PB = 5.6438), suggesting that preferable cleavage sites exist.

In order to get an insight on the source of the effect of diAA on cleavage, we computed the ratio between the observed diAA probability and the expected probability if one assumes that the left and right positions around the cleavage site are independent around cleavage sites, the distribution of hydrophobic amino acids differs from the expected behavior; while negatively charged and aliphatic amino acids behave as expected. In the peptide body, the behavior of the vast majority of amino acids is much closer to the expected (Supp Mat. Figure S2). In order to measure the statistical dispersion of the above distributions over the cleavage site and the peptide body, we used the Average Absolute Deviation (AAD) from the median (see Methods). Over the cleavage site the AAD was 0.595, while over the peptide body the AAD was 0.398. These values indicate a greater dispersion over the cleavage site than in the peptide body and emphasize the importance of the cleavage site to the characterization of a possible classifier.

To summarize, both the single and amino acid usage and the diAA usage are different around the cleavage site, within the epitope and outside the epitope. Moreover, the AA distribution at the two sides of the cleavage sites is not independent. We further checked whether the different AA distribution can be used to learn a classifier for phagosomal and lysosomal cleavage.

3.2. Single amino acid scores

MHC class II epitope precursor cleavage probably takes place in multiple compartments. We here denote it as “endosomal” cleavage as a simplification. The cleavage process itself may be affected by the full structure of the protein, since proteins are not fully denatured before their cleavage in endosomal compartments (in contrast with proteasomal cleavage in the cytosol) (Collins et al., 1991). However, in the current analysis, we are only interested in relatively short peptides (up to 30 amino acids) that were observed to bind MHC class II molecules (Sercarz and Maverakis, 2003). We assume that within such peptides the secondary structure leaves most AAs exposed to cleavage, and we thus use a sequence-based and not a structure-based score. Note that we do not propose a cleavage site prediction algorithm; rather we proposed a combined score to classify full peptides as either the result of endosomal cleavage or not. The score is based on the peptide sequence as well as the C and N termini flanking residues. We used the end-product of the endosomal process (i.e., naturally processed CD4+ T cell epitopes) as a learning set for the cleavage process.

Cleaved peptide must theoretically fulfill three conditions: cleavage at the C and N termini and no cleavage within the peptide body. We thus propose a score based on the residues in the body of the peptide and in the C and N termini flanking residues. The body of the peptide was defined as all positions, except for the first and the last positions in the peptide. The C and N termini flanking residues are the first positions within and outside the peptide surrounding the cleavage site. We propose a scoring method that contains these three elements: the contribution of the C terminal, the contribution of the N terminal, and the contribution of the body of the peptide (Figure 2). Formally, we assume that:

where P(C) and P(N) are the cleavage probabilities at the C and N termini, and P(pep_i) is the cleavage probability at position i of the peptide body. The independence assumption allows for the translation to a sum by taking the log of the two sides of the equation to produce:

The goal of the machine learning algorithms proposed in the following sections is to produce estimates of S_C(C) = ln(P(C)), S_P(pep_i) = ln(1 − P(pep_i)) and S_N(N) = ln(P(N)) (Figure 2). As a first assumption, we test the possibility that the cleavage probability at a given site is the product of the effect of the AA to the left and right of the cleavage site, which leads to:

where S_CL(AA_CL) is the contribution of the position to the left of the C terminal to the cleavage, and all other elements follow the same notation. If we further assume that the behavior at the C and N termini are similar, Eq 3 simplifies to:

where S_C/NL(AA_CL) is the contribution of the position to the left of the C terminal to the cleavage, but the score is similar for the C and N termini.

As is the case for many machine learning problems, a tradeoff must be obtained between the complexity of the score used and the precision obtained during the training stage. We thus tested classifiers with scores of varying complexity. The complexity of the scores is expressed either by the number of features used for the score or in the scoring function itself. The simplest possibility is a score incorporating only the AA distribution, with either different (Eq. 3) or similar (Eq. 4) scores for the C and N termini. The single AA composition exhibits a well-defined distribution for each region (Figure 1), and is a reasonable candidate for a classifier. A Support Vector Machine (SVM) based classifier was built with either linear or quadratic kernel functions. We tested classifiers learnt with different sizes of the negative learning set, and different costs for the classification error of the learning set (C) (see Methods for the precise definition). The best accuracy obtained was 0.75 with sensitivity and specificity of 0.69 and 0.82, respectively. The best sensitivity obtained was 0.81 with specificity and accuracy of 0.64 and 0.73, respectively. The best specificity obtained was 0.99 with sensitivity and accuracy of 0.3 and 0.65, respectively (Table 1).

The main two models of peptide processing in the MHC class II pathway are either that proteins are first cut in their natural form and then the resulting peptides bind MHC class II molecules (‘cut first, bind later’ model) (Sette et al., 1989; Sercarz and Maverakis, 2003), or that proteins are first partially cut. The resulting long peptides bind MHC-II molecules and then are further trimmed (‘bind first, trim later model’/‘MHC-guided processing’) (Sercarz and Maverakis, 2003). If full cleavage precedes binding, then peptides bound to MHC class II molecules are the product of fully degraded proteins. If this is indeed the case, it would be reasonable to assume similar hydrolyses activity for the N and C termini. This would lead to a similar computed cleavage probability for both termini. However, if the MHC-guided processing model is correct, then a stage of binding between denatured proteins and the MHC molecule exists prior to trimming by exopeptidases. Exopetidases are roughly divided into two basic families: Aminopeptidases and Carboxypeptidases, and it is only reasonable to assume that the trimming probabilities differ depending on the termini. While the amino acid distribution around the C and N termini were slightly different, using scores based on different cleavage properties for the N and C termini leads to no significant improvement in the precision. Thus, the observed AA distribution around cleavage sites does not explicitly support a MHC-guided processing model.

If the N and C termini do not have the same effect on the cleavage probability, then the score should contain 40 more parameters for the different contribution of the AA flanking the N and C termini, leading to the 100 parameter classifier in Eq. 3, instead of the classifier of Eq. 4. We have tested such a classifier and obtained a slightly better accuracy than the classifier with N and C termini united (Table 1). We thus do not observe a significant difference between the cleavage properties at the C and N termini.

The three scores produced in the above section highlight important aspects of the cleavage process. When looking at the score separating the C and N cleavage sites, one observes a high correlation between the two sites (r=0.42, p<0.01). Most of the difference is due to positions with a limited contribution to the cleavage, leading again to the conclusion that the cleavage process is approximately similar in C and N.

If one averages the scores of the C and N termini, one obtains again a 60 parameter score, which is very similar to the 60 parameter scores obtained with and without the in vitro cleavage (Figure 3). Since the different scores are quite similar, we averaged them and compared the contribution of each amino acid to the score. The most important contribution is from the cleavage sites, while the body of the peptide has a more limited contribution to the variability of the score, with the main limitation being that a Cysteine in the body of the peptide significantly increases the probability that the peptides should be produced and not cleaved in its center (Figure 4).

At the right of the cleavage site, there is a clear negative contribution of hydrophilic amino acids and a positive contribution of hydrophobic amino acids. This contribution is well demonstrated by correlation of cleavages introduced by individual cathepsins with the more hydrophobic regions within the HEL sequence (Pluger et al., 2002). On the left side of the cleavage, charged amino acids have a negative contribution to the cleavage, while more neutral or aliphatic amino acids increase the cleavage probability.

We further examined the contribution of amino acids at a distance of more than 1 from the cleavage site to the peptide cleavage probability. We repeated the analysis described in the previous section, and extended the score to include flanking regions up to 5 amino acids away from the cleavage site at each side (Supp. Mat., Figure S3). However, no significant improvement in the learning performance was shown (Accuracy of 0.682).

3.3. Amino Acid Pair score

In order to further improve the accuracy of the classifier, we tested scores incorporating not only the single AA distribution, but also the diAA distribution. Note that enlarging the score to include diAA increases the complexity of the problem from 20/40 parameters to 400 variables in each region. The resulting score is:

If one assumes again that the behaviors at the C and N termini are similar, Eq. 5 simplifies to:

3.4. Inclusion of in vitro cleavage experiments

A second element introduced to further improve the precision of the algorithm is endosomal proteases cleavage sites, as determined by in vitro experiments (Rawlings and Barrett, 1999). These in-vitro experiments were merged with the current data (naturally processed CD4+ T cell epitopes). These two data sources are independent and each provide a different type of information. In vitro cleavage experiments provide information about the activity of specific proteases. However, such experiments do not provide information about the relative importance of each protease in the combined cleavage process. On the other hand, naturally processes CD4+ T cell epitopes provide a comprehensive end product of the cleavage process, but they are affected by other processes, mainly the MHC class II binding, and as such may be biased. We created a score based on the average diAA frequency in the peptide body (20X20 parameters) and the flanking regions (unified C and N – 20X20 parameters), yielding an overall dimension of 800. The combination of the two data sources produced a better classifier than is obtained by each separately (Table 1). In this case again the SVM results were better than the ones obtained in the SA analysis (data not shown). The SVM algorithm was used with varying error cost factors: 10⁻⁶≤C≤10⁰, negative training data set size: 1000 and 3000 peptides. The SVM classifier was applied either with or without the in vitro cleavage data and with linear or quadratic kernels. Within the different SVM, the accuracy of linear kernel functions was lower than the quadratic kernel function. These best classification results are summarized in Table 1. The best classifiers are available at http://peptibase.cs.biu.ac.il/PepCleave_cd4/

3.5. Optimal classifier with reduced dimension

The gain obtained in the learning accuracy may be offset by the loss in the generalization resulting from the large number of parameters used. We tested whether an optimal balance between the learning and generalization by testing scores with different numbers of variables. We ranked the features by the value of $\frac{{\mu }_{\mathit{pos}}-{\mu }_{\mathit{neg}}}{\sqrt{{\sigma }_{\mathit{pos}}+{\sigma }_{\mathit{neg}}}}$ in descending order. An SVM based classifier was built with either linear or quadratic kernel functions with a gradually increasing number of features. We tested classifiers learnt with different sizes of the negative learning set, and different costs for the classification error of the learning set (C) (see Methods for a precise definition). Features were taken from the top of the ranking list created by the above criterion. The thresholds for each classifier were set to maximize accuracy. The best accuracy was achieved with 240 out 800 features and a linear classifier with in vitro data included (accuracy of 0.77 with sensitivity and specificity of 0.82 and 0.71 relatively (Table 2); threshold = −0.2023). While this score is less natural from an immunological point of view, since it ignores multiple diAA combinations, given the amount of existing observations, it is one of the best that can be produced. Using data extraction algorithms (see Methods) to further reduce dimensionality didn’t improved the accuracy (Figure S4).

3.6. Effect of MHC class II binding

The predicted difference between the cleaved peptides and the non cleaved was based on MHC class II naturally processed eluted peptides. One could claim in principle that the source of the difference is due to the MHC class II binding. This claim is not very probable since, in contrast with MHC class I, the binding core of MHC class II epitopes is limited to a small part of the peptide. Thus, many different cleaved peptides can prime the appropriate T cells. However, in order to ensure that the results of the algorithm are not due to MHC class II binding, we used the learnt classifier (linear SVM learning on 240 out of 800 features description) to evaluate the score of MHC class II binding peptides. The flanking regions of 6169 peptides were extracted from their source protein, and only peptides with a single source protein were used. The fraction of MHC class II binding peptide positives was not significantly different than the fraction of random peptides ((sensitivity of 0.3197 (Table 2)); Figure 5).

3.7. T Cell response epitopes

An interesting question emerging in the context of a classifier is the predicted score of T cell activating peptides (epitopes). In the case of MHC class I molecules, an externally introduced epitope is considered to be naturally processed, since otherwise there would be no T cell clone responding to it. Thus T cell epitopes indeed have high predicted proteasomal cleavage and TAP binding scores as one can see in the case of HIV and Herpes Viruses (Vider-Shalit et al., 2007) (Vider-Shalit et al., 2009). However, in the case of MHC class II, this correlation may break, since binding is based on the properties of the typically nine-amino acid long. The cleavage, on the other hand, depends on positions that may be far from this nine amino acid binding site. Thus, T cell response (of an externally administered peptide) and cleavage may be independent in the case of MHC class II. To test if this is indeed the case, we repeated the analysis performed for MHC class II binding peptides, using T cell epitopes, and indeed there is no difference between the fraction of T cell epitopes defined to be positive in the cleavage score and the fraction of random peptides defined to be positive (Figure 5).

To conclude, the low sensitivity of positive MHC-II binders and T cell responders indicate that the classifier is indeed based on cleavage features. To validate that this observations are not sensitive to the threshold we used (−0.2023), we compared the sensitivity of the same data sets with different thresholds (Figure S5). For any threshold we used, the sets of MHC-II binders and CTL responders had low sensitivity, similar to the random set.

3.8. Biochemical elements affecting cleavage

In order to understand the cleavage mechanism, we have looked for correlations between the score obtained in the previous section and multiple biochemical properties. We have here focused on the cleavage score that yielded the optimal accuracy (240 out of 800 variables). The propensity scales relate to single AA and the weight vector are based on diAA. We have thus averaged the weight vector over all positions on the left of the diAA, yielding a length 20 vector. A parallel analysis was then performed when averaging on the position to the right of the cleavage site. The resulting vectors all showed a high correlation to structural properties of amino acids, such as the linker propensity, which is determined by the rigidity and hydrophobicity of the amino acids (George and Heringa, 2002). The highest correlations coefficient calculated are: for the left side of cleavage site with Linker propensity from non-helical (annotated by DSSP) dataset (George and Heringa, 2002) and with Linker propensity from 1-linker dataset (George and Heringa, 2002) r=0.803 p<0.005 and r=0.763, p<0.005, for the right side of cleavage site r=−0.627, p<0.005, for the left side of diAAs in peptide body with Linker propensity from 3-linker dataset (George and Heringa, 2002) r=−0.619, p<0.005 and for the left side of diAAs in peptide body with Linker propensity from 3-linker dataset (George and Heringa, 2002) r= −0.694, p<0.005. Multiple other measures were highly correlated with the score, most of which are related to the secondary structure of the peptide.

These correlations highlight the importance of secondary structural aspects and the available surface area for the cleavage process. The structural effect is probably the most important difference between proteasomal cleavage operating on fully denatured peptides and endosomal cleavage affecting natural proteins. The proposed score was fully sequence based. We expect structure based scores to provide an even better accuracy. However, the limited precision of structure prediction algorithms and the number of features describing the structure make such an algorithm more complex, and would probably require significantly larger learning sets.

3.9. Relation between body of peptide and cleavage site

The AA that have the greatest impact on the cleavage score (highest absolute weight in the cleavage score) are Arg, Gly and Lys on the left side, and Ala, Gly and Ser on the right side. Since the actual cleavage site profile is a consequence of the peptidases participating in the cleavage process, the calculated cleavage profile site should express these peptidases substrates.

We have computed the correlation between the contributions of each amino acid to the experimental profile of the two peptidases shown to have a distinct roles in the endosomal presentation in vivo: cathepsins L and S (Villadangos and Ploegh, 2000) and the score vector. We found a significant correlation for the left side of cleavage site, but not for the right side of the cleavage site. We have no clear explanation at this stage for this difference.