Graph pyramids for protein function prediction

Motivation and contributions

Relationships among biological sequences can be effectively represented by building a PPS network. However, in protein families, modularity, local clustering and scale-free topology coexist [19]. Thus use of single graph for modeling them has not been efficient. Here we propose the graph pyramid approach, where multiple graph features are used at different levels for modeling protein families. Along with this, the proposed algorithm using hierarchical voting scheme, tries to blend important characteristics from the PPI network and ISS methods for protein classification. This makes it possible to more objectively and reasonably predict the protein functions with high accuracy.

An edge connecting a pair of nodes with the largest weight in the PPS network represents the strongest match and it accounts for the global feature in similarity space. Only considering the nodes connected through the edges having lower weights, constitutes weak local features and with also considering interactions among weakly interacting nodes boost up the local features. In addition, PPS networks are analyzed hierarchically in the form of the Graph Pyramid (GP). It helps to extract more vital information about proteins from the network topology and looks for stronger global and local features. Unlike most of the existing methods, the proposed approach also shows that 'how the query sequence interacts with each protein in the family'. PPS interaction topology shows 'small world network' [19] properties like the PPI network (see Figure 2(b)). This helps and guides in devising important graph features (discussed in section 'Graph structured features').

This paper is organized as follows. In the section 'Methods', we present the PPS network construction approach and its modeling via graph pyramid. The same section also discusses the important graph topological or graph structured features and the hierarchical protein classification algorithm. Experimental evaluation results of the proposed method, their comparison with existing prediction methods and the techniques of searching for the optimal algorithm parameters are presented in section 'Results'. We conclude our work with discussions about related issues and the direction of our future work.

Graph construction

The proposed algorithm uses the graphical modeling of protein sequences from each protein family/class. Consider the large protein database having M classes. Let the set of class labels be given as ${ℂ}_M=\left\{c_1,c_2,\cdots ,c_M\right\}$. Consider one of the protein classes, $c\in {ℂ}_M$, then let $s_{t_i}^c$ be the i ^th training sequence in that class. Similarly $s_q^{c_q}$ is the new query or test sequence and its class label c _q is what we have to find.

These edges form a set $E^c=\left\{e_{1,1}^c,e_{1,2}^e,\cdots ,e_{N,N}^c\right\}$. Now the graph of a protein class $c\in {ℂ}_M$ is given by G(V ^c , E ^c ). This is a weighted and an undirected graph. An edge weight is nothing but the strength of protein sequence similarity.

To construct the graph of each protein class, we just need to consider interaction (EB-score) among all protein sequences within that class. Number of proteins in a class is far smaller than that of the entire database. So graphs of all protein classes can be easily and independently constructed by using protein-protein BLAST within the corresponding classes.

Graph analysis

In protein similarity graphs, modularity, local clustering and scale-free topology [19] coexist. To explain this phenomenon we need the hierarchy, so graphs are analyzed in the hierarchical manner. At each hierarchical level, the edges with weights lower than a certain threshold are pruned. Now surviving edges are considered to be weightless. So the graph structure changes along hierarchical levels, as well as the graph becomes unweighted and remains undirected. This hierarchical analysis helps to extract different graph features for weakly similar hits (sequence matches) and thus captures the complex relationship between sequence similarity and protein function.

The corresponding graph is given as $G\left(V^c,E_t^c\right)$. For notation simplicity lets represent it as G(V ^c , E ^c , t), instead of $G\left(V^c,E_t^c\right)$ and note that G(V ^c , E ^c ) = G(V ^c , E ^c , 0).

and note that $|{\mathbb{S}}_1\uplus {\mathbb{S}}_2|=|{\mathbb{S}}_1+{\mathbb{S}}_2|$, thus added set contains repetitive elements when ${\mathbb{S}}_1\cap {\mathbb{S}}_2\ne \mathrm{0̸}$. For example, let ${\mathbb{S}}_1=\left\{c_1,c_2\right\}$ and ${\mathbb{S}}_2=\left\{c_2,c_3\right\}$ then ${\mathbb{S}}_1\uplus {\mathbb{S}}_2=\left\{c_1,c_2,c_3,c_2\right\}$. This operation obeys the associative and the commutative laws like numerical addition. As defined earlier, $s_q^{c_q}$ is the query sequence from the unknown class c _q . Now $V_q^c=V^c\uplus s_q^{c_q}$, and edges among the vertices in the set is given by an edge set $E_q^c$. After adding $s^{c_q}$ to the original graph G(V ^c , E ^c ), we will get the new graph $G\left(V_q^c,E_q^c\right)$.

