BACKGROUND
in order to provide rapid and sensitive annotation for protein sequences, including direct links to structural and functional information, the national center for biotechnology information  initiated the conserved domain database   <cit>  —a collection of position-specific scoring matrices   that are derived from protein multiple sequence alignments.
as a result, web-based blast searches now include a search of the cdd, which allows users to visualize multiple sequence alignments and  structures of proteins sharing significant homology to the query and, within those alignments, key catalytic and ligand-binding residues.
thus blast searches linked to the cdd provide additional clues to the function and underlying mechanism of the query protein and are thereby often more informative, faster and more sensitive than searching against millions of individual protein sequences.
the cdd is comprised of domain models either manually curated at the ncbi or imported from other alignment collections such as pfam  <cit> , smart  <cit> , and tigrfam  <cit> .
a central and unique feature of the cdd is that related domains are organized into hierarchies when evidence exists to support that tree as a representation of the molecular evolution of the protein class.
a significant bottleneck in the cdd pipeline is the curation of these hierarchies and the manual annotation of the corresponding profiles for functionally important residues .
in order to begin automating this process, here we describe statistically-rigorous procedures for automated creation of conserved domain  hierarchies and for annotation of the corresponding profile alignments.
these procedures do automatically, based on objective, empirical criteria, what the cdd resource group and similar groups currently do manually, based on classifications that have been established in the published literature, on phylogenetic and structural analysis and, to some degree, on subjective judgments.
our focus here is to obtain heuristically an initial  cd hierarchy starting from a typically very large multiple sequence alignment for an entire protein class whose domain boundaries remain fixed.
to do this we utilize procedures that obtain both subgroup assignments for aligned sequences and corresponding discriminating patterns associated with protein functional-divergence.
the automated annotation of functionally critical residues is an important outcome of these proposed procedures: just as a large enzyme class conserves residues directly involved in catalysis, protein subgroups conserve residues likely involved in subgroup-specific biochemical properties and mechanisms.
our procedures use statistical criteria to glean this biochemical information from patterns of divergent residues among related sequences in a manner similar to the use, by classical geneticists, of statistical criteria to glean information from patterns of divergent traits among related individuals.
.
by mapping various categories of pattern residues to corresponding pssms, blast searches against these improved cd profiles can reveal those residues most likely responsible for the specific biochemical and biophysical properties of a query protein.
this can accelerate the pace of biological discovery by enabling researchers to obtain valuable clues regarding as-yet-unidentified protein properties through routine web-based blast searches.
other methods that may be similarly described as addressing the protein subfamily classification problem find sequence clusters either based on pairwise similarity  <cit>  or by cutting phylogenetic trees  <cit> .
.
here we take a different approach, namely the hierarchical classification of a protein  class based on functionally-divergent residue signatures.
unlike our approach, many existing methods, though not all , generally focus on the narrower problem of identifying orthologs or on the broader problem of clustering a database into unrelated protein classes rather than on constructing a hierarchy of domain profiles for a specific protein class.
an approach, which is, in certain respects, similar to the one described here , is the statistical coupling analysis method of lockless and ranganathan  <cit>  for detecting sets of correlated residues in protein sequences  <cit> .
because our approach identifies residues associated with protein functional divergence, it is also related to "functional subtype" prediction  methods  <cit> , but is distinct inasmuch as these related methods typically predict specific residue functions  that are sufficiently well-understood to allow benchmarking  <cit> .
instead, our approach lets the data itself reveal its most statistically striking properties without making assumptions about the types of residues to be identified.
it is further distinguished from each of these related methods in at least several of the following respects:  it does not require that the input alignment be partitioned into divergent subsets beforehand; this is unlike many  <cit> , though not all  <cit>  fsp methods.
it has a rigorous statistical basis.
.
it is designed for very large input alignments.
for optimization it relies on bayesian sampling, which has a solid scientific basis  <cit> .
.
it separates out unrelated and aberrant sequences automatically.
it can identify multiple categories of co-conserved residues within a given protein.
it addresses concurrently the problems of protein subfamily classification and of identifying residues associated with protein functional divergence.
and  it can accomplish all of this automatically starting from a single multiple sequence alignment.
problem definition and solution strategy
here we address the following biological and algorithmic problem: we are given as input a  multiple sequence alignment corresponding to a particular protein class.
our objective is to partition this alignment into a tree of sub-alignments, termed a cd hierarchy, each subtree of which corresponds to a sub-alignment of sequences sharing a certain pattern that most distinguishes them from those sequences associated with the parent node of the subtree and with any other subtrees attached to that parent node.
we interpret these distinguishing residue signatures as associated with functional divergence of the protein class.
as mentioned above, our focus is on obtaining an initial, suboptimal hierarchy that is comparable to current cdd curated hierarchies and that can serve as a starting point for further optimization using either manual or automated methods.
here we describe statistically-based heuristic procedures that, in conjunction with bayesian sampling, can obtain such an initial hierarchy from a multiple sequence alignment.
bayesian sampling over contrast alignment models
our approach relies on bayesian markov chain monte carlo  sampling  <cit> , which starts with an arbitrary model having a certain  probability that it  could have generated the input sequence alignment data.
then, in successive iterative steps, a number of alternative values for a given parameter of the current model are evaluated .
a value for this parameter is then sampled proportional to its probability.
this iterative process continues until convergence on the most likely models.
here we optimize in this way contrast alignment models, each of which consists of a pattern and a set of labels assigning each sequence to either a foreground partition or a background partition corresponding to sequences that either generally conserve or fail to conserve the pattern, respectively.
a schematic representation of a contrast alignment and the corresponding probability distribution are described in figure  <dig>  to sample alternative contrast alignment models we use a mcmc sampling strategy  <cit> , termed bayesian partitioning with pattern selection   <cit> .
the sampler converges on a model where the pattern best distinguishes the foreground from the background sequences.
multiple category functional divergence models
more recently, a multiple category bpps sampler was developed  <cit>  with a view to optimally assigning aligned sequences to various nodes within a predefined protein domain hierarchy based on functionally-divergent residue signatures.
thus the mcbpps sampler aims to precisely define both the sequences belonging to each subgroup and the patterns most distinctive of each subgroup within a specific protein class.
however, because the mcbpps sampler does not define the hierarchy of contrast alignments, it requires that the user provide  both a functional divergence -table  and seed sequences for each divergent subgroup.
sequences are assigned to subgroups through bayesian inference).
each row of a fd-table corresponds to a distinct functionally divergent subgroup of the input sequences and each column corresponds to a distinct contrast alignment whose foreground and background partitions are specified by the symbols in the table.
such a table is shown in figure  <dig>  which also illustrates the correspondence between a tree representing the hierarchical relationships between functionally divergent subgroups and the fd-table and between a column in the table and the contrast alignment; these are shown above and below the table, respectively.
given the relationships specified by the fd-table, the mcbpps sampler stochastically reassigns aligned sequences to alternative subgroups and alternative patterns to each foreground partition until convergence on an optimal  set of contrast alignments.
modeling the functional divergence of an entire protein class in this way is substantially more powerful than using a single contrast alignment because:  in principle, it can optimally model every  divergent subgroup within an entire protein class concurrently.
it sets up a stringent competition between functionally-divergent categories for pattern residues, thereby defining each pattern and partition much more precisely.
it eliminates problematic sequences, which would otherwise tend to obscure analyses, by modeling them explicitly.
problematic sequences include, for example, related proteins that have further functionally diverged to become outliers, pseudogene products and other non-functional proteins, and unrelated or erroneous sequences.
and , by defining multiple categories of pattern residues within individual proteins it can reveal, in the light of available structural information, functionally important residue interactions.
.
here we describe and apply an automated multiple category bpps program that generates its own fd-table and seed sequences automatically and therefore merely requires a multiple sequence alignment as input.
the number and nature of the partitions and the patterns is completely determined by the program.
when used in conjunction with procedures for viewing structural interactions involving pattern residues, the amcbpps sampler automates and enhances the creation and annotation of cdd hierarchical alignments.
and, when linked into web-based blast searches, this can make previously inaccessible molecular information widely available.
RESULTS
in this section, we lay out the basic amcbpps algorithm, illustrate an implementation of the algorithm as applied to p-loop gtpases, compare its performance against various manually-curated cd hierarchies, further evaluate its performance using both delete-half jackknifing and simulations, and apply it to several large protein classes for which existing hierarchies or alignments are currently unavailable.
algorithm
the amcbpps algorithm aims to identify the hierarchical relationships between functionally-divergent subgroups within an entire protein domain class based on the differentiating patterns present in that class.
it does this by defining:  the number of sequence sets,  the members of each set,  the hierarchical relationships between sets and  the corresponding functionally divergent patterns.
this is accomplished in three steps.
steps  <dig> and  <dig> constitute the novel aspect of the program by providing input to the mcbpps sampler in step 3; these first two steps also speed up convergence in step  <dig> by providing a better starting point for the mcbpps sampler .
conceptual aspects of the algorithm corresponding to steps  <dig> and  <dig> are illustrated in figure  <dig> .
these first two steps are performed heuristically and thus constitute an informed guess on how to best model the protein class .
the sampler then improves upon this model by optimizing sequence and pattern assignments within this hierarchy.
the hierarchy returned by the amcbpps program may be edited  and then optimized again by the mcbpps program.
such editing may, for example, further subcategorize previously identified or miscellaneous subgroups.
identifying simple subgroups
step  <dig> of the amcbpps algorithm  first generates a forest of simple  trees, the leaves of which correspond both to functionally-divergent subgroups within the protein domain class and to sub-alignments of the input alignment.
this is accomplished by obtaining seed sequences and using them to create these trees, as follows:  all closely-related sequences  from the same phylum are clustered into a common set .
all pairs of moderately related sequences  from distinct phyla are stored on a heap  using their pairwise scores as the key.
iteratively remove the top-scoring cross-phylum pair from the heap and merge their two disjoint sets into one set.
.
once a disjoint set contains sequences from a pre-defined minimum number of distinct phyla , the sequences of highest rank from each phylum are used to seed a new subgroup; the disjoint set is then labeled to avoid picking this subgroup repeatedly.
keep generating subgroups in this way until a pre-defined number of seed sequence sets are obtained , at which point a simple fd-table is constructed where each subgroup node is a direct descendent of the root node.
.
and  repeat substeps iii-v until all sequence pairs have been removed from the heap.
to ensure that the subgroups are sufficiently diverse, we require that each seed set consensus sequence share less than a specified level of sequence identity with other seed set consensus sequences .
taken together, these sub-steps favor the identification of the most conserved and phylogenetically diverse subgroups.
for each of these fd-tables  the mcbpps sampler assigns each of the multiply aligned input sequences to a subgroup  and determines the differentiating conserved pattern for each contrast alignment .
to ensure that subgroups at different levels of the hierarchy are identified, the algorithm performs multiple runs using various numbers of leaf nodes and various prior probability settings for p, p and p .
convergence on protein subfamilies is favored by specifying a high number of leaf nodes , by lowering the  prior probability for assigning a sequence to a leaf node, p, , by setting the  prior probability, p, to favor a lower degree of background contamination by assigning more pseudo-observations to pattern matching residues and fewer pseudo-observations to contaminating residues  and by raising the  prior probability that a column corresponds to a pattern position, p .
the rationale for choosing these settings is that, for subfamilies, membership is more exclusive, sequences are more highly conserved and, consequently, conserved patterns more extensive.
in order to maximize the influence of the data on model optimization.)
convergence on a super-family is favored by specifying a single subgroup and by altering these prior parameter settings accordingly .
default settings are based on applications to actual protein sequences, though it should be noted that the influence of these prior settings is minor.
hence these priors primarily function as tuning parameters to help gently guide the sampler into finding a variety of functionally divergent subgroups.
to avoid finding the same subgroup repeatedly, sequences assigned to a subgroup in a previous run are prohibited from being used as seeds in subsequent runs.
subfamilies can also be identified recursively; that is, by rerunning the program on a single subgroup in order to find subgroups within subgroups .
the pseudocode for this step of the amcbpps algorithm is given in methods.
defining a hierarchy for the protein class
once individual subgroup sets are identified in step  <dig> , the program hierarchically arranges these into a more complex tree, from which a fd-table is obtained.
it does this using efficient bitwise set operations  <cit> , standard network algorithms  <cit>  and pattern-based statistical criteria, namely the contrast alignment log-likelihood ratio  used by the mcbpps sampler  <cit>  .
these hierarchically-arranged subgroup sets are ‘fuzzy’ due to the uncertainty associated with set membership .
thus step  <dig> of the algorithm  determines which of the input sets :  are the same set;  are distinct sets; or  are supersets of another set or sets.
step  <dig> is subdivided into three sub-steps:  merge each collection of subgroup sets deemed to be identical into a single set;  cluster related subgroup sets into common supersets; and  create a tree representation of the subgroup hierarchy, which is done using a breadth-first scanning algorithm  <cit>  to find a shortest path tree  <cit> .
step 2c also refines the tree to eliminate inappropriate overlap between sets while also eliminating nodes from the tree that, as a result of this refinement process, are no longer statistically significant.
the pseudocode for step  <dig> is given in methods.
from this tree a fd-table is then generated as input to step  <dig>
the mcbpps sampler  and further refinements
the output from step  <dig> provides a starting point for mcbpps sampling, which optimizes the patterns and partitions corresponding to the fd-table.
the basic statistical and algorithmic aspects of the mcbpps sampler were previously described  <cit> .
to further expand a cd hierarchy the output files obtained from an initial amcbpps analysis can also be used to recursively analyze, in the same way, several of the larger subgroups.
to do this, the output alignment file for a major subgroup is used as an input file for the amcbpps program.
likewise, a cd hierarchy can also be refined by editing the fd-table manually and then applying the mcbpps sampler, as was previously described  <cit> .
to speed up analysis of a given subtree, such manually edited fd-tables  may be designed to expand subgroups within that subtree, while modeling the other branches of the hierarchy only at the highest levels .
in keeping with the automated theme of this article, however, we will not describe in detail nor apply these approaches here.
implementation and testing
the amcbpps algorithm was implemented in c++ , applied to various protein classes and the output compared to manually-curated cdd alignment hierarchies .
a wide range of cdd hierarchies—from preliminary to well-developed releases —were examined in this way.
input multiple alignments were obtained by using the ncbi hierarchy of cd alignments as input to the mapgaps program  <cit> , which detected and aligned related protein sequences within the ncbi nr, env_nr and translated est protein databases.
to obtain input alignments corresponding to large protein classes for which cdd hierarchies are not yet available we used alternative procedures, as described in methods.
illustrative example: p-loop gtpases.
to familiarize the reader we begin by illustrating our approach with an analysis of p-loop gtpases.
using an input alignment of  <dig>  p-loop gtpases, the amcbpps program returned the fd-table shown in figure  <dig>  .
it also returns a corresponding set of contrast alignments, which highlight the pattern residues identified by the sampler; one such alignment is shown in .
note that the sampler will reject heuristically proposed subgroups whose existence is not supported by the data .
further subdivision of the hierarchy in figure  <dig> may be accomplished by recursively applying the amcbpps sampler to a previously-identified subgroup.
illustrates this for the ras-like gtpases by showing an expanded subtree corresponding to the column  <dig> foreground partition in figure  <dig>  by applying the amcbpps sampler recursively in this way, a very extensive hierarchy may be obtained.
‡ after removing identical sequences and sequences that fail to align with at least 75% of the domain.
* numbers in parentheses indicate the nodes retained after insignificant nodes were removed by the mcbpps program.
† the log-likelihood ratio in nats.
§ the time  is for steps  <dig> and  <dig> of the algorithm only; step  <dig> can be parallelized to run in less than 10% of the time shown.
criteria for comparing hierarchies
to assess how well the amcbpps program performs relative to curated cd hierarchies, we compared its output against  <dig> manually curated cdd hierarchies .
before considering this analysis, however, we first need to discuss the criteria used to evaluate and compare hierarchies.
lack of gold standards
cdd hierarchies have been carefully constructed by expert curators and therefore come the closest to a benchmark set for evaluating the amcbpps sampler.
however, as this study reveals, certain aspects of cdd hierarchies lack statistical support or are incomplete or incorrect for various reasons: for example, cdd hierarchies are typically at different stages in an ongoing refinement process, and, for protein domain classes consisting of tens or hundreds of thousands of sequences, the number of possible hierarchies to consider is astronomical, which makes optimization through manual curation extremely difficult.
furthermore, due to the stochastic nature of and the inability to directly observe evolutionary divergence, it is impossible to eliminate the inherent uncertainties associated with protein classification.
hence, for the present study our aim is merely to replicate the current manual curation process by generating hierarchies of comparable quality automatically, thereby dramatically speeding up the current labor-intensive curation process.
comparison criteria for this analysis
despite the absence of a gold standard, the statistical criteria used by the amcbpps program provide a way to compare two hierarchies for the same conserved domain.
it does this by determining objectively whether or not  the sequences in each protein subgroup have diverged from the evolutionarily related subgroups indicated by a specific hierarchy.
this measure is expressed as a log-likelihood ratio , where non-positive values indicate a lack of statistical support for a functionally divergent event within the hierarchy.
such a comparison is performed as follows: we are given two heuristic methods for obtaining a  hierarchy: one manual and one automatic.
to compare the two methods, we first use each hierarchy  as input to the mcbpps sampler, which then optimizes the patterns and sequence partitions associated with that hierarchy and returns an optimized log-likelihood ratio .
because this optimizes the automatically-generated and manually-curated hierarchies in the same way based on the same statistical criteria, the only difference is that the hierarchies and seed alignments were obtained either automatically or through manual curation.
thus, by comparing their optimized llr scores, we can obtain a measure of the relative performance of the two methods.
in addition, we also determine the degree of overlap between the two hierarchies as a qualitative indication of the similarity of the two hierarchies.
the results of such comparisons are summarized in table  <dig> and figure  <dig>
evaluation of the amcbpps program
to evaluate the amcbpps program over a wide variety of input, we chose the  <dig> conserved domains given in table  <dig>  these domains vary in the numbers of members detected in the protein databases , in the lengths of their conserved core , and in the size  and complexity  of their curated hierarchies.
comparisons with manually curated cdd hierarchies
based on the llr statistic the automatically-generated hierarchies  are comparable to the corresponding manually-curated hierarchies  and are, in fact, slightly better on average .
manual and amcbpps hierarchies  were also compared qualitatively by determining the degree to which the sequence sets corresponding to the nodes in one hierarchy overlap with those of the other hierarchy; in additional file 1: figure s <dig> illustrates such a comparison for the pi3kc_like domain hierarchy .
this provides a detailed comparison of two hierarchies by computing how the sequences at each level of one hierarchy are assigned to each level of the other hierarchy and vice versa.
such comparisons indicate that the differences between the curated and amcbpps hierarchies are mainly due to the following:  a node in one hierarchy being modeled as a subtree in the other.
additional child nodes being added to parent nodes in one hierarchy but not the other.
a subtree in one hierarchy being split into unrelated subtrees  in the other due to a failure to join these to a common internal  node.
and  inherently ambiguous sequences that can’t be clearly assigned to a specific node in the hierarchy; such sequences may correspond to pseudogene products or functionally defective members of a protein class that are difficult to categorize because they harbor degenerate subgroup patterns.
of course, the larger and more functionally divergent hierarchies are more challenging.
most of the differences between the manual and automated hierarchies are due to fundamental differences between the two approaches .
curated hierarchies may rely, in part, on information that is  ignored by the amcbpps program, such as subfamily-conserved inserts and 3d structures—though, on the other hand, the amcbpps program utilizes a far greater amount of sequence data that is also up-to-date.
in contrast, some of the cdd heirarchies may be out-of-date or still incomplete.
the amcbpps program also requires each node to  correspond to at least  <dig> sequences  in order to avoid statistical biases due to small sample size.
cdd curators, however, may construct subgroups containing fewer sequences.
likewise, the amcbpps program selects seed sequences from at least three or four distinct phyla in order to avoid sampling biases introduced by orthologous sequences from closely related organisms.
hence, it will fail to identify subgroups that only occur in vertebrates, for instance.
in contrast, cdd curators may choose representative sequences  from more closely related taxa.
due to such restrictions, an amcbpps-generated hierarchy tends to have fewer nodes  because it is prohibited from identifying certain cdd-defined subgroups.
this also tends to lower the llr, which  increases with the number of nodes in the hierarchy.
despite these differences, after taking these considerations into account, we found the amcbpps hierarchies to be persuasively consistent with the corresponding cdd hierarchies .
a perceived “unfair” advantage of the amcbpps algorithm might be that it utilizes the same statistical model to construct a hierarchy that is used to score that hierarchy, whereas curators do not.
however, subsequent  optimization of both hierarchies using the mcbpps sampler should counteract this advantage.
that is, assuming that the curated hierarchy is in fact superior and that our statistical model is biologically meaningful, then optimal partitioning of the sequences and optimal pattern assignment by the sampler for both types of hierarchies should result in a superior llr score for the curated hierarchy.
unsurprisingly, our analysis also indicates that the hierarchies obtained both manually and automatically are typically suboptimal.
for example, manual and amcbpps hierarchies for the s8/s <dig> peptidase domain  had llrs of  <dig>  and  <dig>  nats, respectively, whereas a hybrid hierarchy containing features of both of these has a llr of  <dig>  nats.
figure  <dig> likewise illustrates the construction of a hybrid hierarchy for the class i aars-like core domain , which improves the llr from  <dig>   to  <dig>  nats.
examining llr statistics can suggest other ways in which to improve cdd hierarchies.
for example, within the glo_edi_brp_like hierarchy , the mcbpps sampler rejected an intermediate node  to which it had assigned a llr of − <dig> nats.
an investigation to determine why this occurred revealed that, based as well on the criteria used by the cdd curators, the cd <dig> intermediate node is incompatible with several of its leaf nodes, which are therefore better modeled as direct descendents of the root node.
for these and other domains in table  <dig>  suggested improvements in cdd hierarchies based on llrs were corroborated through manual inspection by the cdd resource group.
delete-half jackknife analyses
a bootstrap or jackknife  <cit>  procedure can be used to estimate confidence levels for evolutionary trees  <cit> .
however, applying this approach to a cd hierarchy is complicated by the potential run-to-run variability both in the number of the leaf nodes and in the associated sequence sets.
thus existing evolutionary tree bootstrap and jackknife procedures, which require that each of the sampled trees use the same set of leaf nodes, cannot be used.
instead, we implemented a delete-half jackknife procedure that—though unable to provide quantitative confidence levels for specific features of a hierarchy—can nevertheless provide a qualitative assessment of the run-to-run variability of amcbpps-generated hierarchies.
this involved running the amcbpps program on each of  <dig> different input alignments  after randomly removing half of the sequences.
given ten such runs for each domain, we compared the consistency of the resultant hierarchies from run to run as follows: for each pair of runs  we determined how those input sequences shared by both hierarchies  were partitioned among the nodes of one hierarchy relative to the other hierarchy.
an example output file is shown in .
for these analyses we found that, among the leaf node sets in one tree that share at least one sequence in common with a leaf node set in the other tree, on average 47% share precisely the same set of sequences  and 74% share more than 90% of their sequences in common.
moreover, in most cases where an identical sequence set is not found, the missing sequences were typically assigned, not to unrelated leaf nodes, but either to a parent node further up the tree or to the rejected sequence set.
among the remaining cases, a node in one hierarchy is either split into multiple nodes or  split between nodes in the other hierarchy.
at times a hierarchy could end up omitting certain nodes due to the delete-half jackknife procedure removing sequences belonging to certain phyla resulting in insufficient phylogenetic diversity to seed the formation of a subgroup.
of course the topologies  of the jackknife trees found by the sampler also differ, which is a common problem associated with evolutionary trees consisting of large numbers of distantly related sequences.
this is presumably due in large part to the amcbpps algorithm failing to find the optimal topology—an issue that, in the future, we will address by sampling over alternative topologies.
of course, both this future sampler and the jackknife procedure applied here will be useful for identifying the most reliable features of a hierarchy.
taken together, these results confirm the observation we made in the previous section, namely that the amcbpps program generally finds a suboptimal hierarchy that, nevertheless, provides a good starting point both for curation and further automation.
output from these jackknife analyses are available at http://www.chain.umaryland.edu/amcbpps/jackknife.txt.
simulations
as an additional check, we implemented a procedure to generate simulated sequences from profile hmms where each such profile corresponds to a node from one of the  <dig> domain hierarchies used in the jackknife analysis.
the rationale for doing this was to determine how well the amcbpps program identifies sequences corresponding to predefined subgroups.
note that this procedure captures sequence features of each subgroup, but not how those subgroups are hierarchically arranged.
for each node of each hierarchy we generated the same number of aligned sequences as were assigned to that node in the original hierarchy.
after running the amcbpps program on each of these simulated alignments, we determined the degree to which each set of related simulated sequences were correctly modeled as belonging to a single subgroup.
an example output file in  illustrates how the amcbpps program correctly categorized nearly all of these sequences given the structure of the inferred hierarchy .
output from simulations for the  <dig> domains is available at http://www.chain.umaryland.edu/amcbpps/simulate.txt.
time complexity
the computationally most intensive routine in step  <dig> of the amcbpps program is an all-versus-all pairwise comparison of pre-aligned sequences .
this has a time complexity of o = o where k is the number of aligned columns, m is the number of sequences and n = k·m is the effective size of the input alignment.
in addition, step  <dig> involves a simpler version of the step  <dig> algorithm that, of course, exhibits the same time complexity as step  <dig>  as well as other operations that perform better than o   <cit> .
hence the time complexity for step  <dig> is o.
the time complexity of steps 2– <dig> is unclear based on the underlying algorithm.
therefore, using a plot of the run times for the amcbpps analyses in table  <dig> versus the size of the input alignments, we estimate that the time required for steps 2– <dig> scales as o .
assuming that the asymptotic time complexity for steps 2– <dig> is indeed o, which admittedly may not be the case given our empirically-based approach, then whether or not o is better than o depends on the ratio of k to m <dig>   step  <dig>  which is o, performs asymptotically worse when k < m <dig> and steps 2– <dig>  which is o, is worse when k > m <dig>  since for essentially all protein domains k < m <dig> the time complexity of the amcbpps program  appears to be o. it is important to note, however, that step  <dig>  required less time than steps 2– <dig> in our analyses—even on the largest input alignments—suggesting that constant factors rather than asymptotics are influencing program performance.
because the run times for steps 2– <dig> are also likely to depend on the size of the hierarchy generated by the program in step  <dig>  figure 7b plots the run times versus the number of aligned residues times the number of nodes in the hierarchy .
this yields a slightly improved, essentially linear dependency.
this observed time complexity is largely due to step  <dig> being more or less independent of the number of nodes, which is achieved by computing conditional posterior probabilities  for each column of the fd-table only when considering the assignment of a sequence to one of two possible new partitions rather than to one of a typically much larger number of rows.
thus the amcbpps program can be applied to very large multiple sequence alignments, which is important given the current rapid increase in sequence data.
analysis of protein domains lacking cd hierarchies
there are a significant number of protein domains for which a cdd hierarchy has not yet been constructed.
in some  cases a single curated alignment is available as a starting point.
to test the performance of the amcbpps program in such cases, we chose  <dig> domains, for which curated alignments were available, and two domains, for which we first constructed an alignment using bayesian multiple alignment methods  <cit>  .
we then applied the mapgaps program to these alignments to obtain much larger multiple alignments as input to the amcbpps program.
the amcbpps-generated hierarchies were then evaluated by mapping each node’s sequence subgroup onto a phylogenetic tree constructed for all sequences in the hierarchy.
such sequence alignment derived phylogenetic trees are used by cdd curators, both to get started on an initial subfamily classification and to iteratively refine that initial hierarchy .
this reveals that the amcbpps-generated hierarchies agree very well with how the cdd curators would subgroup the sequences based on such a tree.
additional file 1: figure s <dig> shows part of such a sequence tree computed from the input sequences used for rna recognition motif domains  with the amcbpps hierarchy mapped onto the tree using a color coding scheme.
this confirms that the amcbpps program can substantially speed up the curation process when starting from scratch.
‡ for these non-cdd curated alignments were used as input.
§ the time  is for steps  <dig> and  <dig> of the algorithm only.
unaligned sequences were aligned using the multiple alignment procedures cited in methods to generate an input alignment for the amcbpps program.
CONCLUSIONS
currently the construction and annotation of cd hierarchies relies on the labor intensive process of manual curation.
this has created a bottleneck hindering the cdd  <cit>  from achieving the goal of comprehensive coverage of the protein domain universe.
the incorporation of the amcbpps program into the cdd curation pipeline can help automate this process while also providing a statistical measure of the quality of cd hierarchies.
likewise, the delete-half jackknife procedure applied here can provide qualitative estimates of the reliability of various features of a given hierarchy.
and, because the amcbpps program models protein domains based on those residue signatures that most distinguish each functionally divergent subgroup within a protein class from other subgroups, it can also accelerate the annotation of domain profiles.
by linking these profiles to the cn3d viewer  <cit>  structural features associated with likely functionally critical residues can be identified within web-based blast searches.
in previous studies  we have mapped key residues identified using the mcbpps sampler to available crystal structures in this way, thereby obtaining insights into biological functions and mechanisms.
such information also facilitates structural evaluation of sequence alignments.
of course, starting from the procedures described here, the cdd pipeline can be further automated and improved in various ways along similar lines.
for example, we have demonstrated that our bayesian alignment methods can be used to generate, for major protein classes , large multiple alignments in the aligned block-based format required by the cdd.
these, in turn, can serve as input alignments for generating protein domain hierarchies.
moreover, these alignment procedures could be refined to utilize information regarding pattern residue 3d structural interactions to identify and correct misaligned regions automatically .
likewise, protein domain hierarchies generated by the amcbpps program could be further optimized by implementing sampling operations to add or remove leaves and branches.
more sophisticated taxonomic schemes could be devised for distinguishing conserved patterns due to functional constraints rather than to recent common descent.
taken together, these enhancements will accelerate the construction of an optimal, comprehensive set of hierarchically arranged cd profiles.
this will free up curators to focus less on the tedious and labor intensive aspects of database construction and more on biological interpretation, a task that computational and statistical procedures cannot perform.
having such a comprehensive set of well annotated, high quality cd profiles will summarize what is known about each type of domain.
through application of the mapgaps program, these cd hierarchies could be used to obtain up-to-date, very large and highly accurate multiple sequence alignments of an entire protein class for in-depth computational analyses.
and by mapping various categories of pattern residues to corresponding structures, blast searches against these improved cd profiles can reveal those residues most likely responsible for the specific biochemical properties of a query protein.
this can accelerate the pace of biological discovery by enabling researchers to obtain valuable clues regarding as-yet-unidentified protein biochemical and biophysical properties.

 in order to provide rapid and sensitive annotation for protein sequences, including direct links to structural and functional information, the national center for biotechnology information initiated the conserved domain database —a collection of position-specific scoring matrices that are derived from protein multiple sequence alignments. as a result, web-based blast searches now include a search of the cdd, which allows users to visualize multiple sequence alignments and structures of proteins sharing significant homology to the query and, within those alignments, key catalytic and ligand-binding residues. thus blast searches linked to the cdd provide additional clues to the function and underlying mechanism of the query protein and are thereby often more informative, faster and more sensitive than searching against millions of individual protein sequences. the cdd is comprised of domain models either manually curated at the ncbi or imported from other alignment collections such as pfam, smart, and tigrfam. a central and unique feature of the cdd is that related domains are organized into hierarchies when evidence exists to support that tree as a representation of the molecular evolution of the protein class. a significant bottleneck in the cdd pipeline is the curation of these hierarchies and the manual annotation of the corresponding profiles for functionally important residues. in order to begin automating this process, here we describe statistically-rigorous procedures for automated creation of conserved domain hierarchies and for annotation of the corresponding profile alignments. these procedures do automatically, based on objective, empirical criteria, what the cdd resource group and similar groups currently do manually, based on classifications that have been established in the published literature, on phylogenetic and structural analysis and, to some degree, on subjective judgments. our focus here is to obtain heuristically an initial cd hierarchy starting from a typically very large multiple sequence alignment for an entire protein class whose domain boundaries remain fixed. to do this we utilize procedures that obtain both subgroup assignments for aligned sequences and corresponding discriminating patterns associated with protein functional-divergence. the automated annotation of functionally critical residues is an important outcome of these proposed procedures: just as a large enzyme class conserves residues directly involved in catalysis, protein subgroups conserve residues likely involved in subgroup-specific biochemical properties and mechanisms. our procedures use statistical criteria to glean this biochemical information from patterns of divergent residues among related sequences in a manner similar to the use, by classical geneticists, of statistical criteria to glean information from patterns of divergent traits among related individuals.. by mapping various categories of pattern residues to corresponding pssms, blast searches against these improved cd profiles can reveal those residues most likely responsible for the specific biochemical and biophysical properties of a query protein. this can accelerate the pace of biological discovery by enabling researchers to obtain valuable clues regarding as-yet-unidentified protein properties through routine web-based blast searches. other methods that may be similarly described as addressing the protein subfamily classification problem find sequence clusters either based on pairwise similarity or by cutting phylogenetic trees.. here we take a different approach, namely the hierarchical classification of a protein class based on functionally-divergent residue signatures. unlike our approach, many existing methods, though not all, generally focus on the narrower problem of identifying orthologs or on the broader problem of clustering a database into unrelated protein classes rather than on constructing a hierarchy of domain profiles for a specific protein class. an approach, which is, in certain respects, similar to the one described here, is the statistical coupling analysis method of lockless and ranganathan for detecting sets of correlated residues in protein sequences. because our approach identifies residues associated with protein functional divergence, it is also related to "functional subtype" prediction methods, but is distinct inasmuch as these related methods typically predict specific residue functions that are sufficiently well-understood to allow benchmarking. instead, our approach lets the data itself reveal its most statistically striking properties without making assumptions about the types of residues to be identified. it is further distinguished from each of these related methods in at least several of the following respects: it does not require that the input alignment be partitioned into divergent subsets beforehand; this is unlike many, though not all fsp methods. it has a rigorous statistical basis.. it is designed for very large input alignments. for optimization it relies on bayesian sampling, which has a solid scientific basis.. it separates out unrelated and aberrant sequences automatically. it can identify multiple categories of co-conserved residues within a given protein. it addresses concurrently the problems of protein subfamily classification and of identifying residues associated with protein functional divergence. and it can accomplish all of this automatically starting from a single multiple sequence alignment. problem definition and solution strategy here we address the following biological and algorithmic problem: we are given as input a multiple sequence alignment corresponding to a particular protein class. our objective is to partition this alignment into a tree of sub-alignments, termed a cd hierarchy, each subtree of which corresponds to a sub-alignment of sequences sharing a certain pattern that most distinguishes them from those sequences associated with the parent node of the subtree and with any other subtrees attached to that parent node. we interpret these distinguishing residue signatures as associated with functional divergence of the protein class. as mentioned above, our focus is on obtaining an initial, suboptimal hierarchy that is comparable to current cdd curated hierarchies and that can serve as a starting point for further optimization using either manual or automated methods. here we describe statistically-based heuristic procedures that, in conjunction with bayesian sampling, can obtain such an initial hierarchy from a multiple sequence alignment. bayesian sampling over contrast alignment models our approach relies on bayesian markov chain monte carlo sampling, which starts with an arbitrary model having a certain probability that it could have generated the input sequence alignment data. then, in successive iterative steps, a number of alternative values for a given parameter of the current model are evaluated. a value for this parameter is then sampled proportional to its probability. this iterative process continues until convergence on the most likely models. here we optimize in this way contrast alignment models, each of which consists of a pattern and a set of labels assigning each sequence to either a foreground partition or a background partition corresponding to sequences that either generally conserve or fail to conserve the pattern, respectively. a schematic representation of a contrast alignment and the corresponding probability distribution are described in figure to sample alternative contrast alignment models we use a mcmc sampling strategy, termed bayesian partitioning with pattern selection. the sampler converges on a model where the pattern best distinguishes the foreground from the background sequences. multiple category functional divergence models more recently, a multiple category bpps sampler was developed with a view to optimally assigning aligned sequences to various nodes within a predefined protein domain hierarchy based on functionally-divergent residue signatures. thus the mcbpps sampler aims to precisely define both the sequences belonging to each subgroup and the patterns most distinctive of each subgroup within a specific protein class. however, because the mcbpps sampler does not define the hierarchy of contrast alignments, it requires that the user provide both a functional divergence -table and seed sequences for each divergent subgroup. sequences are assigned to subgroups through bayesian inference). each row of a fd-table corresponds to a distinct functionally divergent subgroup of the input sequences and each column corresponds to a distinct contrast alignment whose foreground and background partitions are specified by the symbols in the table. such a table is shown in figure which also illustrates the correspondence between a tree representing the hierarchical relationships between functionally divergent subgroups and the fd-table and between a column in the table and the contrast alignment; these are shown above and below the table, respectively. given the relationships specified by the fd-table, the mcbpps sampler stochastically reassigns aligned sequences to alternative subgroups and alternative patterns to each foreground partition until convergence on an optimal set of contrast alignments. modeling the functional divergence of an entire protein class in this way is substantially more powerful than using a single contrast alignment because: in principle, it can optimally model every divergent subgroup within an entire protein class concurrently. it sets up a stringent competition between functionally-divergent categories for pattern residues, thereby defining each pattern and partition much more precisely. it eliminates problematic sequences, which would otherwise tend to obscure analyses, by modeling them explicitly. problematic sequences include, for example, related proteins that have further functionally diverged to become outliers, pseudogene products and other non-functional proteins, and unrelated or erroneous sequences. and, by defining multiple categories of pattern residues within individual proteins it can reveal, in the light of available structural information, functionally important residue interactions.. here we describe and apply an automated multiple category bpps program that generates its own fd-table and seed sequences automatically and therefore merely requires a multiple sequence alignment as input. the number and nature of the partitions and the patterns is completely determined by the program. when used in conjunction with procedures for viewing structural interactions involving pattern residues, the amcbpps sampler automates and enhances the creation and annotation of cdd hierarchical alignments. and, when linked into web-based blast searches, this can make previously inaccessible molecular information widely available.  in this section, we lay out the basic amcbpps algorithm, illustrate an implementation of the algorithm as applied to p-loop gtpases, compare its performance against various manually-curated cd hierarchies, further evaluate its performance using both delete-half jackknifing and simulations, and apply it to several large protein classes for which existing hierarchies or alignments are currently unavailable. algorithm the amcbpps algorithm aims to identify the hierarchical relationships between functionally-divergent subgroups within an entire protein domain class based on the differentiating patterns present in that class. it does this by defining: the number of sequence sets, the members of each set, the hierarchical relationships between sets and the corresponding functionally divergent patterns. this is accomplished in three steps. steps and constitute the novel aspect of the program by providing input to the mcbpps sampler in step 3; these first two steps also speed up convergence in step by providing a better starting point for the mcbpps sampler. conceptual aspects of the algorithm corresponding to steps and are illustrated in figure. these first two steps are performed heuristically and thus constitute an informed guess on how to best model the protein class. the sampler then improves upon this model by optimizing sequence and pattern assignments within this hierarchy. the hierarchy returned by the amcbpps program may be edited and then optimized again by the mcbpps program. such editing may, for example, further subcategorize previously identified or miscellaneous subgroups. identifying simple subgroups step of the amcbpps algorithm first generates a forest of simple trees, the leaves of which correspond both to functionally-divergent subgroups within the protein domain class and to sub-alignments of the input alignment. this is accomplished by obtaining seed sequences and using them to create these trees, as follows: all closely-related sequences from the same phylum are clustered into a common set. all pairs of moderately related sequences from distinct phyla are stored on a heap using their pairwise scores as the key. iteratively remove the top-scoring cross-phylum pair from the heap and merge their two disjoint sets into one set.. once a disjoint set contains sequences from a pre-defined minimum number of distinct phyla, the sequences of highest rank from each phylum are used to seed a new subgroup; the disjoint set is then labeled to avoid picking this subgroup repeatedly. keep generating subgroups in this way until a pre-defined number of seed sequence sets are obtained, at which point a simple fd-table is constructed where each subgroup node is a direct descendent of the root node.. and repeat substeps iii-v until all sequence pairs have been removed from the heap. to ensure that the subgroups are sufficiently diverse, we require that each seed set consensus sequence share less than a specified level of sequence identity with other seed set consensus sequences. taken together, these sub-steps favor the identification of the most conserved and phylogenetically diverse subgroups. for each of these fd-tables the mcbpps sampler assigns each of the multiply aligned input sequences to a subgroup and determines the differentiating conserved pattern for each contrast alignment. to ensure that subgroups at different levels of the hierarchy are identified, the algorithm performs multiple runs using various numbers of leaf nodes and various prior probability settings for p, p and p. convergence on protein subfamilies is favored by specifying a high number of leaf nodes, by lowering the prior probability for assigning a sequence to a leaf node, p,, by setting the prior probability, p, to favor a lower degree of background contamination by assigning more pseudo-observations to pattern matching residues and fewer pseudo-observations to contaminating residues and by raising the prior probability that a column corresponds to a pattern position, p. the rationale for choosing these settings is that, for subfamilies, membership is more exclusive, sequences are more highly conserved and, consequently, conserved patterns more extensive. in order to maximize the influence of the data on model optimization.) convergence on a super-family is favored by specifying a single subgroup and by altering these prior parameter settings accordingly. default settings are based on applications to actual protein sequences, though it should be noted that the influence of these prior settings is minor. hence these priors primarily function as tuning parameters to help gently guide the sampler into finding a variety of functionally divergent subgroups. to avoid finding the same subgroup repeatedly, sequences assigned to a subgroup in a previous run are prohibited from being used as seeds in subsequent runs. subfamilies can also be identified recursively; that is, by rerunning the program on a single subgroup in order to find subgroups within subgroups. the pseudocode for this step of the amcbpps algorithm is given in methods. defining a hierarchy for the protein class once individual subgroup sets are identified in step, the program hierarchically arranges these into a more complex tree, from which a fd-table is obtained. it does this using efficient bitwise set operations, standard network algorithms and pattern-based statistical criteria, namely the contrast alignment log-likelihood ratio used by the mcbpps sampler. these hierarchically-arranged subgroup sets are ‘fuzzy’ due to the uncertainty associated with set membership. thus step of the algorithm determines which of the input sets: are the same set; are distinct sets; or are supersets of another set or sets. step is subdivided into three sub-steps: merge each collection of subgroup sets deemed to be identical into a single set; cluster related subgroup sets into common supersets; and create a tree representation of the subgroup hierarchy, which is done using a breadth-first scanning algorithm to find a shortest path tree. step 2c also refines the tree to eliminate inappropriate overlap between sets while also eliminating nodes from the tree that, as a result of this refinement process, are no longer statistically significant. the pseudocode for step is given in methods. from this tree a fd-table is then generated as input to step the mcbpps sampler and further refinements the output from step provides a starting point for mcbpps sampling, which optimizes the patterns and partitions corresponding to the fd-table. the basic statistical and algorithmic aspects of the mcbpps sampler were previously described. to further expand a cd hierarchy the output files obtained from an initial amcbpps analysis can also be used to recursively analyze, in the same way, several of the larger subgroups. to do this, the output alignment file for a major subgroup is used as an input file for the amcbpps program. likewise, a cd hierarchy can also be refined by editing the fd-table manually and then applying the mcbpps sampler, as was previously described. to speed up analysis of a given subtree, such manually edited fd-tables may be designed to expand subgroups within that subtree, while modeling the other branches of the hierarchy only at the highest levels. in keeping with the automated theme of this article, however, we will not describe in detail nor apply these approaches here. implementation and testing the amcbpps algorithm was implemented in c++, applied to various protein classes and the output compared to manually-curated cdd alignment hierarchies. a wide range of cdd hierarchies—from preliminary to well-developed releases —were examined in this way. input multiple alignments were obtained by using the ncbi hierarchy of cd alignments as input to the mapgaps program, which detected and aligned related protein sequences within the ncbi nr, env_nr and translated est protein databases. to obtain input alignments corresponding to large protein classes for which cdd hierarchies are not yet available we used alternative procedures, as described in methods. illustrative example: p-loop gtpases. to familiarize the reader we begin by illustrating our approach with an analysis of p-loop gtpases. using an input alignment of p-loop gtpases, the amcbpps program returned the fd-table shown in figure. it also returns a corresponding set of contrast alignments, which highlight the pattern residues identified by the sampler; one such alignment is shown in. note that the sampler will reject heuristically proposed subgroups whose existence is not supported by the data. further subdivision of the hierarchy in figure may be accomplished by recursively applying the amcbpps sampler to a previously-identified subgroup. illustrates this for the ras-like gtpases by showing an expanded subtree corresponding to the column foreground partition in figure by applying the amcbpps sampler recursively in this way, a very extensive hierarchy may be obtained. ‡ after removing identical sequences and sequences that fail to align with at least 75% of the domain. * numbers in parentheses indicate the nodes retained after insignificant nodes were removed by the mcbpps program. † the log-likelihood ratio in nats. § the time is for steps and of the algorithm only; step can be parallelized to run in less than 10% of the time shown. criteria for comparing hierarchies to assess how well the amcbpps program performs relative to curated cd hierarchies, we compared its output against manually curated cdd hierarchies. before considering this analysis, however, we first need to discuss the criteria used to evaluate and compare hierarchies. lack of gold standards cdd hierarchies have been carefully constructed by expert curators and therefore come the closest to a benchmark set for evaluating the amcbpps sampler. however, as this study reveals, certain aspects of cdd hierarchies lack statistical support or are incomplete or incorrect for various reasons: for example, cdd hierarchies are typically at different stages in an ongoing refinement process, and, for protein domain classes consisting of tens or hundreds of thousands of sequences, the number of possible hierarchies to consider is astronomical, which makes optimization through manual curation extremely difficult. furthermore, due to the stochastic nature of and the inability to directly observe evolutionary divergence, it is impossible to eliminate the inherent uncertainties associated with protein classification. hence, for the present study our aim is merely to replicate the current manual curation process by generating hierarchies of comparable quality automatically, thereby dramatically speeding up the current labor-intensive curation process. comparison criteria for this analysis despite the absence of a gold standard, the statistical criteria used by the amcbpps program provide a way to compare two hierarchies for the same conserved domain. it does this by determining objectively whether or not the sequences in each protein subgroup have diverged from the evolutionarily related subgroups indicated by a specific hierarchy. this measure is expressed as a log-likelihood ratio, where non-positive values indicate a lack of statistical support for a functionally divergent event within the hierarchy. such a comparison is performed as follows: we are given two heuristic methods for obtaining a hierarchy: one manual and one automatic. to compare the two methods, we first use each hierarchy as input to the mcbpps sampler, which then optimizes the patterns and sequence partitions associated with that hierarchy and returns an optimized log-likelihood ratio. because this optimizes the automatically-generated and manually-curated hierarchies in the same way based on the same statistical criteria, the only difference is that the hierarchies and seed alignments were obtained either automatically or through manual curation. thus, by comparing their optimized llr scores, we can obtain a measure of the relative performance of the two methods. in addition, we also determine the degree of overlap between the two hierarchies as a qualitative indication of the similarity of the two hierarchies. the results of such comparisons are summarized in table and figure evaluation of the amcbpps program to evaluate the amcbpps program over a wide variety of input, we chose the conserved domains given in table these domains vary in the numbers of members detected in the protein databases, in the lengths of their conserved core, and in the size and complexity of their curated hierarchies. comparisons with manually curated cdd hierarchies based on the llr statistic the automatically-generated hierarchies are comparable to the corresponding manually-curated hierarchies and are, in fact, slightly better on average. manual and amcbpps hierarchies were also compared qualitatively by determining the degree to which the sequence sets corresponding to the nodes in one hierarchy overlap with those of the other hierarchy; in additional file 1: figure s illustrates such a comparison for the pi3kc_like domain hierarchy. this provides a detailed comparison of two hierarchies by computing how the sequences at each level of one hierarchy are assigned to each level of the other hierarchy and vice versa. such comparisons indicate that the differences between the curated and amcbpps hierarchies are mainly due to the following: a node in one hierarchy being modeled as a subtree in the other. additional child nodes being added to parent nodes in one hierarchy but not the other. a subtree in one hierarchy being split into unrelated subtrees in the other due to a failure to join these to a common internal node. and inherently ambiguous sequences that can’t be clearly assigned to a specific node in the hierarchy; such sequences may correspond to pseudogene products or functionally defective members of a protein class that are difficult to categorize because they harbor degenerate subgroup patterns. of course, the larger and more functionally divergent hierarchies are more challenging. most of the differences between the manual and automated hierarchies are due to fundamental differences between the two approaches. curated hierarchies may rely, in part, on information that is ignored by the amcbpps program, such as subfamily-conserved inserts and 3d structures—though, on the other hand, the amcbpps program utilizes a far greater amount of sequence data that is also up-to-date. in contrast, some of the cdd heirarchies may be out-of-date or still incomplete. the amcbpps program also requires each node to correspond to at least sequences in order to avoid statistical biases due to small sample size. cdd curators, however, may construct subgroups containing fewer sequences. likewise, the amcbpps program selects seed sequences from at least three or four distinct phyla in order to avoid sampling biases introduced by orthologous sequences from closely related organisms. hence, it will fail to identify subgroups that only occur in vertebrates, for instance. in contrast, cdd curators may choose representative sequences from more closely related taxa. due to such restrictions, an amcbpps-generated hierarchy tends to have fewer nodes because it is prohibited from identifying certain cdd-defined subgroups. this also tends to lower the llr, which increases with the number of nodes in the hierarchy. despite these differences, after taking these considerations into account, we found the amcbpps hierarchies to be persuasively consistent with the corresponding cdd hierarchies. a perceived “unfair” advantage of the amcbpps algorithm might be that it utilizes the same statistical model to construct a hierarchy that is used to score that hierarchy, whereas curators do not. however, subsequent optimization of both hierarchies using the mcbpps sampler should counteract this advantage. that is, assuming that the curated hierarchy is in fact superior and that our statistical model is biologically meaningful, then optimal partitioning of the sequences and optimal pattern assignment by the sampler for both types of hierarchies should result in a superior llr score for the curated hierarchy. unsurprisingly, our analysis also indicates that the hierarchies obtained both manually and automatically are typically suboptimal. for example, manual and amcbpps hierarchies for the s8/s peptidase domain had llrs of and nats, respectively, whereas a hybrid hierarchy containing features of both of these has a llr of nats. figure likewise illustrates the construction of a hybrid hierarchy for the class i aars-like core domain, which improves the llr from to nats. examining llr statistics can suggest other ways in which to improve cdd hierarchies. for example, within the glo_edi_brp_like hierarchy, the mcbpps sampler rejected an intermediate node to which it had assigned a llr of − nats. an investigation to determine why this occurred revealed that, based as well on the criteria used by the cdd curators, the cd intermediate node is incompatible with several of its leaf nodes, which are therefore better modeled as direct descendents of the root node. for these and other domains in table suggested improvements in cdd hierarchies based on llrs were corroborated through manual inspection by the cdd resource group. delete-half jackknife analyses a bootstrap or jackknife procedure can be used to estimate confidence levels for evolutionary trees. however, applying this approach to a cd hierarchy is complicated by the potential run-to-run variability both in the number of the leaf nodes and in the associated sequence sets. thus existing evolutionary tree bootstrap and jackknife procedures, which require that each of the sampled trees use the same set of leaf nodes, cannot be used. instead, we implemented a delete-half jackknife procedure that—though unable to provide quantitative confidence levels for specific features of a hierarchy—can nevertheless provide a qualitative assessment of the run-to-run variability of amcbpps-generated hierarchies. this involved running the amcbpps program on each of different input alignments after randomly removing half of the sequences. given ten such runs for each domain, we compared the consistency of the resultant hierarchies from run to run as follows: for each pair of runs we determined how those input sequences shared by both hierarchies were partitioned among the nodes of one hierarchy relative to the other hierarchy. an example output file is shown in. for these analyses we found that, among the leaf node sets in one tree that share at least one sequence in common with a leaf node set in the other tree, on average 47% share precisely the same set of sequences and 74% share more than 90% of their sequences in common. moreover, in most cases where an identical sequence set is not found, the missing sequences were typically assigned, not to unrelated leaf nodes, but either to a parent node further up the tree or to the rejected sequence set. among the remaining cases, a node in one hierarchy is either split into multiple nodes or split between nodes in the other hierarchy. at times a hierarchy could end up omitting certain nodes due to the delete-half jackknife procedure removing sequences belonging to certain phyla resulting in insufficient phylogenetic diversity to seed the formation of a subgroup. of course the topologies of the jackknife trees found by the sampler also differ, which is a common problem associated with evolutionary trees consisting of large numbers of distantly related sequences. this is presumably due in large part to the amcbpps algorithm failing to find the optimal topology—an issue that, in the future, we will address by sampling over alternative topologies. of course, both this future sampler and the jackknife procedure applied here will be useful for identifying the most reliable features of a hierarchy. taken together, these results confirm the observation we made in the previous section, namely that the amcbpps program generally finds a suboptimal hierarchy that, nevertheless, provides a good starting point both for curation and further automation. output from these jackknife analyses are available at http://www.chain.umaryland.edu/amcbpps/jackknife.txt. simulations as an additional check, we implemented a procedure to generate simulated sequences from profile hmms where each such profile corresponds to a node from one of the domain hierarchies used in the jackknife analysis. the rationale for doing this was to determine how well the amcbpps program identifies sequences corresponding to predefined subgroups. note that this procedure captures sequence features of each subgroup, but not how those subgroups are hierarchically arranged. for each node of each hierarchy we generated the same number of aligned sequences as were assigned to that node in the original hierarchy. after running the amcbpps program on each of these simulated alignments, we determined the degree to which each set of related simulated sequences were correctly modeled as belonging to a single subgroup. an example output file in illustrates how the amcbpps program correctly categorized nearly all of these sequences given the structure of the inferred hierarchy. output from simulations for the domains is available at http://www.chain.umaryland.edu/amcbpps/simulate.txt. time complexity the computationally most intensive routine in step of the amcbpps program is an all-versus-all pairwise comparison of pre-aligned sequences. this has a time complexity of o = o where k is the number of aligned columns, m is the number of sequences and n = k·m is the effective size of the input alignment. in addition, step involves a simpler version of the step algorithm that, of course, exhibits the same time complexity as step as well as other operations that perform better than o. hence the time complexity for step is o. the time complexity of steps 2– is unclear based on the underlying algorithm. therefore, using a plot of the run times for the amcbpps analyses in table versus the size of the input alignments, we estimate that the time required for steps 2– scales as o. assuming that the asymptotic time complexity for steps 2– is indeed o, which admittedly may not be the case given our empirically-based approach, then whether or not o is better than o depends on the ratio of k to m step which is o, performs asymptotically worse when k < m and steps 2– which is o, is worse when k > m since for essentially all protein domains k < m the time complexity of the amcbpps program appears to be o. it is important to note, however, that step required less time than steps 2– in our analyses—even on the largest input alignments—suggesting that constant factors rather than asymptotics are influencing program performance. because the run times for steps 2– are also likely to depend on the size of the hierarchy generated by the program in step figure 7b plots the run times versus the number of aligned residues times the number of nodes in the hierarchy. this yields a slightly improved, essentially linear dependency. this observed time complexity is largely due to step being more or less independent of the number of nodes, which is achieved by computing conditional posterior probabilities for each column of the fd-table only when considering the assignment of a sequence to one of two possible new partitions rather than to one of a typically much larger number of rows. thus the amcbpps program can be applied to very large multiple sequence alignments, which is important given the current rapid increase in sequence data. analysis of protein domains lacking cd hierarchies there are a significant number of protein domains for which a cdd hierarchy has not yet been constructed. in some cases a single curated alignment is available as a starting point. to test the performance of the amcbpps program in such cases, we chose domains, for which curated alignments were available, and two domains, for which we first constructed an alignment using bayesian multiple alignment methods. we then applied the mapgaps program to these alignments to obtain much larger multiple alignments as input to the amcbpps program. the amcbpps-generated hierarchies were then evaluated by mapping each node’s sequence subgroup onto a phylogenetic tree constructed for all sequences in the hierarchy. such sequence alignment derived phylogenetic trees are used by cdd curators, both to get started on an initial subfamily classification and to iteratively refine that initial hierarchy. this reveals that the amcbpps-generated hierarchies agree very well with how the cdd curators would subgroup the sequences based on such a tree. additional file 1: figure s shows part of such a sequence tree computed from the input sequences used for rna recognition motif domains with the amcbpps hierarchy mapped onto the tree using a color coding scheme. this confirms that the amcbpps program can substantially speed up the curation process when starting from scratch. ‡ for these non-cdd curated alignments were used as input. § the time is for steps and of the algorithm only. unaligned sequences were aligned using the multiple alignment procedures cited in methods to generate an input alignment for the amcbpps program.  currently the construction and annotation of cd hierarchies relies on the labor intensive process of manual curation. this has created a bottleneck hindering the cdd from achieving the goal of comprehensive coverage of the protein domain universe. the incorporation of the amcbpps program into the cdd curation pipeline can help automate this process while also providing a statistical measure of the quality of cd hierarchies. likewise, the delete-half jackknife procedure applied here can provide qualitative estimates of the reliability of various features of a given hierarchy. and, because the amcbpps program models protein domains based on those residue signatures that most distinguish each functionally divergent subgroup within a protein class from other subgroups, it can also accelerate the annotation of domain profiles. by linking these profiles to the cn3d viewer structural features associated with likely functionally critical residues can be identified within web-based blast searches. in previous studies we have mapped key residues identified using the mcbpps sampler to available crystal structures in this way, thereby obtaining insights into biological functions and mechanisms. such information also facilitates structural evaluation of sequence alignments. of course, starting from the procedures described here, the cdd pipeline can be further automated and improved in various ways along similar lines. for example, we have demonstrated that our bayesian alignment methods can be used to generate, for major protein classes, large multiple alignments in the aligned block-based format required by the cdd. these, in turn, can serve as input alignments for generating protein domain hierarchies. moreover, these alignment procedures could be refined to utilize information regarding pattern residue 3d structural interactions to identify and correct misaligned regions automatically. likewise, protein domain hierarchies generated by the amcbpps program could be further optimized by implementing sampling operations to add or remove leaves and branches. more sophisticated taxonomic schemes could be devised for distinguishing conserved patterns due to functional constraints rather than to recent common descent. taken together, these enhancements will accelerate the construction of an optimal, comprehensive set of hierarchically arranged cd profiles. this will free up curators to focus less on the tedious and labor intensive aspects of database construction and more on biological interpretation, a task that computational and statistical procedures cannot perform. having such a comprehensive set of well annotated, high quality cd profiles will summarize what is known about each type of domain. through application of the mapgaps program, these cd hierarchies could be used to obtain up-to-date, very large and highly accurate multiple sequence alignments of an entire protein class for in-depth computational analyses. and by mapping various categories of pattern residues to corresponding structures, blast searches against these improved cd profiles can reveal those residues most likely responsible for the specific biochemical properties of a query protein. this can accelerate the pace of biological discovery by enabling researchers to obtain valuable clues regarding as-yet-unidentified protein biochemical and biophysical properties. 
