On the basis of Bayesian probabilistic inference, Gaussian process (GP) is a powerful machine learning method for nonlinear classification and regression, but has only very limited applications in the new areas of computational vaccinology and immunoinformatics. In the current work, we present a paradigmatic study of using GP regression technique to quantitatively model and predict the binding affinities of over 7000 immunodominant peptide epitopes to six types of human major histocompatibility complex (MHC) proteins. In this procedure, the sequence patterns of diverse peptides are characterized quantitatively and the resulting variables are then correlated with the experimentally measured affinities between different MHC and their peptide ligands, by using a linearity- and nonlinearity-hybrid GP approach. We also make systematical comparisons between the GP and two sophisticated modeling methods as partial least square (PLS) regression and support vector machine (SVM) with respect to their fitting ability, predictive power and generalization capability. The results suggest that GP could be a new and effective tool for the modeling and prediction of MHC-peptide interactions and would be promising in the field of computer-aided vaccine design (CAVD).