CRF-based models of protein surfaces improve protein-protein interaction site predictions

Residue-wise score-based prediction methods

Let x _r be the data relevant for a residue r in a given protein chain. These methods then employ a function f(x _r ,λ), where λ are some coefficients which have been learned through the training. The value of f(x _r ,λ) then determines, whether r is rated as an interface or not. The linear regression method [4, 5], the scoring function method [6–11], the neural network method [12–17], and the support vector machine method [18–25] are of this kind.

Probabilistic methods

Let X be the data relevant for a protein chain, where these data are assumed to stem from a random source thus obeying a random distribution. X, which alternatively is called the observation, typically includes the structure. The label sequence of the residues Y that classifies each individual residue either as interface or as non-interface is assumed to be random, too. Typically, probabilistic methods use the conditional probability distribution (Y | X) to determine a classification y* of the residues of maximal posterior probability (y* | x). Naive Bayesian methods [26], Bayesian network methods [27], hidden Markov models (HMMs) [26], and linear-chain Conditional Random Fields taking the backbone as underlying graphical structure [28] fall in this category. Using posterior decoding on the basis of the forward-backward algorithm, both HMMs and CRFs are residue-wise score-based prediction methods, where the binary decision is made by thresholding the posterior probabilities of classifying the residues as interface.

Notations

We use Latin uppercase letters when referring to random aspects of the objects denoted by them. In contrast, lowercase letters denote arbitrarily chosen but fixed objects. In this context boldface letters indicate vectors, the corresponding non-boldface letters their coefficients.

The vast majority of methods use the 3D structure of the target protein chain in form of a PDB file as input [4–13, 15, 17–21, 23–25]. However, a few methods are not requiring a 3D structure and rather use sequences only [14, 16, 22]. We here consider the problem with a given 3D structure of the target protein chain. Sequence-based input may include a multiple sequence alignment of related proteins from which, for example, sequence conservation can be inferred. When the 3D structure of an unbound binding partner is also available, protein-protein docking methods can be applied. This has also been exploited to provide feedback from docking to the more specific problem of interface prediction [29]. We here consider the case where the binding partner’s 3D structure is not given. Nor requires the presented method the sequence of the binding partner. Albeit, we tested on homodimers only as we here rather focus on our new method rather than on features or types of proteins. The protein features used for interface prediction in the literature are reviewed in the Methods section as far as we make use of them in this article.

Most of the current studies for predicting interaction sites of proteins that use a probabilistic method are restricted by treating the residues of the proteins as independent vertices. Li et al. have taken the backbone neighborhood into account thus modeling the protein as a sequence [28] using what can be called a line CRF or linear-chain CRF. The features they define on the label pair of two backbone neighbors have the effect of smoothing the predicted labels along the protein sequence. Decisive is, however, that they were the first who used conditional random fields (CRFs) for interface prediction. CRFs in turn have come into use for solving sequence labeling problems due to Lafferty et al.[30]. See [31] for an overview. From the mathematical point of view they take advantage of the fact that they model the conditional probability (Y | X) rather than the joint probability (Y, X). Recently there has been an explosion of interest in conditional random fields (CRFs) with successful applications. It has been shown that CRFs have the abilities for solving sequence labeling problems like part-of-speech tagging (POST) [32] and natural language processing [33]. Furthermore in the web extraction problem, in which the web-sites are modeled as two dimensional grid graphs, CRFs perform well [34]. One of their outstanding benefits over many other statistical models is that a CRF can easily describe the dependencies of observations.

As proteins are folded into three dimensional structures, spatial relationships create dependencies between residues. For example, we find on the test data described below that the correlation coefficient between spatial neighbors that are not also sequence neighbors (distance ≤3.5 Å) is 0.45. This is only slightly lower than the correlation coefficient between residues that are sequence neighbors (0.49). As there are more than three times as many spatial pairs of neighbors than sequence neighbors at this threshold it is reasonable from a modeling standpoint to use a model that respects all dependencies induced by spatial proximity, not only the dependencies induced by proximity along the backbone.

There are many papers using spatial neighborhood information of residues to predict-protein interaction sites (see e.g. [2, 13, 21, 28]). However, the spatial information of proteins was only integrated into the feature functions, but not represented in the model. For probabilistic models, the difference between the two ways to integrate spacial information is that in previous models the label of the i-th residue Y _i is conditional independent from the labels of other residues given data X and – in the case of linear CRFs or HMMs – given the labels of Y_i−1 and Y_i+1. Even when neighborhood information is only used for spatial smoothing of the labels, the intuitive advantage over, say, an SVM classifier that uses spatial neighborhood in the features but classifies each residue independently, is that not-patch-like candidate labelings are explicitly punished. In contrast, such an independent classifier-approach may have a tendency to predict individual interface residues ‘sprinkled’ around the protein surface [28].