Graph Structured Features (GSF)

Most of the real world and biological (scale-free) networks communicate via a few highly connected nodes known as Hubs. These hubs determine the properties of networks [19]. In real world networks like airline route maps, the important cities form hubs. Proteins with high degrees of connectedness are more likely to be essential for survival than proteins with lesser degrees [22]. Gene duplication leads to growth and preferential attachment in biological networks [23]. This leads to translating the proteins having high similarity. This shows the possibility of hub formation in the protein family graphs, G(V ^c , E ^c ), in the similarity space.

Figure 2(b) shows building of the PPS network with a varying threshold for one of the COG [6] protein families. We can see that as the threshold is lowered, trivially more edges are formed but most of them are associated with only particular nodes (hubs). Thus hubs are getting stronger and becoming more evident in the PPS network. We are not interested in the detailed assessment of whether the network is scale-free (a power-law degree distribution [19]) or not. But the above analysis helps to guide us for finding proper features which take graph structure (i.e. complex relationships among the protein sequences) into account. Also, the different protein families have different characteristics. Thus use of a single graph feature may not be effective. Features are selected such that they could extract different but vital network information.

Average Clustering coefficient (AC)

For a node n, the clustering coefficient C _n , measures the extent to which neighbors of n are also the neighbors of each other [19]. Thus it is nothing but the density of sub-graph induced by the neighborhood of n. Consider the graph G(V, E) with n ∈ V . Let ℕ _n be the number of neighbors of n and ${\mathbb{E}}_n$ be the number of edges between them, then the AC is given by,

$AC\left(G\left(V,E\right)\right)=\frac{1}{\left|V\right|}{\displaystyle \sum _{n\in V}}\left(\frac{2{\mathbb{E}}_n}{{\mathbb{N}}_n\left({\mathbb{N}}_n-1\right)}\right).$

(7)

A clique is a maximal complete sub-graph where all the vertices are connected. C _n quantifies how close the neighbors of a node are, to form a clique among themselves. It represents the potential modularity of a network and C _n of the most real networks is much larger than that of a random networks. AC distribution is found to be effective for an identification of a modular organization of the metabolic networks [24]. Consider the example shown in Figure 3(a). Node c has 4 neighbors $\left({\mathbb{N}}_c\right)$, having 2 connected edges $\left({\mathbb{E}}_c\right)$ among them (a to b and d to e), which forms $C_c=\frac{1}{3}$ and then $AC=\frac{13}{15}$. When $s_q^{c_q}$ is attached to G(V ^c , E ^c ), it may change its AC. For a given the change in AC is given as,

$\text{Δ}AC\left(c,t\right)=AC\left(G\left(V_q^c,E_q^c,t\right)\right)-AC\left(G\left(V^c,E^c,t\right)\right).$

(8)

Algorithm

Graphs are analyzed hierarchically (as discussed in section 'Graph analysis') and threshold plays an important role in making hierarchical graph structure. If the query can interact at the higher layer of GP (see Figure 2(a)), then it means it is a strong interaction and it accounts for the global feature. Because in GP as the level rises, the threshold also increases and at every level, the graph edges can only be formed if their strength is greater than the given threshold. Here interactions are in the sequence similarity space, so a sstrong interaction means the highest similarity between corresponding sequences, which occurs when both sequences match at most of the nucleic acids (i.e. match globally). Hence strong interaction accounts for the global feature and similarly weak interactions account for local features.

Let $\mathbb{T}=\{\left\{t_{AC}\right\},\left\{t_{RC}\right\},\left\{t_{SM}\right\},\left\{t_{TR}\right\},\left\{t_{GE}\right\}\}$, be a set of sets, defining few threshold levels corresponding to each GTF. Consider, t _AC = {t₁, t₂, t₃} and t₁ < t₂ < t₃. Let us consider there are only 3 protein classes i.e. $\left|{ℂ}_M\right|=3$. Figure 4 explains the hierarchical query classification subroutine (algorithm 1) for the feature AC. Each q _i is analyzed first at the highest level (t₃) where we look for class $c\in {ℂ}_M$, having ΔAC(c, t) > T _AC and collect them in ${ℂ}_{AC}\subseteq {ℂ}_M$. For query q₁ we can not find any such classes c at t₃, so we descend the GP to t₂ and discover ${ℂ}_{AC}=\left\{c_2,c_3\right\}$ with threshold t ^∗ = t². The secondary threshold (T _AC ) is necessary, otherwise there will be many spurious classes i.e. false positives (FP), having nonzero change in average clustering coefficient $\left(\text{Δ}AC\left(\cdot \right)\ne 0\right)$.

