This work focuses on the challenges of using regression learning algorithms in the design of highly active peptides. Let A be the set of all amino acids and be the set of all possible peptides. Throughout this paper we will assume that we have a dataset where is the amino acid sequence of the i-th peptide in S and is its bioactivity (which could be its binding affinity to some target protein, its antimicrobial activity, or some other desirable activity for peptides). The regression approach consists of learning a predictor h from the training dataset S. Let h(x) be the estimated bioactivity of x according to the predictor h. There are many ways to represent such a predictor but this work focuses on string kernel based predictors. Many learning algorithm, such as the Support Vector Regression, Ridge Regression, and Gaussian Processes, produce predictors h whose output h(x), on input x, is given by

where is the weight on the i-th training example, is a string kernel, denotes the inner product between feature vectors, and is the feature map associated to K. The weight vector is obtained by minimizing the learning algorithm objective function.

For the discovery of peptide inhibitors and peptides that could act as drug precursors, we are interested in finding the peptide maximizing the predicted bioactivity:

The complexity of this combinatorial problem will obviously depend on the string kernel used to learn h. For few kernels like the Hamming kernel, solving Equation (2) is trivial [1]. However, this kernel is not well suited for peptides [2]. On another hand, the Generic String (GS) kernel [2] is well suited for peptides. It is defined as follows:

where k controls the length of compared k-mers, encodes the physico-chemical properties of k-mers by mapping each of the k amino acids to a real valued vector containing d properties, controls the penalty incurred when the physico-chemical properties of two k-mers differ, and controls the penalty incurred when two k-mers are not sharing the same position in their respective peptides. Depending on the chosen hyper-parameters, this kernel can be specialized to eight known kernels [2], namely the Hamming kernel, the Blended Spectrum [3], the Radial Basis Function (RBF), the Oligo [4], and the Weighted degree [5]. Since the proposed approach uses the GS kernel it is also valid for all of these kernels.

It was recently shown [6] that when K is the Generic String (GS) kernel, and when we restrict peptides to be of length l, the peptide can be found in polynomial time for any predictor h in the form of Equation (1). Their approach maps the combinatorial problem to a directed acyclic graph (DAG) that is basically a de Bruijn graph with weights on the arcs. Then, they show that finding the longest (weighted) path in this graph is equivalent to finding . Finally, the longest path can be found in polynomial time by dynamic programming since the graph is acyclic.

However, with the GS kernel, given two peptides x and of the same length, the Euclidean norms of their feature vectors and can differ substantially, i.e. we can have. Let us first show, with a simple example, why this can be problematic. We will then show how to avoid this problem through kernel normalization. Such normalization does not have any computational impact in the learning phase of the regression algorithm. It does however impact substantially the prediction phase, leading to a harder combinatorial problem.

Hence, let us consider the strings “AAAAA” and “ABCDE” and the Spectrum kernel [7] (also known as the n-gram kernel) with k-mers of size two. The string “AAAAA” has a norm of, while “ABCDE” has a norm of . Hence the Spectrum kernel is sensitive to k-mers repetitions in the string. Note that this problem is shared by many other string kernels. In the case of the GS kernel, it is possible to avoid this problem by fixing the hyper-parameter to 0, which forces a constant norm, i.e. K(x, x) = c for all . However, K(x, x) will vary whenever , as it is the case for the Blended Spectrum and the Oligo kernels, for example. A consequence of a non-constant norm is that h(x) will heavily depend on the norm of . Hence, will be more likely a peptide having a feature vector with a large norm, i.e., a peptide having many repetitions. This is generally an undesired bias. It is easy to overcome this problem by normalizing kernel values, in that case the output function becomes

This optimization problem is also a pre-image problem, but written in a slightly different form. We conjecture that solving Equation (5) is NP-Hard when K is the GS kernel. Given the similarity between the problems of Equation (2) and Equation (5), the difference in their computational complexity is unexpected.

