eCAMBer: efficient support for large-scale comparative analysis of multiple bacterial strains

The closure procedure

The closure procedure is the first step of eCAMBer. The input consists of genome sequences and genome annotations for a set of closely related bacterial strains. In this procedure gene annotations are iteratively transferred among the set of considered strains, until no new ORFs (open reading frames) are identified.More precisely, a gene annotation is transferred to a new location if its BLAST hit extended to the nearest in-frame stop codon is acceptable. Analogous to CAMBer, a BLAST hit extension to the nearest stop codon is acceptableif it satisfies the following conditions:

In this procedure eCAMBer, unlike CAMBer, takes advantage of working with closely related genomes. In contrast to the old approach, in each iteration, instead of using each ORF sequence as a query, it first identifies groups of ORFs with exactly identical sequences. This approach avoids use of the same ORF sequence multiple times as a BLAST query.

The pseudocode for the closure procedure implemented in eCAMber is given in Algorithm ??, which we now describe in more details. First, we start with the set of input annotations $A_s^0$, for each strain s in the set of considered strains S. Each ORF annotation (or simply ORF) is defined by a tuple (start, end, strand, contig, strain). Then, in ith iteration we compute the set of BLAST queries Q ⁱ as the set of distinct ORF sequences among all strains, which have not been used as BLAST queries yet. Next, we calculate in parallel, for each strain, BLAST results for all sequence queries in Q ⁱ . For each strain s∈S, all acceptable BLAST hit extensions $H_s^i$ are added to the strain annotations, defining $A_s^{i+1}\leftarrow A_s^i\cup H_s^i$. Next, the set of newly identified sequences across all genomes H ⁱ is computed, which is then used to update the set of BLAST queries for the next iteration Qⁱ⁺¹←H ⁱ ∖D ⁱ , where D ⁱ denotes the set of all distinct sequences before the ith annotation. The procedure stops when no new ORF sequences are identified, hence Q ⁱ =∅. For each strain s∈S, we denote by $A_s^c$ the set of annotations produced by the closure procedure above. We further denote by A ^c the set of all ORFs produced by the closure procedure.

Here, we also recall the notion of a multigene, introduced in our previous work [11], to account for the situation when multiple ORFs share the same stop codon in the annotations produced during the closure procedure. These ORFs are called multigene elements and represent putative gene translation units. Each multigene is represented by a tuple (e n d,s t r a n d,c o n t i g,s t r a i n,e l t s), where elts is the set of ORFs constituting the multigene. Also, for each strain s∈S, we denote by $M_s^c$ the set of multigenes resulting from the closure procedure.

Figure 2 presents a schematic view of the implementation of the closure procedure in eCAMBer.

The careful reader may also notice two important differences between the closure procedure in CAMBer and eCAMBer. In particular, eCAMBer uses unique ORF sequences, rather than ORF annotations, as queries against all strain genomes and, thus, does not repeat a BLAST query when the same ORF sequence corresponds to multiple ORF annotations. In contrast, firstly, CAMBer uses all ORF sequences as queries and, thus, may repeat a query BLAST several times. Secondly, CAMBer BLASTs a query against all strains’ genomes except the strain from which the query is taken. The second difference may potentially lead to different outcomes generated by these two approaches.

Since BLAST computations are the most time-consuming operation in each iteration of the closure procedure, we express the time complexity of one iteration of the closure procedure by the number of performed BLAST computations. Let k=|S| denote the number of considered strains and let $n={\mathit{\text{max}}}_{s\in S}\left|A_s^i\right|$ be the maximal number of gene annotations per strain, in iteration i. Let, d=|D ⁱ | denote the number of distinct gene sequences among all gene annotations in all considered strains. Then, the time complexity of one iteration of the closure procedure implemented in eCAMBer can be expressed as O(d·k), whereas it is O(n·k²) for CAMBer. Here, it should be noted that, potentially, if every annotated ORF sequence in S is different, then $\left|D^i\right|=\sum _{s\in S}\left|A_s^i\right|=O(n·k)$. However, as our case study experiments show, d is usually much smaller than n·k (see Figure 3).

