BACKGROUND
recently there has been a rapid growth in the number of putative proteins derivable from new genomic and metagenomic data  <cit> .
the extended use of environmental shotgun sequencing to study diverse microbial systems has made metagenomics a vastly growing field leading to a flux of data, calling for development and application of new tools that allow its investigation  <cit> .
conventional tools for predicting the function of a protein from its sequence are based on sequence-similarity  <cit>  or sequence-motifs  <cit> .
here we outline a relatively simple and straight-forward method that is applicable to large numbers of sequences.
its purpose is finding whether each protein in the data is an enzyme and, if so, what its ec classification is.
this data mining of enzymes  is based on the specific peptide  method of  <cit> , and is carried out by comparing the sequences of all proteins with a list of all sps and looking for matches of the latter in the data.
sps are strings of amino-acids, extracted from enzyme sequences using the motif extraction algorithm mex  <cit> .
they are selected for their specificity to levels of the enzyme commission  4-level functional hierarchy.
we have updated the sp set of  <cit>  by extracting it from all swiss-prot enzymes of july 1st,  <dig>  more details are provided in methods.
using sps for prediction of enzymatic function needs some further decisions as to what to do if various sp hits on the same protein have ec assignments that are not consistent with one another.
moreover, one should decide when a single sp hit is sufficient to make a prediction.
the methodology developed here relies on coverage length  of consistent sp hits.
this is further described below, when testing performance on an enzyme test set, and when discussing a ten-organism test-set that contains non-enzymatic as well as enzymatic proteins.
we develop a random model for the latter to assess the effect of accidental sp matches.
the resulting methodology, which we call data mining of enzymes , is being applied to analyze several metagenomes.
methods
the new sp sets
a novel method based on sequence motifs has been proposed by  <cit> , who have studied enzymes in the swiss-prot database.
they have demonstrated that enzyme functions, as represented by the four-level ec hierarchy, can be deduced from the appearance of deterministic short strings of amino-acids, denoted as specific peptides , on these enzymes.
the sps were derived from enzyme sequence data using an unsupervised motif extraction algorithm mex  <cit> , and filtered by the ec so that each sp is specific to a particular ec branch, specifying the ec function that the enzyme performs.
thus, if an extracted motif is found to occur on enzymes belonging to only one ec number , this peptide will be declared to be an sp labeled with this ec number.
if, however, the motif occurs on several ec numbers, all of which share the same 3rd-level hierarchy , the motif is declared as an sp with labeling at the third level of ec hierarchy, etc.
the sps of  <cit>  comprise on average  <dig>  amino-acids , and were shown to compete favorably with a smith-waterman based svm classifier.
usage of the sp methodology is demonstrated by our web-tool http://adios.tau.ac.il/dme.
given the sequence of an enzyme, this tool searches through the set of all sps and finds which of them coincide with substrings of the sequence, indicating where they lie, what is the ec assignment associated with each sp, and provides the ec predicted by the dme method for the protein that is being queried.
kunik et al  <cit>  have investigated  <dig>  enzyme sequences of the  <dig>  swiss-prot release of october  <dig>  we have used the same methodology and applied it to all enzymes in the swiss-prot/enzyme records of july 1st,  <dig>  the number of enzymes that have a single ec assignment is  <dig> .
applying mex and filtering it by ec levels in the same way as  <cit> , we have obtained  <dig>  sps.
this new 1st list of sps serves as the basis for developing and analyzing our methodology.
in making the prediction of an ec number  based on one sp match, or several sp matches that have the same ec number assignment, we require that the total number of amino-acids of the protein matched with these sps be at least seven.
we refer to this number as the coverage-length l. if l at level  <dig>  l <dig>  is less than  <dig>  we check for sp hits that are consistent at level  <dig> of the ec hierarchy, i.e.
have identical first three digits in their assignments.
once again, a prediction is made if l at level  <dig>  l <dig>  is at least  <dig>  in principle, the threshold of l at every ec level can be viewed as a parameter of our method.
reducing l increases recall at the expense of lowering precision, as will be discussed below.
test data were downloaded from swiss-prot release  <dig> on july 1st,  <dig>  we consider two types of test sets.
the "enzyme test set" consists of all enzymes integrated into swiss-prot between july 1st of  <dig> and  <dig>  the " <dig> organisms test set" consists of proteins of e. coli and  <dig> other bacteria  containing enzymes from the same period of  <dig> to  <dig>  and all other proteins incorporated into swiss-prot by july 1st,  <dig>  a summary of all the relevant data is displayed in table  <dig>  it includes also information about precision and recall  that will be further discussed in the first results section.
these values are obtained by determining 3rd level ec assignments, using coverage-length of l <dig> ≥  <dig>  precision values of 100% on the training sets are of course trivial results of specificity.
54% of the proteins in the 1st training set carry swiss-prot annotations of 'active site', 'binding site' or 'metal binding site' at specific locations of single amino-acids.
sps cover these functionally important sites significantly more than other loci on proteins, thus indicating biological significance of sps .
sp matches that overlap such sites are compiled, and the corresponding sps are denoted as annotated sps .
we have thus compiled a list of  <dig>  asps.
all appear at least four times in the training set, and the location of the annotation is consistent in the different appearances.
most asps carry single annotations ,  <dig> asps carry two annotations and  <dig> asps carry all three annotations.
a second set of sps is extracted from swiss-prot data dated july 27th,  <dig>  this training set, consisting of all singly annotated enzymes, contains  <dig>  proteins.
it has led to  <dig>  sps.
their length distribution is presented in additional file  <dig>  figure s <dig>  this set includes  <dig>  sps with labels corresponding to ec levels  <dig> and  <dig> .
only sps with ec labels at levels  <dig> and  <dig> are relevant for the assignment of ec level- <dig> annotations to proteins, and hence for the calculation of recall included in table  <dig>  it should be emphasized that only  <dig>  of the swiss-prot annotated enzymes in the training set carry ec annotations at levels  <dig> and  <dig>  they are the ones on which the ec predictions at level  <dig> are tested, leading to the recall result of 94%.
the 2nd sp set is being used for the analysis of metagenomic data and is incorporated in our web-tool at http://adios.tau.ac.il/dme.
estimate of accidental sp matches
proteins that do not possess enzymatic functions may still have a substring that matches an sp.
such sp matches will be called 'accidentals'.
their occurrences can be modeled by sp hits on random protein sequences.
such random sequences are generated from real data by scrambling the order of the amino-acids in every protein, conserving only first-order statistics.
<dig> such sets were produced in order to measure the expected random hits.
estimates of the probabilities of accidental occurrences of sps are derived below for the  <dig> organism test-set and for sargasso sea data.
recall-precision analysis of ec annotations in enzymes
comparing the results of our method with an expert-method  we face three possible situations when dealing with a collection of enzymes: p|p where the model prediction coincides with that of the expert, p|dp where the expert provides a different ec assignment, and np|p where the model provides no prediction for enzymes whose ec assignments are given by the expert.
accordingly we define the following measures in terms of number of occurrences:
this is a generalization of the common terms used in binary classification problems where p|p, p|dp and np|p are replaced by true-positive, false-positive and false-negative correspondingly.
recall-precision analysis of ec annotations in proteins
extending the previous analysis to a collection of proteins we have to add two more possibilities: p|np, where the new method has an ec prediction whereas the expert does not have one, and np|np where both do not have any ec assignment.
whereas the latter corresponds to true-negative in a binary classifier, the former, p|np could be added to p|dp as 'false-positive'.
since, however, there are many cases where the absence of an ec assignment does not imply that the protein in question is not an enzyme, we opt to define a new measure, putative novelty ratio, as the fraction of such p|np out of all the predictions of the model:
other measures one can define are
they are the analogs of
and
in conventional binary classifications.
RESULTS
analysis of the methodology
analysis of the enzyme test set using the 1st sp set
in making the prediction of an ec number  based on one sp match, or several sp matches that have the same ec number assignments, we require that the total number of amino-acids of the protein matched with these sps be at least seven.
we refer to this number as the coverage-length l. in principle, the threshold of l at every ec level can be viewed as a parameter of our method.
reducing l increases recall at the expense of lowering precision.
this is exemplified in table  <dig>  where we analyze our enzyme test set and show precision and recall at 3rd ec level as function of the l <dig> threshold.
although precision turns out to be quite high, even for low l <dig> values, recall is low when compared to what blast  <cit>  can achieve on this test-set.
using the most significant outcome of a blast search against the 1st training set as its prediction, and limiting the most significant e-value to stay below e- <dig>  we find blast precision of 98% and recall of 95%, to be compared with dme values of  <dig> % and 70% when setting l <dig> ≥  <dig>  thus while precision is similar, dme loses on recall.
there is no direct relation between dme and blast, although high coverage-length l values of dme go usually hand in hand with very low e-values of blast.
differences may occur for low l values of dme, and relatively high e-values in blast.
we refer to kunik et al.
<cit>  for a discussion of such examples .
the advantages of sps in resolving classification problems in situations of remote homology have been discussed and exemplified by  <cit> .
it is worthwhile pointing out that the fact that one can abide by such a small threshold value of l ≥  <dig> is strongly connected to our requirement that the sp matches on the protein's sequence be exact.
if one were to allow for insertions or deletions or replacements, such as the blosum <dig> matrix  <cit> , this would not work.
based on various trials we may state that, whereas reliance on blosum works well for blast searches over large sequences, it ruins predictivity and specificity of sp searches even if only single amino-acid changes are allowed.
analysis of the ten organism test-set
the ten organism test-set contains  <dig>  proteins of e. coli and  <dig> other bacteria listed in additional file  <dig>  table s <dig>  proteins for this dataset were downloaded from swiss-prot on july 1st  <dig> and include all proteins that had no ec annotation in swiss-prot prior to july 1st  <dig>  the intersection between the  <dig> organism test-set and 1st training set used to build the sps is void and allows us to develop and test the sp methodology on general proteomic data rather than on enzymes only.
sp search on this dataset, using our 1st set of  <dig>  sps , leads to the results shown in figure  <dig>  sorted according to the number n of sp matches.
<dig>  proteins have at least one sp match .
some of them may be due to random hits and our task is to resolve which of the hit proteins should be recognized as enzymes and what their ec assignments should be.
as before, we propose to rely on coverage length.
we judge the prediction not by how many sp hits  are observed, but by l <dig>  the number of amino-acids matched by all sp hits whose ec assignment is identical within the first three digits of the ec number.
in order to have some intuition about the expected noise level, we compare in table  <dig> sp hits on real data with random model results for different values of l <dig>  entries of l <dig> =  <dig> refer to either no sp hits, or hits by sps that have labels with ec levels  <dig> and  <dig> but none at ec levels  <dig> and  <dig>  the columns random and stdev refer to the average and standard deviation of three random sets.
noise is the ratio of random/real.
all  <dig> proteins of the ten-organism test-set were included in this search.
we will use l <dig> ≥  <dig> as our threshold criterion, as in the enzyme data-set discussed in the previous section.
we note that predictions based on l <dig> =  <dig> may still have a large uncertainty, however from l <dig> =  <dig> onwards random hits become very small.
our threshold criterion leads to the results displayed in table  <dig>  with precision =  <dig> %, recall =  <dig> %, accuracy =  <dig> % and putative novelty =  <dig> %.
the interest in this exercise is twofold: to see how well our method performs on unassigned proteins, i.e.
true-negatives, and how good our predictions are for putative novelties.
indeed, our accuracy turns out to be high,  <dig> %, which proves that we have correct negative assignments.
seven out of the  <dig> putative novelties  have been annotated by swiss-prot since july  <dig>  six out of which are at levels  <dig> or  <dig>  all observations are consistent with the predictions, as shown in table  <dig>  quoted here are also all coverage lengths on which the predictions were based.
note that also the one based on coverage length  <dig> has been validated.
all this may be viewed as an indication  of the validity of dme predictions.
the first six entries in table  <dig> belong to e. coli, and the last protein belongs to bacillus cereus.
l <dig> to l <dig> are the coverage-lengths at ec levels  <dig> to  <dig> respectively
classification based on annotated sps
it has been noted by  <cit>  and  <cit>  that some of the sps can be demonstrated to play important biological roles since they carry crucial amino-acids known to serve as active sites, binding sites or metal binding sites.
such annotations are available for 54% of the enzymes in the 1st swiss-prot training set.
selecting only sps that carry these annotations we obtain a set of  <dig>  annotated sps , a mere 7% of all sps.
we have tested it on the enzyme test set.
using annotation predictions at the third level of ec we find precision  <dig> % and recall  <dig> %.
the limited recall is due to the fact that asps have been derived from only 54% of the training set.
nonetheless they possess the advantage of being selected due to their demonstrated operational importance to the catalytic function.
because of their limited recall we have not used the asps as the primary tool for large scale analysis; however we list their properties in our web tool http://adios.tau.ac.il/dme.
any queried protein can be analyzed by this tool for sp hits and the expected dme prediction.
the appearance of asps may serve as providing additional credence to the prediction, as well as specifying the positions of expected active or binding sites.
metagenomic analysis
analysis of sargasso-sea data
after verifying dme on the two test-sets we turn to an analysis of the  <dig> , <dig> records in the sargasso sea protein data  <cit> .
the average length of these proteins is  <dig> amino-acids, with sd =  <dig>  for this analysis we employ our 2nd set of sps, updated on july  <dig>  in order to reduce random hits, we have further limited our sp set to include only peptides of length  <dig> amino-acids or more.
using a random set of  <dig>  proteins selected from these data, we generated three randomized protein sets from which we calculated the probabilities of accidental matches.
the results are displayed as function of l <dig> in table  <dig>  the columns random and stdev refer to the average and standard deviation of three random sets.
noise is the ratio of random/real.
similar results are obtained for l <dig>  the results of table  <dig> are slightly better than table  <dig>  the reason is that we have limited ourselves here to sps of individual length  <dig> or more.
once again we choose l =  <dig> as our threshold for dme predictions.
applying dme with this threshold we obtain ec assignments at levels  <dig> and  <dig> for  <dig>  proteins.
all assignments are provided in additional file  <dig>  tables s2-s <dig>
in figure  <dig> we display a histogram of the  <dig> largest ec sub-subclasses  that emerge from our dme analysis.
the category with the largest number of different proteins is  <dig> .
<dig>  corresponding to aminoacyl-trna synthetases .
since there are about  <dig> aars enzymes expected for each organism, this allows us to estimate the content of the metagenome to be of order of  <dig> species or so.
looking at level  <dig> annotations, i.e.
at specific aars enzymes, we find that their numbers vary from  <dig> to  <dig>  these differences may be due both to different occurrences of aars sequences in the sample, and to different efficiencies of the sp methodology for different aarss.
the order of magnitude of  <dig> different species remains a reasonable estimate.
the same order of magnitude can be derived from another source.
venter et al.
<cit>  have provided some information about single copy proteins  in trying to arrive at estimates of the number of species involved.
one such protein is the gyrase subunit b enzyme, gyrb.
the same enzyme has also been proposed by watanabe et al.
<cit>  for the purpose of spanning a database for identification and classification of bacteria.
gyrb is one of several protein families belonging to ec  <dig> .
<dig>  .
checking through the sps belonging to this ec we have found a subset that is specific to gyrb only.
using this subset we estimate the number of gyrb copies in the sargasso-sea data to be  <dig>  which is close to the number of maximal fragment depth of  <dig> quoted in table  <dig> of  <cit> , and is in the same ball-park as the aars estimate.
in addition to  <dig> .
<dig>  enzymes we observe the following leading categories:  <dig> .
<dig> ,  <dig> .
<dig> ,  <dig> .
<dig> , and  <dig> .
<dig> .
there are several ec numbers  that are particularly abundant.
they are presented in table  <dig>  where we list all cases that appear more than  <dig> times in the data.
some of them have already been mentioned above: the dna gyrase, and its role in estimating the number of species, and the two ecs belonging to the susubclass of  <dig> .
<dig> , playing important roles in rna and dna polymerases.
all our predictions for the enzymatic annotations of the sargasso-sea data are presented in additional file  <dig>  tables s2-s <dig>  we wish to point out that some of the enzymes contain two or more ec assignments.
table  <dig> reports some of these occurrences.
included here are the most abundant observations of dual ec assignments, sorted by the numbers of proteins exhibiting the two annotations.
the first and the last entries in table  <dig> have many analogs in currently known doubly-annotated enzymes in swiss-prot.
checking all proteins we find that the sp hits that belong to the two different ec numbers do not overlap on the protein sequences, thus falling comfortably into the categorization of two different catalytic domains.
it is interesting to note that finding multiple domains is easier with sps than it is with blast: we will not miss out on a small domain of a protein that may be overshadowed by sequence similarities with a larger protein domain, and we can immediately check whether the different catalytic regions lie on disjoint sections of the protein.
a full list of the doubly annotated sargasso-sea enzymes is presented in additional file  <dig>  table s <dig>  a further list of triple-enzymatic annotations is presented in additional file  <dig>  table s <dig>
human gut metagenome
gill et al.
<cit>  have analyzed the dna sequences obtained from fecal dna of two healthy adults - 'subject 7' a female aged  <dig> and 'subject 8' a male aged  <dig>  we have analyzed the resulting proteins  with our dme method.
the two proteomes of subjects  <dig> and  <dig> consist of  <dig>  and  <dig>  proteins correspondingly.
we predict enzymatic annotations for  <dig>  proteins of subject  <dig> and  <dig>  proteins of subject  <dig>  these numbers are relatively lower than the enzymatic content of sargasso-sea.
numbers of  <dig> .
<dig> enzymes are predicted to be  <dig> and  <dig> for subjects  <dig> and  <dig> respectively.
thus the number of different species contained in these samples is scaled down by two-orders of magnitude compared to the sargasso-sea data, which is quite reasonable given the size of the databases.
further comparisons between the three metagenomes are offered in the next section.
enzymatic profile
trying to compare different metagenomes with each other one has obviously to resort to some normalization method.
normalizing the results of a histogram like figure  <dig> by the total number of enzymes that we find, we obtain a spectrum characteristic of the genome or metagenome we study, which we will refer to as its enzymatic profile.
in spite of the obvious similarities, there exist differences among the three histograms.
we use the absolute value of the difference of any two distributions as the difference measure .
taking into account all level- <dig> ec predictions we obtain the distances between the different distributions presented in table  <dig>  as expected, the two gut metagenomes are the closest pair.
it has been emphasized by  <cit>  and by  <cit>  that the functional characteristics of a metagenome vary with the environment in which it is being found.
hence we expect the genetic enzymatic profiles to vary accordingly.
our exercise shows that the gross features of microbial communities may be similar, thus more attention will have to be paid to smaller details, in particular emphasizing the cases where the relative differences between ec categories are the largest.
this may become a useful tool in the future.
we wish to close this section by emphasizing that the three metagenomic profiles are different from those derived from the genome of e. coli, and very different from human.
the comparisons are presented in additional file  <dig>  figure s <dig>  drawn according to the top  <dig> categories of e. coli, and additional file  <dig>  figure s <dig>  displaying the top  <dig> categories of human.
it is quite evident that the weights  of different ec categories change considerably from human to e. coli to bacterial metagenomes.
this implies that enzymatic profiles contain information that may be of value in future studies of novel genetic material.
discussion
using sps it seems quite straightforward to perform data-mining of enzymes.
there are however several provisos: a) although a majority of enzymes carry sps, there exists a minority that does not; hence not all enzymes are expected to be discovered in a new dataset.
b) sps were substantiated on a training set, and their generalization carries with it some error, even on a test set composed of enzymes only.
errors may be due to a) changes in the official ec classification of an enzyme, or b) real biological changes such as evolutionary loss of an active site in a protein that resembles a known enzyme but has no catalytic function, or c) random appearance of sps on proteins that have no catalytic activity.
errors due to reclassification of ec numbers cannot be controlled in any a-priori manner.
the question of functionality loss can be partially checked through searching for the absence of annotated sps in cases where such annotations may be expected for the enzyme in question.
this demonstrates the importance of detailed corroboration of each individual prediction of the large-scale method studied here.
the third source of errors, due to random appearance of sps on proteins other than enzymes, has been taken into account by limiting our predictions to consistent sp hits with minimal coverage length of  <dig>  and specifying the l values of our predictions as a measure of their confidence.
dme is based on deterministic motifs only, i.e.
strings with specific sequences of amino-acids.
comparing it with the well-known motif method of prosite patterns  <cit> , by using available information in swiss-prot, we find that the latter has precision of 97% and recall of only 47% on the enzyme test set, thus falling short of dme predictions.
when comparing dme to blast on the enzyme test-set we found that dme had comparable precision  while blast has much better recall .
note that this comparison was based on the 1st sp set of july  <dig>
it should be appreciated that the comparative procedure based on the enzyme test set has some bias in favor of blast, because the latter serves as one of the inputs to swiss-prot assignments.
as a result, cases of remote homology which may be captured by dme could have been missed by blast-based assignments, as was demonstrated by  <cit>  and by  <cit> .
the sp-based search has two other advantages over blast: it is conceptually simpler, relying only on a look-up table, and it points to specific locations on the queried protein which may be relevant to the expected catalytic function of that enzyme.
hence it may have wide practical implications for enzyme research and development.
in spite of all the precautions outlined in the first paragraph, our predictions concerning the  <dig> organism test-set reported in this paper, do extremely well.
moreover, note that the recall quality of sps on their training sets increased dramatically from 85% in  <dig> to 94% in  <dig> .
this means that the minority of enzymes without sp hits diminishes with time.
the reason is quite clear: mex thrives on redundancy of patterns in the data.
therefore, the more proteins of the same family there are in the database, the better mex will perform.
as these lists fill up in the swiss-prot database, they can be better represented by simple sp motifs.
higher recall on the training set will undoubtedly reflect itself also as higher recall on future test sets, thus suggesting that the gap between the recall of blast vs dme will shrink with time.
indeed, carrying out a dme analysis, based on the 2nd sp set, of  <dig>  enzymes that have been added to swiss-prot from july  <dig> to sep  <dig>   <dig>  we find on this novel test set precision of  <dig> % and recall of  <dig> %.
this is a considerable increase over the recall of 70% of the 1st sp set measured on the enzyme test set .
a straightforward peptide characterization of protein families seemed hopeless a decade or two ago, and hence necessitated the development of more sophisticated approaches such as blast, to quantify sequence similarities.
our analysis demonstrates that this has changed with time  so that nowadays the sp approach may be regarded as a useful tool, leading to valuable information.
such information, for three metagenomic data-sets, has been presented here as an example of the power of our novel methodology.
CONCLUSIONS
the requirement that sp occurrences on protein sequences has some minimal coverage length, e.g.
l ≥  <dig> amino-acids in our analyses, leads to the novel tool of dme.
it is applicable to large genomic and metagenomic data, and provides a good indicator for the enzymatic classification of the queried proteins, based on a look-up table only.
a web tool identifying sp  occurrences on any queried protein sequence, and providing the ec prediction of dme, is available online at http://adios.tau.ac.il/dme.
authors' contributions
uw developed the dme criteria, extracted the test-sets and the 2nd sp set, and performed dme analysis of genomic and metagenomic data.
yl extracted the 1st sp set, and carried out preliminary analysis of the gut metagenomic data.
dh conceived the study, participated in its design and coordination and drafted the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary tables.
table s1: list of the ten organisms used as a test-set.
table s2: list of dme predicted single ec annotations of proteins in sargasso-sea data.
table s3: list of dme predicted double ec annotations of proteins in sargasso-sea data.
table s4: list of dme predicted triple ec annotations of proteins in sargasso-sea data.
click here for file
additional file 2
supplementary figures.
figure s1: length histogram of the 2nd sp set.
figure s2: comparison of enzymatic profiles based on the  <dig> leading categories of e. coli.
figure s3: comparison of enzymatic profiles based on the  <dig> leading categories of human.
click here for file
acknowledgements
we thank uri gophna and eytan ruppin for helpful conversations.
this study was supported in part by fellowships granted to uw and yl by the edmond j. safra bioinformatics program at tel-aviv university.

 recently there has been a rapid growth in the number of putative proteins derivable from new genomic and metagenomic data. the extended use of environmental shotgun sequencing to study diverse microbial systems has made metagenomics a vastly growing field leading to a flux of data, calling for development and application of new tools that allow its investigation. conventional tools for predicting the function of a protein from its sequence are based on sequence-similarity or sequence-motifs. here we outline a relatively simple and straight-forward method that is applicable to large numbers of sequences. its purpose is finding whether each protein in the data is an enzyme and, if so, what its ec classification is. this data mining of enzymes is based on the specific peptide method of, and is carried out by comparing the sequences of all proteins with a list of all sps and looking for matches of the latter in the data. sps are strings of amino-acids, extracted from enzyme sequences using the motif extraction algorithm mex. they are selected for their specificity to levels of the enzyme commission 4-level functional hierarchy. we have updated the sp set of by extracting it from all swiss-prot enzymes of july 1st, more details are provided in methods. using sps for prediction of enzymatic function needs some further decisions as to what to do if various sp hits on the same protein have ec assignments that are not consistent with one another. moreover, one should decide when a single sp hit is sufficient to make a prediction. the methodology developed here relies on coverage length of consistent sp hits. this is further described below, when testing performance on an enzyme test set, and when discussing a ten-organism test-set that contains non-enzymatic as well as enzymatic proteins. we develop a random model for the latter to assess the effect of accidental sp matches. the resulting methodology, which we call data mining of enzymes, is being applied to analyze several metagenomes. methods the new sp sets a novel method based on sequence motifs has been proposed by, who have studied enzymes in the swiss-prot database. they have demonstrated that enzyme functions, as represented by the four-level ec hierarchy, can be deduced from the appearance of deterministic short strings of amino-acids, denoted as specific peptides, on these enzymes. the sps were derived from enzyme sequence data using an unsupervised motif extraction algorithm mex, and filtered by the ec so that each sp is specific to a particular ec branch, specifying the ec function that the enzyme performs. thus, if an extracted motif is found to occur on enzymes belonging to only one ec number, this peptide will be declared to be an sp labeled with this ec number. if, however, the motif occurs on several ec numbers, all of which share the same 3rd-level hierarchy, the motif is declared as an sp with labeling at the third level of ec hierarchy, etc. the sps of comprise on average amino-acids, and were shown to compete favorably with a smith-waterman based svm classifier. usage of the sp methodology is demonstrated by our web-tool http://adios.tau.ac.il/dme. given the sequence of an enzyme, this tool searches through the set of all sps and finds which of them coincide with substrings of the sequence, indicating where they lie, what is the ec assignment associated with each sp, and provides the ec predicted by the dme method for the protein that is being queried. kunik et al have investigated enzyme sequences of the swiss-prot release of october we have used the same methodology and applied it to all enzymes in the swiss-prot/enzyme records of july 1st, the number of enzymes that have a single ec assignment is. applying mex and filtering it by ec levels in the same way as, we have obtained sps. this new 1st list of sps serves as the basis for developing and analyzing our methodology. in making the prediction of an ec number based on one sp match, or several sp matches that have the same ec number assignment, we require that the total number of amino-acids of the protein matched with these sps be at least seven. we refer to this number as the coverage-length l. if l at level l is less than we check for sp hits that are consistent at level of the ec hierarchy, i.e. have identical first three digits in their assignments. once again, a prediction is made if l at level l is at least in principle, the threshold of l at every ec level can be viewed as a parameter of our method. reducing l increases recall at the expense of lowering precision, as will be discussed below. test data were downloaded from swiss-prot release on july 1st, we consider two types of test sets. the "enzyme test set" consists of all enzymes integrated into swiss-prot between july 1st of and the " organisms test set" consists of proteins of e. coli and other bacteria containing enzymes from the same period of to and all other proteins incorporated into swiss-prot by july 1st, a summary of all the relevant data is displayed in table it includes also information about precision and recall that will be further discussed in the first results section. these values are obtained by determining 3rd level ec assignments, using coverage-length of l ≥ precision values of 100% on the training sets are of course trivial results of specificity. 54% of the proteins in the 1st training set carry swiss-prot annotations of'active site','binding site' or'metal binding site' at specific locations of single amino-acids. sps cover these functionally important sites significantly more than other loci on proteins, thus indicating biological significance of sps. sp matches that overlap such sites are compiled, and the corresponding sps are denoted as annotated sps. we have thus compiled a list of asps. all appear at least four times in the training set, and the location of the annotation is consistent in the different appearances. most asps carry single annotations, asps carry two annotations and asps carry all three annotations. a second set of sps is extracted from swiss-prot data dated july 27th, this training set, consisting of all singly annotated enzymes, contains proteins. it has led to sps. their length distribution is presented in additional file figure s this set includes sps with labels corresponding to ec levels and. only sps with ec labels at levels and are relevant for the assignment of ec level- annotations to proteins, and hence for the calculation of recall included in table it should be emphasized that only of the swiss-prot annotated enzymes in the training set carry ec annotations at levels and they are the ones on which the ec predictions at level are tested, leading to the recall result of 94%. the 2nd sp set is being used for the analysis of metagenomic data and is incorporated in our web-tool at http://adios.tau.ac.il/dme. estimate of accidental sp matches proteins that do not possess enzymatic functions may still have a substring that matches an sp. such sp matches will be called'accidentals'. their occurrences can be modeled by sp hits on random protein sequences. such random sequences are generated from real data by scrambling the order of the amino-acids in every protein, conserving only first-order statistics. such sets were produced in order to measure the expected random hits. estimates of the probabilities of accidental occurrences of sps are derived below for the organism test-set and for sargasso sea data. recall-precision analysis of ec annotations in enzymes comparing the results of our method with an expert-method we face three possible situations when dealing with a collection of enzymes: p|p where the model prediction coincides with that of the expert, p|dp where the expert provides a different ec assignment, and np|p where the model provides no prediction for enzymes whose ec assignments are given by the expert. accordingly we define the following measures in terms of number of occurrences: this is a generalization of the common terms used in binary classification problems where p|p, p|dp and np|p are replaced by true-positive, false-positive and false-negative correspondingly. recall-precision analysis of ec annotations in proteins extending the previous analysis to a collection of proteins we have to add two more possibilities: p|np, where the new method has an ec prediction whereas the expert does not have one, and np|np where both do not have any ec assignment. whereas the latter corresponds to true-negative in a binary classifier, the former, p|np could be added to p|dp as'false-positive'. since, however, there are many cases where the absence of an ec assignment does not imply that the protein in question is not an enzyme, we opt to define a new measure, putative novelty ratio, as the fraction of such p|np out of all the predictions of the model: other measures one can define are they are the analogs of and in conventional binary classifications.  analysis of the methodology analysis of the enzyme test set using the 1st sp set in making the prediction of an ec number based on one sp match, or several sp matches that have the same ec number assignments, we require that the total number of amino-acids of the protein matched with these sps be at least seven. we refer to this number as the coverage-length l. in principle, the threshold of l at every ec level can be viewed as a parameter of our method. reducing l increases recall at the expense of lowering precision. this is exemplified in table where we analyze our enzyme test set and show precision and recall at 3rd ec level as function of the l threshold. although precision turns out to be quite high, even for low l values, recall is low when compared to what blast can achieve on this test-set. using the most significant outcome of a blast search against the 1st training set as its prediction, and limiting the most significant e-value to stay below e- we find blast precision of 98% and recall of 95%, to be compared with dme values of % and 70% when setting l ≥ thus while precision is similar, dme loses on recall. there is no direct relation between dme and blast, although high coverage-length l values of dme go usually hand in hand with very low e-values of blast. differences may occur for low l values of dme, and relatively high e-values in blast. we refer to kunik et al. for a discussion of such examples. the advantages of sps in resolving classification problems in situations of remote homology have been discussed and exemplified by. it is worthwhile pointing out that the fact that one can abide by such a small threshold value of l ≥ is strongly connected to our requirement that the sp matches on the protein's sequence be exact. if one were to allow for insertions or deletions or replacements, such as the blosum matrix, this would not work. based on various trials we may state that, whereas reliance on blosum works well for blast searches over large sequences, it ruins predictivity and specificity of sp searches even if only single amino-acid changes are allowed. analysis of the ten organism test-set the ten organism test-set contains proteins of e. coli and other bacteria listed in additional file table s proteins for this dataset were downloaded from swiss-prot on july 1st and include all proteins that had no ec annotation in swiss-prot prior to july 1st the intersection between the organism test-set and 1st training set used to build the sps is void and allows us to develop and test the sp methodology on general proteomic data rather than on enzymes only. sp search on this dataset, using our 1st set of sps, leads to the results shown in figure sorted according to the number n of sp matches. proteins have at least one sp match. some of them may be due to random hits and our task is to resolve which of the hit proteins should be recognized as enzymes and what their ec assignments should be. as before, we propose to rely on coverage length. we judge the prediction not by how many sp hits are observed, but by l the number of amino-acids matched by all sp hits whose ec assignment is identical within the first three digits of the ec number. in order to have some intuition about the expected noise level, we compare in table sp hits on real data with random model results for different values of l entries of l = refer to either no sp hits, or hits by sps that have labels with ec levels and but none at ec levels and the columns random and stdev refer to the average and standard deviation of three random sets. noise is the ratio of random/real. all proteins of the ten-organism test-set were included in this search. we will use l ≥ as our threshold criterion, as in the enzyme data-set discussed in the previous section. we note that predictions based on l = may still have a large uncertainty, however from l = onwards random hits become very small. our threshold criterion leads to the results displayed in table with precision = %, recall = %, accuracy = % and putative novelty = %. the interest in this exercise is twofold: to see how well our method performs on unassigned proteins, i.e. true-negatives, and how good our predictions are for putative novelties. indeed, our accuracy turns out to be high, %, which proves that we have correct negative assignments. seven out of the putative novelties have been annotated by swiss-prot since july six out of which are at levels or all observations are consistent with the predictions, as shown in table quoted here are also all coverage lengths on which the predictions were based. note that also the one based on coverage length has been validated. all this may be viewed as an indication of the validity of dme predictions. the first six entries in table belong to e. coli, and the last protein belongs to bacillus cereus. l to l are the coverage-lengths at ec levels to respectively classification based on annotated sps it has been noted by and that some of the sps can be demonstrated to play important biological roles since they carry crucial amino-acids known to serve as active sites, binding sites or metal binding sites. such annotations are available for 54% of the enzymes in the 1st swiss-prot training set. selecting only sps that carry these annotations we obtain a set of annotated sps, a mere 7% of all sps. we have tested it on the enzyme test set. using annotation predictions at the third level of ec we find precision % and recall %. the limited recall is due to the fact that asps have been derived from only 54% of the training set. nonetheless they possess the advantage of being selected due to their demonstrated operational importance to the catalytic function. because of their limited recall we have not used the asps as the primary tool for large scale analysis; however we list their properties in our web tool http://adios.tau.ac.il/dme. any queried protein can be analyzed by this tool for sp hits and the expected dme prediction. the appearance of asps may serve as providing additional credence to the prediction, as well as specifying the positions of expected active or binding sites. metagenomic analysis analysis of sargasso-sea data after verifying dme on the two test-sets we turn to an analysis of the, records in the sargasso sea protein data. the average length of these proteins is amino-acids, with sd = for this analysis we employ our 2nd set of sps, updated on july in order to reduce random hits, we have further limited our sp set to include only peptides of length amino-acids or more. using a random set of proteins selected from these data, we generated three randomized protein sets from which we calculated the probabilities of accidental matches. the results are displayed as function of l in table the columns random and stdev refer to the average and standard deviation of three random sets. noise is the ratio of random/real. similar results are obtained for l the results of table are slightly better than table the reason is that we have limited ourselves here to sps of individual length or more. once again we choose l = as our threshold for dme predictions. applying dme with this threshold we obtain ec assignments at levels and for proteins. all assignments are provided in additional file tables s2-s in figure we display a histogram of the largest ec sub-subclasses that emerge from our dme analysis. the category with the largest number of different proteins is. corresponding to aminoacyl-trna synthetases. since there are about aars enzymes expected for each organism, this allows us to estimate the content of the metagenome to be of order of species or so. looking at level annotations, i.e. at specific aars enzymes, we find that their numbers vary from to these differences may be due both to different occurrences of aars sequences in the sample, and to different efficiencies of the sp methodology for different aarss. the order of magnitude of different species remains a reasonable estimate. the same order of magnitude can be derived from another source. venter et al. have provided some information about single copy proteins in trying to arrive at estimates of the number of species involved. one such protein is the gyrase subunit b enzyme, gyrb. the same enzyme has also been proposed by watanabe et al. for the purpose of spanning a database for identification and classification of bacteria. gyrb is one of several protein families belonging to ec.. checking through the sps belonging to this ec we have found a subset that is specific to gyrb only. using this subset we estimate the number of gyrb copies in the sargasso-sea data to be which is close to the number of maximal fragment depth of quoted in table of, and is in the same ball-park as the aars estimate. in addition to. enzymes we observe the following leading categories:.,.,., and.. there are several ec numbers that are particularly abundant. they are presented in table where we list all cases that appear more than times in the data. some of them have already been mentioned above: the dna gyrase, and its role in estimating the number of species, and the two ecs belonging to the susubclass of., playing important roles in rna and dna polymerases. all our predictions for the enzymatic annotations of the sargasso-sea data are presented in additional file tables s2-s we wish to point out that some of the enzymes contain two or more ec assignments. table reports some of these occurrences. included here are the most abundant observations of dual ec assignments, sorted by the numbers of proteins exhibiting the two annotations. the first and the last entries in table have many analogs in currently known doubly-annotated enzymes in swiss-prot. checking all proteins we find that the sp hits that belong to the two different ec numbers do not overlap on the protein sequences, thus falling comfortably into the categorization of two different catalytic domains. it is interesting to note that finding multiple domains is easier with sps than it is with blast: we will not miss out on a small domain of a protein that may be overshadowed by sequence similarities with a larger protein domain, and we can immediately check whether the different catalytic regions lie on disjoint sections of the protein. a full list of the doubly annotated sargasso-sea enzymes is presented in additional file table s a further list of triple-enzymatic annotations is presented in additional file table s human gut metagenome gill et al. have analyzed the dna sequences obtained from fecal dna of two healthy adults -'subject 7' a female aged and'subject 8' a male aged we have analyzed the resulting proteins with our dme method. the two proteomes of subjects and consist of and proteins correspondingly. we predict enzymatic annotations for proteins of subject and proteins of subject these numbers are relatively lower than the enzymatic content of sargasso-sea. numbers of. enzymes are predicted to be and for subjects and respectively. thus the number of different species contained in these samples is scaled down by two-orders of magnitude compared to the sargasso-sea data, which is quite reasonable given the size of the databases. further comparisons between the three metagenomes are offered in the next section. enzymatic profile trying to compare different metagenomes with each other one has obviously to resort to some normalization method. normalizing the results of a histogram like figure by the total number of enzymes that we find, we obtain a spectrum characteristic of the genome or metagenome we study, which we will refer to as its enzymatic profile. in spite of the obvious similarities, there exist differences among the three histograms. we use the absolute value of the difference of any two distributions as the difference measure. taking into account all level- ec predictions we obtain the distances between the different distributions presented in table as expected, the two gut metagenomes are the closest pair. it has been emphasized by and by that the functional characteristics of a metagenome vary with the environment in which it is being found. hence we expect the genetic enzymatic profiles to vary accordingly. our exercise shows that the gross features of microbial communities may be similar, thus more attention will have to be paid to smaller details, in particular emphasizing the cases where the relative differences between ec categories are the largest. this may become a useful tool in the future. we wish to close this section by emphasizing that the three metagenomic profiles are different from those derived from the genome of e. coli, and very different from human. the comparisons are presented in additional file figure s drawn according to the top categories of e. coli, and additional file figure s displaying the top categories of human. it is quite evident that the weights of different ec categories change considerably from human to e. coli to bacterial metagenomes. this implies that enzymatic profiles contain information that may be of value in future studies of novel genetic material. discussion using sps it seems quite straightforward to perform data-mining of enzymes. there are however several provisos: a) although a majority of enzymes carry sps, there exists a minority that does not; hence not all enzymes are expected to be discovered in a new dataset. b) sps were substantiated on a training set, and their generalization carries with it some error, even on a test set composed of enzymes only. errors may be due to a) changes in the official ec classification of an enzyme, or b) real biological changes such as evolutionary loss of an active site in a protein that resembles a known enzyme but has no catalytic function, or c) random appearance of sps on proteins that have no catalytic activity. errors due to reclassification of ec numbers cannot be controlled in any a-priori manner. the question of functionality loss can be partially checked through searching for the absence of annotated sps in cases where such annotations may be expected for the enzyme in question. this demonstrates the importance of detailed corroboration of each individual prediction of the large-scale method studied here. the third source of errors, due to random appearance of sps on proteins other than enzymes, has been taken into account by limiting our predictions to consistent sp hits with minimal coverage length of and specifying the l values of our predictions as a measure of their confidence. dme is based on deterministic motifs only, i.e. strings with specific sequences of amino-acids. comparing it with the well-known motif method of prosite patterns, by using available information in swiss-prot, we find that the latter has precision of 97% and recall of only 47% on the enzyme test set, thus falling short of dme predictions. when comparing dme to blast on the enzyme test-set we found that dme had comparable precision while blast has much better recall. note that this comparison was based on the 1st sp set of july it should be appreciated that the comparative procedure based on the enzyme test set has some bias in favor of blast, because the latter serves as one of the inputs to swiss-prot assignments. as a result, cases of remote homology which may be captured by dme could have been missed by blast-based assignments, as was demonstrated by and by. the sp-based search has two other advantages over blast: it is conceptually simpler, relying only on a look-up table, and it points to specific locations on the queried protein which may be relevant to the expected catalytic function of that enzyme. hence it may have wide practical implications for enzyme research and development. in spite of all the precautions outlined in the first paragraph, our predictions concerning the organism test-set reported in this paper, do extremely well. moreover, note that the recall quality of sps on their training sets increased dramatically from 85% in to 94% in. this means that the minority of enzymes without sp hits diminishes with time. the reason is quite clear: mex thrives on redundancy of patterns in the data. therefore, the more proteins of the same family there are in the database, the better mex will perform. as these lists fill up in the swiss-prot database, they can be better represented by simple sp motifs. higher recall on the training set will undoubtedly reflect itself also as higher recall on future test sets, thus suggesting that the gap between the recall of blast vs dme will shrink with time. indeed, carrying out a dme analysis, based on the 2nd sp set, of enzymes that have been added to swiss-prot from july to sep we find on this novel test set precision of % and recall of %. this is a considerable increase over the recall of 70% of the 1st sp set measured on the enzyme test set. a straightforward peptide characterization of protein families seemed hopeless a decade or two ago, and hence necessitated the development of more sophisticated approaches such as blast, to quantify sequence similarities. our analysis demonstrates that this has changed with time so that nowadays the sp approach may be regarded as a useful tool, leading to valuable information. such information, for three metagenomic data-sets, has been presented here as an example of the power of our novel methodology.  the requirement that sp occurrences on protein sequences has some minimal coverage length, e.g. l ≥ amino-acids in our analyses, leads to the novel tool of dme. it is applicable to large genomic and metagenomic data, and provides a good indicator for the enzymatic classification of the queried proteins, based on a look-up table only. a web tool identifying sp occurrences on any queried protein sequence, and providing the ec prediction of dme, is available online at http://adios.tau.ac.il/dme. authors' contributions uw developed the dme criteria, extracted the test-sets and the 2nd sp set, and performed dme analysis of genomic and metagenomic data. yl extracted the 1st sp set, and carried out preliminary analysis of the gut metagenomic data. dh conceived the study, participated in its design and coordination and drafted the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 supplementary tables. table s1: list of the ten organisms used as a test-set. table s2: list of dme predicted single ec annotations of proteins in sargasso-sea data. table s3: list of dme predicted double ec annotations of proteins in sargasso-sea data. table s4: list of dme predicted triple ec annotations of proteins in sargasso-sea data. click here for file additional file 2 supplementary figures. figure s1: length histogram of the 2nd sp set. figure s2: comparison of enzymatic profiles based on the leading categories of e. coli. figure s3: comparison of enzymatic profiles based on the leading categories of human. click here for file acknowledgements we thank uri gophna and eytan ruppin for helpful conversations. this study was supported in part by fellowships granted to uw and yl by the edmond j. safra bioinformatics program at tel-aviv university.
