EPMLR: sequence-based linear B-cell epitope prediction method using multiple linear regression

Sliding window size selection

To evaluate the effect of sliding window size n on the prediction performance, we conducted modelling trials on BEOD dataset using different window sizes from 5 to 19, representing the range in which peptides can be synthesized relatively easily for immune experiments. As shown in the Figure 1, the F-measure value of 10-fold cross-validation test achieved its highest value when the window size n was 15. Moreover, at 15 point, the F-measures obtained by the 10-fold cross-validation test and the self-consistency test are very close to each other, which further validates the reliability of the performance using sliding window size of 15. It is generally accepted that the closer the F-measures obtained by the cross-validation and self-consistency tests are, the more reliable the performance of the cross-validation test is. Therefore, in this work, 15 was set as the default window size.

Prediction performance

We performed 300 experiments on 10-fold cross-validation utilizing 300 sub-datasets that are the same in the positive datasets but different in the negative datasets. For each trial, the positive dataset of 4405 epitopes are exactly same with BEOD’s 4405 epitopes while the negative dataset of 4405 non-epitopes are randomly selected from BEOD’s 8467 non-epitopes. The ROC plots for the best and worst performances among the 300 trials are shown in Figure 2. The performances of all 300 trials are summarized in Table 1. As shown in Table 1 and Figure 2, the variance of the 300 results is large, with Sn ranging from 83.5% to 81.7%, P from 77.6% to 55.7%, F-measure from 0.805 to 0.663, and AUC from 0.893 to 0.673. These large discrepancies corroborate our speculation of the noise of non-epitopes even if they are experimentally verified and support our means of randomly constructing many negative sub-datasets and reporting the average result instead of the best result. In conclusion, our sequence-based linear B-cell epitope prediction method achieved an average Sn of 81.8 ± 0.8% (95% CI), P of 64.1 ± 0.2% (95% CI), F-measure of 0.719 ± 0.08 (95% CI), and AUC of 0.728 using 10-fold cross-validation.

Comparison with Other Prediction Methods

We compared our EPMLR method with the methods of ABCpred [10], AAP [11] and BCPred [13] through applying their web servers to the BEOD dataset. The ROC plots for performances of ABCpred, AAP, BCPred and EPMLR are shown in Figure 3. The AUC values for ABCpred, AAP, BCPred and EPMLR are 0.547, 0.582, 0.615 and 0.728, respectively. It is clear from the ROC plots that EPMLR produced better performance in comparison with ABCpred, AAP and BCPred.

Next, we compared our method with SVMTriP method which is a recently published large dataset based method [21]. We performed a 5-fold cross-validation on the SVMTriP dataset. Our method obtained Sn of 80.56% and P of 54.9% which is similar to the performance of SVMTriP method (Sn of 80.1%, P of 55.2%) using 5-fold cross-validation. Our method observed similar Sn (81.8% vs. 80.56%) but a decreased P (64.1% vs. 54.9%) on the BEOD dataset and SVMTriP dataset. The decreased P value could be resulted from the fact that the negative non-epitope dataset of the SVMTriP dataset was from the remaining segments which have not been marked as epitopes in the corresponding antigen sequences.

Similarly, we compared with LBtope method which is the most recently published large dataset developed method [25]. We applied our method to the Lbtope_Fixed_non_redundant dataset (LFNR) whose epitopes and non-epitopes were all experimentally verified. Using the same experimental procedure of LBtope, on the LFNR dataset, our method obtained an AUC of 0.62, which is comparable to the AUCs (0.57 ~ 0.69) obtained by LBtope method by training using 5-fold cross-validation on 90% of the data and testing on the remaining 10% of the data with various features.

Datasets

In this work, we used BEOracle dataset because it is a large dataset and both epitopes and non-epitopes of BEOracle dataset were experimentally verified. Through combining entries from IEDB, BCIPEP and AntiJen databases, Wang and his colleagues constructed the BEOracle dataset [20]. They extended these epitope sequences equally on both sides to get epitopes of a final length of 100 amino acids using the Uniprot identifiers associated with them. Further, we trimmed BEOracle dataset (100-mer) from both ends equally to extract the core 20-mer peptides. Finally, we obtained 4,405 epitopes and 8,467 non-epitopes and we called this unbalanced BEOracle-Derived dataset (BEOD).

To alleviate the problem caused by the noise of negative dataset, we then constructed 300 sub-datasets which were the same in the positive dataset but different in the negative dataset. Each of the 300 sub-datasets contained the whole 4,405 epitopes of BEOD and an equal number of 4,405 non-epitopes randomly selected from the 8,467 BEOD non-epitopes. These 300 sub-datasets were used to perform the 300 experiments using the same algorithm.