Importantly, the number of I/O operations per iteration is also significantly decreased, from O(n·k²) in CAMBer to O(k) in eCAMBer.

Consolidation graphs

Having the closure procedure computed we represent its results in the form of graph structures, called consolidation graphs.

First, we introduce the conceptual representation, called the ORF consolidation graph. In this graph G _O =(V _O ,E _O ), each node o∈V _O represents an ORF annotation in $A_s^c$, for some s∈S. There is an undirected edge {o₁,o₂}∈E _O between a pair of ORFs, if there is an acceptable BLAST hit from the sequence of o₁ to o₂ or from the sequence of o₂ to o₁. We additionally assume, that there are no self-edges, i. e. o₁≠o₂.

Second, we recall the definition of the multigene consolidation graph, introduced in our previous work [11]. In this graph G _M =(V _M ,E _M ) each node m∈V _M represents a multigene in $M_s^c$, for some s∈S. There is an undirected edge {m₁,m₂}∈E _M between a pair of multigenes, if there is a pair of ORFs o₁∈e l t s(m₁) and o₂∈e l t s(m₂), such that there is an edge between them in the ORF consolidation graph (i.e., such that {o₁,o₂}∈E _O ).

Finally, we introduce the sequence consolidation graph, which is the structure used in the implementation of eCAMBer, as it is a compact representation of the information stored in the ORF consolidation graph and the multigene consolidation graph. In this graph G _S =(V _S ,E _S ,E _B ), nodes represent distinct ORF sequences. There are two types of edges, E _B called BLAST-hit edges, and E _S called shared-end edges. There is an undirected shared-end edge {x,y}∈E _S between a pair of sequence nodes if there is a multigene having two elements with these sequences. There is an undirected BLAST-hit edge {x,y}∈E _B between a pair of sequence nodes if there is an acceptable BLAST hit from x to an ORF with sequence y, or if there is an acceptable BLAST from y to an ORF with sequence x.

Figure 4 illustrates the correspondence between the ORF consolidation graph, sequence consolidation graph and the multigene consolidation graph.

Homologous gene clusters

The second step of eCAMBer is to determine homologous gene families as connected components of the multigene consolidation graph G _M . There is a natural one-to-one correspondence between the connected components of the multigene consolidation graph and the connected components of the sequence consolidation graph (the latter connected components are obtained by taking the union of E _S and E _B ). So, in eCAMBer, we do this using connected components of the sequence consolidation graph G _S , because it tends to be smaller for closely related genomes. The obtained set of homologous gene families is represented as a set of disjoint multigene clusters, denoted by C _M .

Refinement procedure

The third step of eCAMBer is the refinement procedure. The goal of the refinement procedure is splitting the homologous gene families, represented by multigene clusters, to obtain anchors. We call a multigene cluster an anchor, if it includes at most one multigene for every strain. Analogously, we call a multigene cluster non-anchor, if there is a strain which includes at least two multigenes in the cluster. Multigenes in the same anchor are potentially orthologous to each other, whereas a non-anchor contains at least two multigenes that are non-orthologous. Following CAMBer, we use genomic context information to decompose non-anchors into smaller multigene clusters that can emerge as anchors, as described below.

The input for the refinement procedure consists of the set of multigene clusters C _M , the sequence consolidation graph G _M , and the multigene annotations $M_s^c$, for each strain s∈S. We start with classifying the set C _M of multigene clusters into two disjoint sets of anchors and non-anchors, denoted C _A and C _N , respectively. We also sort all multigenes within strain contigs by positions of their stop codons. We reconstruct the subgraph of the multigene consolidation graph, called the refinement graph. In this graph G _R =(V _R ,E _R ), nodes V _R are constituted by the subset of multigenes, which belong to non-anchor clusters. There is an edge {m₁,m₂}∈E _R , between a pair of multigenes coming from different strains, if there is an edge {m₁,m₂}∈E _M , and the two multigenes belong to the same multigene cluster. By $E_R^{\{s_1,s_2\}}$ we denote the subset of edges between multigenes from a pair of strains s₁ and s₂. We omit details of the reconstruction of the refinement graph for brevity.

