Protein domain architectures provide a fast, efficient and scalable alternative to sequence-based methods for comparative functional genomics

A functional comparative genome analysis is essential to understand the mechanisms underlying bacterial evolution and adaptation. Detection of functional orthologs using standard global sequence similarity methods faces several problems; the need for defining arbitrary acceptance thresholds for similarity and alignment length, lateral gene acquisition and the high computational cost for finding bi-directional best matches at a large scale. We investigated the use of protein domain architectures for large scale functional comparative analysis as an alternative method. The performance of both approaches was assessed through functional comparison of 446 bacterial genomes sampled at different taxonomic levels. We show that protein domain architectures provide a fast and efficient alternative to methods based on sequence similarity to identify groups of functionally equivalent proteins within and across taxonomic boundaries, and it is suitable for large scale comparative analysis. Running both methods in parallel pinpoints potential functional adaptations that may add to bacterial fitness.


Introduction
Comparative analysis of genome sequences has been pivotal to unravel mechanisms shaping bacterial evolution like gene duplication, loss and acquisition 1,2 , and helped in shedding light on pathogenesis and genotype-phenotype associations 3,4 .
Comparative analysis relies on the identification of sets of orthologous and paralogous genes and subsequent transfer of function to the encoding proteins. Technically orthologs are defined as best bi-directional hits (BBH) obtained via pairwise sequence comparison among multiple species and thus exploits sequence similarity for functional grouping. Sequence similarity-based (SB) methods present a number of shortcomings. First, a generalized minimal alignment length and similarity cut-off need to be arbitrarily selected for all, which may hamper proper functional grouping. Second, sequence and function might differ across evolutionary scales. Protein sequences change faster than protein structure and proteins with same function but with low sequence similarity have been identified 5,6 . SB methods may fail to group them hampering a functional comparison. This limitation becomes even more critical when comparing either phylogenetically distant genomes or gene sequences that were acquired with horizontal gene transfer events. Recent technological advancements are resulting in thousands of organisms and billions of proteins being sequenced 7 which increases the need of methods able to perform comparisons at the larger scales.
To overcome these bottlenecks, protein domains have been suggested as an alternative for defining groups of functionally equivalent proteins 8-10 and have been used to perform comparative analyses of Escherichia coli 9 , Pseudomonas 10 , Streptococcus 11 and for protein functional annotation 12,13 . A protein domain architecture describes the arrangement of domains contained in a protein and is exemplified in Figure 1. As protein domains capture key structural and functional features, protein domain architectures may be considered to be better proxies to describe functional equivalence than a global sequence similarity 14 . The concept of using the domain architecture to precisely describe the extent of functional

Amendments from Version 2
We have amended the manuscript as suggested by the reviewer. Specifically: • The Abstract and Introduction no longer state that the requirements of the SB approach, time and memory, need to scale quadratically with the number of genomes.
• We have modified the Discussion to further emphasize that DAB is similar to SB methods which extend existing groups into new genomes.
• We have also rephrased the reviewers' comment regarding the extensive use of DAB to define domain families, as we think it might further clarify the text.
• The sentence "Our aim was to investigate whether using HMMs instead of sequence similarity would yield similar results" has been modified as suggested, to: "Our aim was to investigate whether using domain architectures instead of sequence similarity alone would yield similar results."  equivalence is exemplified in Figure 2. Moreover, once the probabilistic domain models have been defined, mining large sets of individual genome sequences for their occurrences is a considerably less demanding computational task than an exploration of all possible bi-directional hits between them 15,16 .
Domain architectures have been shown to be preserved at large phylogenetic distances both in prokaryotes and eukaryotes 17,18 . This lead to the use of protein domain architectures to classify and identify evolutionarily related proteins and to detect homologs even across evolutionarily distant species [19][20][21][22] . Structural information encoded in domain architectures has also been deployed to accelerate sequence search methods and to provide better homology detection. Examples are CDART 23 which finds homologous proteins across significant evolutionary distances using domain profiles rather than direct sequence similarity, or DeltaBlast 24 where a database of pre-constructed position-specific score matrix is queried before searching a protein-sequence database. Considering protein domain

Retrieval of domain architecture
The positions (start and end on the protein sequence) of domains having Pfam 28 , TIGRFAMs 29 and InterPro 30 identifiers were extracted through SPARQL querying of the graph database and domain architectures were retrieved for each protein individually. InterPro aggregates protein domain signatures from different databases. Here no pruning for redundancies has been done. Identification of domains was done using the intrinsic InterPro cut-off that represents in each case the e-values and the scoring systems of the member databases 30 . The domain starting position was used to assess relative position in the case of overlapping domains; alphabetic ordering was used to order domains with the same starting position or when the distance between the starting position of overlapping domains was < 3 amino acids.
Labels indicating N-C terminal order of identified domains were assigned to each protein using the starting position of the domains: the same labels were assigned to proteins sharing the same domain architecture.
Sequence similarity based clustering To make a direct comparison possible only protein sequences containing at least one protein domain signature were considered for analysis. BBH were obtained using Blastp (2.2.28+) with an E-value cutoff of 10 −5 and -max_target_seqs of 10 5 . OrthaGogue (1.0.3) 35 combined with MCL (14-137) 36 was used to identify protein clusters on the base of sequence similarity.
Domain architecture based clustering Domain architecture based clusters were built by clustering proteins with the same labels using bash terminal commands (sort, awk). The number of proteins sharing a given domain architecture in each genome was stored in a 446 × 21054 (genomes × domain architectures) matrix; from this a binarized presence-absence matrix was obtained and used solely for principal component analysis.
Heaps' law fitting and pan-genome openness assessment A Heaps' law model was fit to the abundance matrices using 5 × 10 3 random genome ordering permutations and the micropan R package 37 .
Software SAPP, a Semantic Annotation Pipeline with Provenance which stores results in a graph database 10 , used for genome handling and annotation, is available at http://semantics.systemsbiology.nl. Matrix manipulations and multivariate analysis were performed using the R software (3.2.2).

Results
SB and DAB approaches were compared by considering eight sets of genome sequences sampled at different taxonomic levels, from species to order, preserving phylogenetic diversity (see Table 1). Each set contained 60 genome sequences, except for Listeria monocytogenes for which only 26 complete genomes were publicly available. To facilitate the comparison between DAB and SB clusters only protein sequences that contained at least one domain were considered. On average, 85% of the protein sequences contain at least one domain from the InterPro database (see Table 1). Values content, order, recurrence and position has been shown to increase the accuracy of protein function prediction 25 and has led to the development of tools for protein functional annotation, such as UniProt-DAAC 26 which uses domain architecture comparison and classification for the automatic functional annotation of large protein sets. The systematic assessment and use of domain architectures is enabled by databases containing protein domain information such as UniProt 27 , Pfam 28 , TIGRFAMs 29 , InterPro 30 , SMART 31 and PROSITE 32 , that also provide graphical view of domain architectures.
Building on these observations we aim at exploring the potential of domain architecture-based (DAB) methods for large scale functional comparative analysis by comparing functionally equivalent sets of proteins, defined using domain architectures, with standard clusters of orthogonal proteins obtained with SB methods. We compared the SB and DAB approach by analysing i) the retrieved number of singletons (i.e. clusters containing only one protein) and ii) the characteristics of the inferred pan-and core-genome size considering a selection of bacterial genomes (both gram positive and negative) sampled at different taxonomic levels (species, genus, family, order and phylum). We show that the DAB approach provides a fast and efficient alternative to SB methods to identify groups of functionally equivalent/related proteins for comparative genome analysis and that the functional pan-genome is more closed in comparison to the sequence based pan-genome. DAB approaches can complement standardly applied sequence similarity methods and can pinpoint potential functional adaptations.

