Estimating probabilistic context-free grammars for proteins using contact map constraints

Basic notations

Let Σ be a non-empty finite set of atomic symbols (representing for instance amino acid species). The set of all finite strings over this alphabet is denoted by Σ^∗. Let |x| denote the length of a string x. The set of all strings of length n is denoted by Σⁿ = {x ∈ Σ^∗:|x| = n}. Let x = x₁…x_n be a sequence in Σⁿ.

Contact constraints

Most protein sequences fold into complex spatial structures. Two amino acids at positions i and j in the sequence x are said to be in contact if distance between their coordinates in spatial structure d(i, j) is below a given threshold τ. A full contact map for a protein of length n is a binary symmetric matrix m^full = (m_i,j)_n×n such that m_i,j = [d(i, j) < τ], where [x] is the Iverson bracket (see Fig. 1C). Usually only a subset of the contacts is considered (see ‘Protein contact constraints’). A (partial) contact map for a protein of length n is a binary symmetric matrix m = (m_i,j)_n×n such that m_i,j = 1⟹d(i, j) < τ. Let d_u(i, j) be the length of the shortest path from ith to jth leaf in UST u for x. Given a threshold δ, UST u is said to be consistent with a contact map m of length n if (1)${\displaystyle m_{i,j}=1⟹d_u\left(i,j\right)<\delta .}$

For a contact map m of length n, let ${\mathcal{U}}_n^{\mathsf{m}}$ denotes the subset of ${\mathcal{U}}_n$ consistent with m, and ${\mathcal{U}}_x^{\mathsf{m}}$ denotes the subset of ${\mathcal{U}}_x$ consistent with m. Note that ${\mathcal{U}}_x^{\mathsf{m}}={\mathcal{U}}_n^{\mathsf{m}}\cap {\mathcal{U}}_x$. Analogous notations apply to parse trees.

Estimation

Learning grammar $\mathcal{G}=〈\Sigma ,V,v_0,R,\theta 〉$ can be seen as inferring the unfixed components of $\mathcal{G}$ with the aim of shifting the probability mass from the entire space of unlabeled syntactic trees ${\mathcal{U}}_{\ast }$ to the set of unlabeled syntactic trees for the target population 𝔘_target. In practice, only a sample of the target population can be used for learning, hence estimation is performed on 𝔘_sample⊆𝔘_target. Note that even in the most general case the set of terminal symbols Σ is implicitly determined by the sample; moreover the start symbol v₀ is typically also fixed. A common special case considered in this work confines learning grammar $\mathcal{G}$ to estimating θ for a fixed quadruple of non-probabilistic parameters 〈Σ, V, v₀, R〉 (which fully determine the non-probabilistic grammar G underlying $\mathcal{G}$). Given inferred grammar ${\mathcal{G}}_{\ast }$ and a query set of unlabeled syntactic trees ${\mathcal{U}}_{\mathrm{query}}$, probability $prob\left({\mathcal{U}}_{\mathrm{query}}\mid {\mathcal{G}}_{\ast }\right)$ is an estimator of the likelihood that ${\mathcal{U}}_{\mathrm{query}}$ belongs to population 𝔘_target.

Definitions

Let $\ddot{\mathcal{G}}=〈\Sigma ,V,v_0,R,\theta 〉$ be a probabilistic context-free grammar such that V = V_T⊎V_N, R = R_a⊎R_b⊎R_c, and ${\displaystyle \begin{array}{c}{R}_a=\left\{r_i:V_T\to \Sigma \right\},\hfill \\ R_b=\left\{r_j:V_N\to \left(V_N\cup V_T\right)\left(V_N\cup V_T\right)\right\},\hfill \\ R_c=\left\{r_k:V_N\to V_T{V}_N{V}_T\right\}.\hfill \end{array}}$ Subsets R_a, R_b and R_c are referred to as lexical, branching, and contact rules, respectively. Joint subset R_b∪R_c is referred to as structural rules. Grammars which satisfy these conditions are hereby defined to be in the Chomsky Form with Contacts (CFC). It happens that the toy grammar in Fig. 1B is in CFC. When a CFC grammar satisfies R_c = ∅, it is in the Chomsky Normal Form (CNF).