In the next section, we will present a low complexity upper bound on Equation (4) that makes it a good candidate for a branch and bound search to solve Equation (5).

Since the number of peptides grows exponentially with its length, it becomes impossible to evaluate all solutions for large peptides, we propose a branch and bound scheme to guide this search. A branch and bound algorithm starts by dividing the search space into disjoint subspaces. For example, one subspace could be all peptides ending with the string “DE”. For a maximization problem, an upper bound on the best achievable solution is computed for each of these subspaces. Then, the subspace with the highest upper bound is further divided into smaller subspaces. Finally, the search stops when a subspace can no longer be divided (a leaf is reached in the search tree), or when the upper bound value is lower than the value of an already achieved solution (i.e., an already reached leaf in the search tree). A branch and bound approach can thus avoid exploring a large part of the search space.

Algorithm 1 gives the specifics of the branch and bound algorithm applied to our case. The search algorithm alternates between a greedy phase and a branch and bound phase. The greedy phase is important to ensure that leaves of the search tree are quickly visited. This allows good but sub-optimal solutions to be returned by the algorithm if the allowed computational time expires. Whenever a node is visited, the bound is computed for all its children and they are added to the priority queue accordingly. This greedy process is repeated until a leaf is reached. Then, the node with the largest bound is visited and the greedy process starts again. At all time, the best solution found so far is kept and the search stops when the bound of the node on top of the priority queue is smaller than the value of the best solution.

Let be the set of all possible strings of length l that end with , in other words, is a set of strings having their last p characters fixed. Our goal is to have a function F that upper bounds for every . In other words:

To do so, let , where

are, respectively, a lower and an upper bound.

When K is the GS kernel with hyper-parameters and , the lower bound f can be obtained as follows:

where

and

Note that the lower bound f is not attained since it under-estimates the value of the string minimizing K(x, x).

For , note that is what [6] proposed to maximize. Their approach uses a dynamic program- ming table to compute the longest path in a graph. It is relatively easy to modify their algorithm to return the table instead of the longest path. In that way, given a string with suffix , it is possible to determine in constant time, by accessing the dynamic programming table, the prefix from maximizing . The computation of g is thus very efficient, the algorithm of [6] only needs to be done once before the branch and bound search, then can be computed in constant time for any . Finally, g is the least upper bound (or suprema) since there is always a string with . In other words, there are no tighter bound.

We followed the protocol suggested by [6] and, as a proof of concept, we used the same dataset they used: 101 cationic antimicrobial pentadecapeptides (CAMPs) from the synthetic antibiotic peptides database [8]. Peptide antibacterial activities are expressed as the logarithm of bactericidal potency. As in [6], we used kernel ridge regression as the learning algorithm. Except when stated otherwise, all hyper-parameters for the GS kernel () and the kernel ridge regression () were chosen by standard cross-validation. We learned two predictors of antimicrobial potency, one uses unnormalized kernel values (thus, the same predictor used in [6]), the other was trained using normalized kernel values. We refer to these predictors respectively as h and .

The method presented in [6] was used to identify the peptide of maximal predicted bioactivity according to h and the branch and bound was used to identify , the peptide maximizing . Both approaches found the same peptide, which is “WWKWWKRLRRLFLLV”.

Note that both methods are able to output more than one sequence. The method of [6] uses a k-longest path algorithm to obtain the k peptides of maximal bioactivity. It is also possible for the branch and bound by stopping the search only when the bound on top of the priority queue is lower than the k-th peptide found. These peptides can be used, for example, for motif generation, for multiple peptide synthesis or in combinatorial chemistry. To highlight the differences between the methods, they were used to list the top 1000 peptides and the lists were compared: 70.5% of the 1000 peptides were found by both methods. Then, the Pearson correlation coefficient (PCC) was computed between the rank of the 679 peptides present in both lists: a PCC of 0.56 was obtained. Hence, the ranking of peptides found by h and by differ significantly.