Maximum value of H(·) is 1, so all the arguments (classes c) are assigned to ${ℂ}_{AC}$ whenever it produces output 1. For reducing FP, the subroutine for SM and TR slightly changes (see algorithm 2). Here we look for k maximally influenced classes from ${ℂ}_{RC}$ by the query. Thus classes from only ${ℂ}_{RC}$ are assessed (voted) again by the features SM and TR. Subroutine for GE finds the class from ${ℂ}_{SM}\cup {ℂ}_{TR}$ for which ΔGE(·) is maximum. Thus this produces a hierarchical class voting scheme which helps to improve the classification accuracy and to reduce the computational load.

Each GSF has an ability to extract different information from different levels of the protein family GP. However, applying each GSF to entire GP, is computationally inefficient when dealing with large numbers of protein families. In addition, it may add up the FP, when a decision is being made in the much lower GP level than the level defined by t ^∗ . To avoid these problems, the GP based hierarchical voting scheme is necessary. And the rational behind placing different GSF at different GP levels, is explained in the next 'Results' section.

Algorithm 1 Graph pyramid search subroutine

input: $s_q^{c_q}$; secondary threshold T _AC ; primary threshold $t_{AC}=\left\{t_1,t_2,\cdots ,t_n\right\}$, where t _i >t _j for i >j

output: ${ℂ}_{AC},t^{*}$

1: ${ℂ}_{AC}=\varnothing$

2: for t ∈ t _AC from t _n to t₁ do

3:   if ${ℂ}_{AC}=\varnothing$ then

4:      for all classes $c\in {ℂ}_M$do

5:         I _AC = H (ΔAC(c,t)−(T _AC ))

6:         ${ℂ}_{AC}={ℂ}_{AC}\uplus \delta \left(c,I_{AC}\right)$

7:         t* = t

8:      end for

9:   end if

10: end for

Dataset and evaluation details

Proposed method is evaluated on entire COG database [6]. It is the protein database of Clusters of Orthologous Groups (COG). It consists of 4,873 COG (protein families), having total 138,458 proteins from 66 different genomes. Approximately 10% sequences from each COG, are selected randomly, which has produced 14,086 test sequences. This procedure is repeated 5 times further to get average performance. First, each GSF is tested independently, for various thresholds without using hierarchical voting scheme. Here, the class which produces maximum change in the GSF for a given $s_q^{c_q}$, is selected as the output. These output labels are produced with either correct decision (cd), wrong decision (wd) or no decision (nd, when $\left|c_q\right|\ne 1$). Then let the performance measure be defined as, $precision=\frac{cd}{cd+wd}$.

Rational behind hierarchical voting

Figure 5(a) shows the plots of precision Vs threshold for all GSF. High precision indicates low wrong decisions. AC and RC produce high precision as the threshold rises, so these GSF are appointed to work at higher levels of GP. So at a high threshold, it is more likely that ${ℂ}_{AC}$ and ${ℂ}_{RC}$ contain true output class. Thus it is sufficient to apply other GSF, to the GP generated from either ${ℂ}_{AC}$ or ${ℂ}_{RC}$. On the other hand, GE produces high precision for low thresholds. One of the possible reasons behind this is that, the higher the threshold, the sparser will be the graph. So eigen-decomposition of the adjacency matrix to calculate GE will not give any information as ΔGE is close to zero for all classes. While at low threshold, the original graph becomes dense and query also interacts with almost all nodes in the true output class graph. Immensity in the interaction at lower threshold, helps GE to detect the true class easily and correctly. This forces GE to work at lower levels of GP, with the smallest search space as ${ℂ}_{SM}\cup {ℂ}_{TR}$. With similar arguments, SM and T R are placed at intermediate GP levels, allowing them to look for c _q in ${ℂ}_{RC}$. So in the algorithm 2, hierarchy as well as input search space for GSF are organized carefully. This helps to reduce the search space dramatically for other GSF and thus speeds up the algorithm (Figure 6 for speed comparison).

Algorithm 2 Classification by hierarchical voting

training: all G(V^c,E^c) are constructed $\forall c\in {ℂ}_M$

