BACKGROUND
the identification of functionally relevant features in protein sequences is an important task for gaining insight into their molecular and biological activities.
commonly used feature classifications systems focus on protein regions of different lengths ranging from single residues in active site representations and relatively short sequence motifs to much longer protein domains.
the identification of these functional modules is often of immediate importance for guiding molecular and evolutionary studies of genes and genomes, such as experimental or computational discoveries of drug targets, catalytic residues and ligand binding sites  <cit> .
due to the greater evolutionary constraints, the functionally important regions in proteins tend to be more conserved among related sequences than their less relevant regions.
as a result of this basic similarity-function principle, one can predict the functional features in proteins relatively reliably by identifying their conserved regions  <cit> .
the same information is often useful to predict differences of the catalytic and substrate specificities within subgroups of protein families by identifying their specificity determining residues  <cit> .
profile hidden markov models  provide the basis of very efficient approaches for modeling longer conserved regions in protein families, which are referred to as protein domains  <cit> .
these domain models usually co-align well with longer functional and structural units of proteins, such as protein folds  <cit> .
the genome regions coding for protein domains, rather than entire genes, are often considered the functional base units of protein evolution.
because domain models are relatively complex by covering longer conserved sequence areas, the identification of essential sub-features within protein domains can greatly facilitate their functional characterization.
well known examples are the conserved protein motifs from the prosite database  <cit> .
these much shorter patterns frequently map to residues within protein domains that are directly involved in the core functions of many proteins, such as the coordination of the catalytic centers of enzymes.
the most specific and functionally insightful information about known or predicted active sites is provided by protein structure-based resources, such as the catalytic site atlas , castp, actsitepred, consurf and pdbsite  <cit> .
the utility spectrum of these structure-based resources is typically restricted to proteins that share sequence similarity with proteins of known 3d structure.
this requirement of structure information makes these methods less suited for functional site predictions of many membrane proteins or other difficult to crystallize protein classes.
thus, it is important to develop additional tools that can be used for the prediction of functionally relevant features of all protein classes.
conservation analyses are widely used alternatives for this purpose  <cit> .
typically, these methods aim to identify conserved residues in multiple sequence alignments of related proteins.
based on the above principle, these conserved sites tend to be functionally more important than more variable ones.
more recently developed approaches incorporate additional information with conservation data, such as secondary structure predictions, solvent accessibility data and other parameters  <cit> .
in addition, mistry et al.
<cit>  have developed a set of strict rules that allows the transfer of experimentally validated active site information to other sequences within the same enzyme family.
a disadvantage of most conserved residue approaches is the difficulty of using their data sets without major modifications for search applications in order to identify novel proteins containing these features.
the more information rich motif and domain models are usually more effective in this regard.
this is also facilitated by the availability of many efficient motif or domain search algorithms in this area.
much of the information available in conserved sequence databases is the direct result of mining the available protein space with existing feature prediction tools.
this includes very established databases on protein motif or domain information, such as prosite, interpro and pfam  <cit> .
however, the annotation and curation process of the conserved features provided by these databases is still a very time consuming and largely manual curation processes by many experts in the field.
therefore, the development of additional functional feature prediction methods, that can facilitate the automation of various steps in this laborious annotation process, will be of great importance for the field.
here we propose an automated method for identifying conserved protein motifs by creating sub-hmms from custom or existing profile hmm data sets, such as pfam.
the method builds on existing profile hmm domain models and expands their utility spectrum to motif discovery.
the approach has many applications for studying protein functions.
first, it is useful for predicting the most highly conserved and functionally relevant sequence motifs in protein families.
second, it provides an effective alternative for profile-based similarity searches to detect sequences with short similarities in any order.
finally, it can be used for the characterization of domains of unknown function by associating them with sub-hmms from functionally characterized domains.
the most closely related method for modeling protein families by a fragment-based approach was proposed by plotz and fink  <cit> .
their goal was to minimize the number of parameters used by the model in order to improve its performance on small training sets.
to achieve this, the authors started with a signal-like protein sequence representation  <cit>  and trained a new model on this data set.
their model consisted of sub-protein units  that were concatenated in an order learned from the data set.
each spu of this method is an hmm by itself.
in contrast to this, our method uses pre-calculated profile hmms to discover functionally relevant motifs in protein domains.
in addition, our method has the ability to allow any combination of sub-hmms to occur in any order.
another related method is meta-meme  <cit> .
this method also minimizes the number of model parameters.
it accomplishes this by concatenating short pssms  instead of hmms, which are generated by its sister tool meme  <cit> .
this approach is similar to the blockmaker program  <cit> , which also models conserved regions with un-gapped pssms.
our method differs from these approaches significantly by retaining full hmms of the most highly conserved sub-regions within protein domain families.
this allows us to model more complex consensus regions containing gaps.
the method developed by sun and buhler  <cit>  attempts to speed up searching with profile hmms by extracting un-gapped subsections  of hmms and then modifying the match distributions in each position to make each block as sensitive as possible.
these blocks are then used as pre-filters to eliminate sequences which would not match the whole hmm well.
our proposed protein sub-hmm method starts with a profile hmm that has been trained on the multiple sequence alignment of a protein family.
we then extract the most conserved sub-hmms from the original hmm.
a robust scoring method is used to predict the presence of the sub-hmms in any protein sequence of interest.
the hmms required for this approach can be easily generated from unaligned protein sequences of interest by aligning them with a multiple sequence alignment program and then generating an hmm for them with tools like hmmer  <cit>  or sam  <cit> .
alternatively, one can use existing protein family hmms from databases like pfam  <cit> .
the latter approach is taken in this paper for benchmarking the proposed protein sub-hmm method.
RESULTS
a profile hidden markov model of a sequence family is a statistical model over sequences whose structure consists of a number of states and transitions between states.
for each state z there is a distribution, p over a set of observations, x ∈ x. in our case, x is the set of amino acids.
a transition matrix t defines the probability of transitioning from state z <dig> to state z <dig>  we can view this transition matrix as a graph in which a link exists from z <dig> to z <dig> if t >  <dig>  figure  <dig> shows the structure used for aligning protein sequences  <cit> .
for each nominal position i there are three possible states: a match state mi, an insert state ii, and a delete state di.
p  is a distribution over amino acids occurring at position i. p is a background distribution, which is the probability of each amino acid occurring given no other information.
this state is used to model noise sections in the input sequence.
the delete state does not have a real observation distribution; it requires that nothing be observed .
this is used to model sections of the input sequence which have been lost.
the parameters of an hmm can be learned using the expectation maximization   <cit>  algorithm given a set of observed protein sequences , producing a model tuned to this set of protein sequences.
once the model has been trained, we can take another protein sequence, s, and ask what is the most likely sequence of hmm states to generate s, and what is the probability of that combination of states and observations.
this is done with the viterbi algorithm  <cit> .
to rank the results, it is common to calculate the log-odds:
in this equation, pback is the probability of s, assuming each amino acid has been drawn independently from the background distribution, while phmm is the probability that the hmm would generate the state sequence z and the observed sequence s. a positive score means that s is more likely to be derived from the hmm than randomly generated from the background distribution.
a more detailed description of profile hmms can be found in  <cit> .
extraction of sub-hmms
our sub-hmm method is built on top of the well-established profile hmm framework described above.
the algorithm consists of a simple but effective two step procedure for extracting the most highly conserved regions from profile hmms .
first, the kullback-leibler divergence is calculated for all columns of a profile hmm  <cit> .
second, after a series of normalization and smoothing steps , the most information rich hmm regions are excised from the original profile hmms.
the resulting sub-hmms have the same structure as the original profile hmms, but they are usually much shorter.
typically, the method will extract several non-overlapping sub-hmms from a single domain model, especially when its most conserved regions are highly localized and discontinuous.
a more detailed outline of the algorithm for extracting sub-hmms and using them for scoring their presence in protein sequences is described in the methods section.
in the following outline we first describe our sub-hmm experiments and provide several performance comparisons to related tools.
subsequently, we use our tool to find sequences that share short sequence features encoded in our sub-hmms.
properties of sub-hmms
sub-hmms were extracted from pfam domain families using hmmer <dig> and hmmer <dig> models  <cit> .
pfam  <dig>  was used for all experiments, whereas pfam  <dig>  was mainly used in the performance comparisons with hmmer <dig>  this is because pfam has adopted hmmer <dig> models only very recently, and at this point many of its families have not been as rigorously tested and curated by experts in the field as in the earlier hmmer2-based releases.
using our new sub-hmm method, we extracted  <dig>  sub-hmms  from the pfam  <dig>  database .
this database consisted of  <dig>  domain profile hmms with  <dig> , <dig> unique protein sequences associated with at least one domain.
due to the presence of multiple domains in many sequences, the data set contained a total of  <dig> , <dig> family memberships.
the length distributions of the original pfam hmms and our sub-hmms for all families are shown in figure  <dig>  as expected the sub-hmms are much shorter than the original pfam hmms, with an average length of  <dig> residues compared to  <dig> residues, respectively.
this has several advantages for the goals of this study.
first, the sub-hmms have a length distribution similar to the size of many known functional motifs, which is essential for predicting features with related properties  <cit> .
second, their shorter length reduces the computation time for scoring a protein.
finally, it reduces the number of parameters, which should improve the accuracy of the detector.
the table provides the sizes of the different data sets used and generated by this study using pfam  <dig> .
subsequently, we performed several benchmark tests to determine the performance of the new sub-hmm method in identifying functionally relevant sequence features and searching for sequences sharing them.
for this, we determined the presence of each pfam hmm and our sub-hmms in all protein sequences from the pfam database by applying the scoring system described in the materials section.
we found that the processing time of our method is comparable to hmmer <dig>  the slightly better time performance of our method by a factor  <dig>  is most likely due to the lower complexity of its sub-hmm models.
the sub-hmm method showed comparable time improvements when using it with the hmmer <dig> software.
cross-validation with prosite and csa
next, we determined how well the sub-hmm method performed in identifying known motifs that are likely to be of functional relevance.
this was addressed by comparing the extracted sub-hmms from the pfam  <dig>  database with the hand curated conserved protein motifs from the prosite database.
if the sub-hmms are enriched in functionally relevant candidates, then one would expect a high degree of overlap with the motifs from the prosite database.
this should be the case because the prosite motifs are derived from a comparable protein knowledge space as the sub-hmms generated by this study.
the overlaps were determined by comparing the matching positions of the two fragment data sets in their corresponding protein family sequences.
for counting overlaps, we used relatively conservative filtering criteria: the two fragment models had to have 50% of their matching protein sequences in common and the overlaps had to occur in least 95% of the common protein members.
in addition, we consider a sub-hmm to match only if it has a score of  <dig> or higher.
furthermore, we compute the probability of this event happening by chance and require that it be less than  <dig> .
according to these comparisons,  <dig>  of the  <dig>  sub-hmms overlapped with  <dig> of the  <dig>   prosite motifs by at least 10% of the length of the shortest fragment.
the probability of finding ≥ <dig> matches just by chance was estimated to be <  <dig>  * 10- <dig> .
of these  <dig>  prosite motifs,  <dig> were associated by pfam with one or more of its protein families.
the number of matching families for varying percent overlaps is shown in table  <dig>  an example of a matching pair is shown in figure  <dig>  the full result set is available in additional file 1:prosite-comp.tar.
the numbers of sub-hmms are listed that overlapped with prosite motifs.
the first column provides the relative overlap among the two feature types.
the second and third columns contain the number of overlapping sub-hmms and prosite motifs, respectively.
the details of the filter settings used in these comparisons is given in the result and discussion section.
the column tp contains the number of true positives that we identified out of the  <dig> prosite families annotated by pfam  <dig> .
the last column tpr gives the corresponding true positive rate.
a similar test was performed for the catalytic residue annotations from the catalytic site atlas   <cit> .
this is a database of active site residues from enzymes represented in the protein data bank .
due to their functional importance, most of these residues are highly conserved within protein families.
in our tests, we considered only those sites which are supported by the literature and also mapped to protein domain regions in the pfam data set.
this left us with  <dig> sites mapping to  <dig> proteins.
subsequently, we counted how many sub-hmms overlapped with these sites and found that  <dig> sub-hmms overlapped with csa residues.
these corresponded to  <dig> active sites from  <dig> proteins.
thus, our sub-hmm data set contained 70% of these active sites.
the probability of observing ≥ <dig> overlaps among the two data sets just by chance is <  <dig>  * 10- <dig>  the complete result set of this analysis is available in additional file 2:csa-comp.
the considerable agreement of our method with the prosite and csa data sets indicates that the sub-hmm method is efficient in identifying many of the known functionally important residues in protein families.
therefore, it is reasonable to assume that the novel conserved regions, identified by this study, are a useful resource for characterizing the functional hotspots in protein sequences of known or unknown function in the future.
search performance comparisons
to compare the sensitivity and selectivity performance of the sub-hmm method with the widely used hmmer <dig> software, we tested how well each method could recover the members of each domain family from all proteins in the entire pfam  <dig>  database.
we used the scores computed for each protein to generate an roc  curve for each method .
this allowed us to compare the methods without choosing a fixed threshold, which is usually hard to define a priori.
in this preliminary test, we used the original pfam hmms for the hmmer <dig> method, and the sub-hmms extracted by our method from the same pfam hmms.
as a test sample, all proteins in pfam were used.
this experimental design gives a slight advantage to both methods, because the pfam hmms are trained on a representative subset of proteins that overlaps with the total protein set in each family.
despite this limitation, the difference in performance is still meaningful due to the identical starting conditions for both methods.
figure  <dig> shows the resulting roc curves for assembling all  <dig>  families.
the results show that the hmmer <dig> method has a higher sensitivity at false positive rates less than  <dig> , but the sub-hmm method performs slightly better at higher false positive rates.
due to the much shorter profiles used by our method, it is expected to have a higher false positive rate when it is benchmarked against a test data set that is based on the family assignments of complete domain models.
we also performed more rigorous comparisons of our method against hmmer <dig>  hmmer <dig>  sam and psi-blast  <cit> .
additionally, we tested our sub-hmm method with hmmer <dig> profile hmms.
in this case the sub-hmms where excised from hmmer <dig> models and the hmmer <dig> search tool was used to map and score the individual sub-hmms to the sequences.
we then combined the scores as described in the methods section.
in the following text of this section, the sub-hmm experiments performed with hmmer <dig> and hmmer <dig> are referred to sub-hmm-hmmer <dig> and sub-hmm-hmmer <dig>  respectively.
in all tests we trained the models ourselves by randomly selecting 20% of the members from each protein family, but the training data were not included in the test data sets.
hmmer <dig>  hmmer <dig> and sam use a multiple sequence alignment for the model building step.
since it was not our goal to test the alignment quality, we used the curated domain alignments provided by pfam as input to all methods.
although sam can create its own alignments, we forced it to use the alignments we provided to make this method more comparable to hmmer <dig> and hmmer <dig>  for psi-blast, we first created multiple sequence alignments for all the training data sets using clustalw.
subsequently, we built pssms to search the test data set with psi-blast.
for all methods, we compared how well they could recover the remaining 80% in each protein family from the combined set of all test sequences.
due to computational resource constraints, it was not possible to test these methods on all pfam families.
instead we created two smaller subsets of families, one composed of smaller families and one composed of larger families.
the small family set contained  <dig> families randomly selected from pfam  <dig>  with of  <dig> to  <dig> members, while the large set contained  <dig> families with more than  <dig> members.
in addition, we tested the different methods on the hmmer3-based pfam  <dig>  data set.
to maximize the comparability of the results, we selected only families that were available in both pfam releases and fell into the same size categories.
for the small set, we found  <dig> families in pfam  <dig>  but only  <dig> of them had less than  <dig> members.
for the large set,  <dig> families were also available in pfam  <dig>  and all of them contained more than  <dig> members.
the roc plots for all comparisons are shown in figures  <dig> and  <dig>  for the experiments with pfam  <dig> , the results indicate that the sub-hmm-hmmer <dig>  sub-hmm-hmmer <dig> and psi-blast methods perform better on the small family set than on the large one, while hmmer <dig>  hmmer <dig> and sam show an opposite performance trend.
when comparing the six methods, both sub-hmm methods perform at least as well as hmmer <dig>  whereas sam, hmmer <dig> and psi-blast show the best performance in assembling the families from both family size categories.
direct comparisons of the pfam  <dig>  and pfam  <dig>  results indicate that hmmer <dig>  psi-blast, sam, sub-hmm-hmmer <dig> and sub-hmm-hmmer <dig> perform very similarly on the small family set, while hmmer <dig> improves slightly.
these trends are almost identical for the large family set, except that sub-hmm-hmmer <dig> performs better on this data set.
since our method is designed to find short sequence similarities, it is expected to have a lower selectivity  than the other methods when reassembling family relationships that are based on longer domain similarities.
in fact, such a performance characteristics on known family data sets is required for discovering novel conserved fragments in sequences that do not necessarily belong to the same domain family.
the latter is the main utility feature of the sub-hmm method.
discovery of conserved fragments in protein families with sub-hmms
to evaluate the utility spectrum of sub-hmms for conserved feature discovery, we determined for each sub-hmm excised from pfam  <dig>  its matching profile against different domain families in the same pfam release.
to define a match, we required a sub-hmm to match at least 50% of the sequences in each pfam family with a log-odds score of  <dig> or higher.
table  <dig> shows how many sub-hmms from pfam domains of unknown function  matched pfam families of known function  and vice versa.
a sub-duf is defined as a sub-hmm that was extracted from a duf, whereas a sub-dkf was extracted from a dkf.
interestingly, the sub-dkfs shows considerable overlaps with the prosite data set, whereas the sub-dufs do not overlap with prosite at all .
the latter is expected because prosite focuses on motifs from functionally characterized proteins.
this also indicates that our sub-duf data sets contains many novel conserved and potentially functional motifs that are not represented in prosite.
the table lists the numbers of sub-dkfs and sub-dufs which matched in addition to their source families other dkf and duf families.
a sub-hmm is considered to have matched a pfam  <dig>  family if it scores greater than  <dig> on more than 50% of its members.
the last column contains the counts of sub-hmms that also overlapped with prosite motifs.
a similar approach was used for constructing networks of pfam  <dig>  families by their common sub-hmm matches.
the obtained clusters in this network showed many similarities to the clusters from the pfam clan database, but also significant differences  <cit> .
the variation of information  coefficient  <cit>  for the two network sets was  <dig> .
this score has a range from  <dig> to log =  <dig> , with lower scores indicating more similar clusterings.
two small sub-graphs of the sub-hmm based domain network are shown in figures  <dig> and  <dig>  the box in figure  <dig> encloses those families which are part of a clan according to the pfam database.
in this case the sub-hmm-based grouping of families agrees almost perfectly with the corresponding pfam clan assignment.
in contrast to this, figure  <dig> gives an example of a new cluster of domains predicted by our method.
such differences in the results of the two methods are expected, because the pfam clans are assembled with a profile hmm to profile hmm alignment method  <cit>  that is fundamentally different from our sub-hmm method.
the large number of sub-hmms matching different pfam domains indicates the usefulness of our sub-hmm approach for discovering short sequence features that are conserved among different protein domains.
due to their high conservation, an important functional role for many of these features can be expected.
many of the sub-dkfs will be useful for assigning potential functions to dufs.
a much more comprehensive study on applying our sub-hmm approach to biologically relevant questions will be published in an experimental journal.
CONCLUSIONS
we have developed a simple but effective method for identifying the most highly conserved residues in protein sequences in a fully automated manner.
its design strategy is highly practical and versatile by making efficient use of a well-established bioinformatic infrastructure, such as existing domain databases and profile hmm search tools.
in addition, the conserved patterns, identified by this study, are useful for characterizing proteins of unknown function by associating them with those of known function by their common sub-hmms.
furthermore, the sub-hmm search method appears to be a very effective tool for finding sequences that share only very short sequence similarities with a sensitivity performance similar to hmmer <dig>  the possibility to ignore the order of different sub-hmm matches in sequences is another advantage, which will allow the identification of more complex similarity arrangements among otherwise unrelated sequences.

 the identification of functionally relevant features in protein sequences is an important task for gaining insight into their molecular and biological activities. commonly used feature classifications systems focus on protein regions of different lengths ranging from single residues in active site representations and relatively short sequence motifs to much longer protein domains. the identification of these functional modules is often of immediate importance for guiding molecular and evolutionary studies of genes and genomes, such as experimental or computational discoveries of drug targets, catalytic residues and ligand binding sites. due to the greater evolutionary constraints, the functionally important regions in proteins tend to be more conserved among related sequences than their less relevant regions. as a result of this basic similarity-function principle, one can predict the functional features in proteins relatively reliably by identifying their conserved regions. the same information is often useful to predict differences of the catalytic and substrate specificities within subgroups of protein families by identifying their specificity determining residues. profile hidden markov models provide the basis of very efficient approaches for modeling longer conserved regions in protein families, which are referred to as protein domains. these domain models usually co-align well with longer functional and structural units of proteins, such as protein folds. the genome regions coding for protein domains, rather than entire genes, are often considered the functional base units of protein evolution. because domain models are relatively complex by covering longer conserved sequence areas, the identification of essential sub-features within protein domains can greatly facilitate their functional characterization. well known examples are the conserved protein motifs from the prosite database. these much shorter patterns frequently map to residues within protein domains that are directly involved in the core functions of many proteins, such as the coordination of the catalytic centers of enzymes. the most specific and functionally insightful information about known or predicted active sites is provided by protein structure-based resources, such as the catalytic site atlas, castp, actsitepred, consurf and pdbsite. the utility spectrum of these structure-based resources is typically restricted to proteins that share sequence similarity with proteins of known 3d structure. this requirement of structure information makes these methods less suited for functional site predictions of many membrane proteins or other difficult to crystallize protein classes. thus, it is important to develop additional tools that can be used for the prediction of functionally relevant features of all protein classes. conservation analyses are widely used alternatives for this purpose. typically, these methods aim to identify conserved residues in multiple sequence alignments of related proteins. based on the above principle, these conserved sites tend to be functionally more important than more variable ones. more recently developed approaches incorporate additional information with conservation data, such as secondary structure predictions, solvent accessibility data and other parameters. in addition, mistry et al. have developed a set of strict rules that allows the transfer of experimentally validated active site information to other sequences within the same enzyme family. a disadvantage of most conserved residue approaches is the difficulty of using their data sets without major modifications for search applications in order to identify novel proteins containing these features. the more information rich motif and domain models are usually more effective in this regard. this is also facilitated by the availability of many efficient motif or domain search algorithms in this area. much of the information available in conserved sequence databases is the direct result of mining the available protein space with existing feature prediction tools. this includes very established databases on protein motif or domain information, such as prosite, interpro and pfam. however, the annotation and curation process of the conserved features provided by these databases is still a very time consuming and largely manual curation processes by many experts in the field. therefore, the development of additional functional feature prediction methods, that can facilitate the automation of various steps in this laborious annotation process, will be of great importance for the field. here we propose an automated method for identifying conserved protein motifs by creating sub-hmms from custom or existing profile hmm data sets, such as pfam. the method builds on existing profile hmm domain models and expands their utility spectrum to motif discovery. the approach has many applications for studying protein functions. first, it is useful for predicting the most highly conserved and functionally relevant sequence motifs in protein families. second, it provides an effective alternative for profile-based similarity searches to detect sequences with short similarities in any order. finally, it can be used for the characterization of domains of unknown function by associating them with sub-hmms from functionally characterized domains. the most closely related method for modeling protein families by a fragment-based approach was proposed by plotz and fink. their goal was to minimize the number of parameters used by the model in order to improve its performance on small training sets. to achieve this, the authors started with a signal-like protein sequence representation and trained a new model on this data set. their model consisted of sub-protein units that were concatenated in an order learned from the data set. each spu of this method is an hmm by itself. in contrast to this, our method uses pre-calculated profile hmms to discover functionally relevant motifs in protein domains. in addition, our method has the ability to allow any combination of sub-hmms to occur in any order. another related method is meta-meme. this method also minimizes the number of model parameters. it accomplishes this by concatenating short pssms instead of hmms, which are generated by its sister tool meme. this approach is similar to the blockmaker program, which also models conserved regions with un-gapped pssms. our method differs from these approaches significantly by retaining full hmms of the most highly conserved sub-regions within protein domain families. this allows us to model more complex consensus regions containing gaps. the method developed by sun and buhler attempts to speed up searching with profile hmms by extracting un-gapped subsections of hmms and then modifying the match distributions in each position to make each block as sensitive as possible. these blocks are then used as pre-filters to eliminate sequences which would not match the whole hmm well. our proposed protein sub-hmm method starts with a profile hmm that has been trained on the multiple sequence alignment of a protein family. we then extract the most conserved sub-hmms from the original hmm. a robust scoring method is used to predict the presence of the sub-hmms in any protein sequence of interest. the hmms required for this approach can be easily generated from unaligned protein sequences of interest by aligning them with a multiple sequence alignment program and then generating an hmm for them with tools like hmmer or sam. alternatively, one can use existing protein family hmms from databases like pfam. the latter approach is taken in this paper for benchmarking the proposed protein sub-hmm method.  a profile hidden markov model of a sequence family is a statistical model over sequences whose structure consists of a number of states and transitions between states. for each state z there is a distribution, p over a set of observations, x ∈ x. in our case, x is the set of amino acids. a transition matrix t defines the probability of transitioning from state z to state z we can view this transition matrix as a graph in which a link exists from z to z if t > figure shows the structure used for aligning protein sequences. for each nominal position i there are three possible states: a match state mi, an insert state ii, and a delete state di. p is a distribution over amino acids occurring at position i. p is a background distribution, which is the probability of each amino acid occurring given no other information. this state is used to model noise sections in the input sequence. the delete state does not have a real observation distribution; it requires that nothing be observed. this is used to model sections of the input sequence which have been lost. the parameters of an hmm can be learned using the expectation maximization algorithm given a set of observed protein sequences, producing a model tuned to this set of protein sequences. once the model has been trained, we can take another protein sequence, s, and ask what is the most likely sequence of hmm states to generate s, and what is the probability of that combination of states and observations. this is done with the viterbi algorithm. to rank the results, it is common to calculate the log-odds: in this equation, pback is the probability of s, assuming each amino acid has been drawn independently from the background distribution, while phmm is the probability that the hmm would generate the state sequence z and the observed sequence s. a positive score means that s is more likely to be derived from the hmm than randomly generated from the background distribution. a more detailed description of profile hmms can be found in. extraction of sub-hmms our sub-hmm method is built on top of the well-established profile hmm framework described above. the algorithm consists of a simple but effective two step procedure for extracting the most highly conserved regions from profile hmms. first, the kullback-leibler divergence is calculated for all columns of a profile hmm. second, after a series of normalization and smoothing steps, the most information rich hmm regions are excised from the original profile hmms. the resulting sub-hmms have the same structure as the original profile hmms, but they are usually much shorter. typically, the method will extract several non-overlapping sub-hmms from a single domain model, especially when its most conserved regions are highly localized and discontinuous. a more detailed outline of the algorithm for extracting sub-hmms and using them for scoring their presence in protein sequences is described in the methods section. in the following outline we first describe our sub-hmm experiments and provide several performance comparisons to related tools. subsequently, we use our tool to find sequences that share short sequence features encoded in our sub-hmms. properties of sub-hmms sub-hmms were extracted from pfam domain families using hmmer and hmmer models. pfam was used for all experiments, whereas pfam was mainly used in the performance comparisons with hmmer this is because pfam has adopted hmmer models only very recently, and at this point many of its families have not been as rigorously tested and curated by experts in the field as in the earlier hmmer2-based releases. using our new sub-hmm method, we extracted sub-hmms from the pfam database. this database consisted of domain profile hmms with, unique protein sequences associated with at least one domain. due to the presence of multiple domains in many sequences, the data set contained a total of, family memberships. the length distributions of the original pfam hmms and our sub-hmms for all families are shown in figure as expected the sub-hmms are much shorter than the original pfam hmms, with an average length of residues compared to residues, respectively. this has several advantages for the goals of this study. first, the sub-hmms have a length distribution similar to the size of many known functional motifs, which is essential for predicting features with related properties. second, their shorter length reduces the computation time for scoring a protein. finally, it reduces the number of parameters, which should improve the accuracy of the detector. the table provides the sizes of the different data sets used and generated by this study using pfam. subsequently, we performed several benchmark tests to determine the performance of the new sub-hmm method in identifying functionally relevant sequence features and searching for sequences sharing them. for this, we determined the presence of each pfam hmm and our sub-hmms in all protein sequences from the pfam database by applying the scoring system described in the materials section. we found that the processing time of our method is comparable to hmmer the slightly better time performance of our method by a factor is most likely due to the lower complexity of its sub-hmm models. the sub-hmm method showed comparable time improvements when using it with the hmmer software. cross-validation with prosite and csa next, we determined how well the sub-hmm method performed in identifying known motifs that are likely to be of functional relevance. this was addressed by comparing the extracted sub-hmms from the pfam database with the hand curated conserved protein motifs from the prosite database. if the sub-hmms are enriched in functionally relevant candidates, then one would expect a high degree of overlap with the motifs from the prosite database. this should be the case because the prosite motifs are derived from a comparable protein knowledge space as the sub-hmms generated by this study. the overlaps were determined by comparing the matching positions of the two fragment data sets in their corresponding protein family sequences. for counting overlaps, we used relatively conservative filtering criteria: the two fragment models had to have 50% of their matching protein sequences in common and the overlaps had to occur in least 95% of the common protein members. in addition, we consider a sub-hmm to match only if it has a score of or higher. furthermore, we compute the probability of this event happening by chance and require that it be less than. according to these comparisons, of the sub-hmms overlapped with of the prosite motifs by at least 10% of the length of the shortest fragment. the probability of finding ≥ matches just by chance was estimated to be < * 10-. of these prosite motifs, were associated by pfam with one or more of its protein families. the number of matching families for varying percent overlaps is shown in table an example of a matching pair is shown in figure the full result set is available in additional file 1:prosite-comp.tar. the numbers of sub-hmms are listed that overlapped with prosite motifs. the first column provides the relative overlap among the two feature types. the second and third columns contain the number of overlapping sub-hmms and prosite motifs, respectively. the details of the filter settings used in these comparisons is given in the result and discussion section. the column tp contains the number of true positives that we identified out of the prosite families annotated by pfam. the last column tpr gives the corresponding true positive rate. a similar test was performed for the catalytic residue annotations from the catalytic site atlas. this is a database of active site residues from enzymes represented in the protein data bank. due to their functional importance, most of these residues are highly conserved within protein families. in our tests, we considered only those sites which are supported by the literature and also mapped to protein domain regions in the pfam data set. this left us with sites mapping to proteins. subsequently, we counted how many sub-hmms overlapped with these sites and found that sub-hmms overlapped with csa residues. these corresponded to active sites from proteins. thus, our sub-hmm data set contained 70% of these active sites. the probability of observing ≥ overlaps among the two data sets just by chance is < * 10- the complete result set of this analysis is available in additional file 2:csa-comp. the considerable agreement of our method with the prosite and csa data sets indicates that the sub-hmm method is efficient in identifying many of the known functionally important residues in protein families. therefore, it is reasonable to assume that the novel conserved regions, identified by this study, are a useful resource for characterizing the functional hotspots in protein sequences of known or unknown function in the future. search performance comparisons to compare the sensitivity and selectivity performance of the sub-hmm method with the widely used hmmer software, we tested how well each method could recover the members of each domain family from all proteins in the entire pfam database. we used the scores computed for each protein to generate an roc curve for each method. this allowed us to compare the methods without choosing a fixed threshold, which is usually hard to define a priori. in this preliminary test, we used the original pfam hmms for the hmmer method, and the sub-hmms extracted by our method from the same pfam hmms. as a test sample, all proteins in pfam were used. this experimental design gives a slight advantage to both methods, because the pfam hmms are trained on a representative subset of proteins that overlaps with the total protein set in each family. despite this limitation, the difference in performance is still meaningful due to the identical starting conditions for both methods. figure shows the resulting roc curves for assembling all families. the results show that the hmmer method has a higher sensitivity at false positive rates less than, but the sub-hmm method performs slightly better at higher false positive rates. due to the much shorter profiles used by our method, it is expected to have a higher false positive rate when it is benchmarked against a test data set that is based on the family assignments of complete domain models. we also performed more rigorous comparisons of our method against hmmer hmmer sam and psi-blast. additionally, we tested our sub-hmm method with hmmer profile hmms. in this case the sub-hmms where excised from hmmer models and the hmmer search tool was used to map and score the individual sub-hmms to the sequences. we then combined the scores as described in the methods section. in the following text of this section, the sub-hmm experiments performed with hmmer and hmmer are referred to sub-hmm-hmmer and sub-hmm-hmmer respectively. in all tests we trained the models ourselves by randomly selecting 20% of the members from each protein family, but the training data were not included in the test data sets. hmmer hmmer and sam use a multiple sequence alignment for the model building step. since it was not our goal to test the alignment quality, we used the curated domain alignments provided by pfam as input to all methods. although sam can create its own alignments, we forced it to use the alignments we provided to make this method more comparable to hmmer and hmmer for psi-blast, we first created multiple sequence alignments for all the training data sets using clustalw. subsequently, we built pssms to search the test data set with psi-blast. for all methods, we compared how well they could recover the remaining 80% in each protein family from the combined set of all test sequences. due to computational resource constraints, it was not possible to test these methods on all pfam families. instead we created two smaller subsets of families, one composed of smaller families and one composed of larger families. the small family set contained families randomly selected from pfam with of to members, while the large set contained families with more than members. in addition, we tested the different methods on the hmmer3-based pfam data set. to maximize the comparability of the results, we selected only families that were available in both pfam releases and fell into the same size categories. for the small set, we found families in pfam but only of them had less than members. for the large set, families were also available in pfam and all of them contained more than members. the roc plots for all comparisons are shown in figures and for the experiments with pfam, the results indicate that the sub-hmm-hmmer sub-hmm-hmmer and psi-blast methods perform better on the small family set than on the large one, while hmmer hmmer and sam show an opposite performance trend. when comparing the six methods, both sub-hmm methods perform at least as well as hmmer whereas sam, hmmer and psi-blast show the best performance in assembling the families from both family size categories. direct comparisons of the pfam and pfam results indicate that hmmer psi-blast, sam, sub-hmm-hmmer and sub-hmm-hmmer perform very similarly on the small family set, while hmmer improves slightly. these trends are almost identical for the large family set, except that sub-hmm-hmmer performs better on this data set. since our method is designed to find short sequence similarities, it is expected to have a lower selectivity than the other methods when reassembling family relationships that are based on longer domain similarities. in fact, such a performance characteristics on known family data sets is required for discovering novel conserved fragments in sequences that do not necessarily belong to the same domain family. the latter is the main utility feature of the sub-hmm method. discovery of conserved fragments in protein families with sub-hmms to evaluate the utility spectrum of sub-hmms for conserved feature discovery, we determined for each sub-hmm excised from pfam its matching profile against different domain families in the same pfam release. to define a match, we required a sub-hmm to match at least 50% of the sequences in each pfam family with a log-odds score of or higher. table shows how many sub-hmms from pfam domains of unknown function matched pfam families of known function and vice versa. a sub-duf is defined as a sub-hmm that was extracted from a duf, whereas a sub-dkf was extracted from a dkf. interestingly, the sub-dkfs shows considerable overlaps with the prosite data set, whereas the sub-dufs do not overlap with prosite at all. the latter is expected because prosite focuses on motifs from functionally characterized proteins. this also indicates that our sub-duf data sets contains many novel conserved and potentially functional motifs that are not represented in prosite. the table lists the numbers of sub-dkfs and sub-dufs which matched in addition to their source families other dkf and duf families. a sub-hmm is considered to have matched a pfam family if it scores greater than on more than 50% of its members. the last column contains the counts of sub-hmms that also overlapped with prosite motifs. a similar approach was used for constructing networks of pfam families by their common sub-hmm matches. the obtained clusters in this network showed many similarities to the clusters from the pfam clan database, but also significant differences. the variation of information coefficient for the two network sets was. this score has a range from to log =, with lower scores indicating more similar clusterings. two small sub-graphs of the sub-hmm based domain network are shown in figures and the box in figure encloses those families which are part of a clan according to the pfam database. in this case the sub-hmm-based grouping of families agrees almost perfectly with the corresponding pfam clan assignment. in contrast to this, figure gives an example of a new cluster of domains predicted by our method. such differences in the results of the two methods are expected, because the pfam clans are assembled with a profile hmm to profile hmm alignment method that is fundamentally different from our sub-hmm method. the large number of sub-hmms matching different pfam domains indicates the usefulness of our sub-hmm approach for discovering short sequence features that are conserved among different protein domains. due to their high conservation, an important functional role for many of these features can be expected. many of the sub-dkfs will be useful for assigning potential functions to dufs. a much more comprehensive study on applying our sub-hmm approach to biologically relevant questions will be published in an experimental journal.  we have developed a simple but effective method for identifying the most highly conserved residues in protein sequences in a fully automated manner. its design strategy is highly practical and versatile by making efficient use of a well-established bioinformatic infrastructure, such as existing domain databases and profile hmm search tools. in addition, the conserved patterns, identified by this study, are useful for characterizing proteins of unknown function by associating them with those of known function by their common sub-hmms. furthermore, the sub-hmm search method appears to be a very effective tool for finding sequences that share only very short sequence similarities with a sensitivity performance similar to hmmer the possibility to ignore the order of different sub-hmm matches in sequences is another advantage, which will allow the identification of more complex similarity arrangements among otherwise unrelated sequences. 