Non-terminal symbols in V_T, which can be rewritten only into terminal symbols are referred to as lexical non-terminals, while non-terminal symbols in V_N are referred to as structural non-terminals. Comparing the CFC grammar with the profile HMM, each match state of the latter can be identified with a unique lexical non-terminal, and emissions from a given state—with a set of lexical rules rewriting the non-terminal corresponding to the state.

Let m be a contact matrix compatible with the context-free grammar, i.e., no pair of positions in contact overlaps nor crosses boundaries of other pairs in contact (though pairs can be nested one in another): ${\displaystyle \forall \left(i,j\right)m_{i,j}=1\wedge \left(i\le k\le j\oplus i\le l\le j\right)\Rightarrow m_{k,l}=0,}$ where ⊕ denotes the exclusive disjunction, and positions in contact are separated from each other by at least 2: ${\displaystyle \forall \left(i,j\right)i<j+2.}$

Let distance threshold in tree δ = 4. Then a complete parse tree y generated by $\ddot{\mathcal{G}}$ is consistent with m only if for all m_i,j = 1 derivation ${\displaystyle {\alpha }_{1,i-1}{v}_k{\alpha }_{j+1,n}\stackrel{\ast }{\Rightarrow }{\alpha }_{1,i-1}{x}_i{v}_l{x}_j{\alpha }_{j+1,n}}$ is performed with a string of production rules ${\displaystyle \left[v_k\to v_t{v}_l{v}_u\right]\left[v_t\to x_i\right]\left[v_t\to x_j\right],}$ where α_i,j ∈ (Σ∪V)^j−i+1, v_k, v_l ∈ V_N and v_t, v_u ∈ V_T.

According to this definition, the left-hand (right-hand) side parse tree in Fig. 1E is consistent (not consistent) with the contact map in Fig. 1C.

Parsing

Given an input sequence x of length n and a grammar in the CFC form $\ddot{\mathcal{G}}$, $prob\left(x\mid \ddot{\mathcal{G}}\right)\equiv prob\left({\mathcal{Y}}_x\mid \ddot{\mathcal{G}}\right)={\sum }_{y\in {\mathcal{Y}}_x}prob\left(y\mid \ddot{\mathcal{G}}\right)$ can be calculated in O(n³) by a slightly modified probabilistic Cocke-Kasami-Younger bottom-up chart parser (Cocke, 1969; Kasami, 1965; Younger, 1967). Indeed, productions in R_a⊎R_b conforms to the Chomsky Normal Form (Chomsky, 1959), while it is easy to see that productions in R_c requires only O(n²). The algorithm computes $prob\left(x\mid \ddot{\mathcal{G}}\right)=prob\left({\mathcal{Y}}_x\mid \ddot{\mathcal{G}}\right)$ in chart table P of dimensions n × n × |V|, which effectively sums up probabilities of all possible parse trees ${\mathcal{Y}}_x$. In the first step, probabilities of assigning lexical non-terminals V_T for each terminal in the sequence x are stored in the bottom matrix P₁ = P[1, :, :]. Then, the table P is iteratively filled upwards with probabilities $\mathsf{P}\left[j,i,v\right]=prob\left(v\stackrel{\ast }{\Rightarrow }{x}_i\dots x_{i+j-1}\mid v\in V,\ddot{\mathcal{G}}\right)$. Finally, $prob\left({\mathcal{Y}}_x^{\mathsf{m}}\mid \ddot{\mathcal{G}}\right)=\mathsf{P}\left[n,1,v_0\right]$.