input: $s_q^{c_q},r,p,k$, primary threshold set $\mathbb{T}$, secondary thresholds T _AC , T _RC

output: c _q

1: ${ℂ}_{AC}\leftarrow \underset{c\in {ℂ}_M}{\text{arg max}}H\left(\text{Δ}AC\left(c,t^{*}\right)-T_{AC}\right)$

2:${ℂ}_{RC}\leftarrow \underset{c\in {ℂ}_{AC}}{\text{arg max}}H\left(\text{Δ}RC\left(c,t^{*},r\right)-T_{RC}\right)$

3: ${ℂ}_{SM}\leftarrow \underset{c\in {ℂ}_{RC}}{\text{arg k-max}}\text{ Δ}SM\left(c,t^{*},p\right)$

4: ${ℂ}_{TR}\leftarrow \underset{c\in {ℂ}_{RC}}{\text{arg k-max}}\text{ Δ}TR\left(c,t^{*}\right)$

5: ${ℂ}_{GE}\leftarrow \underset{c\in {ℂ}_{SM}\cup {ℂ}_{TR}}{\text{arg k-max}}\text{Δ}GE\left(c,t^{*}\right)$

6: \∗ get the class label having the highest frequency ∗\

7: ψ _q = mode $\left({ℂ}_{AC}\uplus {ℂ}_{RC}\uplus {ℂ}_{SM}\uplus {ℂ}_{TR}\uplus {ℂ}_{GE}\right)$

8: \∗ resolve the indecisive case step by step ∗\

9: if |ψ _q | ≥ 2 then

10:   ${ℂ}_{SM}\leftarrow \underset{c\in {ℂ}_{RC}}{\text{arg max}}\text{ Δ}SM\left(c,t^{*}\right)$

11:   ${ℂ}_{TR}\leftarrow \underset{c\in {ℂ}_{RC}}{\text{arg max}}\text{ Δ}TR\left(c,t^{*}\right)$

12:   ${\psi }_q^{\prime }=\text{mode(}{ℂ}_{SM}\uplus {ℂ}_{TR}\uplus {ℂ}_{GE}\text{)}$

13:    if ${\psi }_q^{\prime }\ge 2$ then

14:      $c_q={ℂ}_{GE}$

15:   else

16:      $c_q={\psi }_q^{\prime }$

17:   end if

18: else

19:   c _q = ψ _q

20: end if

Sometimes, for the query $s_q^{c_q}$, having subtle interactions with many classes, it is difficult for all GSF to come up with an unique agreement about c _q . When it happens, the threshold would have already hit the bottom of its range. Thus, with earlier reasoning, the solution would be to rely only on GE to find c _q , and 14 ^th step in the algorithm 2 does the same.

Deciding algorithmic parameters

Figure 5(b) shows the precision Vs threshold plots for SM and RC with different parameters. This helps to decide optimal parameters (p ^∗ , r ^∗ ) for them at the particular threshold. In the implementation, p ^∗ = 2 for all GP levels, while r ^∗ is set 3 for lower and 5 for higher GP levels. Other thresholds are set using the validation set by analyzing the maximum values for each GSF. And each set in the$\mathbb{T}$ has the uniformly quantized numbers from 0 to maximum GSF value.

Performance analysis

Figure 6 shows the normalized time taken by different classifiers for testing 14,086 sequences. In Majority voting scheme, first all GSF classify each sequence from the large pool of testing sequences, and then the voting begins. This slows down the scheme. On the other hand, in the proposed GP hierarchical scheme, the testing pool is gradually shrunken down. So the subsequent GSF have to investigate only small set of sequences, which likely to contain the true protein class. Which in turn speeds up the proposed algorithm along with maintaining high accuracy.

GP based modeling of protein families provides an extra advantage of fast incremental learning. In this procedure, $s_q^{c_q}$ is added back to cq after its classification. This is an instance based learning and it only requires slight modification to the graph of protein family c _q , like replacing G(V ^{c q} , E ^{c q} ) with $G\left(V_q^{c_q},E_q^{c_q}\right)$. Experiment consists of randomly selecting various (10 to 90) percent data from each $c\in {ℂ}_M$ for training purpose. This procedure is repeated 5 times and the Figure 7(b) shows the average classification performance comparison between incremental and batch learning. For small amounts of training data, batch learning performs poorly. On the other hand, incremental learning takes advantage of each correctly classified $s_q^{c_q}$ and produces better performance even during the scarcity of the training data.