For this reason, a general CRF seems to be more suitable for the task. However, inference for general CRFs is intractable. In this paper, pairwise conditional random fields (pCRFs) are utilized. Specializing general CRFs, only node cliques and edge cliques are taken into consideration in pCRFs. A pCRF retains most spatial information of proteins, can be specified with the same number of parameter as a line CRF and approximate inference remains feasible with the generalization of the Viterbi algorithm introduced here. Taking pattern from piecewise training methods [35], we disentangled the labels of nodes and edges to train the model.

In order to take advantage of a residue-wise score-based predictor, we model the protein surface by means of a pCRF, where the observation is solely a sequence of surface residue scores between 0 and 1 output by the predictor. We then utilize a generalized Viterbi algorithm and piecewise training. The resulting tool tries to enhance the predictor chosen on the surface of the protein under study. It is the aim of this paper to demonstrate effectiveness of this approach provided that the interface residue scores and the non-interface residue scores are appropriately distributed.

Using conditional random fields to model protein surfaces

For every protein under study that has n surface residues, a pair of random vectors (X, Y) is considered. The vector X is the observation that represents the knowledge about this protein that is utilized in the prediction, e.g. the 3D structure of the target protein and a multiple sequence alignment together with homologs.

The vector Y is a random sequence of length n over the alphabet {I,N} that labels the index set {1,2,…,n}, which in turn is called the set of positions (of the surface residues). The label I represents interface residues, whereas the label N represents non-interface residues. {I,N} ⁿ is the set of all label sequences of length n over {I,N}. We will also call them assignments as the term ‘label sequence’ may lead to confusion when applied below to subsets of {1,2,…,n} that are not contiguous sequences.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be the neighborhood graph, where $\mathcal{V}=\{1,2,\dots ,n\}$ is the set of positions, is the set of edges that typically results from an atom-distance-based neighborhood definition for positions. We assume for convenience in notation that has no isolated nodes. Cases with isolated nodes could trivially be reduced to cases without isolated nodes. Let be the set of ’s cliques, which we refer to as node cliques. For a node clique $c\in \mathcal{C}$ and an assignment y we denote by y _c the restriction of y to the positions belonging to the node clique c. For c={i} and c={i,j} we write y _i and $\left(y_i,y_j\right)$ rather than y_{i} and y_{i,j}.

The preceding notation is also used in the slightly more general case of partial label assignments to arbitrarily chosen subsets of the set of positions . Formally, let ${\mathbf{y}}_{\mathcal{S}}$ denote $\left\{(i,y_i)|i\in \mathcal{S},y_i\in \{\text{I},\text{N}\}\right\}$. Given two partial assignments ${\mathbf{y}}_{{\mathcal{S}}_1}$ and ${\mathbf{y}}_{{\mathcal{S}}_2}$ are identical on ${\mathcal{S}}_1\cap {\mathcal{S}}_2$, the union ${\mathbf{y}}_{{\mathcal{S}}_1}\cup {\mathbf{y}}_{{\mathcal{S}}_2}$ is well-defined.

Let us call $ln\left(Z\right(\mathbf{x}\left)\mathbb{P}\right(\mathbf{y}\mid \mathbf{x}\left)\right)$ the score of the label sequence y given the observation x.

A CRF is called a pairwise CRF (pCRF) if Φ _c ≡0, for all node cliques c larger than two. The remaining features Φ _i and Φ_i,j are referred to as node features and edge features, respectively. Thus, every position $i\in \mathcal{V}$ and every edge $(i,j)\in \mathcal{E}$ is represented by the pair (Φ _i (N,x),Φ _i (I,x)) and by the quadruplet (Φ_{i,j} (N,N,x),Φ_{i,j}(I,N,x),Φ_{i,j}(N,I,x),Φ_{i,j}(I,I,x)).

need to be calculated in a training phase.

In the most general sense, protein characteristics are real-valued evaluations of positions and pairs of adjacent positions (edges of the neighborhood graph), respectively, that are correlated with our position labeling problem. We use a standard step function technique to obtain base features from protein characteristics, rather than taking the raw values of the characteristics. To make our paper self-contained, let us describe this technique for short.

A protein characteristic depends on the observation and either a node or an edge. Each protein characteristic, such as e.g. the relative solvent-accessible surface area of a residue, is transformed into several binary features by binning, i.e. we distinguish only a few different cases rather than the whole range of the characteristic. Assuming the common case of real-valued characteristics, the bins are a partition of the reals into intervals. The use of this discretization allows to approximate any shape of dependency of the labels on the characteristics, rather than assuming a fixed shape such as linear or logarithmic.