profile hmms provide effective methods for modeling the conserved regions of protein families. to identify these conserved motifs efficiently, we propose a method for extracting the most information-rich regions in protein families from their profile hmms. sub-hmms extend the application spectrum of profile hmms to motif discovery. finally, the sub-hmm method showed very promising results as an alternative search method for identifying proteins that share only short sequence similarities. cross-validations with the prosite and csa databases confirmed the efficiency of the method in predicting most of the known functionally relevant motifs and residues. the data set also allowed us to link at least pfam domains of known and unknown function by their common sub-hmms. their most interesting utility is the identification of the functionally relevant residues in proteins of known and unknown function. a limitation of the resulting domain models is the difficulty to pinpoint their much shorter functional sub-features, such as catalytically relevant sequence motifs in enzymes or ligand binding signatures of receptor proteins. the motif data generated by this study is a valuable knowledge resource for characterizing protein functions in the future.
 profile hmms provide effective methods for modeling the conserved regions of protein families. a limitation of the resulting domain models is the difficulty to pinpoint their much shorter functional sub-features, such as catalytically relevant sequence motifs in enzymes or ligand binding signatures of receptor proteins.  to identify these conserved motifs efficiently, we propose a method for extracting the most information-rich regions in protein families from their profile hmms. the method was used here to predict a comprehensive set of sub-hmms from the pfam domain database. cross-validations with the prosite and csa databases confirmed the efficiency of the method in predicting most of the known functionally relevant motifs and residues. at the same time, novel conserved regions could be predicted. the data set also allowed us to link at least pfam domains of known and unknown function by their common sub-hmms. finally, the sub-hmm method showed very promising results as an alternative search method for identifying proteins that share only short sequence similarities.  sub-hmms extend the application spectrum of profile hmms to motif discovery. their most interesting utility is the identification of the functionally relevant residues in proteins of known and unknown function. additionally, sub-hmms can be used for highly localized sequence similarity searches that focus on shorter conserved features rather than entire domains or global similarities. the motif data generated by this study is a valuable knowledge resource for characterizing protein functions in the future. 