Genome sequence retrieval
Bacterial species were chosen on the basis of the availability of fully sequenced genomes in the public domain: two species (Listeria monocytogenes and Helicobacter pylori), three genera (Streptococcus, Pseudomonas, Bacillus), one family (Enterobacteriaceae), one order (Corynebacteriales), and one phylum (Cyanobacteria) were selected. For each, 60 genome sequences were considered, except for L. monocytogenes for which only 26 complete genome sequences were available. Maximal diversity among genome sequences was ensured by sampling divergent species (when possible) at each taxonomic level. Genome sequences were retrieved from the European Nucleotide Archive database (www.ebi.ac.uk/ena). A full list of genomes analyzed is available in the Data availability section.
De novo genome annotation To avoid bias due to different algorithms used for the annotation of the original deposited genome sequences, all genomes were de novo re-annotated using the SAPP framework (1.0.0) 10 . In particular, the FASTA2RDF, GeneCaller (implementing Prodigal (2.6.2) 33 ) and InterPro (implementing interproscan-5.17-56.0) 34 ) modules were used to handle, re-annotate the genome sequences and store the results in the RDF data model. This resulted in 446 annotated genomes (7 × 60 genomes + 1 × 26 genomes). For each annotation step the provenance information (E-value cut off, score, originating tool or database) was stored together with annotation information in a graph database (RDF-model) and can be reproduced through the SAPP framework (http://semantics.systemsbiology.nl). Cluster formation based on sequence similarity A standard BBH workflow was used to obtain SB protein clusters for the eight sets. We started by calculating the total number of clusters, corresponding to the pan-genome size, as shown in Table 1. Then we considered protein cluster persistence, that is the number of genomes where at least one member of the cluster is present, divided by the total number of genomes considered. Results are shown in Figure 3.
The ratio between the size of the core-genome (clusters with persistence of 1, i.e. present in all genomes) and the number of singletons decreased with evolutionary distance (see Table 1). It ranged from 3.51 and 3.07 at species level (H. pylori and L. monocytogenes respectively) to 0.05 and 0.06 when considering members of the same order (Corynebacteriales) and phylum (Cyanobacteria) respectively. A similar pattern is observed when directly comparing the sizes of the pan-and core-genomes of the sampled genomes. Within the gram negative bacteria this ratio ranges from 0.69 for members of the same species (H. pylori) to 0.05 for members of the same phylum (Cyanobacteria) with intermediate values (0.12) for sequences from the same genus (Pseudomonas).

Cluster formation based on domain architectures
Domain architectures directly rely on the definition of protein domain models: those were retrieved from Pfam, InterPro and TIGRFAMs databases. However, TIGRFAMs results were not further considered because of a lower coverage. As shown in Table 1, as expected partly overlapping results were obtained when different domain databases were used. The number of singletons was larger when using InterPro rather than Pfam and for the latter we also observed larger core-genome size. These discrepancies can be due to the fact InterPro aggregates different resources (including Pfam and TIGRFAMs) and domain signatures arising from different databases are integrated with different identifiers in InterPro. In light of this we focused on results obtained using Pfam whose current release (30.0) contains hidden Markov models for over 16300 domain families. Size and persistence of groups of functionally equivalent proteins obtained using Pfam domains are presented in Figure 4.
Similar to what has been observed in the SB case we observed a decrease of the ratio between the size of the core genome and the number of singletons when higher taxonomic levels are considered. For organisms of the same species (H. pylori and L. monocytogenes) the ratio was 5.09 and 4.30, respectively, while for member of the same order (Corynebacteriales) and phylum (Cyanobacteria) it was 0.55 and 0.009 respectively. Similarly, also the ratio between the size of the core-and pan-genome decreases as higher taxonomic levels are considered, ranging from 0.54 for H. pylori to 0.04 for Cyanobacteria.

Comparison of DAB and SB clusters
We compared the clusters obtained using both approaches and the proteins assigned to them. The number of one-to-one relationships (indicating a complete agreement) between SB and DAB clusters is indicated in Table 2 and ranges from 648 (for H. pylori) to 1680 (in Pseudomonas) corresponding to 50% and 25% of the pan-genome. This indicates that results of SB and DAB clustering tend to be more similar when working at closer phylogenetic distances. However, more complicated cases occur when proteins in a single SB cluster are assigned to various DAB clusters including singletons and vice versa. An overview of the possible mismatches between SB and DAB clusters is in Figure 5. The observed frequency of the different types of cluster mismatches are given in Figure 6. We observed that single domain architectures predominated the one-to-one clusters as shown in Table 3.
For L. monocytogenes we found 378 1d → 1s DAB cluster mismatches, ( Figure 5, panel A, top case) meaning that in those cases sequences in a DAB cluster are a subset of the sequences in the corresponding SB cluster. This lower number of sequences in the DAB cluster could be due to, for instance an insertion or expansion of a domain, leading to SB clustered sequences with partly overlapping but distinct domain architectures as is depicted in Figure 1. Similarly, there are 399 1s → 1d clusters. Each of these cases represent a sequence cluster where all the sequences share the same domain architecture, but other sequences exist with the same architecture that have not been included in the cluster due to a too low similarity score. The low similarity between sequences with the same domain architecture could be due to a horizontal acquisition of the gene or to a fast protein evolution at the sequence level. Genes acquired from high phylogenetic distances could greatly vary in sequence while presenting the same domain architecture.
Proteins contained in a single DAB cluster but assigned to multiple SB clusters contain mostly ABC transporters-like (PF00005) or Major Facilitator Superfamily (MFS, PF07690) domains. This is not surprising considering that such generic functions are usually associated with a high sequence diversity. Conversely, ABC transporters are found in multiple DAB clusters. However, many of them are grouped into a single SB cluster with ATPase domain containing proteins (1s → Nd case).
We observed distinct architectures with one of two very similar domains, the GDSL-like Lipase/Acylhydrolase and the GDSL-like Lipase/Acylhydrolase family domain (PF00657 and PF13472 respectively) and those architectures were often seen clustered using a SB approach. However, architectures containing both domains were also identified, pointing to a degree of functional difference as a result of convergent or divergent evolution. Still, the corresponding sequences remain similar enough as to be indistinguishable when a SB approach is used. For SB clustering we also observed the case of identical protein sequences not clustered together, probably because of the tie breaking implementation when BBH are scored. In all cases we found the size of both the pan-and the coregenome to be larger when a SB approach is used to identify gene clusters and SB approaches lead to a larger number of singletons than DAB ones. This indicates that DAB clusters are assigned to several SB clusters, many of them consisting of just one protein.
When going from species to phylum level, the ratio between the number of DAB and SB singletons changes from 0.48 and 0.41 (for H. pylori and L . monocytogenes respectively) to 0.19 and 0.40 when considering organisms of a higher taxonomic level (Corynebacteriales and Cyanobacteria respectively).
We investigated the predicted size of the pan-genome upon addition of new sequences. Heaps' law regression can be used to estimate whether the pan-genome is open or closed 38 through the fitting of the decay parameter α ; α < 1 indicates openness of the pan-genome (indicating that possibly many clusters remain to be identified within the considered set of sequences), while α > 1 indicates a closed one; the α values are given in Table 4. In all cases the pan-genome is predicted to be open; however, α values obtained using DAB clusters (α DAB ) are systematically closer to one than the α SB obtained with the standard sequence similarity approach.
The α DAB value retrieved for L. monocytogenes is strikingly low. Heaps law regression relies on the selected genomes providing a uniform sampling of selected taxon, here species. Analysis of the domain content of the selected genomes shows a divergent behaviour of strain LA111 (genome id GCA_000382925-1). This behaviour is clear in Figure 7, where GCA_000382925-1 appears as an outlier of the L. monocytogenes group. Removal of these outlier leads to α DAB = 1.04 and α SB = 0.64, which emphasizes the need for uniform sampling prior to Heaps regression analysis.  DAB comparison across multiple taxa DAB clusters can be labelled by their domain architecture and since this is a formal description of functional equivalence, results of independently obtained analyses can be combined. Figure 7 shows the results of a principal component analysis of the combined DAB clusters for selected genomes from eight taxa. The first two components account for a relatively low explained variance (29%) still grouping of genomes from the same taxa is apparent. High functional similarity among genomes of the same species (H. pylori and L . monocytogenes) is reflected by the compact clustering, while phylogenetically more distant genomes appear scattered in the functional space defined by the principal components.

Discussion
We have shown that domain architecture-based methods can be used as an effective approach to identify clusters of functionally  can be stored and re-used for further analysis. This provides an effective scalable approach for large scale functional comparisons which by and large is independent of phylogenetic distances between species.
The chosen set of domain models and the database used as a reference greatly impact the results. InterPro aggregates protein domain signatures from different databases, which leads to redundancy of the domain models. This redundancy causes overlaps between the entries and an increase of the granularity of the clusters retrieved: this can bias downwards the size of the pan-genome and upwards the size of the core-genome, as shown in Table 1.
In InterPro this redundancy is taken into account by implementing a hierarchy of protein families and domains. The entries at the top of these hierarchies correspond to broad families or domains that share higher level structure and/or function; the entries at the bottom correspond to specific functional subfamilies or structural /functional subclasses of domains 30 . Using Inter-Pro for DAB clustering would require taking into account the hierarchy of protein families and domains: however, this would pose challenges of its own and would require discrimination of the functional equivalence of different signatures within the same hierarchy.
Another source of redundancy are functionally equivalent domains from distantly related sequences. Pfam represents this through related families, termed clans, where relationships may be defined by similarity of sequence, structure or profile-HMM. Clans might Figure 6. Comparison between DAB and SB clusters. On the left DAB is used as a reference and each bar represents the relative frequency of one DAB cluster containing sequences assigned to {1, 2, ... , 5} and 6 or more SB clusters and one-to-one represents the relative frequency of identical cluster. Similarly, on the right SB is used as a reference. Axis labels follow notation in Figure 5.  Members of a clan have diverging sequences and very often SB approaches would recognize the evolutionary distance between the sequences and group them in different clusters. If we were to assume that members of a clan are functionally equivalent and collect them in the same DA cluster, we will have a higher number of cases where a single DA cluster is split in multiple sequence clusters 1d→Ns. Also there would be higher number of cases of sequence clusters with the same DA but no exactly matching the DA clusters (1s→1d cases).
In many cases a one-to-one correspondence could be established between DAB and SB clusters indicating that often the sequence can be used as a proxy for function. At first this may seem a trivial result but it has a profound implication: domain model databases (in this case Pfam) contain enough information, encoded by known domain models, to represent the quasi totality of biological function encoded in the bacterial sequences analyzed here. However, it is important to stress that the comparisons have been performed considering sequences with known domains, representing currently around 85% of the genome coding content, a number that will only increase in the future.
A significant advantage of the DAB method over the SB method is that the domain architecture captured within a cluster can be used as a formal description of the function. Currently, more than 20% of all separable domains in the Pfam database, are so-called domains of unknown function (DUFs). Despite this, in bacterial species they are often essential 42 . With the DAB method they are formally included and often semantically linked to one or more domains of known function.
The starting position of the domains was used to generate labels indicating N-C terminal order of identified domains. The labels were used only for clustering as proteins sharing the same labels were assigned to the same clusters. Choosing instead the mid-point or the C-terminal position could affect the labeling but not the obtained clusters.
A content-wise formal labeling of DAB clusters makes a seamless integration of multiple independently performed DAB analysis possible. This allows for a comparison of potential functionomes across taxonomic boundaries, as presented in Figure 7, while new genomes can be added at a computational cost O(n), with n the number of genomes to be analyzed. SB methods that create orthologous groups require more memory and time as they come at an O(n 2 ) computational cost. Other SB approaches, such as COGNITOR, reduce the computational costs to O(n) by using precomputed databases. In this respect, the DAB approach is similar to the approach implemented in COGNITOR, by searching against existing databases of domains architectures. In this way the DAB approach leverages the extensive amount of work already put into defining domain families.
The bimodal shape of the distributions presented in Figure 3 and Figure 4 indicates the relative role of horizontal gene transfer and vertical descent when shaping bacterial genomes: the first peak accounts for sequences (or functions) only present in a small number of genome sequences which have been a likely acquired by horizontal gene transfer. The second peak accounts for high persistence genetic regions representing genes (or functions) belonging to the taxon core which have been likely acquired by vertical descent.
A measure of the impact of vertical descent and horizontal gene transfer is provided by the ratio between the core-and pan-genome sizes. The number of singletons provides a measure of the number of genes horizontally acquired from species outside the considered group.
Two of the most prominent differences between the two approaches are the number of retrieved singletons and the coreto pan-genome size ratio. Multiple members of the same taxon might acquire the same function through horizontal gene transfer 43 . This is likely to occur given that they would have similar physiological characteristics, hence they would tend to occupy a similar niche or, at least, more similar than when comparing species from different taxa. As the origin of the horizontally acquired genes may vary for each organism, an SB approach will correctly recognize the heterologous origin of the corresponding sequences and those will be assigned to singletons. However, the probabilistic hidden Markov models used for domain recognition are better at recognizing the functional similarity of the considered sequences and clusters them together.
Another indication of the relative impact of horizontal and vertical gene acquisition events is provided by the openness or closedness of the genome. Values for the decay parameter α in Table 3 indicate a relatively large impact of horizontal gene transfer.
Within the considered taxa we observed α DAB > α SB , meaning that the sequence diversity is larger than the functional diversity: upon addition of new genomes to the sample the rate of addition of new sequence clusters appears higher than the rate of addition of new functions.

Limitations of DAB approaches
We have shown that domain architecture-based methods can be used as an effective approach to identify clusters of functionally equivalent proteins, leading to results similar to those obtained by classical methods based on sequence similarity. However, whether DAB methods are more accurate than SB methods to assess functional equivalence will require further analysis. In this light, results of functional conservation for both approaches could be compared in terms of GO similarity and/or EC number 44,45 . Partial domain hits might arise as a result of alignment, annotation and sequence assembly artefacts. To reduce the number of partial domain hits additional pruning could be implemented to distinguish these cases. However, this is an open problem that requires caution as it could influence the functional capacity of an organism and clustering approaches using DA.
The performance of DAB methods may be sub-optimal when dealing with newly sequenced genomes that are not yet well-characterized enough to have all of their domains present in domain databases, since DAB methods will be unable to handle unknown architectural types. Around 15% of the genome coding content corresponds to sequences with no identified protein domains. DAB approaches can be complemented with SB methods to consider these sequences or even protein sequences with low domain coverage, possible indicating the location of protein domains yet to be identified. Since DAB methods rely on the constant upgrading of public resources like UniProt and Pfam databases, an initial assessment of domain coverage appears as a sine qua non condition for application of these methods. DAB approaches could be used to assess the consistency of existing orthologous groups in terms of their domain architectures, at least when domain architectures are expected to be completely known in advance (for instance in the case of micro-evolutionary variations within a species where mutational events may disrupt a protein's function). For other purposes, such as the discovery of a new phyla of cellular life that contains radically different domain architectures, global similarity methods may be preferred 45 .

Conclusions
As protein domain databases have evolved to the point where DAB and SB approaches produce similar results in closely related organisms, the DAB approach provides a fast and efficient alternative to SB methods to identify groups of functionally equivalent/ related proteins for comparative genome analysis. The lower computational cost of DAB approaches makes them the better choice for large scale comparisons involving hundreds of genomes.
Highly redundant databases, such as InterPro, are best suited for domain based protein annotation, but are not effective for DAB clustering if the goal is to identify clusters of functionally equivalent proteins. To enable DAB approaches for highly structured databases, such as Inter-Pro, the hierarchy of protein families and domains within has to be explicitly considered. Currently Pfam is for this task a better alternative.
Differences between DAB and SB approaches increase when the goal is to study bacterial groups spanning wider evolutionary distances. The functional pan-genome is more closed in comparison to the sequence based pan-genome. Both methods have a distinct approach towards horizontally transferred genes, and the DAB approach has the potential to detect functional equivalence even when sequence similarities are low.
Complementing the standardly applied sequence similarity methods with a DAB approach pinpoints potential functional protein adaptations that may add to the overall fitness.

Data availability
List of genomes used for the analysis at different phylogenetic levels. The genomes are grouped per taxonomic lineage used in this study.

Competing interests
No competing interests were disclosed.

Grant information
This work was partly supported by the European Union's Horizon 2020 research and innovation programme (EmPowerPutida, Contract No. 635536, granted to Vitor A P Martins dos Santos).
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Figure S1. Persistence of Sequence Based (SB) clusters. Cluster persistence is defined as the relative number of genomes with at least one protein assigned to the cluster. The plots show frequency of SB clusters according to their persistence. Publicly available and complete genome sequences assigned to each taxon were selected so that phylogenetic diversity within the taxon was preserved, as described in materials and methods. 60 distinct genome sequences were considered for each taxon shown.

Figure S3. Comparison between DAB and SB clusters.
On the left DAB is used as a reference and each bar represents the relative frequency of one DAB cluster containing sequences assigned to {1, 2, . . . , 5} and 6 or more SB clusters and one-to-one represents the relative frequency of identical cluster. Similarly, on the right SB is used as a reference. Axis labels follow notation in Figure 5. The manuscript is overall of much higher quality than it was previously. However, not all of the comments were addressed fully. The most egregious of these oversights is that it still gives the impression that a "straw man" argument is being set up to easily fall when set up against the authors' favored approach. Fortunately this is not quite the case, and therefore this is likely an entirely unintentional effect of not clearly explaining the method that is being used or its comparison to other methods, but even so this issue is quite important since it can easily mislead an unwary reader. e.g., statements about the time and memory requirements of the SB approach having to scale quadratically with the number of genomes to be compared are untrue. In fact, while SB methods to construct orthologous groups do often take advantage of a full all-against-all comparison (and therefore these methods require a quadratic scale), perhaps the more proper comparison of DAB is not to the set of SB methods that construct orthologous groups, but rather to those that extend existing groups into new genomes (much like DAB does not construct domain families, but merely extends these existing families into new genomes, taking into account their architectures while this is being done). This fact is now acknowledged in the manuscript, but is buried deeply in the middle of the Discussion section, and yet the confusing description of the comparison of the quadratic to linear scales also remains at several places in the manuscript (such as in the abstract and the second paragraph of the introduction). Perhaps this was merely an oversight, but in any case this issue should be made much more clear than it currently is. If I understand things correctly, the overall summary seems to be that: the DAB approach, much like several existing SB approaches, leverages the extensive amount of work (much of it done with manual curation) already put into defining domain families, and attempts to extend these families to identify new members of orthologous groups in newly discovered genomes -which it is able to do more accurately than similar SB methods due to taking into account domain architectures. Both of these approaches -DAB and SB -scale linearly with the number of new genomes to be compared. In contrast, there is also a different class of SB methods (such as OrthaGogue, and COGs) that create orthologous groups de novo -these methods require more memory and time since they scale quadratically with the number of new genomes, although this class of method provides the advantage of being able to work even in the absence of domain family information, which DAB is not able to do.

Open Peer Review
Another oversight occurs in the sentence that "Our aim was to investigate whether using HMMs instead of sequence similarity would yield similar results", where instead of HMM I think the authors meant domain architectures? (and since domain architecture comparisons also rely on sequence similarity, perhaps also add the word "alone" after "sequence similarity" to distinguish the use of the latter alone vs. in combination with domain architectures) No competing interests were disclosed.

Competing Interests:
No competing interests were disclosed.

Competing Interests:
Referee Expertise: Orthology identification, comparative genomics, bioinformatics I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Jasper Koehorst
Thanks, we have amended the manuscript as suggested by the reviewer.

Specifically:
We have deleted from the abstract and introduction the statements about the time and memory requirements of the SB approach having to scale quadratically with the number of genomes.
We have modified the paragraph in the discussion to further emphasizes that DAB is similar to SB methods that extend existing groups into new genomes. We have also rephase the reviewers comment on the extensive use DAB does on the amount of work put on defining domain families as we think it might further clarify the text.
The sentence "Our aim was to investigate whether using HMMs instead of sequence similarity would yield similar results" has been modified as suggested to: "Our aim was to investigate whether using domain architectures instead of sequence similarity alone would yield similar results." No competing interests were disclosed. The limitations of global sequence similarity based methods to identify proteins that perform similar functions are well-known. Thus, the approach described in this manuscript of using domain-based clustering of orthologous groups (DAB) represents an exciting and very welcome addition to the field. Or at least it will when it is fully developed, although this manuscript has not convinced me that it outperforms other methods at its current level of development, and I have several substantial reservations about some of its content: As the first reviewer also mentioned, methods such as CDART and DELTA-BLAST (published in 2002 and 2012, respectively) have been available for many years. The latter even seems to aim to perform the exact same function as DAB, by considering domain architectures. How is DAB perform the exact same function as DAB, by considering domain architectures. How is DAB different or better? I suspect that DAB may have greater accuracy since it uses HMMs rather than PSSMs, but this remains to be shown, and DELTA-BLAST is far easier for a user to run, since it is available as a webserver.
The comparison performed in this manuscript appears to fall prey to a straw man argument. In some cases, but not all, re-writing the relevant sections of the manuscript would help to avoid any misconceptions in this regard. a) The issue of replacing a cost with a one upon addition of a new genome was dealt with over O(n ) O(n) 15 years ago, so the statement "On the other hand, addition of a new genome using an SB approach require a new set of all-against-all sequence comparisons which come at a computational cost" is O(n ) false -at least as it is currently written. It is true that building groups of orthologs do require an initial O(n ) computational cost, but once those orthologous groups are formed, methods such as COGNITOR (first published in the year 2000) work extremely quickly and efficiently to assign genes in newly-sequenced genomes to existing groups. In fact, COGNITOR works in the exact same manner in which DAB uses pre-computed domain databases to achieve the much lower cost, although in COGNITOR's case it O(n) searches against a pre-computed database of orthologous groups (of which there are far fewer than domains, so with a smaller "n" it would actually be faster than DAB).
In should be noted that despite DAB's somewhat higher cost, it has the theoretical potential to achieve better accuracy than COGNITOR (at least in some cases) since as a global sequence similarity approach, the latter does not explicitly consider domain architecture. At least not in an automated fashion -doing so would require manual curation of its results, which is often done by careful researchers, but is not a process that is scalable to handle the ever-decreasing cost and ever-increasing amounts of genomic data. Although since a comparison with COGNITOR was not included in the manuscript, either in terms of speed or accuracy, it is unknown how much more useful DAB would be in practice. b) Even the initial cost does not have to be terribly burdensome, since the SIMAP method O(n ) pre-computes and stores BLAST results between all pairs of sequenced genomes anyway, and then uses efficient database retrieval methods to report the stored results. When a new genome is added, new O(n) comparisons have to be made -for a total accumulated cost of , although with the work spread out O(n ) over many years -and these in turn are useful for many other purposes, thus mitigating the construction costs. For instance, the EGGNOG database uses this method to build groups of orthologs. c) Why was only a single SB method chosen to be a representative for this entire class of approaches? Multiple forms of DAP were tested, whereas the only SB method used for comparison was one that uses a strict e-value cutoff of 1e-5, in the form of OrthaGogue and the OrthoMCL method. Also, why was the latter chosen to be this single representative? The latter approach was designed (nearly a decade and a half ago) for eukaryotic organisms, and while it has been applied more recently to bacteria as well, it is by no means the only -or even necessarily the best -approach for prokaryotic genomes. One advantage that it has is that it is completely automated, and thus is "easy" for people to use (even if, as this manuscript points out, horribly slow due to the procedure that it uses). On the other hand, methods like CDART O(n ) and COGTRIANGLES are all also automated (the latter of which uses no arbitrary e-value cutoff -that is, the results are robust to e-values over an immense range such as 1e-5, 1, 10, or even well beyond that on up to 100, or even 1000), and some pre-computed databases (such as COGs, representing the protein families present in the last common ancestor of all cellular life several billions of years ago) even take advantage of further manual validation, and from which pre-computed groups can be identified in newly-sequenced genomes by the fully automated and even easier approaches such as DELTA-BLAST 2 2 2 2 2 2 newly-sequenced genomes by the fully automated and even easier approaches such as DELTA-BLAST and COGNITOR. Is it at least possible that the poorer performance of SB methods in comparison to DAB as shown in the current manuscript is due to the choice of this particular SB method? I for one would have loved to see a comparison against the new release of the COGs database last year, since due to its being manually curated it acts as a sort of "Gold Standard" that can be tested against, with perhaps the EGGNOG groups being used as a more realistic measure of what a purely automated method can do without human supervision. Likely, DAB would fall somewhere in-between, and which would benefit the community of researchers who want to do comparative genomics of prokaryotic organisms to have a fully automated method that was demonstrated to surpass the existing fully automated methods. As it now stands though, DAB has only been shown to surpass OrthoMCL, which is not hard to do at all. Indeed, as seventh paragraph of the Discussion section (starting "Two of the most prominent...") states, unlike DAB, the SB methods were not able to cluster together the proteins with functional similarity but little sequence identity, especially across wider taxonomic ranges -which of course is what would be expected from a SB method that uses an e-value cutoff of 1e-5. d) Above and beyond the choice of SB method, it also seems that there may have been a bug in its implementation. The statement "For SB clustering we also observed the case of identical protein sequences not clustered together, probably because of the tie breaking implementation when BBH are scored." However, this was not supposed to happen, due to the within-species reciprocal BBH procedure that is used. In contrast, the tie breaking refers to between-species comparisons, but as shown in Figure 1 of the OrthoMCL paper (http://www.ncbi.nlm.nih.gov/pubmed/12952885), these two sources of information were supposed to have been combined together to form the final orthologous groups. If the proteins were highly similar (e.g., 99%) then perhaps a tie-breaking could be explained, but for 100% identical proteins -e.g., produced by a tanden duplication event -then they should have been collected into the group. One possibility is that this particular SB method simply was not designed to handle the large numbers of extremely closely-related genome assemblies that are available today, since at the time, very few instances of multiple genomic assemblies were available for the same species. If this explanation was demonstrated to be the reason why these identical proteins were not clustered together, that would be another reason for a user to choose to use DAB over this particular SB method. In any case (bug, design flaw, or something else), this event could greatly contribute to explaining some of the results that were observed whereby this single SB method found so many more singletons than DAB with Pfami.e., fixing the bug, or using some other SB method, may move many of those singletons into clusters. Although it would not explain why DAB with InterPro found even more singletons than this SB method? DAB has a lot of potential, but its limitations need to be made more clear: a) Why and how is the matrix of domain architecture binarized? Specifically, what if multiple copies of a domain are present? And does order matter -such as the architectures shown in Figure 2 of "A+B" and "B+A"? So, would "B+A+A" be a different architecture? And, as another reviewer also pointed out, what about "complicated" domain topologies where domains are interrupted by the insertion of another domain? Another major aspect of partial topologies is if DAB only recognizes some but not all of a newly-discovered architecture. E.g., a protein with architecture A+B+C+D, where A is known but B, C and D domains are not yet known. How would this be handled by DAB? Would it be reduced to appear merely as a single-domain "A" architecture? If so, how could that be distinguished from an architecture such as A+Z, which would also be reduced to appear just as a single-domain A? It seems like global sequence similarity methods might be more useful in those particular scenarios? i.e., if all the above domains were the same length, and a coverage threshold was used, then A+B+C+D could not be put into the same group as A+Z and A. Therefore, DAB seems primarily useful to either quickly extend known architectural types into a newly sequenced genome, but at the cost of not being able to work with unknown types. b) For newly sequenced genomes that are not yet well-characterized enough to have all of their domains 1.

2.
b) For newly sequenced genomes that are not yet well-characterized enough to have all of their domains present in the domain databases, DAB can be severely handicapped in comparison to global sequence similarity methods that do not have this limitation. In particular, Table 1 shows that up to nearly a fifth of the and Cornebacteriales genomes are not able to be assigned to domain families. Even these H. pylori numbers are merely lower-bound estimates, since brand-new architectures are expected to be discovered constantly, and yet these may incorporate at least one element that is known -such as the aforementioned A+B+C+D architecture, where only the A domain is represented in Pfam, but B and C and D are unknown. And yet it seems likely that even the fact that these domains are unknown would go unrecognized by the DAB approach -unless a factor is added to look for large segments of a gene that do not have matches in the databases of known domains. Therefore, the cost of DAB not being able to work with unknown architectural types might be quite high indeed. Worse, the exact value of that cost is also likewise unknown, and yet it would seem to be the single crucial piece of information that is most sorely needed in order to answer the question: does the benefits of DAB outweigh its costs?
If the goal is to bring together groups of proteins that have functional equivalence, then why was the only comparison that was done performed against the presence/absence membership of SB orthology approaches? Would it not have been better to actually measure the functional consistency observed within the SB groups, and within the DAB groups, in order to show that the latter was higher than the former? Many other methods that purport to improve upon the state-of-the-art orthology prediction process do just that -for instance Figure 4 of http://www.ncbi.nlm.nih.gov/pubmed/19148271 shows several comparisons with similarity of GO terms, enzyme nomenclature (EC), gene expression, and syntenic local neighborhood tests, with 12 different methods of orthology prediction. While neighborhood conservation is irrelevant for the issue of functional equivalence, the former three (or at least GO terms) would help to answer whether DAB is truly better than SB at the task of measuring functional equivalence. It would also help to answer whether this improved functional equivalence would be outweighed by the costs of being unable to handle unknown domain architectures, especially for highly divergent new genomes. If not, DAB may still be useful to check the consistency of existing orthologous groups in terms of their architecture, at least when domain architectures are expected to be completely known in advance -e.g., microevolutionary variations within a species where mutational events may disrupt a protein's function -but for other tasks such as the discovery of a new phyla of cellular life that contains radically different domain architectures, global similarity methods may be preferable instead.
Finally, some minor points concerning Figure 2: the vertical arrows seem to be pointing the wrong direction -a gene sequence undoubtebly contains more information content than a mere functional description. e.g., if I were to give you a GO code for molecular function, or biological process, then I could not tell you whether the original gene sequence is closer to one type of bacteria vs another type; but if I had the original gene sequence, then I could answer that question as well as many more.
I did not see a description of how amino acid coordinates are used anywhere else in the manuscript, either in DAB itself or in the comparison? In short, what does "Structure" have to do with anything, other than the general theoretical flow of "sequence begets structure which begets function"? If the purpose of Figure 2 is to describe the flowchart of DAB specifically though, it should focus only on the relevant elements. I suppose Structure could have meant how the sequence alignment was made, but if that were true, then DAB would only work for domain families for which a structure is available, instead of those for which only genomic or individual gene sequence has been provided. 3.