Then, for each unordered pair of strains {s₁,s₂} we perform the following procedure in parallel. First, for each multigene m we identify a pair of its nearest neighbours which belong to anchors with a multigene element present on the opposite strain. Such left and right neighbours of m are denoted as $l_m^{\{s_1,s_2\}}$ and $r_m^{\{s_1,s_2\}}$, respectively. Then, for each edge $\{m_1,m_2\}\in E_R^{\{s_1,s_2\}}$ we check whether it is supported in the sense that it satisfies one of the following conditions: (i) it connects multigenes belonging to a cluster, such that m₁ and m₂ are its only elements in strains s₁ and s₂; (ii) the corresponding pairs $(l_{m_1}^{\{s_1,s_2\}},l_{m_2}^{\{s_1,s_2\}})$ and $(r_{m_1}^{\{s_1,s_2\}},r_{m_2}^{\{s_1,s_2\}})$ belong to the same anchor; (iii) the corresponding pairs $(l_{m_1}^{\{s_1,s_2\}},r_{m_2}^{\{s_1,s_2\}})$ and $(r_{m_1}^{\{s_1,s_2\}},l_{m_2}^{\{s_1,s_2\}})$ belong to the same anchor. If any of the four neighbours does not exist we substitute it with a dummy node, which virtually belongs to any anchor.

Finally, we obtain the refined graph$G_R^{\ast }$ by removal of unsupported edges from G _R . Then, the set of connected components C _R of $G_R^{\ast }$ defines the set of multigene clusters after the split. Finally, we update the set of multigene clusters as $C_M^{\ast }\leftarrow (C_M\setminus C_N)\cup C_R$.

The careful reader may also notice the differences between the refinement procedures implemented in CAMBer and eCAMBer. First, the refinement procedure in CAMBer performs in iterations until no multigene clusters can be split. In eCAMBer the refinement procedure consists of only one iteration. However, since the input and output for the procedure are of the same type, it can be used multiple times, until no new clusters are split. Second, the condition for an edge to be supported in eCAMBer is more relaxed than that in CAMBer. Both approaches, for a pair of multigenes on different strains, identify pairs of their nearest left and right neighbour multigenes (belonging to anchor clusters with elements on both strains). However, CAMBer checks the actual presence of edges between the neighbours, whereas eCAMBer only checks if the identified neighbours match the same pair of clusters. This approach allows eCAMBer to avoid a costly reconstruction of the whole multigene consolidation graph.

TIS voting procedure

The fourth step of eCAMBer is the TIS voting procedure. The goal of the TIS voting procedure is to select the most reliable TIS for each multigene. To do this we implement an approach based on the concept of majority voting. This strategy has also been used to improve genome annotation accuracy in several recent studies [24, 31].

In this procedure, for each multigene m in each multigene cluster $c\in C_M^{\ast }$, we try to find a TIS (originally annotated or transferred) that belongs to a connected component of the ORF consolidation graph, where the connected component satisfies the following two conditions: (i) it has TISs (originally annotated or transferred) present in at least 80% of the multigenes in c; and (ii) it has TISs originally annotated in at least 50% of the multigenes in c, or it has TISs originally annotated in at least twice the number of multigenes in c than all other connected components in c. If such a TIS is found, it is selected as the TIS for m. If such a TIS is not found, but m has an originally annotated TIS, then the originally annotated TIS is selected as the TIS for m. If both of these two cases cannot be applied, the TIS corresponding to the longest ORF in the multigene m is selected. After the TIS voting procedure, every multigene has exactly one TIS selected. Thus, we obtain unambiguous TIS annotation for every gene.

Note that the connected components of the sequence consolidation graph—after shared-end edges have been removed— are in a natural one-to-one correspondence with the connected components in the ORF consolidation graph. So in eCAMBer, we implement the TIS voting procedure using the sequence consolidation graph, as it tends to be smaller for closely related genomes.