using a training set of july and test sets of july we find that the predictive power of sps, both for true-positives and true-negatives, depends on the coverage length of all sp matches. we employ the novel motif method that consists of specific peptides that are unique to specific branches of the enzyme commission functional classification. moreover, sifting through our predictions we find that some of them have been substantiated by swiss-prot annotations by july finally we extract, for production purposes, a novel sp set trained on all swiss-prot enzymes as of july this new set increases considerably the recall of dme. comparing the two on an enzyme test set of july we find that dme has lower recall. in our analysis, l ≥ has led to highly accurate results. the outcome of these analyses can be characterized by the enzymatic profile of the metagenomes, describing the relative numbers of enzymes observed for different ec categories. the new sp set is being applied to three metagenomes: sargasso sea with over, proteins, producing predictions of over enzymes, and two human gut metagenomes.
 predicting the function of a protein from its sequence is a long-standing challenge of bioinformatic research, typically addressed using either sequence-similarity or sequence-motifs. we employ the novel motif method that consists of specific peptides that are unique to specific branches of the enzyme commission functional classification. we devise the data mining of enzymes methodology that allows for searching sps on arbitrary proteins, determining from its sequence whether a protein is an enzyme and what the enzyme's ec classification is.  we extract novel sp sets from swiss-prot enzyme data. using a training set of july and test sets of july we find that the predictive power of sps, both for true-positives and true-negatives, depends on the coverage length of all sp matches. dme is quite different from blast. comparing the two on an enzyme test set of july we find that dme has lower recall. on the other hand, dme can provide predictions for proteins regarded by blast as having low homologies with known enzymes, thus supplying complementary information. we test our method on a set of proteins belonging to bacteria, dated july establishing the usefulness of the coverage-length cutoff to determine true-negatives. moreover, sifting through our predictions we find that some of them have been substantiated by swiss-prot annotations by july finally we extract, for production purposes, a novel sp set trained on all swiss-prot enzymes as of july this new set increases considerably the recall of dme. the new sp set is being applied to three metagenomes: sargasso sea with over, proteins, producing predictions of over enzymes, and two human gut metagenomes. the outcome of these analyses can be characterized by the enzymatic profile of the metagenomes, describing the relative numbers of enzymes observed for different ec categories.  employing sps for predicting enzymatic activity of proteins works well once one utilizes coverage-length criteria. in our analysis, l ≥ has led to highly accurate results. 