1.
sequence has been provided. Also, the last paragraph of the Discussion uses the word "closeness", but I think "closedness" was intended.
No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Jasper Koehorst
As the first reviewer also mentioned, methods such as CDART and DELTA-BLAST (published in 2002 and 2012, respectively) have been available for many years. The latter even seems to aim to perform the exact same function as DAB, by considering domain architectures. How is DAB different or better? I suspect that DAB may have greater accuracy since it uses HMMs rather than PSSMs, but this remains to be shown, and DELTA-BLAST is far easier for a user to run, since it is available as a webserver.

Following the suggestions made by the other reviewers we have added a paragraph in the Introduction regarding domain architectures, comparison of domain architectures and their use for sequence search. We have also discussed on how these have been included in domain databases and on the preservation of domain architectures at high phylogenetic distances.
We agree that most likely HMMs outperform PSSMs, however as the reviewer says, that is a topic that would required a dedicated investigation. Here our goal was to used domain architectures for functional comparative genomics, and we agree that a similar approach could be implemented using PSSM.
Regarding usability, we have used SAPP (semantic annotation platform with provenance) for genome analysis and annotation. SAPP is able to store the results in the RDF data model, that can be then queried using SPARQL. This tool is available with a web interface and is available at 2.

available with a web interface and is available at http://semantics.systemsbiology.nl/
The comparison performed in this manuscript appears to fall prey to a straw man argument. In some cases, but not all, re-writing the relevant sections of the manuscript would help to avoid any misconceptions in this regard.
The issue of replacing a O(n2) cost with a O(n) one upon addition of a new genome was dealt with over 15 years ago, so the statement "On the other hand, addition of a new genome using an SB approach require a new set of all-against-all sequence comparisons which come at a O(n2) computational cost" is false -at least as it is currently written.
We have amended the above mentioned sentence to: On the other hand, addition of a new genome using an SB approach require a new set of all-against-all sequence comparisons which come at a ( ) computational O n cost. However, approaches has been proposed to overcome this shortcomings of SB methods, such as COGNITOR which reduces the computational to ( ) by using O n pre-computed databases.
It is true that building groups of orthologs do require an initial O(n2) computational cost, but once those orthologous groups are formed, methods such as COGNITOR (first published in the year 2000) work extremely quickly and efficiently to assign genes in newly-sequenced genomes to existing groups. In fact, COGNITOR works in the exact same manner in which DAB uses pre-computed domain databases to achieve the much lower O(n) cost, although in COGNITOR's case it searches against a pre-computed database of orthologous groups (of which there are far fewer than domains, so with a smaller "n" it would actually be faster than DAB).
In should be noted that despite DAB's somewhat higher cost, it has the theoretical potential to achieve better accuracy than COGNITOR (at least in some cases) since as a global sequence similarity approach, the latter does not explicitly consider domain architecture. At least not in an automated fashion -doing so would require manual curation of its results, which is often done by careful researchers, but is not a process that is scalable to handle the ever-decreasing cost and ever-increasing amounts of genomic data.
We have commented on the analogy between DAB and COGNITOR.