Clean up procedure

The last step of eCAMBer is the clean up procedure, which is designed to filter out multigene clusters which are likely due to gene annotation errors propagated during the closure procedure.

The input for this procedure consists of the set of multigene clusters $C_M^{\ast }$ and multigene annotations $M_s^c$, for each strain s∈S. For each multigene cluster $c\in C_M^{\ast }$ we compute the following features: (i) l, the median multigene length in c; (ii) p, the ratio of the number of strains with at least one element from c to the total number of strains; (iii) r, the ratio of the number of strains with at least one originally annotated multigene to the total number of strains with at least one element from c; (iv) v, the ratio of the number of multigenes in the cluster that are overlapped by a longer multigene to the total number of multigenes in the cluster.

Then, we update the set of multigene clusters $C_M^{\ast }$, by removing of multigene clusters for which ($p<\frac{1}{3}$ or $r<\frac{1}{3}$) and (l<150 or v>0.5).

Other features

In order to make eCAMBer more user friendly we have added a functionality for downloading genome sequences and genome annotations from the PATRIC database, for the set of selected strains within a species. The downloaded data is automatically formatted as input for eCAMBer. Additionally, eCAMBer integrates Prodigal to generate input gene annotations.

Furthermore, eCAMBer generates output compatible with CAMberVis [27], a tool for simultaneous visualization of multiple genome annotations of bacterial strains. CAMBerVis also handles visualization of genome annotation inconsistencies.

Comparison of running times

The above results also suggest that eCAMBer scales well to larger datasets.

Large case studies

We examine the scalability of eCAMBer to large datasets by running it on 10 datasets for the 10 species with the highest number of sequenced strains in the PATRIC database [2], in the 16 March 2013 release. All datasets consist of genome sequences and annotations for the sets of strains within the same species. Experiments for all of these datasets were conducted on a machine with 24 processor cores, out of which 20 were used.

In order to further investigate the scalability of eCAMBer, we check how the number of distinct gene sequences increases, when more strains are included. For this experiment, we chose the largest dataset of 569 strains of E. coli. Next, we sorted all genomes from the smallest to the largest. The plots (Figure 3) present the number of annotated genes and the number of gene sequences in a cumulative manner. We observe that the total number of distinct sequences grows much slower than the total number of gene annotations, suggesting sub-linear growth of the number of distinct gene sequences. Thus, according to our theoretical considerations, the algorithm implemented in eCAMBer for computing the closure procedure is sub-quadratic with respect to the number of strains included.

This experiment also shows that the strategy applied in eCAMBer to work with unique ORF sequences, rather than ORF annotations, leads to a sequence consolidation graph that is significantly smaller than the corresponding ORF consolidation graph. For example, in the largest dataset for 569 strains of E. coli, there is about 12.4mln nodes (ORF annotations) and 2.8bln edges in the ORF consolidation graph, whereas there are only about 1.6mln nodes (unique ORF sequences), 1.3mln shared-end edges, and 55.9mln BLAST-hit edges in the sequence consolidation graph.

Annotation consistency

We also investigate ability of eCAMBer to identify annotation inconsistencies and to improve the consistency of annotations. As a case study, we use the set of 20 E. coli strains with manually curated annotations, deposited in the ColiScope database [5], available through the web-based interface MaGe [32]. Pseudogenes were excluded from the analysis. On this dataset we run the closure procedure, followed by: the refinement procedure, the TIS voting procedure, and the clean up procedure. For comparison we also include annotations for the same set of strains, but downloaded from the PATRIC database [2].

In order to assess the improvement of annotation consistency, after running eCAMBer, we calculated the mean absolute difference in the number of annotated multigenes between two neighbour strains. It is 311 for the original annotations from ColiScope vs. 159 after applying eCAMBer. Analogous statistics on the dataset from PATRIC are 409 for the original annotations and 311 after applying eCAMBer.