where y=N,I, and ι=0,1,…,γ−1 and s₀:=−∞,s _γ :=∞.

In both cases we set γ=5.

Devising a generalized Viterbi algorithm for pCRFs

The problem of finding a most probable label sequence y^∗ given an observation x is NP-hard for general pCRFs [31]. In this subsection we present a heuristic that approximately solves this problem.

To this end, we first devise an algorithm, which we call generalized Viterbi algorithm. It computes an optimal label sequence, where the posterior probability of y^∗ given x is maximized. Unfortunately, its run-time is in too many cases not acceptable. That is why we transform it in a second step into a feasible, time-bounded approximation algorithm.

The generalized Viterbi algorithm

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be the neighborhood graph underlying the protein under study. For any assignment (label sequence) y and any subset ${\mathcal{V}}^{\prime }$ of , let ${\mathbf{y}}_{{\mathcal{V}}^{\prime }}$ denote the partial assignment of y with respect to ${\mathcal{V}}^{\prime }$. (This is in line with the notation y _c (c a position clique) introduced earlier in this study).

If ${\mathcal{V}}_1,{\mathcal{V}}_2,\dots ,{\mathcal{V}}_r$ are pairwise disjoint position sets, the assignment for ${\mathcal{V}}_1\cup {\mathcal{V}}_2\cup \dots \cup {\mathcal{V}}_r$ canonically resulting from assignments ${\mathbf{y}}_{{\mathcal{V}}_1},{\mathbf{y}}_{{\mathcal{V}}_2},\dots ,{\mathbf{y}}_{{\mathcal{V}}_r}$ is denoted by ${\mathbf{y}}_{{\mathcal{V}}_1}\cup {\mathbf{y}}_{{\mathcal{V}}_2}\cup \dots \cup {\mathbf{y}}_{{\mathcal{V}}_r}$. For ${\mathcal{V}}^{\prime }\subset \mathcal{V}$, the score $\mathrm{s}{\mathcal{V}}^{\prime }({\mathbf{y}}_{{\mathcal{V}}^{\prime }}\mid \mathbf{x})$ is defined by

$\begin{array}{l}{\mathrm{s}}_{{\mathcal{V}}^{\prime }}({\mathbf{y}}_{{\mathcal{V}}^{\prime }}\mid \mathbf{x}):=\sum _{i\in {\mathcal{V}}^{\prime }}{\Phi }_i\left(y_i,\mathbf{x}\right)+\sum _{\begin{array}{c}i,j\in {\mathcal{V}}^{\prime }\\ \{i,j\}\in \mathcal{E}\end{array}}{\Phi }_{i,j}\left(y_i,y_j,\mathbf{x}\right).\end{array}$

Then the problem of determining a most probable label sequence y^∗ given an observation x can be reformulated as

$\begin{array}{l}{\mathbf{y}}^{\ast }=\underset{\mathbf{y}}{\text{argmax}}{\mathrm{s}}_{\mathcal{V}}(\mathbf{y}\mid \mathbf{x}).\end{array}$

This is the case, because it suffices to consider the score.

To put this into practice, we devised an algorithm we call generalized Viterbi. On the one hand, it is analogous to the classical Viterbi algorithm. On the other hand, there is a major difference. In our case there is no canonical order in which the positions of are traversed. Having explained our algorithm for any order, we show how to calculate a fairly effective one. In what follows, we assume that the positions not yet touched are held in a dynamic queue. Those positions having already left the queue form the history set$\mathcal{\mathscr{H}}\subseteq \mathcal{V}$.

Assume that the subgraph of induced by has connected components ${\mathcal{\mathscr{H}}}_1,{\mathcal{\mathscr{H}}}_2$, …, ${\mathcal{\mathscr{H}}}_m$. For μ=1,2,…,m, let ${\mathcal{ℬ}}_{\mu }\subseteq {\mathcal{\mathscr{H}}}_{\mu }$ be the so-called boundary component associated with ${\mathcal{\mathscr{H}}}_{\mu }$ defined by ${\mathcal{ℬ}}_{\mu }:=\left\{i\in {\mathcal{\mathscr{H}}}_{\mu }\left|\exists j\notin \mathcal{\mathscr{H}},\{i,j\}\in \mathcal{E}\right.\right\}$. The complement ${\mathcal{\mathscr{H}}}_{\mu }\setminus {\mathcal{ℬ}}_{\mu }$ is the interior of the μ-th history component. See Figure 1 for an example.

