Performance of Hidden Markov Models in Recovering the Standard Classification of Glycoside Hydrolases

Introduction

Glycoside hydrolases (GHs) are Carbohydrate-Active enZymes (CAZymes) that assist the breakdown of glycoside bonds of complex sugars into carbohydrates. ¹ They are extensively studied because of their role in the degradation of cellulose, hemicellulose, and lignin. ² Over the past 2 decades, studies have been conducted aiming to organize GH diversity^3–7 based on their mechanism of action, ⁸ enzyme structure,^7–10 type of substrates, and the domains of life in which GHs occur.^11,12 Glycoside hydrolases were first classified based on substrate specificity and the type of catalytic reaction. ¹³ However, these classifications failed to account for both structural and evolutionary features. Therefore, a new classification scheme was proposed by Henrissat, ³ which is currently widely used and may be regarded as the standard GH classification. The Henrissat classification assigned GH to families based on the similarity of amino acid sequences and, consequently, on their secondary structure, combining methods of alignment and hydrophobic cluster analysis.^3,14,15 Glycoside hydrolase families were later compiled in the CAZy database according to the Henrissat classification, which currently encompasses >130 families. ⁶ Since then, this classification of CAZymes has been used in studies dealing with a diverse array of subjects, from the structure, specificity, and efficiency of GHs to the analysis of metagenomes.^10,16–18

Currently, an automatic method of assigning sequences to GH groups is required ¹⁹ because of the huge generation of new data from genome projects via high-throughput sequencing technologies. The CAZy database employs a semi-automatic annotation using Pfam hidden Markov model (HMM) profiles ²⁰ and BLAST (Basic Local Alignment Search Tool), ²¹ which are posteriorly curated by experts manually. ¹ Although this procedure has generated a widely implemented GH classification, it is difficult to reproduce automatically. Several CAZy GH families are absent in databases such as Pfam, which makes difficult a cross-correspondence between protein databases. ²⁰

Recently, statistical methods based on HMMs have shown good performance in homology inference of amino acid sequences and have been widely used in evolutionary studies of protein families.^22–24 Thus, HMM may result in a viable alternative to the semi-automated approach of CAZy. As a consequence of the HMM method, a profile consisting of amino acid variation probabilities along alignment sites is inferred, indicating possible insertions and deletions. ²⁵ In addition, HMM profiles indicate conserved and nonconserved positions within members of the same family. ²² An extensive collection of annotated HMM profiles is currently available in protein databases, such as Pfam, ²⁰ which use this statistical approach for improvement of the classification criteria used to group protein families into clans. ²⁰ Families belonging to the same clan share similarities in amino acid sequences, structure, and HMM profile. ²⁰

An automatic classification of GH families based on the CAZy database was proposed by Yin et al ²⁶ and is currently available in the dbCAN database. For each CAZy GH family, dbCAN identifies signature domains and builds an HMM model to represent them. ²⁶ However, dbCAN does not contain all GH families, and no overall evaluation has ever been conducted of the performance of the HMM method in recovering the standard GH classification implemented in CAZy.

In this study, we conducted the first comparative meta-analysis of the degree of correspondence between HMM-inferred profiles and their corresponding CAZy database classifications, which is largely regarded as the standard GH classification. When HMM profiles correctly recovered the GH family classification, we considered this family to be automated using this statistical approach. We showed that HMM profiles recovered the CAZy family classification well. Although mismatches between HMM and CAZy were found, they were not random, and errors could be systematically associated with specific GH families. This indicates that further developments in the automatic assignment of amino acid sequences to CAZy protein families may be targeted at those particular cases, allowing for a consistent classification of GH sequences among protein databases.

Methods

We downloaded, from GenBank, approximately 300 000 GH sequences classified by the CAZy database as of January 2017. This data set consisted of all available GH sequences from cellular organisms, including 5469 sequences from Archaea, 240 650 sequences from Bacteria, and 55 080 sequences from Eukarya. A total of 127 out of the 135 GH families listed in CAZy were studied. Eight families were excluded from the analysis because they were absent in cellular organisms. For the sake of comparison, we also analyzed the performance of HMM models readily available in the dbCAN database. ²⁶