The SVMTriP dataset, which was introduced by Yao B et al. [21], consists of 4925 epitopes and 4925 non-epitopes. Originally, total of 65,456 B-cell linear epitopes were downloaded from IEDB (version June 11th, 2012) and the identical epitopes and those possibly related to T-cell are removed. Next, truncation and extension technique was applied to get fixed length pattern. Finally, 4925 non-redundant epitope sequences were obtained after > = 30% similarity process by BLAST [28]. For the negative dataset, the same number of equal-length sub-sequences were extracted from the non-epitopic segments in the corresponding antigen sequences.

The Lbtope_Fixed_non_redundant dataset (LFNR), which was introduced by Singh H et al. [25], consists of 7824 B-cell epitopes and 7853 non-epitopes. Originally, total of experimentally validated 49694 B-cell epitopes and 50324 non B-cell epitopes were obtained from the IEDB in Jan 2012. After truncation and extension, sequences with fixed length were created. Then identical epitopes and common patterns in both types of patterns were removed. Finally, after 80% non-redundant process by CD-HIT [29], 7824 B-cell epitopes and 7853 non-epitopes were kept. This non-redundant and fixed length dataset was named Lbtope_Fixed_non_redundant.

Algorithm

In this study, we constructed an epitope prediction model based on primary sequence information. The modeling trial was performed as follows.

Here, subscripts j and k denote position j and k in the window. R _j is a 19-D vector with the component for the residue at position j as 1 and the others as 0. α(1, 2 … 19|ω) is the coefficient vector for 19 amino acids (with one omitted). ${\displaystyle \sum _{j=1}^n\alpha \left(1,2\dots 19|\omega \right)R_j}$ represent the features of occurrence of amino acid type from the first position to the last position for an n-mer window sequence. B _j and B _k are the normalized hydrophilicity values of residues at positions j and k, while β _j,k(ω) is the coefficient combining the residue pair. ${\displaystyle \sum _{j=1}^{n-1}}{\displaystyle \sum _{k=j+1}^n}{\beta }_{j,k}\left(\omega \right)B_j{B}_k$ represent the features of autocorrelation of the hydrophobicity index of residue pair (residue R _j at position j and residue R _k at position k) for an n-mer window sequence. Similarly, S _j and S _k are the normalized side chain mass values of residues at positions j and k, while γ _j,k(ω) is the coefficient combining the residue pair. ${\displaystyle \sum _{j=1}^{n-1}}{\displaystyle \sum _{k=j+1}^n}{\gamma }_{j,k}\left(\omega \right)S_j{S}_k$ represent the features of autocorrelation of the side chain mass of residue pair for an n-mer window sequence. V(R _j R _k ) is a 500-D vector whose components refer to 500 most important position specific residue pairs R _j R _k , while δ _j,k(ω) is the coefficient combining the residue pair. In model training, we compared the 500 R _j R _k with all R _j R _k (n × (n − 1)/2 in total) existed in a window, the value of a component of V(R _j R _k ) is set as 1 if the R _j R _k to which the component referred exists in the window, otherwise as 0. ${\displaystyle \sum _{j=1}^{n-1}}{\displaystyle \sum _{k=j+1}^n}{\delta }_{j,k}\left(\omega \right)V\left(R_j{R}_k\right)$ represent the feature of occurrence of selected residue pairs in an n-mer window sequence.

where f _t (R _j R _k ) represents the occurrence frequency of a R _j R _k derived from the epitope (t = 1) and non-epitope (t = 0) in the training dataset, respectively. P _t represents the naturally occurring probability of a R _j R _k based on the relative sizes of the epitope and non-epitope datasets in the training dataset (for example, P ₁ = P ₂ = 0.5 if the size of the epitope dataset is equal to the size of the non-epitope dataset). All D(R _j R _k ) values were ranked by the descending orders. Finally, 500 R _j R _k with the largest value of D(R _j R _k ) were selected. Here, we selected 500 components because the curve of all D(R _j R _k ) values by descending order shows as exponential decay and the point of inflection is about 500 ( Additional file 1).

On the training dataset, all the fitting coefficients in Equation (1) were determined by the MLR method [30]. Once the coefficient matrix is obtained, we adopted the same sliding window procedure with the 20-mer peptides on the testing dataset. Each of the n-sized window ω _i of the 20-mer peptide was predicted to be an epitope or not with an epitope propensity score Q(ω _i ). For any 20-mer peptide, there are 21 − n windows and the epitope propensity score of the 20-mer peptide was calculated by taking the average of all 21 − n Q(ω _i ) scores. In this representation, every 20-mer peptide in the testing dataset is scored for its propensity to be an epitope or a non-epitope.

Performance Measures

In 10-fold cross-validation test, the original dataset is randomly partitioned into 10 equal size subsets. Of the 10 subsets, a single subset is retained as the validation data for testing the model, and the remaining 10-1 subsets are used as training data. The cross-validation process is then repeated 10 times, with each of the 10 subsets used exactly once as the validation data. The 10 results can then be averaged to produce a single estimation.

where TP, TN, FP, and FN represent the number of true positive, true negative, false positive, and false negative cases, respectively.