For assignments ${\mathbf{y}}_{{\mathcal{ℬ}}_1},{\mathbf{y}}_{{\mathcal{ℬ}}_2},\dots ,{\mathbf{y}}_{{\mathcal{ℬ}}_m}$ of the boundary components ${\mathcal{ℬ}}_1,{\mathcal{ℬ}}_2,\dots ,{\mathcal{ℬ}}_m$, the Viterbi variables ${\text{vit}}_{{\mathcal{\mathscr{H}}}_1}\left({\mathbf{y}}_{{\mathcal{ℬ}}_1}\right),{\text{vit}}_{{\mathcal{\mathscr{H}}}_2}\left({\mathbf{y}}_{{\mathcal{ℬ}}_2}\right),\dots ,{\text{vit}}_{{\mathcal{\mathscr{H}}}_m}\left({\mathbf{y}}_{{\mathcal{ℬ}}_m}\right)$ are defined as

$\begin{array}{ll}{\text{vit}}_{{\mathcal{\mathscr{H}}}_{\mu }}\left({\mathbf{y}}_{{\mathcal{ℬ}}_{\mu }}\right)& :=\underset{{\mathbf{y}}_{{\mathcal{\mathscr{H}}}_{\mu }\setminus {\mathcal{ℬ}}_{\mu }}}{max}{\mathrm{s}}_{{\mathcal{\mathscr{H}}}_{\mu }}\left({\mathbf{y}}_{{\mathcal{\mathscr{H}}}_{\mu }\setminus {\mathcal{ℬ}}_{\mu }}\cup {\mathbf{y}}_{{\mathcal{ℬ}}_{\mu }}\mid \mathbf{x}\right)\\ \times (\mu =1,2,\dots ,m).\end{array}$

(6)

The Viterbi variables can be represented as a set of tables, one table of size $2^{\left|{\mathcal{ℬ}}_{\mu }\right|}$ for each boundary component ${\mathcal{ℬ}}_{\mu }$. In the case where a boundary component is empty the table reduces to a single number.

At any stage, the algorithm stores the connected components ${\mathcal{\mathscr{H}}}_1,{\mathcal{\mathscr{H}}}_2$, …, ${\mathcal{\mathscr{H}}}_m$ of the current history set , the corresponding boundary components ${\mathcal{ℬ}}_1,{\mathcal{ℬ}}_2$, …, ${\mathcal{ℬ}}_m$, and Viterbi variable values ${\text{vit}}_{{\mathcal{\mathscr{H}}}_1}\left({\mathbf{y}}_{{\mathcal{ℬ}}_1}\right),{\text{vit}}_{{\mathcal{\mathscr{H}}}_2}\left({\mathbf{y}}_{{\mathcal{ℬ}}_2}\right),\dots ,{\text{vit}}_{{\mathcal{\mathscr{H}}}_m}\left({\mathbf{y}}_{{\mathcal{ℬ}}_m}\right),$ where ${\mathbf{y}}_{{\mathcal{ℬ}}_1},{\mathbf{y}}_{{\mathcal{ℬ}}_2},\dots ,{\mathbf{y}}_{{\mathcal{ℬ}}_m}$ range over all possible assignments of corresponding boundary component. We store for every assignment on the boundary, a maximizing interior assignment. This assignment is the argmax of (6) but is determined with the dynamic programming recursions defined below. Let us call these data the current state of the algorithm. It mainly consists of record sets indexed by the boundary labelings.

At the very beginning the queue contains all positions, the history set and the corresponding boundary component are empty. As long as the position queue is not empty, the top element v is extracted and the state is updated as follows.

Adjoining v to the history set , there are two cases to distinguish. Either position v is not adjacent to any other position of any old boundary component (see Figure 2) or adjoining position v to results in adding it to some connected component of the old history set or even merging together two or more of them (see Figure 3).

In the first case we simply have to take over all the old connected components, boundary sets and Viterbi variables. Moreover, we perform the instructions

$\begin{array}{ll}{\mathcal{\mathscr{H}}}_{m+1}& \leftarrow {\mathcal{ℬ}}_{m+1}\leftarrow \left\{v\right\},{\text{vit}}_{{\mathcal{\mathscr{H}}}_{m+1}}\left(\text{N}\right)\leftarrow {\mathrm{s}}_{{\mathcal{\mathscr{H}}}_{m+1}}(\text{N}\mid \mathbf{x}),\\ {\text{vit}}_{{\mathcal{\mathscr{H}}}_{m+1}}\left(\text{I}\right)\leftarrow {\mathrm{s}}_{{\mathcal{\mathscr{H}}}_{m+1}}(\text{I}\mid \mathbf{x}).\end{array}$