In the dataset of 20 E. coli strains from ColiScope database, after the closure procedure, eCAMBer identifies 73 gene families which have the following property: each family has a member in every strain, and for each family exactly one strain has a missing original annotation in that family. The top three strains with the highest number of missing gene annotations of that type are: Sd197 (13), 2a 2457T (8) and 536 (7). The most well-studied strain K-12 MG1655 has four missing annotations of the above described type. These annotations were added by eCAMBer during the closure procedure.

Based on this case study, we also investigate how eCAMBer improves consistency of TISs. There are 8038 pairs of originally annotated genes with different TISs, but with identical sequence (including 100bp. upstream region from the TIS of the longer annotation). This number was reduced to 4230 after applying the TIS majority voting procedure and the clean up procedure.

This case study also shows that inconsistencies, which come from annotation errors, are present even for a very well-studied bacterial organism like E. coli. Note also that the discussed annotation inconsistencies were identified among strains with annotations curated by the same laboratory.

Comparison of features of eCAMBer and other tools

Different tools also aim in solving different annotation problems. For example, the GMV pipeline only identifies and solves TIS annotation inconsistencies, whereas Mugsy-Annotator also tries to identify missing genes. Our new tool, eCAMBer, is capable of resolving TIS inconsistencies, as well as removal of overannotated genes and addition of missing genes (Table 4). Our previous tool only identifies annotation inconsistencies, but it does not propose corrections.

Support for multithreading is a valuable feature for computationally demanding problems. Thus, it should be noted that eCAMBer has the most comprehensive support for multithreading among the tools considered. It allows the use of multiple threads for each of its steps. The GMV pipeline and CAMBer support multithreading only for BLAST computations. Mugsy-Annotator does not support it (Table 4).

Evaluation of annotation accuracy

In order to evaluate accuracy of annotations produced by eCAMBer, Mugsy-Annotator and the GMV pipeline, we apply the tools to annotations produced by the automatic annotation pipeline in PATRIC [2] for the set of 20 E. coli strains with manually curated annotations in the ColiScope database [5]. As an alternative dataset of input annotations for the same set of strains we use annotations generated using Prodigal [16].

In all our comparative experiments we run Mugsy-Annotator and the GMV pipeline with default parameters. It should also be mentioned that both Mugsy-Annotator and the GMV pipeline output lists of suggestions of changes to input annotations, rather than actually output the corrected annotations. We post-processed these proposed lists of changes to generate the output annotations used for the comparative experiments.

First we assess the correctness of the changes introduced to the input annotations based on the dataset of gene annotations with experimental support available in the EcoGene 3 database [31]. This dataset consist of 922 gene annotations for the K-12 MG1655 strain. From this set we excluded four genes: fdhF, prfB, rph’, insN’; since their sequences corresponding to the annotated coordinates are disrupted (the length of the sequence from the start codon to the stop codon is not divisible by 3). Additionally, we ran one iteration of the eCAMBer closure procedure to transfer the set of 918 gene annotations on the remaining 19 strains. The transferred gene annotations share at least 80% of sequence identity with original annotations for strain K-12 MG1655.

Since the extended dataset of annotations from Ecogene 3 constitutes only about 20% of all genes in the 20 strains of E. coli it is not sufficient for direct assessment of overall quality of changes introduced by eCAMBer and other tools. In particular we cannot conclude if a gene annotation is correct or not based on its absence in this dataset (so that there is no gene annotations in the dataset sharing the same stop codon). Thus, we perform further assessment of the quality of changes introduced relying on manually curated annotations for the set of 20 E. coli strains in the ColiScope dataset [5]. It is a reasonable choice as a gold standard, since many of the annotations have experimental support. In particular, the annotation for the strain K-12 MG1655 contains 901 out of 918 gene annotations present in the dataset described previously. For comparison, for this strain, there are only 841 and 883 such gene annotations for PATRIC and Prodigal, respectively.