In this respect, the DAB approach is similar in to the approach implemented in COGNITOR, by searching against existing databases of domains architectures.
Although since a comparison with COGNITOR was not included in the manuscript, either in terms of speed or accuracy, it is unknown how much more useful DAB would be in practice.
The focus of the paper was not to propose a comparative analysis of different methods but rather to present and contextualize the use of domain architecture for comparative genomics. However, we want to stress that we are not claiming that DA methods are superior to SB but that are an efficient and scalable alternative.
Even the initial O(n2) cost does not have to be terribly burdensome, since the SIMAP 2 5.

6.
Even the initial O(n2) cost does not have to be terribly burdensome, since the SIMAP method pre-computes and stores BLAST results between all pairs of sequenced genomes anyway, and then uses efficient database retrieval methods to report the stored results. When a new genome is added, O(n) new comparisons have to be made -for a total accumulated cost of O(n2), although with the work spread out over many years -and these in turn are useful for many other purposes, thus mitigating the construction costs. For instance, the EGGNOG database uses this method to build groups of orthologs.
Why was only a single SB method chosen to be a representative for this entire class of approaches? Multiple forms of DAB were tested, whereas the only SB method used for comparison was one that uses a strict e-value cutoff of 1e-5, in the form of OrthaGogue and the OrthoMCL method. Also, why was the latter chosen to be this single representative?
We have added the following to the discussion: To asses whether DAB results were consistent with those of SB methods we chosen. OrthaGogue as a representative of the latter class. Several tools such as COGNITOR and MultiPARANOID are available that implement different algorithm solutions to the task of identifying homologous sequences; however, despite different implementation, they all rely on sequence similarity as a proxy for functional equivalence. Here we considered SB methods as a golden standard for functional comparative genomics, especially when organisms within close evolutionary proximity are considered. Our aim was to investigate whether using HMMs instead of sequence similarity would yield similar results, thereby justifying their use for large scale functional genome comparisons. Regarding domain architectures, we have explored different alternatives, as we have seen that the chosen database or set of reference domains plays a critical role, an example is the low coverage of TIGRFAM preventing obtention of reasonable clusters.
The latter approach was designed (nearly a decade and a half ago) for eukaryotic organisms, and while it has been applied more recently to bacteria as well, it is by no means the only -or even necessarily the best -approach for prokaryotic genomes.
One advantage that it has is that it is completely automated, and thus is "easy" for people to use (even if, as this manuscript points out, horribly slow due to the O(n2) procedure that it uses).
On the other hand, methods like CDART and COGTRIANGLES are all also automated (the latter of which uses no arbitrary e-value cutoff -that is, the results are robust to e-values over an immense range such as 1e-5, 1, 10, or even well beyond that on up to 100, or even 1000), and some pre-computed databases (such as COGs, representing the protein families present in the last common ancestor of all cellular life several billions of years ago) even take advantage of further manual validation, and from which pre-computed groups can be identified in newly-sequenced genomes by the fully automated and even easier approaches such as DELTA-BLAST and COGNITOR. Is it at least possible that the poorer performance of SB methods in comparison to DAB as shown in the current manuscript is due to the choice of this particular SB method? I for one would have loved to see a comparison against the new release of the COGs database last year, since due to its being manually curated it acts as a sort of "Gold Standard" that can be tested against, with perhaps the EGGNOG groups being used as a more realistic measure of what a purely 7.