New extended version of the algorithm (Fig. 2) computes $prob\left({\mathcal{Y}}_x^{\mathsf{m}}\mid \ddot{\mathcal{G}}\right)$, i.e., it considers only parse trees ${\mathcal{Y}}_x^{\mathsf{m}}$ which are consistent with m. To this goal it uses an additional table C of dimensions ∑(m)∕2 × n × |V_T|. After completing P₁ (lines 10–12), probabilities of assigning lexical non-terminals V_T at positions involved in contacts are moved from P₁ to C (lines 13–21) such that each matrix C_p = C[p, :, :] corresponds to p-th contact in m. In the subsequent steps C can only be used to complete productions in R_c; moreover both lexical non-terminals have to come either from P₁ or C, they can never be mixed (lines 35–40). The computational complexity of the extended algorithm is still O(n³) as processing of productions in R_c has to be multiplied by iterating over the number of contact pairs in m, which is O(n) since the cross-serial dependencies are not allowed.

Calculating $prob\left({\mathcal{U}}_n^{\mathsf{m}}\mid \ddot{\mathcal{G}}\right)$

This section shows effective computing $prob\left({\mathcal{U}}_n^{\mathsf{m}}\mid \ddot{\mathcal{G}}\right)$, which is the denominator for the contrastive estimation of ${\mathcal{G}}_{\mathrm{CE}\left(\mathsf{m}\right)}$ (cf. ‘Estimation’). Given a sequence x of length n, a corresponding matrix m of size n × n and a grammar $\ddot{\mathcal{G}}$, the probability of the set of trees over any sequence of length n consistent with m is ${\displaystyle prob\left({\mathcal{U}}_n^{\mathsf{m}}\mid \ddot{\mathcal{G}}\right)\equiv \sum _{x\in {\Sigma }^n}prob\left({\mathcal{U}}_x^{\mathsf{m}}\mid \ddot{\mathcal{G}}\right)=\sum _{x\in {\Sigma }^n}\sum _{y\in {\mathcal{Y}}_x^{\mathsf{m}}}prob\left(y\mid \ddot{\mathcal{G}}\right).}$ Given grammar $\ddot{\mathcal{G}}$, any complete derivation r is a composition $r=\dot{r}\circ \tilde{r}$, where $\dot{r}\in {\left(R_a\right)}^{\ast }$ and $\tilde{r}\in {\left(R_b\cup R_c\right)}^{\ast }$. Let y be the parse tree corresponding to derivation r, and let $\tilde{y}$ be an incomplete parse tree corresponding to derivation $\tilde{r}$. Note that for any y corresponding to $r=\dot{r}\circ \tilde{r}$ there exists one and only one $\tilde{y}$ corresponding to $\tilde{r}$. Let ${\tilde{\mathcal{Y}}}_x^{\mathsf{m}}$ denote the set of such incomplete trees $\tilde{y}$. Note that labels of the leaf nodes of $\tilde{y}$ are lexical non-terminals ∀(i)α_i,i ∈ V_T, and that $\dot{r}$ represents the unique left-most derivation $yield\left(\tilde{y}\right)\stackrel{\ast }{\Rightarrow }x$. Thus, ${\displaystyle \sum _{x\in {\Sigma }^n}\sum _{y\in {\mathcal{Y}}_x^{\mathsf{m}}}prob\left(y\mid \ddot{\mathcal{G}}\right)=\sum _{x\in {\Sigma }^n}\sum _{\tilde{y}\in {\tilde{\mathcal{Y}}}_x^{\mathsf{m}}}prob\left(\tilde{y}\mid \ddot{\mathcal{G}}\right)\cdot prob\left(yield\left(\tilde{y}\right)\stackrel{\ast }{\Rightarrow }x\mid \ddot{\mathcal{G}}\right).}$ Note that value of the expression will not change if the second summation is over $\tilde{y}\in {\tilde{\mathcal{Y}}}_n^{\mathsf{m}}$ since $\forall \left(\tilde{y}⁄\in {\tilde{\mathcal{Y}}}_x^{\mathsf{m}}\right)prob\left(yield\left(\tilde{y}\right)\stackrel{\ast }{\Rightarrow }x\mid \ddot{\mathcal{G}}\right)=0$. Combining with observation that $prob\left(\tilde{y}\mid \ddot{\mathcal{G}}\right)$ does not depend on x, the expression can be therefore rewritten as: ${\displaystyle \sum _{x\in {\Sigma }^n}\sum _{y\in {\mathcal{Y}}_x^{\mathsf{m}}}prob\left(y\mid \ddot{\mathcal{G}}\right)=\sum _{\tilde{y}\in {\tilde{\mathcal{Y}}}_n^{\mathsf{m}}}prob\left(\tilde{y}\mid \ddot{\mathcal{G}}\right)\cdot \sum _{x\in {\Sigma }^n}prob\left(yield\left(\tilde{y}\right)\stackrel{\ast }{\Rightarrow }x\mid \ddot{\mathcal{G}}\right).}$ However, if $\ddot{\mathcal{G}}$ is proper, then $\forall \left(\tilde{y}\in {\tilde{\mathcal{Y}}}_n^{\mathsf{m}}\right){\sum }_{x\in {\Sigma }^n}prob\left(yield\left(\tilde{y}\right)\stackrel{\ast }{\Rightarrow }x\mid \ddot{\mathcal{G}}\right)=1$, as: ${\displaystyle \begin{array}{cc}\hfill & \sum _{x\in {\Sigma }^n}prob\left(yield\left(\tilde{y}\right)\stackrel{\ast }{\Rightarrow }x\mid \ddot{\mathcal{G}}\right)=\sum _{x\in {\Sigma }^n}\prod _{i=1}^n\theta \left({\alpha }_{i,i}\to x_i\right)=\hfill \\ \hfill & \sum _{x\in {\Sigma }^n}\theta \left({\alpha }_{1,1}\to x_1\right)\cdot \cdots \cdot \theta \left({\alpha }_{n,n}\to x_n\right)=\hfill \\ \hfill & \theta \left({\alpha }_{1,1}\to a_1\right)\cdot \theta \left({\alpha }_{2,2}\to a_1\right)\cdot \cdots \cdot \theta \left({\alpha }_{n-1,n-1}\to a_1\right)\cdot \theta \left({\alpha }_{n,n}\to a_1\right)+\hfill \\ \hfill & \theta \left({\alpha }_{1,1}\to a_1\right)\cdot \theta \left({\alpha }_{2,2}\to a_1\right)\cdot \cdots \cdot \theta \left({\alpha }_{n-1,n-1}\to a_1\right)\cdot \theta \left({\alpha }_{n,n}\to a_2\right)+\hfill \\ \hfill & ⋮\hfill \\ \hfill & \theta \left({\alpha }_{1,1}\to a_{\left|\Sigma \right|}\right)\cdot \theta \left({\alpha }_{2,2}\to a_{\left|\Sigma \right|}\right)\cdot \cdots \cdot \theta \left({\alpha }_{n-1,n-1}\to a_{\left|\Sigma \right|}\right)\cdot \theta \left({\alpha }_{n,n}\to a_{\left|\Sigma \right|}\right)=\hfill \end{array}}$ ${\displaystyle \left(\begin{array}{cc}\hfill & \theta \left({\alpha }_{1,1}\to a_1\right)\cdot \theta \left({\alpha }_{2,2}\to a_1\right)\cdot \cdots \cdot \theta \left({\alpha }_{n-1,n-1}\to a_1\right)+\hfill \\ \hfill & \theta \left({\alpha }_{1,1}\to a_1\right)\cdot \theta \left({\alpha }_{2,2}\to a_1\right)\cdot \cdots \cdot \theta \left({\alpha }_{n-1,n-1}\to a_2\right)+\hfill \\ \hfill & ⋮\hfill \\ \hfill & \theta \left({\alpha }_{1,1}\to a_{\left|\Sigma \right|}\right)\cdot \theta \left({\alpha }_{2,2}\to a_{\left|\Sigma \right|}\right)\cdot \cdots \cdot \theta \left({\alpha }_{n-1,n-1}\to a_{\left|\Sigma \right|}\right)\hfill \end{array}\right)\cdot \sum _{s=1}^{\left|\Sigma \right|}\theta \left({\alpha }_{n,n}\to a_s\right),}$ where a_s ∈ Σ. Since $\ddot{\mathcal{G}}$ is proper then $\forall \left(v\in V_T\right){\sum }_{s=1}^{\left|\Sigma \right|}\theta \left(v\to a_s\right)=1$ and therefore the entire formula evaluates to 1, which can be easily shown by iterative regrouping. This leads to the final formula: ${\displaystyle prob\left({\mathcal{U}}_n^{\mathsf{m}}\mid \ddot{\mathcal{G}}\right)=\sum _{\tilde{y}\in {\tilde{\mathcal{Y}}}_n^{\mathsf{m}}}prob\left(\tilde{y}\mid \ddot{\mathcal{G}}\right).}$ Technically, ${\sum }_{\tilde{y}\in {\tilde{\mathcal{Y}}}_n^{\mathsf{m}}}prob\left(\tilde{y}\mid \ddot{\mathcal{G}}\right)$ can be readily calculated by the bottom-up chart parser by setting ∀(r_k ∈ R_a)θ(r_k) = 1.