Next, Figure 5 presents the assessment of TIS changes introduced during the TIS voting procedure based on the ColiScope dataset. It shows the assessment of the TIS changes introduced to the input PATRIC annotations, with respect to each of the 20 E. coli strains. Statistic presented in this figure distinguishes three different scenarios: (i) a correct TIS annotation is changed to an incorrect one (orange); (ii) an incorrect TIS annotation is changed to another incorrect TIS (yellow); (iii) an incorrect TIS is changed to the correct one (green). Since for each gene, there is only one TIS annotation considered as correct, there is no possible change from one correct TIS to another one. For each strain the majority of TIS changes introduced by eCAMBer is correct. Additional file 2 shows analogous panel figures for results of eCAMBer, Mugsy-Annotator and GMV on both PATRIC and Prodigal annotations. Rows 5 to 8 of Table 6 summarize the overall impact of eCAMBer and Mugsy-Annotator on TIS annotations. Remarkably, 70% (1591 out of 2260) of TIS changes introduced by eCAMBer to PATRIC annotations were correct. For comparison, only 43% of the TIS changes introduced by Mugsy-Annotator were correct.

Figure 6 presents the assessment of gene additions and removals introduced during the closure and the clean up procedures, respectively. It shows the assessment of the changes introduced to the input PATRIC annotations, with respect to each of the 20 E. coli strains. Statistic presented in this figure distinguishes four different scenarios: (i) a missing genome annotation is correctly added during the closure procedure (blue); (ii) a wrong gene annotation is correctly removed during the clean up procedure (green); (iii) a wrong gene annotation is incorrectly added during the closure procedure (red); and (iv) a correct gene annotation is incorrectly removed during the clean up procedure (orange). It can be seen that, for each strain, the majority of changes introduced by eCAMBer is correct. Additional file 3 shows analogous panel figures for results of Mugsy-Annotator and eCAMBer on both PATRIC and Prodigal annotations. The first four rows of Table 6 summarize the overall impact of eCAMBer and Mugsy-Annotator on gene presence. The results show that eCAMBer outperforms Mugsy-Annotator in this aspect. For example, 70% of the changes introduced by eCAMBer to PATRIC annotations were correct, whereas it was only 26% for Mugsy-Annotator.

Finally, we investigate how the whole pipelines implemented in eCAMBer, Mugsy-Annotator and GMV improve the overall annotation accuracy. Here, the accuracy is measured by f₁ statistic, defined as $2·\frac{\mathit{\text{precision}}·\mathit{\text{recall}}}{\mathit{\text{precision}}+\mathit{\text{recall}}}$, where $\mathit{\text{precision}}=\frac{\mathit{\text{TP}}}{\mathit{\text{TP}}+\mathit{\text{FP}}}$ and $\mathit{\text{recall}}=\frac{\mathit{\text{TP}}}{\mathit{\text{TP}}+\mathit{\text{FN}}}$. Here, TP, FP and FN denote true positive, false positive and false negative prediction, respectively. Since a pair of gene annotations may have the same stop codon, but different TISs, we keep track on the results for both stop codon predictions and for the TIS predictions.

Results of eCAMBer on PATRIC annotations in this experiment are presented in Figure 7. Note that each correctly identified TIS determines also its correctly identified stop codon, but not the other way round. Thus, the accuracy for the TIS prediction is lower than for the stop codons. As the figure shows, eCAMBer improves annotation accuracy, for each strain, both in terms of TIS annotations and stop codon annotations. Additional file 4 shows analogous panel figures for results of eCAMBer, Mugsy-Annotator and GMV on both PATRIC and Prodigal annotations. Rows 9 and 12 of Table 6 summarize the change in accuracy when running different tools on PATRIC and Prodigal annotations. It is clear from this table that eCAMBer outperforms other tools. For example, eCAMBer increased the f1 statistic of initial annotations of Prodigal (for gene starts) from 0.764 to 0.775, whereas the application of GMV improved it only by 0.001 and the application of Mugsy-Annotator decreased it by 0.027. In the case of PATRIC annotations, application of Mugsy-Annotator improved the accuracy from 0.680 to 0.682. However, the accuracy of annotations after eCAMBer increased to 0.703.