The overlap in the lists is attributed to the value of that was chosen during cross-validation for h. As explained earlier, when is 0, the unnormalized predictor h will not suffer from differences in the norm of . In the case of antimicrobial peptides, was found to be the optimal value for h. This suggests that the method of [6] offers some resistance to variation in norm when is small. Indeed, for small , all examples have about the same norm.

To further highlight the difference between the methods, we intentionally fixed to infinity and cross-validated all other parameters for h and and compared the best peptides found. Situations were would have to be set to infinity are not at all unlikely. For example, cyclic peptides have no N-terminus and no C-terminus. For that reason, there is no origin from which we can express the positions of k-mers in these peptides. The method of [6] found the peptide “FKKIFKKIFKKIFKF” using the predictor h and the branch and bound approach found the peptide “WKKIFKKIWKFRVFK” using the predictor . The peptide identified by h shares almost no similarity with peptides of the training set and is basically composed of repetition of the k-mer “FKK”. In contrast, the peptide identified by shares many substructures with the most bioactive peptides of the training set. This tends to point out that in situation where example norms vary a lot, h clearly favors examples having a large norm. However, this bias is unjustified and is not related to a biological reality.

We proposed a bound for maximizing the inference function of kernel methods that uses normalized string kernels. Moreover, the bound is also valid for solving the pre-image of a variety of string kernels. Empirical results show that the method proposed by [6] can suffer from a dominance of the norm for certain strings, a problem which is present with many string kernels. In these situations, the proposed method was shown to overcome this problem. Tighter bounds should take advantage of the proposed framework and allow the discovery of novel peptide inhibitors. Finally, applications for drugs based on cyclic peptides and other structured output applications are expected.

The authors would like to thank Prudencio Tossu for his help in the experimentations. Computations were performed on the Colosse supercomputer at Universit´e Laval (resource allocation project: avt-710-aa), under the auspices of Calcul Qu´ebec and Compute Canada. AR is recipient of a Master’s Scholarship from the Fonds de recherche du Qu´ebec - Nature et technologies (FRQNT). This work was supported by the FRQNT (FL & MM; 2013-PR-166708).

[1] S´ebastien Gigu`ere, Franc¸ois Laviolette, Mario Marchand, and Khadidja Sylla. Risk bounds and learning algorithms for the regression approach to structured output prediction. In International Conference on Machine Learning (ICML), 2013.

[2] S´ebastien Gigu`ere, Mario Marchand, Franc¸ois Laviolette, Alexandre Drouin, and Jacques Corbeil. Learning a peptide-protein binding affinity predictor with kernel ridge regression. BMC Bioinformatics, 14, 2013.

[3] John Shawe-Taylor and Nello Cristianini. Kernel methods for pattern analysis. Cambridge university press, 2004.

[4] P. Meinicke, M. Tech, B. Morgenstern, and R. Merkl. Oligo kernels for datamining on biological sequences: A case study on prokaryotic translation initiation sites. BMC Bioinformatics, 5, 2004.

[5] Gunnar R¨atsch and S¨oren Sonnenburg. Accurate Splice Site Detection for Caenorhabditis elegans. In B and J. P. Vert, editors, Kernel Methods in Computational Biology, pages 277–298. MIT Press, 2004.

[6] S´ebastien Gigu`ere, Franc¸ois Laviolette, Mario Marchand, Denise Tremblay, Sylvain Moineau, ´Eric Biron, and Jacques Corbeil. Machine learning assisted design of highly active peptides for drug discovery. Under review, Submitted to MLCB, 2014.

[7] Christina S Leslie, Eleazar Eskin, and William Stafford Noble. The spectrum kernel: A string kernel for svm protein classification. In Pacific symposium on biocomputing, volume 7, pages 566–575. World Scientific, 2002.

[8] David Wade and Jukka Englund. Synthetic antibiotic peptides database. Protein and peptide letters, 9(1):53–57, 2002.