In the second case position v is adjacent to some boundary components, say ${\mathcal{ℬ}}_{m^{\prime }},{\mathcal{ℬ}}_{m^{\prime }+1},\dots ,{\mathcal{ℬ}}_m$. Then the old history components ${\mathcal{\mathscr{H}}}_{m^{\prime }},{\mathcal{\mathscr{H}}}_{m^{\prime }+1},\dots ,{\mathcal{\mathscr{H}}}_m$ and the current position v are merged together:

${\mathcal{\mathscr{H}}}_{\mathit{\text{tmp}}}\leftarrow {\mathcal{\mathscr{H}}}_{m^{\prime }}\cup {\mathcal{\mathscr{H}}}_{m^{\prime }+1}\cup \dots \cup {\mathcal{\mathscr{H}}}_m\cup \left\{v\right\}.$

The other history set components and corresponding Viterbi variables are not affected.

For μ=m^′,m^′+1,…,m, let ${\mathcal{R}}_{\mu }\subseteq {\mathcal{ℬ}}_{\mu }$ be the set of all positions out of ${\mathcal{ℬ}}_{\mu }$ that are no longer boundary nodes after having adjoined v to the history set. The nodes in ${\mathcal{R}}_{\mu }$ are removed from the boundary ${\mathcal{ℬ}}_{\mu }$ after the iteration. Let ${\widetilde{\mathcal{ℬ}}}_{\mu }$ be the complement of ${\mathcal{R}}_{\mu }$ in ${\mathcal{ℬ}}_{\mu }$. By inspecting the edges incident to the current position v, all these sets can be computed in linear time.

The new boundary set ${\mathcal{ℬ}}_{\mathit{\text{tmp}}}$ is then either ${\widetilde{\mathcal{ℬ}}}_{m^{\prime }}\cup {\widetilde{\mathcal{ℬ}}}_{m^{\prime }+1}\cup \dots \cup {\widetilde{\mathcal{ℬ}}}_m$ or ${\widetilde{\mathcal{ℬ}}}_{m^{\prime }}\cup {\widetilde{\mathcal{ℬ}}}_{m^{\prime }+1}\cup \dots \cup {\widetilde{\mathcal{ℬ}}}_m\cup \left\{v\right\}$, where it can be checked in linear time whether or not v is a new boundary position.

We are now in a position to calculate the new Viterbi variables ${\text{vit}}_{{\mathcal{\mathscr{H}}}_{\mathit{\text{tmp}}}}\left({\mathbf{y}}_{{\mathcal{ℬ}}_{\mathit{\text{tmp}}}}\right)$, where ${\mathbf{y}}_{{\mathcal{ℬ}}_{\mathit{\text{tmp}}}}$ ranges over all assignments of the new boundary set ${\mathcal{ℬ}}_{\mathit{\text{tmp}}}$.

If $v\notin {\mathcal{ℬ}}_{\mathit{\text{tmp}}}$ then

$\begin{array}{ll}{\text{vit}}_{{\mathcal{\mathscr{H}}}_{\mathit{\text{tmp}}}}\left({\mathbf{y}}_{{\mathcal{ℬ}}_{\mathit{\text{tmp}}}}\right)& \leftarrow \underset{y_v}{max}\left\{{\Phi }_v\left(y_v,\mathbf{x}\right)+\sum _{\mu =m^{\prime }}^m\underset{{\mathbf{y}}_{{\mathcal{R}}_{\mu }}}{max}\right.\\ \left.\times \left({\text{vit}}_{{\mathcal{\mathscr{H}}}_{\mu }}{\mathbf{y}}_{{\widetilde{\mathcal{ℬ}}}_{\mu }}\cup {\mathbf{y}}_{{\mathcal{R}}_{\mu }}+\sum _{\begin{array}{c}w\in {\mathcal{ℬ}}_{\mu }\\ \{v,w\}\in \mathcal{E}\end{array}}{\Phi }_{v,w}\left(y_v,y_w,\mathbf{x}\right)\right)\right\}.\end{array}$

Here, any assignment of a node set is assumed to implicitly define assignments for any subset thereof. Figure 4 illustrates this case of the recursion step. If, however, $v\in {\mathcal{ℬ}}_{\mathit{\text{tmp}}}$, then