Learning

Our evolutionary learning framework used the genetic algorithm where each individual represented a whole grammar, the approach known as the Pittsburgh style (Smith, 1980). For a given underlying non-probabilistic CFG $\ddot{G}$ and the positive training sample, the framework estimated probabilities θ of the corresponding PCFG $\ddot{\mathcal{G}}=〈\ddot{G},\theta 〉$. Unlike previous applications of the framework in which probabilities of the lexical rules were fixed according to representative physicochemical properties of amino acids (Dyrka & Nebel, 2009; Dyrka, Nebel & Kotulska, 2013), in this research probabilities of all rules were subject to evolution. The objective functions were implemented for the maximum-likelihood estimator ${\ddot{\mathcal{G}}}_{\text{ML}}$, and for the constrastive estimators ${\ddot{\mathcal{G}}}_{\mathrm{CE}\left(X\right)}$ and ${\ddot{\mathcal{G}}}_{\mathrm{CE}\left(\mathsf{m}\right)}$. Besides, the setup of the genetic algorithm closely followed that of Dyrka & Nebel (2009).

Performance measures

Performance of grammars was evaluated using a variant of the 8-fold Cross-Validation scheme in which 6 parts are used for training, 1 part is used for validation and parameter selection, and 1 part is used for final testing and reporting results (the total of 56 combinations). The negative set was not used in the training phase. For testing, protein sequences were scored against the null model (a unigram), which assumed global average frequencies of amino acids, no contact information, and the length of query sequence. The amino acid frequencies were obtained using the online ProtScale tool for the UniProtKB/Swiss-Prot database (Gasteiger et al., 2005).