8.
perhaps the EGGNOG groups being used as a more realistic measure of what a purely automated method can do without human supervision. Likely, DAB would fall somewhere in-between, and which would benefit the community of researchers who want to do comparative genomics of prokaryotic organisms to have a fully automated method that was demonstrated to surpass the existing fully automated methods. As it now stands though, DAB has only been shown to surpass OrthoMCL, which is not hard to do at all. Indeed, as seventh paragraph of the Discussion section (starting "Two of the most prominent...") states, unlike DAB, the SB methods were not able to cluster together the proteins with functional similarity but little sequence identity, especially across wider taxonomic rangeswhich of course is what would be expected from a SB method that uses an e-value cutoff of 1e-5.
Above and beyond the choice of SB method, it also seems that there may have been a bug in its implementation. The statement "For SB clustering we also observed the case of identical protein sequences not clustered together, probably because of the tie breaking implementation when BBH are scored." However, this was not supposed to happen, due to the within-species reciprocal BBH procedure that is used. In contrast, the tie breaking refers to between-species comparisons, but as shown in Figure 1 of the OrthoMCL paper (http://www.ncbi.nlm.nih.gov/pubmed/12952885), these two sources of information were supposed to have been combined together to form the final orthologous groups. If the proteins were highly similar (e.g., 99%) then perhaps a tie-breaking could be explained, but for 100% identical proteins -e.g., produced by a tanden duplication event -then they should have been collected into the group. One possibility is that this particular SB method simply was not designed to handle the large numbers of extremely closely-related genome assemblies that are available today, since at the time, very few instances of multiple genomic assemblies were available for the same species. If this explanation was demonstrated to be the reason why these identical proteins were not clustered together, that would be another reason for a user to choose to use DAB over this particular SB method. In any case (bug, design flaw, or something else), this event could greatly contribute to explaining some of the results that were observed whereby this single SB method found so many more singletons than DAB with Pfam -i.e., fixing the bug, or using some other SB method, may move many of those singletons into clusters. Although it would not explain why DAB with InterPro found even more singletons than this SB method?
We have added a paragraph in the discussion regarding why the InterPro hierarchy has to be taken into account, also we mention this in the conclusion section. The hierarchical structure produces an increase in the domain multiplicity as many are related to each other. As a results an artificial variability in the DA is introduced leading to a higher number of singletons.
DAB has a lot of potential, but its limitations need to be made more clear.

We have added a new section to the Discussion: Limitations of DAB approaches
Why and how is the matrix of domain architecture binarized? Specifically, what if multiple copies of a domain are present?
We understand that our phrasing may have caused some confusion and we apologize for unclarity. The matrix of domain architectures is only binarized (presence/absence) to compute the PCA shown in Fig. 8, not to compare DAB and 8.

9.
(presence/absence) to compute the PCA shown in Fig. 8, not to compare DAB and SB clustering.
We have rephrased this in the Materials and Methods section: ...a binarized presence-absence matrix was obtained and used solely for principal component analysis.
[..] does order matter -such as the architectures shown in Figure 2 of "A+B" and "B+A"? So, would "B+A+A" be a different architecture? And, as another reviewer also pointed out, what about "complicated" domain topologies where domains are interrupted by the insertion of another domain? Another major aspect of partial topologies is if DAB only recognizes some but not all of a newly-discovered architecture. E.g., a protein with architecture A+B+C+D, where A is known but B, C and D domains are not yet known. How would this be handled by DAB? Would it be reduced to appear merely as a single-domain "A" architecture? If so, how could that be distinguished from an architecture such as A+Z, which would also be reduced to appear just as a single-domain A? It seems like global sequence similarity methods might be more useful in those particular scenarios? i.e., if all the above domains were the same length, and a coverage threshold was used, then A+B+C+D could not be put into the same group as A+Z and A. Therefore, DAB seems primarily useful to either quickly extend known architectural types into a newly sequenced genome, but at the cost of not being able to work with unknown types.
For newly sequenced genomes that are not yet well-characterized enough to have all of their domains present in the domain databases, DAB can be severely handicapped in comparison to global sequence similarity methods that do not have this limitation. In particular, Table 1 shows that up to nearly a fifth of the H. pylori and Cornebacteriales genomes are not able to be assigned to domain families. Even these numbers are merely lower-bound estimates, since brand-new architectures are expected to be discovered constantly, and yet these may incorporate at least one element that is known -such as the aforementioned A+B+C+D architecture, where only the A domain is represented in Pfam, but B and C and D are unknown. And yet it seems likely that even the fact that these domains are unknown would go unrecognized by the DAB approach -unless a factor is added to look for large segments of a gene that do not have matches in the databases of known domains. Therefore, the cost of DAB not being able to work with unknown architectural types might be quite high indeed. Worse, the exact value of that cost is also likewise unknown, and yet it would seem to be the single crucial piece of information that is most sorely needed in order to answer the question: does the benefits of DAB outweigh its costs?
The reviewer raises a very interesting point regarding how extensive available knowledge on protein domains is. The high agreement between the results of DAB and SB methods is only possible because databases of protein domains have enough information. Still, we believe many domains remain to be identified and in the scenarios the reviewer mentions DAB methods will be limited. We have added the following to the Discussion section, under the "Limitations of DAB approaches" header.
Still around 15% of the genome coding content corresponds to sequences with no identified protein domains. DAB approaches can be complemented with SB methods to consider these sequences or even protein sequences with low domain 10. methods to consider these sequences or even protein sequences with low domain coverage, possible indicating the location of protein domains yet to be identified.

We have extended the paragraph in the Materials and Methods where domain architectures are defined to further emphasize that N-C-terminal domain order is an inherent part of domain architecture definition.
Labels indicating N-C terminal order of identified domains were assigned to each protein using the starting position of the domains: the same labels were assigned to proteins sharing the same domain architecture.
In the Introduction we have added a paragraph regarding the use of protein domain architecture in protein annotations and we have included references to previous works showing that domain order is often key for the function of the protein and that domain duplications/insertions can also alter the function of the protein.

Moreover, a similar point on how domain architectures were defined and how the hierarchical relationships between protein domains, families and clans has been raised by R. Finn and a paragraph has been added in the Discussion (see answer to R. Finn's comments).
If the goal is to bring together groups of proteins that have functional equivalence, then why was the only comparison that was done performed against the presence/absence membership of SB orthology approaches? Would it not have been better to actually measure the functional consistency observed within the SB groups, and within the DAB groups, in order to show that the latter was higher than the former? Many other methods that purport to improve upon the state-of-the-art orthology prediction process do just that -for instance Figure 4 of http://www.ncbi.nlm.nih.gov/pubmed/19148271 shows several comparisons with similarity of GO terms, enzyme nomenclature (EC), gene expression, and syntenic local neighborhood tests, with 12 different methods of orthology prediction. While neighborhood conservation is irrelevant for the issue of functional equivalence, the former three (or at least GO terms) would help to answer whether DAB is truly better than SB at the task of measuring functional equivalence. It would also help to answer whether this improved functional equivalence would be outweighed by the costs of being unable to handle unknown domain architectures, especially for highly divergent new genomes. If not, DAB may still be useful to check the consistency of existing orthologous groups in terms of their architecture, at least when domain architectures are expected to be completely known in advance -e.g., microevolutionary variations within a species where mutational events may disrupt a protein's function -but for other tasks such as the discovery of a new phyla of cellular life that contains radically different domain architectures, global similarity methods may be preferable instead.

We have added the following section dedicated to limitations of DAB methods:
We have shown that domain architecture-based methods can be used as an effective approach to identify clusters of functionally equivalent proteins, leading to results similar to those obtained by classical methods based on sequence similarity. However, whether DAB methods are more accurate than SB methods to assess functional equivalence will require further analysis. In this light, results of 11. assess functional equivalence will require further analysis. In this light, results of functional conservation for both approaches could be compared in terms of GO similarity and/or EC number.
The performance of DAB methods may be sub-optimal when dealing with newly sequenced genomes that are not yet well-characterized enough to have all of their domains present in domain databases, since DAB methods will be unable to handle unknown architectural types.
Around 15% of the genome coding content corresponds to sequences with no identified protein domains. DAB approaches can be complemented with SB methods to consider these sequences or even protein sequences with low domain coverage, possible indicating the location of protein domains yet to be identified. Since DAB methods rely on the constant upgrading of public resources like UniProt and Pfam databases, an initial assessment of domain coverage appears as a sine qua non condition for application of these methods.

DAB approaches could be used to assess the consistency of existing orthologous groups in terms of their domain architectures, at least when domain architectures are expected to be completely known in advance (for instance in the case of micro-evolutionary variations within a species where mutational events may disrupt a protein's function). For other purposes, such as the discovery of a new phyla of cellular life that contains radically different domain architectures, global similarity methods may be preferred.
Finally, some minor points concerning Figure 2: The vertical arrows seem to be pointing the wrong direction -a gene sequence undoubtedly contains more information content than a mere functional description. e.g., if I were to give you a GO code for molecular function, or biological process, then I could not tell you whether the original gene sequence is closer to one type of bacteria vs another type; but if I had the original gene sequence, then I could answer that question as well as many more.
I did not see a description of how amino acid coordinates are used anywhere else in the manuscript, either in DAB itself or in the comparison? In short, what does "Structure" have to do with anything, other than the general theoretical flow of "sequence begets structure which begets function"? If the purpose of Figure 2 is to describe the flowchart of DAB specifically though, it should focus only on the relevant elements. I suppose Structure could have meant how the sequence alignment was made, but if that were true, then DAB would only work for domain families for which a structure is available, instead of those for which only genomic or individual gene sequence has been provided.
The ordering also seems unclear -wouldn't BBHs inform HMM domains, which then in turn inform domain architectures? Or if starting with BBHs, then how could architectures possibly be known prior to knowing the domains themselves? Or if it should be read from top to bottom as shown, how exactly does one start with Function (e.g., a GO term) and then, somehow via Structure, thereby arrive at a Sequence alignment? Specifically, is a Pfam entry a "Function", from which the Sequence alignment is downloaded? Or are Function and the Sequence alignment both part of the starting Pfam entry (and then again, what does any of that have to do with Structure)? From which domains are found (but aren't Pfam entries domains to begin with?), and then BBHs are made from the domain architectures? (an 12. 1.
domains to begin with?), and then BBHs are made from the domain architectures? (an extremely different way of doing the BBH procedure, which is normally done via Sequence alignments). In any case, as pointed out by other reviewers, the methodology used by DAB is not clearly explained in this figure, nor in the manuscript text.
We have edited the Figure for clarity incorporating the reviewer's suggestions.
Also, the last paragraph of the Discussion uses the word "closeness", but I think "closedness" was intended.

The typo has been amended.
No competing interests were disclosed. sequences for the purpose of performing comparative genomics. The principle behind the two approaches, sequenced based clustering and domain based clustering, is described well in the introduction. The motivation of the article is clear and well founded. However, the details provided about how domain assignments are actually performed and handled throughout the experiment generated so many questions, these have clouded the validity of any conclusions.
How was InterPro used to assign a domain architecture? As the database presents a hierarchy of protein families and domains, unlike Pfam and TIGRFAM, there are numerous overlaps between the entries. Some of these are trivial C-terminal to N-terminal overlaps, while others are complex arrangements that cannot be simply represented as described. If three overlapping domains from InterPro are in the same hierarchy, which domain is used? If all member databases are used, this will account for the explosion of clusters in the InterPro based-clustering seen in Table 1. If InterPro accessions are used (e.g. as seen in the condensed view of a sequence on the InterPro website) then numbers are surprising.
How were Family vs Domain "types" handled from InterPro or Pfam? In InterPro, type families tend to be near full length protein families. In Pfam, they represent a more heterogeneous bag of entries that are yet to be established as a 'domain'.
Pfam has a notion of related families, termed clans. Here the entries may not be intended to represent functionally distinct domains, but rather can represent a collection of families representing a continuum of evolution. How are entries belonging to a clan handled? How would the results differ if entries in one clan were treated as a single entity, for example, all P-loop NTPases as CL0023? How does this influence the sequence cluster to domain architecture relationships (schematicly shown in Figure 5).

1.
while others are complex arrangements that cannot be simply represented as described. If three overlapping domains from InterPro are in the same hierarchy, which domain is used? If all member databases are used, this will account for the explosion of clusters in the InterPro based-clustering seen in Table 1. If InterPro accessions are used (e.g. as seen in the condensed view of a sequence on the InterPro website) then numbers are surprising.
All member databases in InterPro were used. We partly took into account trivial Nterminal overlaps by alphabetically ordering the domains when distances between starting position were <3 amino acids. After analysing the results, we agree that this was not enough and this is the most likely cause of the explosion of this of clusters. As the reviewer suggests, taking the the full hierarchy of protein families and domains within InterPro would be required for comparative genome analysis based on domain architectures.
We have now better explained the selection criteria in the Materials and Methods section: The positions (start and end on the protein sequence) of domains having Pfam, TIGRFAMs and InterPro identifiers were extracted through SPARQL querying of the graph database and domain architectures were retrieved for each protein individually. InterPro aggregates protein domain signatures from different databases. Here no pruning for redundancies has been done. Identification of domains was done using the intrinsic InterPro cut-off that represents in each case the e-values and the scoring systems of the member databases. The domain starting position was used to assess relative position in the case of overlapping domains; alphabetic ordering was used to order domains with the same starting position or when the distance between the starting position of overlapping domains was <3 amino acids. Labels indicating N-C terminal order of identified domains were assigned to each protein in such a way that the same labels were assigned to proteins sharing the same domain architecture. We have commented on this point in the discussion, where the use of InterPro is addressed.
This paragraph now reads: The chosen set of domain models and the database used as a reference greatly impact the results. InterPro aggregates protein domain signatures from different databases, which leads to redundancy of the domain models. This redundancy causes overlaps between the entries and an increase of the granularity of the clusters retrieved: this can bias downwards the size of the pan-genome and upwards the size of the core-genome, as shown in Table 1. In InterPro this redundancy is taken into account by implementing a hierarchy of protein families and domains. The entries at the top of these hierarchies correspond to broad families or domains that share higher level structure and/or function; the entries at the bottom correspond to specific functional subfamilies or structural/functional subclasses of domains \cite{mitchell_interpro_2015}. Using InterPro for DAB clustering would require taking into account the hierarchy of protein families and domains: however, this would pose challenges of its own and would require discrimination of the functional equivalence of different signatures within the same hierarchy.

hierarchy.
We have also added the following to the conclusion To enable DAB approaches for highly structured databases, such as InterPro, the hierarchy of protein families and domains within has to be explicitly considered.
How were Family vs Domain "types" handled from InterPro or Pfam? In InterPro, type families tend to be near full length protein families. In Pfam, they represent a more heterogeneous bag of entries that are yet to be established as a 'domain'.
No distinction has been introduced as there don't seem to be general rules that apply to all cases. In the discussion section a paragraph has been added on the effects of the structure of the databases.
Pfam has a notion of related families, termed clans. Here the entries may not be intended to represent functionally distinct domains, but rather can represent a collection of families representing a continuum of evolution. How are entries belonging to a clan handled? How would the results differ if entries in one clan were treated as a single entity, for example, all P-loop NTPases as CL0023? How does this influence the sequence cluster to domain architecture relationships (schematicly shown in Figure 5).
The reviewer raises here an interesting point that we have now discussed. The following has been added to the first paragraph of the discussion section.
Another source of redundancy are functionally equivalent domains from distantly related sequences. Pfam represents this notion through related families, termed clans, where relationship may be defined by similarity of sequence, structure or profile-HMM. Clans might contain functionally equivalent domains, however it is not clear whether this is always the case as the criteria for clan definition includes functional similarity but not functional equivalence. Members of a clan have diverging sequences and very often SB approaches would recognize the evolutionary distance between the sequences and group them in different clusters. If we were to assume that members of a clan are functionally equivalent and collect them in the same DA cluster, we will have a higher number of cases where a single DA cluster is split in multiple sequence clusters 1d→Ns. Also there would be higher number of cases of sequence clusters with the same DA but no exactly matching the DA clusters (1s→1d cases).
Why was the N-terminal starting position used to assess position of the domain?
The following line has been rewritten in the Methods section Labels indicating N-C terminal order of identified domains were assigned to each protein using the starting position of the domains: the same labels were assigned to proteins sharing the same domain architecture.
What is the effect of choosing the mid-point?
We have commented on this in Results and Discussion. The following paragraph has been added: 5. 6.

7.
has been added: The starting position of the domains was used to generate labels indicating N-C terminal order of identified domains. The labels were used only for clustering as proteins sharing the same labels were assigned to the same clusters. Choosing instead the mid-point or the C-terminal position could affect the labeling but it not the obtained clusters.
Both Pfam and TIGRFAM use HMMER version 3, which uses local-local alignment algorithm. How are partial hits to an HMM handled? Would two partial domain matches that occur due to an insertion between two halves of a domain be treated differently (see Triant and Pearson, 2015)?
In the discussion we have added a subsection on the limitations on DAB approaches. There we have added the following:

Partial domain hits might arise as a result of alignment, annotation and sequence assembly artifacts (cite Triant
). To reduce the number of partial domain hits et al. additional pruning could be implemented to distinguish these cases. However, this is an open problem that requires caution as it could influence the functional capacity of an organism and clustering approaches using DA.
The use of domain architectures as an approach for accelerating sequence searching is not that novel, for example, CD-ART has been available for many years. Domain architecture views have been present in most domain databases (e.g. Pfam, SMART, Prosite) for over a decade, and used in genomic contexts. A more extensive overview of the use of domain architectures in the field is desirable.
We have added the paragraph in the introduction regarding domain architectures, comparison of domain architectures and their use for sequence search. We have also discussed on how these have been included in domain databases and, as also suggested by the first reviewer, on the preservation of domain architectures at high phylogenetic distances.
The following paragraph has been added to the introduction: The composite graphs presented in Figures 6, 7 and supplementary figures use different scales, so make the graphs hard to compare.

Figures 6 and 7 have been combined (also supplementary figures).
When the domain based clusters are compared to the sequence based clusters, it would be interesting to understand whether the number of domains that makes up the domain architecture influences the correlations to the sequence based clusters. Do single domain architectures predominated the 1:1 clusters?
We have looked into this and single domain architectures predominated the 1:1 clusters. A table has been added to the text (Table 3).
Many readers may be unaware of the thresholds employed in InterProScan relate to the individual databases, so greater clarity is required.
This point was also raised by A. Rosato. We have further explained the selected thresholds in the material and methods. Identification of domains was done using the intrinsic InterPro cut-off that represents in each case the e-values and the scoring systems of the member databases.
Why is the versioned InterProScan described as a semantic wrapper?
This line has been re-written, now it is explained that the versioned InterProScan stores the output in the RDF data model.

11.
information encoded in domain architectures has also been deployed to accelerate sequence search methods and to provide better homology detection. Examples are CDART (Geer 2002) which finds homologous proteins across significant evolutionary distances using domain profiles rather than direct sequence similarity, or DeltaBlast (Boratyn 2012) where a database of pre-constructed a position-specific score matrix is queried before searching a protein-sequence database. Considering protein domain content, order, recurrence and position has been shown to increase the accuracy of protein function prediction ( Table 1 shows that InterPro domains provide pangenomes that are not only always larger than the pangenomes obtained from Pfam domains but sometimes even larger than SB-derived pangenomes (e.g. H. pylori or Cyanobacteria). How is this possible?
InterPro aggregates protein domain signatures from different databases, which leads to redundancy of the domain models. This redundancy causes overlaps between the entries and an increase of the granularity of the clusters retrieved: this can bias downwards the size of the pan-genome and upwards the size of the core-genome, as shown in Table 1.
The low value of alpha in the Heaps regression for L. monocytogenes afforded by the DAB is striking and should be analyzed in more detail We thank the reviewer for this very interesting observation. We have investigated the low value of alpha in this case and the following paragraph has been added The alpha DAB value retrieved for is strikingly low. Heaps law L. monocytogenes regression relies on the selected genomes providing a uniform sampling of selected taxon, here species. Analysis of the domain content of the selected genomes shows a divergent behaviour of strain LA111 (genome id GCA\_000382925-1). This behaviour is clear in Figure 7 (PCA), where GCA\_000382925-1 appears as an outlier of the group. Removal L.monocytogenes of these outlier leads to alpha DAB=1.04 and alpha SB=0.64, which emphasizes the need for uniform sampling prior to Heaps regression analysis.
The line break after "transfer events" in the second paragraph of the introduction is not needed 11.

12.
The line break has been removed.

This typo has been fixed.
No competing interests were disclosed. Competing Interests: