Enhancing HMM-based protein profile-profile alignment with structural features and evolutionary coupling information

Background Protein sequence profile-profile alignment is an important approach to recognizing remote homologs and generating accurate pairwise alignments. It plays an important role in protein sequence database search, protein structure prediction, protein function prediction, and phylogenetic analysis. Results In this work, we integrate predicted solvent accessibility, torsion angles and evolutionary residue coupling information with the pairwise Hidden Markov Model (HMM) based profile alignment method to improve profile-profile alignments. The evaluation results demonstrate that adding predicted relative solvent accessibility and torsion angle information improves the accuracy of profile-profile alignments. The evolutionary residue coupling information is helpful in some cases, but its contribution to the improvement is not consistent. Conclusion Incorporating the new structural information such as predicted solvent accessibility and torsion angles into the profile-profile alignment is a useful way to improve pairwise profile-profile alignment methods.


Background
Pairwise protein sequence alignment methods have been essential tools for many important bioinformatics tasks, such as sequence database search, homology recognition, protein structure prediction and protein function prediction [1][2][3][4][5]. Following the development of global and local alignment methods of aligning two single sequences [6][7][8], profile-sequence alignment or profile-profile alignment methods such as PSI-BLAST, SAM [9], HMMer [10], HHsearch, HHsuite [4][5][6], which enrich two single sequences with their homologous sequences, has substantially improved both the sensitivity of recognizing remote homologs and the accuracy of aligning two protein sequences.
Due to their relatively high sensitivity in recognizing remote protein homologs, profile-profile alignment methods have become the default structural template identification method for many template-based protein structure modeling methods and servers [11][12][13][14]. For instance, HHsearch, one of top profile-profile alignment tools based on comparing the profile hidden Markov models (HMM) of two proteins, was used by almost all the template-based protein structure prediction methods tested during the last two Critical Assessment of Techniques for Protein Structure Prediction (CASP) [15,16]. The open source package HHsuite contains both the latest implementation of HHSearch that supports a full HMM-HMM alignment-based search on a HMM profile database and a very fast search tool HHblits [5] that reduces the number of unnecessary full HMM pairwise alignment in order to drastically improve its search speed. Moreover, the maximum accuracy (MAC) alignment algorithm is applied in HHsuite, but not in HHsearch. In this work, we aim to introduce new sources of information to improve profile-profile alignments with respect to both the original HHsearch package and the open source HHsuite package, In order to more accurately align the structurally equivalent residues in a target protein and a template protein together, secondary structure information was incorporated into profile-profile sequence alignment methods, yielding the better sensitivity and accuracy [4,17]. Aiming to find the new source of information to further improve the sensitivity and accuracy of pairwise profile-profile alignment, we examine the effectiveness of incorporating into profileprofile alignment methods some new features that have not been used in profile-profile alignments before, including protein solvent accessibility, torsion angles, and the evolutionary residue coupling information [18,19].
Specifically, we add the additional scoring terms for solvent accessibility, torsion angles, and evolutionary residue coupling information into the scoring function of HHsuite [5] in order to enhance the alignment process. According to our evaluation, adding solvent accessibility and torsion angles can improve the alignment accuracy, but incorporating the evolutionary residue coupling information is only useful in some cases.

Methods
We extended an existing profile-profile alignment method within the standard five-step alignment framework of HHsuite [5] shown in Figure 1, including discretization of profile columns, removal of very short or very dissimilar sequences, execution of Viterbi alignment and calculation of E-value and probability, realignment based on the maximum accuracy (MAC) algorithm, and retrieval of alignments by tracing-back. Different from HHsuite, our method applies solvent accessibility and torsion angle information to both the Viterbi alignment and the maximum accuracy alignment, and traces back with the aid of the evolutionary residue coupling information. In the following sections, we focus on describing how to incorporate the new features into the profile-profile method (i.e., HHsuite), while briefly introducing the necessary technical background.

Adding solvent accessibilities and torsion angles into the viterbi alignment
The score of aligning two columns in two protein profiles (namely a query profile q and a template profile t) in HHsuite was calculated according to Equation (1).
in which q i (a) and t j (a) denote the probability of amino acid at position i in the query profile and at position j the template profile, respectively, and f(a) is the background frequency of residue a (a ∈ {1, 2,...., 20}, representing 20 types of amino acids). The best alignment between two profile HMMs was obtained by maximizing the log-sum-odds score S LSO according to Equation (2).
where k denotes the index of columns that query HMM q aligned to template HMM t, i(k) and j(k) are the respective columns in q and t, P tr is the product of all transition probabilities for the path through q and t. The latest version of HHsuite has included the secondary structure information into the calculation of the score. In this work, we further augment the calculation of the score by adding the terms to account for the solvent accessibility, and torsion angles. The Viterbi dynamic program algorithm used five matrices S AB (i.e., AB ∈ {MM, MI, IM, DG, GD}) representing matching different states (M: match, I: insertion, Figure 1 The workflow of the HMM-based profile-profile pairwise alignment. D: deletion; G: Gap [4]) in two HMMs to maximize the augmented log-sum-of-odds score S LSO . They are recursively calculated as: w ss ; w sa ; w tors ∈ 0; 1 ð Þ S IM (i, j) and S GD (i, j) are calculated similarly as S MI (i, j) and S DG (i, j).
The difference between Equation (3) above and the default one in HHsuite is that two new terms (S sa , S tors ) were added to utilize the solvent accessibility and torsion angle information. In Equation (3), S ss (q i , t j ) is the secondary structure score between column i in query HMM (q i ) and column j in template HMM (t j ), which was the same as the one originally used in HHsuite. S sa (q i , t j ) is the solvent accessibility score between q i and t j , and S tors (q i , t j ) is the torsion angle score between q i and t j , which are the new terms introduced in this work. w ss , w sa , and w tors are weights for the secondary structure score, solvent accessibility score and torsion angle score respectively. S shift is the score offset for match-match states. Three weights w ss , w sa , w tors and shift score S shift are set to 0.11, 0.72, 0.4 and −0.03 by default, and can be adjusted by users as well. q i − 1 (M, M) is the transition probability from state M at column i-1 to next state M of in the query HMM, and t j − 1 (M, M) is the transition probability from state M at column j-1 to next state M in the template HMM.
Here we denote this extension of the HHsuite method as HMMsato. HMMsato allows for scoring predicted (or known) solvent accessibilities of one protein against predicted (or known) ones of another protein. DSSP [20] is used to parse the true solvent accessibility of a protein if its tertiary structure is known. PSpro 2.0 [21] is used to predict the solvent accessibility of a protein. The solvent accessibility information can be automatically parsed or predicted in HMMsato, or alternatively provided by a user. The two types of solvent accessibilities (e: exposed, > = 25% of the maximum area of a residue is exposed; b: buried, < 25% of the maximum area of a residue is exposed) are employed. Assuming the predicted or true solvent accessibility states of the i th residue (x i ) of the query protein and the j th residue (y j ) of the template protein are sa(x i ) and sa (y j ), the solvent accessibility score between the two residues S sa (q i , t j ) is defined as: The score is calculated by the kronecker-delta function δ(a, b), which equals to 1 if a = b, 0 otherwise.
Similarly as the solvent accessibility, the torsion angles including both phi angle (φ) and psi angle (ψ) can be automatically predicted by SPINE-X [22,23] or provided by a user. The range of both φ and ψ is (−180,180). Given the query sequence X and template sequence Y, the predicted phi angle and psi angle of the i-th residue x i in the query are denoted as φ(x i ) and ψ(x i ), and those of the j-th residue y j in the template as φ(y j ) and ψ(y j ). The torsion angle score S tors (q i , t j ) between the two residues is calculated as: Realign the profiles by maximum accuracy alignment combining solvent accessibility and torsion angles It has been shown that maximum accuracy (MAC) algorithm can generally create a more accurate alignment than the Viterbi algorithm, while the latter can generate better alignment scores, e-values and probabilities [5,24]. Consequently, the Viterbi algorithm is applied to compute e-values and scores, and the MAC algorithm is chosen to generate the final HMM-HMM pairwise alignment in HHsato by default.
The maximum accuracy algorithm [5,24] creates the local alignment that maximizes the sum of probabilities for each residue pair to be aligned minus a penalty (mact) (i.e., argmax( represents the posterior probability of the match state i in HMM q aligned to the match state j in HMM t. With the parameter mact, users can control the alignment greediness, from nearly global, long alignment (mact = 0) to very precise, short local alignments (mact ≈ 1). The default value of mact is set to 0.3501 in HMMsato as in HHsuite. To find the best MAC alignment path, an optimal sub-alignment score matrix AS is calculated recursively using the posterior probability P q M i e t M j Þ as substitution scores: Here, the Forward-Backward algorithm in local or global mode is applied to calculate the posterior probabilities . The Forward partition function F MM (i, j) and Backward partition function B MM (i, j) are introduced to calculate the posterior probability for pair state ( q M i ; t M j ) according to Equation (9): Five dynamic programming matrices F AB are used to compute the Forward partition function F MM , and AB ∈ {MM, MI, IM, DG, GD}. The top row and left column of the F MM matrix were initialized to 0, and all the matrices were filled recursively: where p min controls the alignment model (0: global alignment mode, 1: local alignment mode). F IM (i, j) and F GD (i, j) are calculated similarly as F MI (i, j) and F DG (i, j). Solvent accessibility score S sa (q i , t j ) and torsion angle score S tors (q i , t j ) are calculated as in the Viterbi alignment.
In analogy to the Forward partition function, the Backward partition function matrix B MM are calculated recursively as follows: Figure 2 Tracing back from the AS matrix by integrating the evolutionary coupling information. In query q, the coupled position of i is k q (i) , and that of i-1 is k q (i − 1). In template t, the coupled position of j is k t (j) , and that of j-1 is k t (j − 1). M q (i) is the corresponding position in template t matched to position i in q during the original tracing-back. M t (j) is the corresponding position in query q matched to position j in t during the original tracing-back. Additional EC scores are added into the corresponding elements in the AS matrix as shown in the figure so that the correct tracing back is performed.  ) and B GD (i, j) are calculated similarly as B MI (i, j) and B DG (i, j).

Trace back maximum accuracy alignments with the evolutionary residue coupling information
The Evolutionary Coupling (EC) stands for the correlation between two positions or columns in a multiple protein sequence alignment or a protein profile [19,20]. It has recently been employed to predict residue-residue contacts [18,19]. In order to improve profile-profile alignment with the evolutionary coupling information, we calculate the mutual information (MI) (one way of calculating EC value) for any two columns (i, j) of each profile according to Equation (12).
N is 21, standing for 20 amino acids plus gap. The joint probability of two residues X i and X j (F ij (X i , X j )) and the probability of residue X i (F i (X i )) are calculated in the same way as in [10]. However, EC ij is calculated as the mutual information (MI) instead of the direct information (DI) based on the global probability model [19] in order to achieve the higher time efficiency. A higher EC value corresponds to a stronger correlation between two columns in the given profile.
Based on the calculated EC value matrices for both the query and template profiles, top highly correlated position pairs with higher EC values for each profile are selected. The evolutionary residue coupling information is then applied to check the counterpart pairs during the process of tracing back through the sub-alignment score matrix AS (see Equation (8)) of the MAC alignment. Specifically, we denote the evolutionary coupled position for position i in query q as k q (i), and the coupled position of position j in template t as k t (j). Moreover, M q (i) denotes the position in template t matched with position i in query q when tracing back the original AS matrix, M t (j) denotes the position in query q matched with  position j in template t when tracing back the original AS matrix, and w ec is the weight for the evolutionary coupling information. The new AS ' matrix integrating the evolutionary coupling information is recalculated as follows during the track back process. Figure 2 illustrates an exampling of taking into account the evolutionary coupling information during the tracing back process to generate the final alignment.

Evaluation data set and metric
We evaluated HMMsato along with HHSearch [4] and HHsuite on the alignments between 106 targets (queries) of the 9 th Critical Assessment of Techniques for Protein Structure Prediction (CASP9) [15,16] and their homologous template proteins (templates) released at the CASP9's web site. The alignment data set has 2,621 pairs of query and template proteins. 1,483 pairs associated with 60 CASP9 targets were used as optimization data set to Bold denotes the two best scores, and an extra superscript of star denotes the highest score. Figure 3 The plot of the SP and TC scores against different values of the weight of solvent accessibility (w sa ).
optimize the parameters of HMMsato, and 1,138 pairs associated with the remaining 46 CASP9 targets were used to test the methods. The reference (presumably true) pairwise alignments of a query-template protein pair was generated by using TMalign [25] to align the tertiary (3D) structures of the two proteins together. The alignments generated by HMMsato and other tools were evaluated by three metrics, including sum-of-pairs (SP) score, true column (TC) score, and the quality of the tertiary structural models of the query proteins built from the alignments. The SP and TC scores are the two standard metrics for evaluating sequence alignment quality [26]. The quality of tertiary structural models indirectly assesses the quality of sequence alignments according to their effectiveness in guiding the construction of protein structural models. The SP score is the number of correctly aligned pairs of residue in the predicted alignment divided by the total number of aligned pairs of residues in the core blocks (i.e., sequence alignment regions precisely determined by structural alignment of structurally equivalent residues in the structures of two proteins) of the true alignment [23]. The TC score is the number of correctly aligned Bold denotes the two best scores, and an extra superscript of star denotes the highest score. columns in the core blocks of the true alignment [27]. The 3D model of a query protein was produced by MODELLER [28] based on both the pairwise alignment generated by an alignment method and the known structure of the template protein in the alignment. We used TM-Score [29] to align a 3D model of a query protein against its true structure to generate TM-scores and GDT-TS scores [30] for the model in order to measure the quality of the alignment used to generate the model, assuming better alignments lead to better 3D models with higher TM-scores and GDT-TS scores. Both TM-score and GDT-TS score are in the range [0, 1] [31].

Optimization of weights for the solvent accessibility, torsion angles and evolutionary coupling information
We estimated the weights of the solvent accessibility, torsion angles and evolutionary residue coupling information on the training alignments step by step. Firstly, we found the best weight value (w sa = 0.72) for solvent accessibility. Then, we identified the best weight value (w tors = 0.4) for torsion angles while keeping the weight for solvent accessibility fixed. Finally, we found the best parameter value (w ec = 0.1) for the evolutionary residue coupling information by keeping w sa and w tors at their optimum values. HHsearch and HHsuite were both evaluated with and without secondary structure information. The default parameter values were used with HHsearch and HHsuite.
Comparison of HMMsato, HHSearch, and HHsuite on the test data set The mean SP and TC scores for the pairwise alignment results generated by HMMsato, HHSearch and HHsuite for 1,138 protein pairs are reported in Table 1. The mean SP score and the mean TC score of HMMsato are 50.39 and 50.02 respectively, higher than HHsearch and HHsuite with or without secondary structure information. The average TM-scores and GDT-TS scores of the 3D models successfully generated from 1,127 out of 1,138 alignments by MODELLER were listed in Table 2. The average TM-score and GDT-TS score of the models generated from the HMMsato alignments are 0.555 and 0.483, respectively, better than those of HHSearch and HHsuite. Furthermore, we carried out the Wilcoxon matched-pair signed-rank test on both SP and TC scores of the three methods on the Figure 5 The plot of the SP score differences between HMMsato and HHsearch with secondary structure (HMMsearch-SS) for all the 1138 testing pairs. X-axis represents the index of the testing pair (1-1138), and y-ray represents the SP score difference (the SP score of HMMsatothe SP score of HHSearch-SS) for all the testing pairs. test data set. The p-values of alignment score differences between HMMsato and the other methods calculated by the Wilcoxon matched-pair signed-rank test are reported in Table 3.
Impact of solvent accessibility, torsion angles and evolutionary coupling information on the alignment accuracy We studied the effect of the solvent accessibility information by solely adjusting the value of its weight w sa . The SP scores and TC scores of the alignments generated by HMMsato with different w sa values on the training data set are shown in Table 4. The results show that incorporating the solvent accessibility information always improves alignment accuracy in comparison with the baseline not using solvent accessibility information (w sa = 0). The highest accuracy is achieved when w sa is set to 0.72. Figure 3 shows the plot of SP scores/TC scores against the different values of w sa. Red curve represents the SP scores and blue represents the TC scores.
We studied the effect of torsion angles on alignments by solely adjusting the value of w tors (weight for torsion angle information) while keeping w sa as 0.72. The SP scores and TC scores of the alignments generated by HMMsato with different w tors values on the training data set are shown in Table 5. The results show that incorporating the torsion angle information also helps improve alignment accuracy. The highest accuracy is achieved when w tors is set to 0.4. Figures 4 shows the TM-scores and GDT-TS scores of the 3D models constructed from the alignments generated by HMMsato with both torsion angles and solvent accessibility with respect to different w tors values.
The effect of evolutionary residue coupling information on alignment accuracy We studied the effect of the evolutionary residue coupling information on alignment accuracy in a similar way. HMMsato worked the best when w ec was 0.1. However, the evolutionary coupling information did not improve the overall alignment accuracy on the training data set, probably due to lack of a large number of diverse sequences in many cases required by the evolutionary coupling calculation to obtain the sufficient discriminative power. Specifically speaking, the alignment quality increased in 57 alignments, stayed the same in 1363 alignments, but decreased in 61 alignments. Similarly, on the test data set, the alignment quality increased in 59 alignments, stayed the same in 1024 alignments, but decreased in 55 alignments. Generally speaking, the evolutionary coupling information contributed to the improvement of alignment accuracy in some cases, but its effect was rather inconsistent. Figure 6 The plot of the average SP score difference between HMMsato and HHSearch-SS for the 46 testing protein targets. X-axis represents the index of the testing targets (1-46), and y-axis represents the score difference.

Comparison of HMMsato and HHSearch with secondary structure information on the test data set
We studied the SP score differences between HMMsato and HHSearch with secondary structure for all the 1138 testing pairs. The plot of the SP score difference (SP score of HMMsato minus SP score of HHSearch) for these pairs is shown in Figure 5. Similarly, the plot of the average SP score difference between HMMsato and HHSearch-SS for the 46 testing protein targets is shown in Figure 6. X-axis represents the index of the testing targets (1-46), and y-axis represents the score difference. Specifically, the alignment quality increased for 24 targets, stayed the same for 2 targets, but decreased for 20 targets. We found that HMMsato often improved the alignment quality for proteins of length ranging from 70 to 450 residues.

Conclusion
We designed a method to incorporate relative solvent accessibility, torsion angles and evolutionary residue coupling information into HMM-based pairwise profile-profile protein alignments. Our experiments on the large CASP9 alignment data set showed that utilizing solvent accessibility and torsion angles improved the accuracy of HMMbased pairwise profile-profile alignments. However, the effect of the evolutionary residue coupling information on alignments is less consistent according to our current experimental setting, even though it may still be a valuable source of information to explore in the future. Particularly, we will use the latest method (i.e., direct information) of calculating evolutionary coupling information to guide the profile alignment process. Furthermore, we will carry out more extensive search of optimal weights for solvent accessibility, torsion angle, secondary structure, and evolutionary coupling information to improve alignment accuracy.