Materials

Probabilistic grammars were estimated for three samples of protein fragments related to functionally relevant gapless motifs (Sigrist et al., 2002; Bailey & Elkan, 1994). Within each sample, all sequences shared the same length, which avoided sequence length effects on grammar scores (this could be resolved by an appropriate null model). For each sample, one experimentally solved spatial structure in the Protein Data Bank (PDB) (Berman et al., 2000) was selected as a representative. The three samples included amino acid sequences of two small ligand binding sites (already analyzed in Dyrka & Nebel (2009) and a functional amyloid (Table 1):

• CaMn: a Calcium and Manganese binding site found in the legume lectins (Sharon & Lis, 1990). Sequences were collected according to the PROSITE PS00307 pattern (Sigrist et al., 2013) true positive and false negative hits. Original boundaries of the pattern were extended to cover the entire binding site, similarly to Dyrka & Nebel (2009). The motif folds into a stem-like structure with multiple contacts, many of them forming nested dependencies, which stabilize anti-parallel beta-sheet made of two ends of the motif (Fig. 3A based on pdb:2zbj (De Oliveira et al., 2008));

• NAP: the Nicotinamide Adenine dinucleotide Phosphate binding site fragment found in an aldo/keto reductase family (Bohren et al., 1989). Sequences were collected according to the PS00063 pattern true positive and false negative hits (four least consistent sequences were excluded). The motif is only a part of the binding site of the relatively large ligand. Intra-motif contacts seem to be insufficient for defining the fold, which depends also on interactions with amino acids outside the motif (Fig. 3B based on pdb:1mrq (Couture et al., 2003));

• HET-s: the HET-s-related motifs r1 and r2 involved in the prion-like signal transduction in fungi identified in a recent study (Daskalov, Dyrka & Saupe, 2015). The largest subset of motif sequences with length of 21 amino acids was used to avoid length effects on grammar scores. When interacting with a related motif r0 from a cooperating protein, motifs r1 and r2 adopt the beta-hairpin-like folds which stack together. While stacking of multiple motifs from several proteins is essential for stability of the structure, interactions between hydrophobic amino acids within a single hairpin are also important. In addition, correlation analysis revealed strong dependency between positions 17 and 21 (Daskalov, Dyrka & Saupe, 2015) (corresponding to L276 and E280 in Fig. 3C based on Van Melckebeke et al. (2010)).

Negative samples were designed to roughly approximate the entire space of protein sequences. They were based on the negative set from (Dyrka & Nebel, 2009), which consisted of 829 single chain sequences of 300–500 residues retrieved from the Protein Data Bank (Berman et al., 2000) at identity of 30% (accessed on 12th December 2006). For each positive sample, the corresponding negative sample was obtained by cutting the basic negative set into overlapping subsequences of the length of positive sequences.

All samples were made non-redundant at level of sequence similarity around 70% using cd-hit (Li & Godzik, 2006), which significantly reduced their cardinalities. The threshold balanced the size of positive samples, distribution of their variability, and inter-fold diversity. Overall diversity of samples ranged from the most homogeneous CaMn (average identity of 49%) to the most diverse HET-s, which consisted of 5 subfamilies (Daskalov, Dyrka & Saupe, 2015) (average identity of 21%). The ratio between negative and positive samples was high and varied from 1190:1 for CaMn to 207:1 for HET-s. Contact pairings were assigned manually and collectively to all sequences in each set based on a selected representative spatial structure in the PDB database (Fig. 3).

Performance

The implementation of the framework for learning PCFGs for protein sequences using contact constraints, presented in ‘Application to contact grammars’ and ‘Evaluation’, is evaluated with reference to learning without the constraints. For grammars with the contact rules (CFC), probabilities of rules θ were estimated either using training samples made of sequences coupled with a contact map, or using sequences alone. For grammars without the contact rules (CNF), probabilities of rules were estimated using sequences alone, since these grammars cannot generate parse trees consistent with contact maps for the distance threshold δ = 4.

Descriptive power.

For evaluation of the descriptive power of the PCFG-CM approach, the rule probabilities were estimated using the maximum-likelihood estimator (denoted ML) and the contrastive estimator with regard to the sequence set (denoted CE(X)). Descriptive value of the most probable parse trees generated using the resulting probabilistic grammars for test sequences without contact information is presented in Table 3. Efficiency of the learning was measured on the basis of the recall at the distance threshold δ = 4 with regard to the context-free compatible contact map m used in the training. Consistency of the most likely parse tree with the protein structure was measured on the basis of the precision of contact prediction at the distance threshold δ = 4 with regard to all contacts in the reference spatial structure with separation in sequence of at least 3. Both measures are not suitable for assessing grammars without contact rules. Therefore, average precision over all thresholds δ was used as a complementary measure of consistency of the most likely trees with the protein structure. Note that the AP scores achievable for a context-free parse tree are reduced by overlapping of pairings.

Table 3:

Descriptive quality of the most likely parse trees derived from sequences alone.

In terms of recall at the distance threshold δ = 4 w.r.t. the training contact map m, and precision at δ = 4 (and AP over thresholds δ) w.r.t. the full contact map of the reference pdb structure for sequence separation 3 +. Note that the shortest length of any path between leaves in the most likely parse trees of the CNF grammar equals 5, which makes measures using δ = 4 unutile.

Grammar	CNF	CFC		CFC		CFC
Estimation	ML	ML		ML		CE(X)
Train w/contacts	n/a	no		yes		yes
Reference	pdb	m	pdb	m	pdb	m	pdb
CaMn	(0.24)	0.45	0.69 (0.53)	0.92	0.87 (0.66)	0.98	0.84 (0.66)
NAP	(0.16)	0.00	0.14 (0.12)	0.96	0.64 (0.29)	0.96	0.64 (0.29)
HET-s	(0.08)	0.02	0.13 (0.14)	0.79	0.52 (0.24)	0.97	0.57 (0.27)

DOI: 10.7717/peerj.6559/table-3

The baseline is the result for grammars with the contact rules estimated without contact constraints. The most likely parse trees generated using these grammars conveyed practically no information about contacts for NAP and HET-s (recall w.r.t. contact map m close to zero) and limited information about contacts for CaMn (moderate recall of 0.45), see Fig. 4. Learning with the contact constraints resulted in increase of the recall to 0.79–0.98, which testified efficiency of the process.

Importantly, consistency of the most likely parse trees with the protein structure measured by the precision followed a similar pattern and increased from 0.13 for HET-s, 0.14 for NAP, and 0.69 for CaMn when grammars with the contact rules were estimated without a contact map, to 0.52–0.57, 0.64, and 0.84–0.87, respectively, when grammars were estimated with a contact map. Accordingly, evaluation in terms of the average precision over distance thresholds indicated that distances in the most likely parse trees better reflected the protein structure if grammars were trained with the contact constraints, as illustrated in Fig. 4.

Searching for related motifs

In this section probabilistic grammars for HET-s r1 and r2 motifs, learned in the proposed estimation scheme, are applied to solving a practical problem of searching for related r0 motifs in a limited-size dataset (around 1,000–5,000 sequences) based on (Dyrka et al., 2014; Daskalov, Dyrka & Saupe, 2015).

Making generalized descriptors

In this section the generalizing potential of PCFG descriptors is illustrated by learning a single grammar for two non-homologous but functionally related Calcium-binding motifs.

Grammatical descriptors.

Due to presumed higher complexity of the model, several variants of grammar rules were again used for training. The best fitting to the training sample was achieved with the same variant as in the previous example. Also in this case, learning with the contact constraints significantly improved fitness to the training data (probability mass distributed over the training set increased roughly 20 times on average over six runs).

The diagram showing the 36 most significant rules (all with probability of at least 0.05) and dependencies between structural non-terminals (possible derivations) of the single best grammar are shown in Fig. 5A. Of note is a pair of structural non-terminal symbols u and v (orange), which can be used to generate paired stretches of hydrophobic (u → ava) and other residues (v → buc). The feature was used to model the pair of beta-strands in the stem part of CaMn (Figs. 5B, 5C). By extending the cooperation between u and v with the derivation path through the structural non-terminal t (pink, v → atb, t → •u•), the grammar generates hydrophobic residues with periodicity of 3, typical to helices, as used in modeling the pair of alpha-helices of the EF hand (Figs. 5D, 5E). To finish a derivation, it is typically necessary to use the structural non-terminal w (green), which is likely to generate lexical non-terminals b and c which emit amino acids with high propensity to binding Calcium (aspartic and glutamic acids, aspargine, serine, and threonine (Bindreither & Lackner, 2009).

Clearly, the grammar has its limitations. The number of only three lexical non-terminals is likely insufficient, as suggested by the unusual merging of hydrophobic alanine with the charged amino acids in one group emitted through symbol b. Also detailed analysis of parse trees reveal inaccuracies possibly resulting from over-generalization. Most notably, the beta-hairpin generating rules (orange) were used to model a part of the binding loop of CaMn (Fig. 5B). Moreover, the residues directly involved in the Calcium binding in 1gsl, according to Ligplot (D130, W133, N135 and D140), were not generated with the non-terminal w. Finally, contact rules used to model the loop of the EF hand did not generate pairs of residues which are actually in contact. Yet, the overall topologies of the trees were rather consistent with the structures.