To evaluate the recovery rate of the standard CAZy GH classification using HMM profiles, multiple alignments were conducted for amino acid sequences from each CAZy GH family, using the Clustal Omega algorithm. ²⁷ The quality of GH alignments was evaluated using the transitive consistency score (TCS) statistics available in T-Coffee platform. ²⁸ As TCS assigns a score for each site in an alignment, it allows verification of the overall alignment’s quality. We were then able to check whether poor-quality alignments affected the inference of HMM profiles. The TCS of more than 500 was regarded as well supported. ²⁸ To evaluate whether the frequency of gaps affected the inference of HMM profiles, the proportion of gaps in each alignment was calculated using Perl scripts. The TCS and gap proportions were estimated for both the alignments of complete sequences from CAZy GH families and the alignments of family domains from dbCAN.

The resulting 127 individual CAZy GH family alignments were then loaded into the HMMER 3.1 platform ²⁹ to build their respective HMM profiles using the hmmbuild software. The inferred CAZy-based HMM profiles were posteriorly combined into a single HMM database using the hmmpress software. This database of HMM profiles was hereafter referred to as HMM-CAZy. Finally, all GH amino acid sequences from cellular organisms available in CAZy, regardless of their CAZy-based classification, were queried against the complete HMM-CAZy database using the hmmscan software. We expected, thus, a one-to-one match between the CAZy GH classifications and their correspondent HMM profiles. Statistically, this was a measure of the power of the HMM method in recovering the CAZy classification. We also downloaded all available individual HMM profiles from dbCAN and built a single database of profiles, hereafter referred to as HMM-dbCAN, using the same strategy employed with the HMM-CAZy. We thus inferred the power of the HMM profiles built solely from protein domains in recovering the CAZy classification. In all steps, default HMMER parameters were used.

If an amino acid sequence assigned by CAZy to a given GH family matched the corresponding HMM profile, built using sequences from the same CAZy-defined GH family, it was counted as a correct match for this GH family because the HMM assignment recovered the standard GH classification. However, a mismatch was computed when an amino acid sequence, defined by CAZy as a member of a given GH family, was not assigned to the respective HMM profile estimated for that same family. If a mismatch was found, we registered the HMM profile with the best hit for that sequence. We summarized the frequency of matches as well as mismatches for each family independently. Information on the steps of the HMM analysis is detailed in Figure 1.

OPEN IN VIEWER

Figure 1.

Steps used to evaluate the performance of retrieving CAZy GH families’ classifications using HMM. GH indicates glycoside hydrolase; HMM, hidden Markov model.

Results and Discussion

When GH amino acid sequences were queried against the HMM-CAZy and HMM-dbCAN databases, positive matches between CAZy classifications and HMM profiles were recovered for 65% and 61% of GH families, respectively. Mismatches between the standard classification and HMM profiles were thus found in 35% and 39% of CAZy families when queried against the profiles of the HMM-CAZy and HMM-dbCAN databases, respectively (Figures 2 and 3, Supplementary Tables S1 and S2). Furthermore, in 24% of cases, identical mismatches were found when analyzing both HMM databases (HMM-CAZy and HMM-dbCAN) (Supplementary Table S3).

OPEN IN VIEWER

Figure 2.

Frequency of errors, as shown in the left pie chart, and discrimination of GH families that presented incongruences between CAZy GH classification and HMM profiles inferred in this study. For each CAZy family, horizontal bars indicate the percentage of mismatches between the CAZy-defined family and the HMM profile. GH indicates glycoside hydrolase; HMM, hidden Markov model.

OPEN IN VIEWER

Figure 3.

Frequency of errors, as shown in the left pie chart, and discrimination of GH families that presented incongruences between CAZy GH classification and HMM profiles downloaded from the dbCAN database. For each CAZy family, horizontal bars indicate the percentage of mismatches between the CAZy-defined family and the HMM profile from the dbCAN database. GH indicates glycoside hydrolase; HMM, hidden Markov model.