$\begin{array}{ll}{\text{vit}}_{{\mathcal{\mathscr{H}}}_{\mathit{\text{tmp}}}}\left({\mathbf{y}}_{{\mathcal{ℬ}}_{\mathit{\text{tmp}}}}\right)\leftarrow {\Phi }_v\left(y_v,\mathbf{x}\right)& +\sum _{\begin{array}{c}w\in {\widetilde{\mathcal{ℬ}}}_{{\mu }^{\prime }}\cup \cdots \cup {\widetilde{\mathcal{ℬ}}}_{\mu }\\ (v,w)\in \mathcal{E}\end{array}}{\Phi }_{v,w}\left(y_v,y_w,\mathbf{x}\right)\\ +\sum _{\mu =m^{\prime }}^m\left\{\underset{{\mathbf{y}}_{{\mathcal{R}}_{\mu }}}{max}\left({\text{vit}}_{{\mathcal{\mathscr{H}}}_{\mu }}\left({\mathbf{y}}_{{\widetilde{\mathcal{ℬ}}}_{\mu }}\cup {\mathbf{y}}_{{\mathcal{R}}_{\mu }}\right)\right.\right.\\ +\left.\left.\sum _{\begin{array}{c}w\in {\mathcal{R}}_{\mu }\\ (v,w)\in \mathcal{E}\end{array}}{\Phi }_{v,w}\left(y_v,y_w,\mathbf{x}\right)\right)\right\}.\end{array}$

Finally, the interior labeling is stored, where the maximum is attained. The algorithm terminates after the last node v from has been processed. In the typical case, where the graph is connected, at termination $m=1,{\mathcal{\mathscr{H}}}_1=\mathcal{V},{\mathcal{ℬ}}_1=\varnothing$.

The running time of the algorithm is $\mathcal{O}\left(n{2}^b\right)$, where b is the size of the largest boundary set and n is the number of surface residues. We call this algorithm generalized Viterbi algorithm as for the case of a graph that is a linear chain 1−2−3−⋯−n of nodes using the node order 1,2,…,n the Viterbi variables we define are the same as in the standard Viterbi algorithm for HMMs. In the case of a graph that is a tree, this algorithm specializes to the Fitch algorithm or an argmax-version of Felsenstein’s pruning algorithm when a leaf-to-root node order is chosen after rooting the tree at an arbitrary node. In both special cases the boundary sets always have size at most 1. The tree example also motivate the use of several history sets at the same time: using a single history set only, one would not be able to achieve a linear running time on trees.

Piecewise training for pCRFs

be the independent identically distributed training sample. For every μ=1,2,…,m, let ${\mathcal{V}}_{\mu }$ and ${\mathcal{E}}_{\mu }$ be the set of positions and edges in the neighborhood graph associated with x _μ , let $n_{\mu }=\left|{\mathcal{V}}_{\mu }\right|$ be the number of positions of the μ-th training example and let ${\{\text{I},\text{N}\}}^{n_{\mu }}$ be the set of all possible label sequences of this graph.

This data set is unbalanced as there are many more non-interface positions as interface positions. As customary for other machine learning approaches such as support vector machines and artificial neural networks [28], we here manipulated the ratio of positive and negative example positions for training in order to obtain reasonable results.

We have amplified the influence of the positive examples, rather than selecting various sets of training data by deleting negative ones as done in [28].

Let ν _I , ν _N , ν _II and ν _NN be the number of interface positions, the number of non-interface positions, the number of interface-interface edges, and the number of non-interface-non-interface edges in , respectively. Then we define the following two amplifier functions for all positions i and for all edges {i,j} of the m neighborhood graphs resulting from the training data .