this approach automates the rapid creation of protein domain hierarchies and thus will eliminate one of the most time consuming aspects of conserved domain database curation. the ncbi conserved domain database consists of a collection of multiple sequence alignments of protein domains that are at various stages of being manually curated into evolutionary hierarchies based on conserved and divergent sequence and structural features. at the same time, it also facilitates protein domain annotation by identifying those pattern residues that most distinguish each protein domain subgroup from other related subgroups. these statistical criteria also provide a means to objectively assess the quality of cd hierarchies, a non-trivial task considering that the protein subgroups are often very distantly related—a situation in which standard phylogenetic methods can be unreliable. a plot of runtimes for the most time-intensive part of the algorithm indicates a nearly linear time complexity so that, even for the extremely large rossmann fold protein class, results were obtained in about a day. here we automate the generation of conserved domain hierarchies using a combination of heuristic and markov chain monte carlo sampling procedures and starting from a multiple sequence alignment. this procedure relies on statistical criteria to define each hierarchy based on the conserved and divergent sequence patterns associated with protein functional-specialization. our aim here is to automatically generate hierarchies that, based on statistical criteria and visual comparisons, are comparable to manually curated hierarchies; this serves as the first step toward the ultimate goal of obtaining optimal hierarchical classifications. these domain models are annotated to provide insights into the relationships between sequence, structure and function via web-based blast searches.
 the ncbi conserved domain database consists of a collection of multiple sequence alignments of protein domains that are at various stages of being manually curated into evolutionary hierarchies based on conserved and divergent sequence and structural features. these domain models are annotated to provide insights into the relationships between sequence, structure and function via web-based blast searches.  here we automate the generation of conserved domain hierarchies using a combination of heuristic and markov chain monte carlo sampling procedures and starting from a multiple sequence alignment. this procedure relies on statistical criteria to define each hierarchy based on the conserved and divergent sequence patterns associated with protein functional-specialization. at the same time this facilitates the sequence and structural annotation of residues that are functionally important. these statistical criteria also provide a means to objectively assess the quality of cd hierarchies, a non-trivial task considering that the protein subgroups are often very distantly related—a situation in which standard phylogenetic methods can be unreliable. our aim here is to automatically generate hierarchies that, based on statistical criteria and visual comparisons, are comparable to manually curated hierarchies; this serves as the first step toward the ultimate goal of obtaining optimal hierarchical classifications. a plot of runtimes for the most time-intensive part of the algorithm indicates a nearly linear time complexity so that, even for the extremely large rossmann fold protein class, results were obtained in about a day.  this approach automates the rapid creation of protein domain hierarchies and thus will eliminate one of the most time consuming aspects of conserved domain database curation. at the same time, it also facilitates protein domain annotation by identifying those pattern residues that most distinguish each protein domain subgroup from other related subgroups. 