Mismatches between the CAZy protein family classification and the retrieved HMM profiles consisted of 2 major types (Figure 4). In the first type, several amino acid sequences from a given CAZy GH family were consistently assigned to an HMM profile estimated for a different CAZy GH family (Figures 2 and 3). For these cases, the greatest discrepancies were observed in GH families 73 and 81, in which the queries against the HMM-CAZy database resulted in a high frequency of incorrect matches with profiles inferred for families 70 and 80, respectively (Figure 2). When we analyzed the HMM-dbCAN database, a high frequency of inconsistencies were recovered for GH families 80 and 134, which resulted in incorrect matches with profiles inferred for families 81 and 135, respectively (Figure 3). We named this type of mismatch Error 1, which accounted for 38% of the mismatches found when using the HMM-CAZy database and 41% when using the HMM-dbCAN database (Figures 2 to 4). This error is equivalent to a measure of statistical bias because amino acid sequences were consistently assigned to an incorrect HMM profile.

OPEN IN VIEWER

Figure 4.

Types of errors found between CAZy GH family classification and HMM profiles. GH indicates glycoside hydrolases; HMM, hidden Markov model.

We also noticed another recurrent type of mismatch, named Error 2, which was recovered when amino acid sequences from a given CAZy GH family were assigned, without any consistent association, to several HMM profiles that were all different from the CAZy classification for those sequences (Figure 4). The frequency of occurrence of Error 2 was higher when using the HMM-CAZy database compared with HMM-dbCAN (62% and 59%, respectively) (Figures 2 and 3). The CAZy families with the highest frequency of Error 2 in both HMM databases were families 5, 16, 18, and 43 (Figures 2 and 3). Amino acid sequences from these families all matched at least 9 HMM profiles that were different to their original CAZy classification (Figures 2 and 3).

Mismatches found in our study were possibly caused by the multimodular structure of GH sequences, which are frequently composed of several protein domains. ¹ Because we conducted HMM analyses based on complete GH sequences, an impact on the family assignment was expected when using HMM profiles built from CAZy families. This issue was not attenuated when we adopted the HMM-dbCAN profile database, in which independent HMM profiles were created using only alignments of GH domains. ²⁶ Estimation of HMM profiles from alignments of complete amino acid sequences resulted in a higher frequency of matches between the HMM search and the CAZy GH classification.

A better performance of our HMM profiles, built from complete sequence alignments, was recovered, despite their lower alignment scores when compared with GH domain sequences downloaded from dbCAN. The low alignment scores obtained for the HMMs from complete sequences were probably due to the higher percentage of gaps in the alignments (circa 1.6×) in comparison with dbCAN domain alignments. In general, the mismatches observed in both HMM data sets were obtained for GH families with the lowest alignment scores (Supplementary Table S4). We hypothesize that the classification criteria of sequence similarity, adopted by CAZy when defining GH families, generate sequence groups that are not evolutionarily related and, consequently, low-quality alignments are obtained. As dbCAN alignments were inferred from sequence domains only, they presented higher alignment scores and lower gap proportions (average of 2%). ²⁶ Nevertheless, the performance of dbCAN HMM profiles in recovering the CAZy classification was poorer.

We conclude that the automated strategy adopted here; that is, using HMM profiles to recover the standard GH classification implemented in CAZy was effective for most families. Moreover, although mismatches were found, they did not occur randomly, as error types were recurrent in specific families. This important finding means that a focus may thus be directed to these problematic GH families, aimed at obtaining a complete correspondence between CAZy and HMM-based methods. Our classificatory approach using HMM was a useful tool to automatize the CAZy classification of GHs; in principle, sequence assignments to GH families can be made objectively using CAZy-derived HMMs built from complete sequences. This is a significant step toward a higher congruence between the classifications adopted by distinct protein databases, which directly affects future studies on biochemical and evolutionary aspects of GHs.