$\begin{array}{ll}\hfill {\eta }_1\left(i\right):& =\left\{\begin{array}{c}\frac{{\nu }_N}{{\nu }_I}-1\text{if}{y}_i=\text{I};\\ 0\text{if}{y}_i=\text{N}.\end{array}\right.\\ \hfill {\eta }_2(i,j)& :=\left\{\begin{array}{c}\frac{{\nu }_{\mathit{\text{NN}}}}{{\nu }_{\mathit{\text{II}}}}-1\text{if}{y}_i=y_j=\text{I};\\ \frac{{\nu }_N}{{\nu }_I}-1\text{if}{y}_i\ne y_j;\\ 0\text{if}{y}_i=y_j=\text{N}.\end{array}\right.\end{array}$

To uniformly govern the influence of the amplifiers, we introduce an amplifier control parameter η₃∈ [ 0,1].

is the training-instance-specific normalization factor.

Unfortunately, maximizing this objective function in general is algorithmically intractable. Taking pattern from Sutton et al.[35] who introduced what they called piecewise training, we deal with this problem by disentangling the labels of nodes and edges. For μ=1,2,…,m, a non-coherent labeling$\mathbf{y}\in {\{\text{I},\text{N}\}}^{{\mathcal{V}}_{\mu }\times {\mathcal{E}}_{\mu }}$ of the neighborhood graph $\left({\mathcal{V}}_{\mu },{\mathcal{E}}_{\mu }\right)$ is any mapping that assigns to every position $v\in {\mathcal{V}}_{\mu }$ and every edge $e\in {\mathcal{E}}_{\mu }$ a label y _v ∈{I,N} and a pair of labels y _e ∈{I,N}², respectively.

as normalization factor. This makes the optimization problem computationally feasible.

The L-BFGS method [38] is used to solve it. That way we obtain the coefficient vectors α and β (see Equations 3), which depend on the amplifier control parameter η₃∈ [ 0,1].

To mitigate the negative consequences of disentanglement, we use a correction factor δ≥1. For any characteristics D and ι=0,1,…,γ−1, the weights of the bases edge features ${\varphi }_{\text{I},\text{I},\iota }^{\left(D\right)}$ and ${\varphi }_{\text{N},\text{N},\iota }^{\left(D\right)}$ (see Equation 5) are all multiplied by δ. Thus a change in classification along an edge is additionally penalized. The correction factor δ is set best between 1.15 and 1.25.

For our implementation of the training, we used the Java CRF package from Sunita Sarawagi at http://crf.sourceforge.net/.

Simulating unimodal scores of various precisions

We randomly chose 120 instances under the uniform distribution from the data set published by Keskin et al.[39] to perform our experiments. Let us refer to this set as KL-subset in the sequel. (It is accessible at http://ppicrf.informatik.uni-goettingen.de/index.html).

Zellner et al.[25] used the following determinations. A residue is defined to be part of the protein surface, if its relative solvent-accessible surface area is at least 5%[17]. A surface residue is said to constitute an inter-facial contact, if there exists at least one atom of this residue which has a van-der-Waals-sphere at a distance of at most 0.5 Å from the van-der-Waals sphere to any atom from a partner chain residue [39].

Based on [3, 12, 15, 20, 41], Li et al.[28] assume an inter-facial contact of a residue on a chain is assumed to be there, if any heavy atom of this residue is at distance of at most 5 Å from any heavy atom from a partner chain. The relative solvent-accessible surface area of surface residue is at least 15%.

Table 1 and Table 2 justify the following conclusion. Enhancing the threshold prediction by our pCRF works provided that the distributions of the interface scores as well as the non-interface scores are unimodal. The enhancement for the data set PlaneDimers is larger than for the KL-subset. This might be caused by the plain interface geometry of the complexes taken from PlaneDimers.

Utilizing the new data set due to Cukuroglu [40]

A main finding of Cukuroglu [40] relevant to protein-protein interface prediction is, that the average interface RASA value is greater than 40%. Since our method is designed to improve performance of a given residue-wise predictor, using this result is not in the scope of this paper. However, a CRF-based predictor integrating features for cliques of size greater than 2 is not beyond the range of current algorithmic capabilities. In such a model a feature set that discretizes the mean RASA value of cliques is promising.

Enhancing the PresCont server prediction on PlaneDimers

For the sake of completeness, we shortly review the residue characteristics used by PresCont.

Scores of local neighborhoods

They are evaluated by means of log-odd ratios of neighboring residue pair frequencies in interfaces as opposed to residue pair frequencies on complementary protein surface areas. The resulting scores are averaged both over the neighborhood of the positions under study and the rows of the MSA associated with the protein.

On the basis of Figure 5 we enhanced PresCont for thresholds θ∈ [ 0.500,0.625]. The decisive factor for this choice is that the PresCont score distributions for interface sites as well as non-interface positions above θ are “sufficiently close to” unimodal distributions. For every such θ, we set all scores less than or equal to θ to zero and then left the classification of all surface residues to the pCRF modified as follows. The residues of score zero are not taken into account when it comes to discretizing the protein characteristics (see Equations 4 and 5). Let us call this enhancing above θ.

To evaluate improvements we proceeded as when compiling Table 2. For every threshold θ under consideration another threshold θ^′ was chosen such that thresholding at θ^′ has the same specificity as enhancing above θ. The results are displayed in Table 4 and visualized for an individual protein by Figure 6. According to Table 4 the increase in sensitivity ranges from 4% to 7%. The true-positive predictions on the surface of the protein with PDB-Entry 1QM4 are compared in Figure 6, where again the specificity of the two classifiers is the same.

	tp	tn	fp	fn	Spec.	Sen.	MCC
Enhancing above 0.500	2181	23182	4145	1414	0.848	0.607	0.362
PresCont	2100	23197	4130	1495	0.849	0.584	0.346
Enhancing above 0.525	2303	22917	4410	1292	0.839	0.641	0.373
PresCont	2206	22912	4415	1389	0.838	0.614	0.353
Enhancing above 0.550	2507	22103	5224	1088	0.809	0.697	0.375
PresCont	2419	22102	5225	1176	0.809	0.673	0.358
Enhancing above 0.575	2560	21992	5335	1035	0.805	0.712	0.380
PresCont	2463	21915	5412	1132	0.802	0.685	0.358
Enhancing above 0.600	2379	22685	4642	1216	0.830	0.662	0.376
PresCont	2253	22780	4547	1342	0.834	0.627	0.356
Enhancing above 0.625	2287	23044	4283	1308	0.843	0.636	0.376
PresCont	2136	23049	4278	1459	0.843	0.594	0.346

The sensitivity increased that way by 4%-7%. For every pair of experiments, the number of true negatives (tn), false negatives (fn), false positives (fp) and true positives (tp) are displayed.

Discussing posterior decoding

As in the case of linear-chain CRFs, the generalized Viterbi algorithm can be transformed into a variant form of the forward algorithm. It might be the case that the following additional problem arises.

Let v₁,v₂,…,v _n be the ordering in which the positions of are traversed by the algorithm, and let $\hat{I}$ denote the set of position indices i<n such that v _i is not an element of the boundary ${\mathcal{ℬ}}^{\left(i\right)}$ of the history set ${\mathcal{\mathscr{H}}}^{\left(i\right)}$ at stage i. If $\hat{I}$ is not empty, we encounter an obstacle when it comes to sampling label sequences. For $i\in \hat{I}$, position v _i is not labeled in the course of the sampling procedure. That is why we augment the neighborhood graph so that those positions no longer exist, all predictions remain unchanged, and the order of magnitude of the running time is not increased. To this end, we complement the ordering v₁,v₂,…,v _n as follows. For every $i\in \hat{I}$, we insert a new node ${\hat{v}}_i$ between v _i and v_i+1. Having extended the neighborhood graph by these nodes not being associated with residue positions of the protein under study and by new edges $\{v_i,{\hat{v}}_i\}(i\in \hat{I})$, where for $i\in \hat{I}$ and $y_0.y_1,y_2\in \{\text{N},\text{I}\}{\Phi }_{{\hat{v}}_i}(y_0,\mathbf{x})={\Phi }_{\{v_i,{\hat{v}}_i\}}(y_1,y_2,\mathbf{x})=0,$ the above mentioned obstacle is eliminated without any influence on the prediction and the order of magnitude of the running time.

Proceeding now in a way analogous to the classical case, in every formula that is a building block of the generalized Viterbi algorithm the following two steps of replacement need to be performed.

First, for every position $i\in \mathcal{V}$, every edge $(i,j)\in \mathcal{E}$, every label y₀∈{I,N}, and every label pair (y₁,y₂)∈{I,N}², we replace Φ _i (y₀,x) with exp(Φ _i (y₀,x)), and Φ_{i,j}(y₁,y₂,x) with exp(Φ_{i,j}(y₁,y₂,x)).

Second, we replace sums with products and then maxima with sums.

Thus we obtain as analogues of the Viterbi variables ${\text{vit}}_{{\mathcal{\mathscr{H}}}_{\mu }}\left({\mathbf{y}}_{{\mathcal{ℬ}}_{\mu }}\right)$ defined by Equation 6 what we call component forward variables ${\text{cf}}_{{\mathcal{\mathscr{H}}}_{\mu }}\left({\mathbf{y}}_{{\mathcal{ℬ}}_{\mu }}\right)$.

For any assignment ${\mathbf{y}}_{{\mathcal{ℬ}}^{\left(i\right)}}$ (i>1), the forward variable ${\mathrm{f}}_i\left({\mathbf{y}}_{{\mathcal{ℬ}}^{\left(i\right)}}\right)$ is a nontrivial linear combination of forward variables ${\mathrm{f}}_{i-1}\left({\mathbf{y}}_{{\mathcal{ℬ}}^{(i-1)}}\right)$, where ${\mathbf{y}}_{{\mathcal{ℬ}}^{(i-1)}}$ ranges over some assignments of the boundary set ${\mathcal{ℬ}}^{(i-1)}$ at stage i−1. Analogous to the linear-chain case, a random backward walk through a state graph, with all possible assignments ${\mathbf{y}}_{{\mathcal{ℬ}}^{\left(i\right)}}$ (i=n,n−1,…,1) being the set of nodes, results in a random labeling of the positions, where each labeling is drawn with its posterior probability.

This sampling technique allows the efficient calculation of posterior probabilities at nodes and edges in a straightforward manner.