Post-calibration finding 1: On average, HMMER3 E-values need to be more stringent than HMMER2 glocal-mode E-values to exhibit the same FPR
The function annotation task via the HMMER algorithm is innately coupled to the domain libraries. In this work, the Pfam library (release 29) with a total of 16295 domain models was examined. Furthermore, to capture the same search space of glocal-mode HMMER2 via the local-mode HMMER3, a calibration (or normalization) procedure is needed to relate the different E-values generated by either tool for the same sequences and domain models (See Methods on “Calibration procedure ..”).
Briefly, this entails the objective comparison of the receiver operating characteristic (ROC) curves generated by the two HMMER builds (i.e., HMMER2 and HMMER3) for each domain model. The query sequence set is specific to every domain model and it is composed of a positive set of α protein sequences (from the seed sequences of the Pfam domain model of interest, actually the potentially “true hits”) and the consensus sequences of the remaining 16294 domain models (if these sequences are hit, they are regarded as “false hits”).
Figure 1 depicts the false-positive rate (FPR) value against the median E-value of negative domain hits in Pfam library (release 29) for the two HMMER builds in a doubly logarithmic plot. The HMMER2 (in solid red) and HMMER3 (in solid blue) plots are shown with their interquartile ranges (IQRs) in dotted red and blue lines respectively. In a nutshell, the FPRs for each domain model are computed over 16294 negative domain consensus sequences of the Pfam library. For each domain, a list of i negative E-values are sorted in ascending order and labelled with its order from 1 to i. The order is in fact the FP value and the FPR range for this domain is, hence, between 1/16294 and i/16294.
Over the full body of the Pfam library, we found in our test calculations that the false-positive (FP) values span between 1 and 626 for HMMER2 and between 1 and 2143 for HMMER3. In turn, this translates to the sampled FPR ranges of between 6.14e-5 (1 out of 16294) and 3.84-2 (626 out of 16294) for HMMER2 and between 6.14e-5 and 1.32e-1 (2143 out of 16294) for HMMER3. Then, at each sample FPR or FP points, the median E-values of the negative domain hits of the Pfam library were taken to give the median E-value ranges of between 2.13e-9 to 2.40e + 1 for HMMER2 and between 1.35e-20 and 7.1e-2 for HMMER3. A reproducible feature to note for HMMER2 software suite is the occurrence of a breakpoint due to the switch in statistical model from the extreme value distribution (EVD) to the logistic model when the similarity score grows larger. In our previous work, the median E-value of Pfam domain (release 23) for this switch occurs at 1.0e-7 [10] which falls within the current (release 29) E-value breakpoint range of between 10−3 and 10−8.
Separately, the upper FPR limit of HMMER2 is lower than that of HMMER3 since it detects much less false-positive hits than HMMER3. This is unsurprising since HMMER2 enforces complete sequence-to-domain alignment when operating in glocal-mode, while HMMER3 does not (and the occurrence of a sequence similar to a fraction of the domain suffices). And for the intended purpose of capturing the search space by HMMER2, the common upper limit of FPR at 5.60e-3 (see label corresponding to 91 false-positives in Fig. 1) is theoretically sufficient for both HMMER variants, albeit at different E-value of 0.1 for HMMER2 and at 4.15e-7 for HMMER3.
In any case, an overall comparison of the two negative domain hits plots reveals that the median E-value needs to be more stringent in HMMER3 than in HMMER2 given a fixed FPR. The magnitude difference is about (106 ~ 1011) across the sampled range. This is consistent with the HMMER manual’s recommended E-value cutoff of < <1 for HMMER3 and <0.1 (between 0.1 and 1 requires manual checking) for HMMER2 to gather trusted hits. As an example, when one imposes the recommended trusted domain E-value cutoff of 0.1 for HMMER2, the corresponding FPR will be at 5.60e-3 (see dotted vertical line). The latter translates into 91 false-positive hits. Meanwhile, the equivalent HMMER3 E-value will be 4.15e-7 at this FPR setting. If extreme stringency is applied to achieve zero FP (i.e., FPR <6.14e-5 or FP <1 given our search database size), then the required HMMER2 and HMMER3 E-values will be at 2.13e-9 and at 1.35e-20 respectively (see leftmost dotted vertical line in Fig. 1).
However, complicated situations can arise as the HMMER2 and HMMER3 E-values enter the lower FPR of below 5.6e-3 (FP = 91) where the IQR of the two plots start to overlap. Essentially, this implies that there can be cases where the HMMER2 E-values can exhibit higher stringency than its HMMER3 equivalent, therefore bucking the general trend of HMMER3’s stringency over HMMER2. As such, a direct comparison of E-values from different HMMER builds even for comparable sequence-to-domain alignments is not always a straightforward exercise.
Post-calibration finding 2: About 20% of HMMER2 domain models have better sensitivity than the HMMER3 counterpart over all the FPR range
After the domain library calibration procedure, the HMMER builds with the better sensitivity can be evaluated for the domain library. Briefly, the area under curve (AUC; see Eq. 1) of the HMMER2 and HMMER3 models are first computed for a fixed FP value, x. These evaluated FP values, which correspond to respective FPRs (in parentheses), are x
1 = 1 (6.14e-5), x
2 = 5 (3.07e-4), x
3 = 10 (6.14e-4), x
4 = 20 (1.20e-3), x
5 = 30 (1.80e-3), x
6 = 40 (2.50e-3), x
7 = 50 (3.10e-3), x
8 = 60 (3.70e-3), x
9 = 70 (4.30e-3), x
10 = 80 (4.90e-3), x
11 = 90 (5.50e-3) and x
12 = 100 (6.14e-3). Furthermore, to enhance the stability of the model selection, a predefined set of FP values A (see Eq. 2) and its associated model counts (Eqs. 3 and 4) are evaluated and taken as the overall criteria for model selection (See Methods on “ Criteria for domain model build…” for details).
Figure 2a depicts the total number (see Eqs. 5–7) of HMMER2 (in red), HMMER3 (in blue) and undetermined (in green) model builds selected by the various sets of FP. The FP set starts with all the evaluated FP values of 1 to 100. As the plot moves along the horizontal axis from left to right, the smaller FP values are dropped progressively. Based on the HMMER2 plot (in red), a stringent set of 3066 HMMER2 models were found to be more sensitive than their HMMER3 counterparts for all evaluated FP sets. As the smaller FP values are dropped, more HMMER2 domains were added until a maximum of 3502 models have been reached. This is about 21.5% of the evaluated Pfam domain library (release 29). Meanwhile, a large proportion of HMMER3 domain models remain more sensitive than their HMMER2 counterparts. Quantitatively speaking, this is between 68.5% (11158/16295) and 78.3% (12757/16295) of the domain library (blue plot). A sizeable percentage of domain models (up to 12.7% of the Pfam library) has shown to be fluctuating between HMMER2 and HMMER3 depending on the evaluated FP value set (green plot) and hence, not deterministic for selection.
Figure 2b shows the normalized AUC difference (see Eq. 1) being expressed as percentage and averaged over the sets of selected HMMER2/HMMER3/undetermined domain models corresponding to Fig. 2a. Basically, a large difference implies higher sensitivity of either the HMMER2 or HMMER3 models. Specifically for the HMMER2 sets (red plot), the average AUC difference ranges between 1.61 and 3.59% while this is held steadily at about 0.3% (blue plot) for the HMMER3 sets. For the sets of undetermined models, the average AUC difference stabilizes to 0% as more low FP values were excluded (green plot) and are hence deemed to perform closer to the HMMER3 models.
Overall, the domain sensitivity evaluation revealed that about 20% (between 3066 and 3502; see Additional file 1 from xHMMER3x2 website) of the HMMER3 models were found unsuitable for the initial domain scanning in the xHMMER3x2 workflow. Interestingly, these 3502 models were biased towards shorter domain models with a median of 171 alignment positions (IQR = 170). The latter is suggestive of alignment quality issues when coupled to the local-mode HMMER3. In any case, these HMMER3 models cannot capture the HMMER2 search space sufficiently and hence, their HMMER2 model counterparts need to be used despite a longer computational speed. Fortunately, the slowdown caused by this 20% of the Pfam library is not detrimental; the slowest xHMMER3x2 will still outperform the standalone glocal-mode HMMER2 at least by an order of magnitude (see below, Results section 5).
In hindsight, larger numbers of domain hits (which can include false-positives at the same time) can generally be expected from local-style search mode (as in HMMER3) compared with the glocal mode regime (as in HMMER2) since alignments over the full domain model are not enforced in the first case. Fortunately, HMMER3 provides larger number of hits and, as a trend, better sensitivity than HMMER2. However, on the flip side, HMMER3’s tendency towards partial domain coverage is insufficient to justify for function annotation transfer since a full domain necessarily implies a unit of function. The extent of HMMER3’s partial coverage over various FPR values is quantified in the next section.
Post-calibration finding 3: Average domain coverage by HMMER3 positive hits is about 88% of the domain model length
In a strict sense, partial sequence-to-domain alignments generally do not justify for a sequence segment having the function of the domain hit since the completeness of a 3D structural unit (as implied by the domain unit) is a precondition for functional annotation transfer. The latter is especially relevant to HMMER3 since it can only operate in local-mode.
To study the completeness of domain coverage by HMMER3, we first computed the average “coverage per domain model” by taking the total coverage of detected positive hits divided by the total detected positive hits at a FPR for a given domain model. Finally, an overall average is taken over all 16295 HMMER3-build Pfam models. Figure 3a depicts across the sampled FPR range of between 6.14e-5 and 6.14e-3 (in logarithmic scale) together with the standard error bars at each evaluated FPR value. Briefly, the “coverage per domain model” gives an idea of how much the sequence-to-domain alignment hits cover the domain model on average. If the domain coverage of a model is complete in all cases, then the average is expected to be one. Across the FPR range, the average domain coverage of the Pfam domain models with HMMER3 straddles at about 0.88 (see horizontal dotted line) with the lower limit (median-1.5IQR) of the average domain coverage at about 0.85.
Although the average “coverage per domain model” remains relatively stable across all FPR values, incomplete domain coverage is expected to have different impact especially in cases of relatively long and of the shorter domain models. Figure 3b gives the distribution of domain length of the Pfam release 29 library. There are a total of 1801 domains with length of 485 alignment positions or more (i.e., >median + 1.5IQR; see vertical line in Fig. 3b). Therefore, the lack of coverage constitutes at least 58 ((1–0.88) × 484.5) missing alignment positions. The latter is about the size of the smallest known domain (e.g., zinc finger models); thus, completeness of the respective domain range in the query sequence is at least doubtful. On the other hand, there are 1429 with domains length of 58 alignment positions and less, where annotation transfers are even more unjustified under incomplete coverage. Taken together, these 3230 (1801 + 1429) implicated models make up almost 20% of the Pfam library. Therefore, despite HMMER3’s generally superior sensitivity, it suffices only for initial fast domain detection. Ultimately, HMMER2 will still be required for re-alignment over the full domain model, especially for the shorter ones.
The xHMMER3x2 workflow and webserver implementation
Figure 4 gives an overview of the three-staged xHMMER3x2 workflow. Briefly, xHMMER3x2 first invokes HMMER3 to scan for initial domain hits over a set of query sequences that matches the search space of HMMER2. The hits are then carried over to HMMER2 for refinement into full sequence-to-domain alignments.
In the first stage of xHMMER3x2, N (3066 ~ 3502) HMMER2 and M (13229 ~ 12793) HMMER3 models can be utilized for the search (where N < <M) based on the sensitivity evaluation between HMMER2 and HMMER3 models in Fig. 2a. The upper E-value limit is set to 0.1 for HMMER3 and 24 for HMMER2 which corresponds to the maximum sampled FP values of 2143 and 626 respectively according to Fig. 1. Essentially, the smaller search space of HMMER2 should be sufficiently captured.
Generally speaking, since the total number of utilized HMMER3 models far exceeds that of HMMER2 (80% versus 20%), the speed of xHMMER3x2 is expected to be closer to the faster HMMER3 than the slower glocal-mode HMMER2. As such, reasonable computational speed can be achieved without sacrificing for detection sensitivity but with the caveat that the HMMER3 hits might be fragmented to various extents. Upon the selection of the appropriate model builds, the respective hmmpfam (HMMER2) and hmmscan (HMMER3) runs are executed with the respective sub-libraries of HMMER2 and HMMER3 domain models.
In the second stage, the two runs from the dominant HMMER3 and the minority HMMER2 will result in K local-mode and L glocal-mode sequence-to-domain alignments respectively for a given query sequence. The resulting K local-mode alignments will be searched against the same K domain models but in glocal-mode HMMER2. The latter is considered as the full domain re-alignment step as intended by glocal-mode HMMER2 originally. Since the number of models used in the realignment step is expected to be much less than the total number of HMMER3 models used (i.e., K < <N < <M), this is an area of major speed gain by xHMMER3x2 over pure HMMER2. At the end of stage two, xHMMER3x2 would have produced K pairs of HMMER2/HMMER3 and L HMMER2 sequence-to-domain hits over K + L domain models for every query sequence.
The third stage of xHMMER3x2 estimates the FPR (and FPs out of 16294) values of the domain hits based on their E-values. For each domain model, the previously calibrated HMMER2/HMMER3 ROCs are purposed to assign the particular domain hits to the appropriate FPR interval based on their E-values. The intervals are bounded by the evaluated FPR (FP) points: 6.14e-5 (1), 3.07e-4 (5), 6.14e-4 (10), 9.21e-4 (15), 1.20e-3 (20), 1.50e-3 (25), 1.80e-3 (30), 2.50e-3 (40), 3.10e-3 (50), 3.70e-3 (60), 4.30e-3 (70), 4.90e-3 (80), 5.50e-3 (90) and 6.14e-3 (100) as depicted in Fig. 4. The estimated FPR (FP) values of the domain hits are simply the nearest margin of their assigned FPR interval.
The overall xHMMER3x2 workflow is implemented as a software suite in PERL language. It is freely available for download at http://xhmmer3x2.bii.a-star.edu.sg for local installation. Since xHMMER3x2 requires the carefully-selected domain models for its proper performance, the list of calibrated models for the recent Pfam domain libraries (release 27, 28 and 29) are also provided at the same site. Alternatively, an online webserver (for up to 1000 FASTA-formatted sequences) exists at the download site as well. When multiple sequences are inputted, the result is returned via email in tabulated form as electronically readable file containing query sequences, Pfam model hits together with their alignments E-values and standardized FPRs. Upon input of a single sequence, the HMMER2 and HMMER3 results are presented graphically-enhanced for easy comparison.
xHMMER3x2 improves on computational speed over HMMER2 and domain-detection sensitivity over HMMER3
To evaluate the real-world performance of xHMMER3x2, the sequences and Pfam domain annotations of the human UniProt/SwissProt sequences were evaluated. The source file was extracted from the public resource at ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.dat.gz.
To create a reference set for evaluation purpose, the sequences were re-annotated via the glocal-mode HMMER2 against the Pfam domain library (release 29). Only the re-annotated domains that overlap with the existing UniProt/SwissProt Pfam annotations were retained. The latter aims to exclude UniProt annotations that cannot be reproduced by this computation as a result of direct inheritance from previous releases, percolation errors, etc. Overall, our re-annotation resulted in 27878 (out of 27988 reported annotations) Pfam domain hits which span across 17517 sequences over 5719 unique domain models. Since the complete human proteome is a moving scientific target [23], the re-annotated set used here is only for benchmarking purpose only. The list of originally-reported and re-computed domain annotations for the human UniProt/SwissProt sequences is available at the xHMMER3x2 website (see Additional file 2). Furthermore, a cutoff E-value of 1e-9 (corresponds to the negative domain E-value for FP < 1 in Fig. 1) was applied to derive a stringent set of 24696 positives domains which now covers 16175 sequences over 5509 unique domain models.
To examine the trade-offs between computational speed and domain-detection sensitivity, these 16175 sequences were subjected to hmmscan/hmmpfam searches via HMMER3 (as reference), glocal-mode HMMER2 and 5 configurations of xHMMER3x2 (where N = 3502, N = top 1500/1000/500 and N = 0) against the Pfam library (release 29). The top most sensitive HMMER2 domain models were ranked according to the AUC evaluation via Eq. 1.
Figure 5 illustrates the domain-detection sensitivity and speed differences against HMMER3 (in log scale) for both HMMER2 and xHMMER3x2. The horizontal dotted line across the x-axis sets the reference where there is no speed difference from HMMER3. The leftmost boxplot depicts the most time-consuming HMMER2 at about 221 times (i.e., M
h2/M
h3 = 102.34; see left-hand side of Eq. 12) slower than HMMER3, but with the improved sensitivity of 4.3% (100% versus 95.7%) over HMMER3. This translates to capturing 1074 more hits than HMMER3. In hindsight, this finding corroborates with our earlier post-calibration finding that some HMMER3 domain models have lower sensitivity than HMMER2.
For the various setting of xHMMER3x2, the fastest is achieved by running all HMMER3 models (i.e., N = 0; see rightmost boxplot) which is only 1.1 times (i.e., 100.04) slower than HMMER3 (due to the necessary full domain re-alignment of glocal-mode HMMER2 in xHMMER3x2; see Fig. 4 Stage 2) but about 201 times (221÷1.1) faster than HMMER2. In hindsight, the re-alignment step of xHMMER3x2 does not impede its overall speed and the condition where K < <N < <M is generally true. However this is not without a caveat that its sensitivity has been compromised to 95.7% (same as HMMER3).
Meanwhile, the trade-offs between speed and sensitivity can be observed when the various sets of selected HMMER2 models were invoked by xHMMER3x2. Based on Fig. 5, as N (i.e., the number of default HMMER2 models) decreases from 3502,1500,1000 to 500 (see second to fifth boxplots from left), xHMMER3x2’s sensitivity decreases from 99.8, 99.3, 99.1 to 98.6% respectively. In contrast, its speed increases accordingly from 32× (M
h2/(M
h3 + K
h3) = 101.50; see right-hand side of Eq. 12), 17× (M
h2/(M
h3 + K
h3) = 101.23), 12× (M
h2/(M
h3 + K
h3) = 101.07) to 6× (M
h2/(M
h3 + K
h3) = 100.80) slower than HMMER3. Alternatively, this translates to 7 × (222÷32), 13 × (221÷17), 18.5 × (221÷12) and 37 × (221÷6.2) faster than HMMER2 respectively. Overall, the best compromise is achieved with the top 1000 most sensitive HMMER2 domains while the remaining as HMMER3 models. For most occasions, this configuration is the recommended and default mode in the xHMMER3x2 webserver since it achieves the best balance between domain-sensitivity (at 99.1%) and speed (12× slower than HMMER3 but 18.5× faster than HMMER2). However, when fast annotation is the main consideration, xHMMER3x2 offers the “glocal-mode HMMER3” as its fastest mode that is almost as fast as the native local-mode HMMER3.
In summary, the xHMMER3x2 framework will always produce the glocal-mode alignments and it can be configured to achieve a balance between speed and domain-detection sensitivity. On one extreme, it can be configured for full speed when all HMMER3 models are selected, thus mimicking a ‘glocal-mode’ HMMER3. On the other end, when all HMMER2 models are deployed, xHMMER3x2 reverts back to HMMER2 with the highest level of domain-detection sensitivity.
False-positive quantification of calibrated domain hits allows for direct comparison between HMMER variants to evaluate for annotation consistency
At the end stage of xHMMER3x2, the FPR (FP) values for all the sequence-to-domain hits would have been estimated accordingly. To evaluate the consistency of domain annotation, the 23622 common domain hits to both HMMER2 and HMMER3 from the preceding analyzed set of 16175 human UniProt/SwissProt sequences were further examined.
We define annotation consistency as the concordance in statistical assessment of same domain hits between HMMER2 and HMMER3 for a particular sequence segment. Figure 6a depicts the 23622 pairs of HMMER2/HMMER3 E-values (as log-log plot) that now covers 5492 unique domain models. The set of blue points denotes the case where the estimated FP values of the paired HMMER2/HMMER3 hits matched exactly, while the red points denote the unmatched cases.
Based on the linear regression analysis, the mathematical relationship of HMMER2 against HMMER3 E-values was found to be y = 0.72x, R
2
= 0.91 (matched cases) and y = 0.77x, R
2
= 0.22 (unmatched cases). In other words, the paired HMMER2/HMMER3 E-values are not equivalent, and non-linear to each other particularly for the unmatched cases given the bad R
2 (coefficient of determination) values. Additionally, about 86.5% of the pairs show higher stringency of HMMER2 E-values over HMMER3, and vice versa for the remaining pairs. To reiterate, direct comparison using E-values to determine concordance proves to be difficult. It is also important to remember that the evaluated Pfam annotations in UniProt are generally constrained to one domain hit per sequence segment. In more realistic settings, the same sequence segment can have multiple significant domain hits. In the absence of FPR estimation, the actual error rates associated to the overlapping domain hits will simply be masked behind a set of inseparable E-values.
Utilizing the estimated FP values of the paired HMMER2/HMMER3 hits, the difference (i.e., FPH2-FPH3) values were taken. The histogram of the 23622 pairs of FP difference values were shown in Fig. 6b. In contrast to the E-value plot, 95% (22428) of the data points converges to the singular point at FPH2 = FPH3 or 0 (see blue dot). When the small difference of FPH2-FPH3 = ±4 (see vertical dotted lines) are also considered since FP values are approximated to their nearest values, the overall statistical concordance of the HMMER pairs makes up about 98% (addition of another 766) of the data points. As such, FP quantification by xHMMER3x2 offers an easy and quick way to check for general concordance of the paired HMMER2/HMMER3 hits.
Unique to xHMMER3x2, its FPR quantification helps to highlight an annotation inconsistency based on an observed left-skewness in the histogram. More specifically, a handful ([FPH2-FPH3 < −4] = 389) data pairs exhibited higher HMMER3 FPR than that of HMMER2 despite similar sequence-to-domain alignments. These 389 outliers, which span across 140 unique Pfam domain models (see Additional file 3 from xHMMER3x2 website). Once again, the problematic domain models belong to the shorter ones; this set of 140 domains has a median length of 166.5 alignment positions (IQR = 204).
To resolve the annotation discrepancies, the sequences containing the 389 domain hits were subjected to a dissectHMMER [13–15, 24] analysis where pairs of HMMER E-values are dissected into their fold-critical and remnant contributions. Briefly, it is the fold-critical parts (the supposedly 3D structural part) of a sequence-to-domain alignment that argues for a similar overall fold, hence similar function whereas the remnant part represents disordered, fibrillary or other non-globular regions that might be not obligatory for the domain fold [13, 14, 24]. As such, the E-values of the fold-critical contributions has to be more significant than the remnant parts to justify for annotation transfer.
Figure 7 shows the boxplots of the total, fold-critical and remnant HMMER2 and HMMER3 E-values of the 389 dissected domain hits. As a baseline for significance, the dotted horizontal line (at -1) denotes the HMMER2’s recommended E-value of 0.1. The first two boxplots on the left shows that the total HMMER2 E-values of the 389 unmatched cases displays higher significance than that of HMMER3. Meanwhile, the fold-critical E-value boxplots display higher significance than the remnant boxplots whether for the case of HMMER2 or HMMER3. Additionally, both the fold-critical boxplots were considered significant when compared to their corresponding remnant counterparts. Taken together, the dissectHMMER analysis is strongly suggestive that these 389 HMMER2/HMMER3-based Pfam annotations were justifiable by their fold-critical plots; the annotation discrepancy is resolved.
Scrutiny of the remnant boxplots further revealed that the remnant E-values of the HMMER3 were so heavily penalized that they exceeded the recommended cutoff of 0.1 (−1 in log10 scale). Therefore, when the scores were added up as part of the total HMMER3 E-values, the latter becomes much less significant than that of HMMER2 and hence their larger FPRs. It is likely that HMMER3 cannot justify any independent non-globular alignment segments of the domain model without significant score compensation from other alignment pieces. Unfortunately, it is beyond the scope of this work to examine why the remnant segments of these particular HMMER3 models appeared to be so heavily penalized though this sheds some light on why certain HMMER3 models were less sensitive than their HMMER2 counterparts.