PanOCT: automated clustering of orthologs using conserved gene neighborhood for pan-genomic analysis of bacterial strains and closely related species

System and software requirements

PanOCT was written in PERL (http://www.perl.org (1 August 2012, date last accessed)) and tested using PERL version 5.10 on Linux CentOS and Mac OS X 10.6 operating systems. To perform BLAST searches, either National Center for Biotechnology Information (NCBI) BLASTALL (23) version 2.2.10 or later (ftp://ftp.ncbi.nih.gov/blast/executables/release/ (1 August 2012, date last accessed)) or WUBLAST 2.0 (now called AB-BLAST available at http://blast.advbiocomp.com (1 August 2012, date last accessed)) are required.

Input requirements

PanOCT requires four input files. The first is either a NCBI (options -m 8 or 9) or WUBLASTP (using the btab program) tabular output file consisting of all-versus-all BLASTP searches of all predicted proteins in each genome to be analyzed. The second input requirement is a text file containing unique genome identifiers, one identifier per line, to determine which genome is to be treated as the reference genome in the output files and which genomes to include in the analysis. The genome identifier can be associated with specific proteins in two ways: (i) by placing the genome identifier after the protein identifier (e.g. NT08AB0001-GENOME_IDENTIFIER) or (ii) in the gene attribute file. The gene attribute file is the third required input file. It is a tab-delimited file containing the following data: contig id, protein identifier (e.g. locus), 5′-coordinate, 3′-coordinate, annotation and genome identifier. The final input file requirement is the protein fasta file used in the all-versus-all BLASTP searches. The protein fasta file is used by PanOCT to calculate the length of each protein, which is necessary in order to compute the BSR.

\In addition to the input file requirements mentioned above, PanOCT has a number of configurable command-line options: amino acid percent identity cutoff (default is 35%), BLAST E-value cutoff (default is 10⁻⁵), minimum percent match length of subject and query (default is 1%), frame-shift overlap parameter (default is 1.33; can be disabled), the number of amino acids at the beginning or end of a match that can be missing and still be considered a full length match between 0 and 100 (default is 20) and the number of blast matches needed to confirm a protein fragment/frame-shift (default is 1).

Output

PanOCT produces seven tab-delimited text output files plus a runtime report file. The most informative output files are the match_table and match_table_id files, which contain the ortholog clusters, one cluster per row. Each column contains protein identifiers with the first column beginning with the reference genome, followed by subsequent genomes ordered as instructed in the genome identifier file. The match_table_id file lists the percent identity for each protein to the reference in addition to the protein identifiers. The third output file has the following fields: locus identifier, annotation and each subsequent column containing the BSR (17). The fourth file is the frame-shift report, showing for each protein fragment the identity of retained and ignored fragments that are considered part of the same frame-shifted gene/pseudogene. Additional optional output files are available and are documented in the README.txt packaged with the distributed tarball. A detailed description of each output file is included with the source code.

PanOCT overview

All ortholog clustering methods start with a measure of protein similarity/homology based on alignment scoring between the proteins. The selection of orthologs is clear when only one protein per genome is homologous to each other (although even in this case compensating gene losses in each organism could cause non-orthologs (i.e. out-paralogs) to be considered orthologous). However, the choice is more complicated when multiple proteins per genome are homologous. When the species are closely related, CGN can be used as a contextual clue to discriminate between orthologs and paralogs. For example, two proteins from the same genome, P and P′, are homologs and have homology to proteins C and C′ from another genome. Based on homology alone, if P is more similar to C than to C′, and C is more similar to P than to P′, C might be assigned the ortholog of P based on RBH. However, by using CGN to distinguish orthologs when the proteins are nearly identical in addition to RBH, C′ can be assigned as the ortholog of P and C the ortholog of P′ (Figure 1).

There is circularity in the above example where CGN evidence for a choice of ortholog is based on prior identification of the orthologs for the neighboring genes. To address this problem, a measure of likely orthology using only homology information must be generated before determining the final orthology using CGN information. This measure of orthology is then combined within a window surrounding each potential ortholog pair (POP) to assign a weighted score (WS), including both homology and CGN measures.

Validation

The complete genome sequences of four Acinetobacter baumannii strains were downloaded from GenBank: ATCC 17978 (CP000521-523), AYE (CU459137-141), SDF (CU468230-233) and AB0057 (CP001182-1183). To obtain consistent structural annotation predictions across genomes, each genome was run through the JCVI annotation pipeline (24,25). These gene predictions were used to generate files containing the combined protein sequences of all four genomes and gene attributes for PanOCT (Supplementary Figures S1 and S2, respectively). NCBI BLASTP results were generated (E-value cutoff 1 × 10⁻⁵, filter off) and used by each of the four clustering methods. Each clustering method was run using default parameters with the exception of Sybil, where the Jaccard-clustering P-value was set to –1 to eliminate the inclusion of paralogs, which reduces Sybil to a simple RBH algorithm. An increase in the inflation parameter of OrthoMCL failed to exclude paralogs, so the default of 1.0 was maintained (data not shown). Cluster results from each of the clustering methods were compared. Clusters where all four methods agreed were assumed to be largely correct and uninformative for comparing relative performance. All other clusters were examined using ClustalW multiple sequence alignments, neighbor-joining trees, high-quality functional annotation evidence (HMM matches to equivalog level PFAM/TIGRFAM families (17)) and genomic context as evidence of orthology. This resulted in only 85 manually curated clusters containing 328 proteins that we felt were clustered with high confidence. The use of genomic context (CGN) as supporting evidence certainly biases the results toward PanOCT, which is the only method of the primary four to make use of CGN. That is why IONs and GOB, which also use CGN, were included in comparisons. These 85 manually curated clusters became the reference set used to evaluate the performance of each clustering method.

The clusters from each of the four primary ortholog clustering methods were given as input to IONs for four different runs. The set of clusters where all four primary methods agreed were used as the high confidence clusters or ‘pillars’ to be input to GOB.

The performance of each orthology detection method was evaluated using sensitivity (TP/(TP + FN)) and the positive predictive value (PPV; TP/(TP + FP)), where TP is the number of true positives, FN is the number of false negatives and FP is the number of false positives. The TP, FN and FP values were calculated over the set of pairwise orthology assertions inherent in each reference cluster. For example, if a reference cluster contained four orthologs (A, B, C and D), there would be six pairwise orthology assertions (A-B, A-C, A-D, B-C, B-D and C-D). For a clustering method, a TP is when the clustering method’s clusters inherently assert a pairwise orthology which is in the reference set, a FP is when the clustering method’s clusters inherently assert a pairwise orthology that is not in the reference set but at least one of the two proteins in the pairwise orthology is present in the reference set, and a FN is when the clustering method’s clusters do not inherently assert a pairwise orthology, which is present in the reference set (Supplementary Figure S3). Clustering methods were given credit for including a protein in a cluster when they included at least one protein fragment from that protein for frame-shifted genes/pseudogenes. Clusters were not penalized if they included multiple protein fragments from the same frame-shifted gene/pseudogene.

To determine the hardware demands of PanOCT compared with the other methods used in this study, resource statistics for each clustering method were tracked by submitting independent jobs to an Oracle Grid Engine (Oracle Corporation, Redwood Shores, CA, USA) grid, using the qacct command to gather information for each run. A total of 25 runs per method were launched. The input for the first runs contained three Escherichia coli genomes, whereas the final runs contained 25 genomes. Each run contained genomes used in the previous run plus one additional genome. All four methods were run on the same order of genomes to eliminate artifacts due to order of addition.

Comparing cluster membership

Except when all four methods agree, it is hard to directly compare clusters. This is because members of a single cluster from one clustering method could correspond to multiple clusters from another method, which may in turn correspond to different clusters from the original method. Therefore, instead of comparing clusters to evaluate the results of each clustering method, the cluster membership for each protein was evaluated. For each protein, two methods agreed if the protein was included in clusters with identical membership and disagreed otherwise. Of 6710 total non-redundant clusters containing 15 180 proteins, all four methods agreed for 86% of proteins (13 041) in 69% of the clusters (4631; Figure 4A). Three methods agreed and one disagreed: PanOCT, InParanoid and Sybil agreed for 4% of proteins; PanOCT, OrthoMCL and Sybil agreed for 3%; InParanoid, OrthoMCL and Sybil agreed for 1% and PanOCT, InParanoid and OrthoMCL agreed for <1% of proteins (Figure 4A).

Comparing to reference clusters

To directly compare and evaluate the behavior of each method, a reference set of 85 manually curated clusters was generated out of the set of clusters shown to disagree among the methods. Clustering method results were compared back with the reference set and cluster agreement or disagreement was scored (Figure 4B). Of the 85 reference clusters, PanOCT agreed with 85 (100%), InParanoid agreed with 19 (22%), OrthoMCL agreed with 38 (45%) and Sybil agreed with 65 (76%). These results can be split into 15 (2⁴–1) possible combined cluster prediction patterns of the four methods (Figure 4B). PanOCT, OrthoMCL and Sybil made the same cluster prediction for the largest number (36 or 42%) of clusters in agreement with the reference set. The next three major patterns of agreement with the reference clusters were PanOCT, InParanoid and Sybil (19 or 22%), PanOCT alone (18 or 21%), and finally PanOCT and Sybil (10 or 12%) (Figure 4B).

Performance

The performance of each method was further compared with the 85 reference clusters, using both sensitivity and PPV. PanOCT had the highest PPV of 1.000, while Sybil, InParanoid and OrthoMCL had PPVs of 0.958, 0.766 and 0.699, respectively (Table 2). PanOCT also had the greatest sensitivity (1.000), followed by OrthoMCL with 0.922, then Sybil and InParanoid at 0.863 and 0.692, respectively. Consistent with the findings in Table 1 where OrthoMCL had the largest number of paralogs in clusters, OrthoMCL had the highest FP rate of all four methods. InParanoid was the worst performer against the reference clusters, having the lowest number of TP, and the second highest FP rate. It is a little misleading to report sensitivity and PPV for only the 85 reference clusters, representing clusters that differed among the four clustering methods, since this is only a fraction of the number of clusters where all four methods agreed. If we presume that the clusters where all four methods agreed are overwhelming correct; hence, treating them as true positives, the PPV values for an expanded set of 4716 clusters (85 + 4631) were 1.000, 0.999, 0.995 and 0.990 for PanOCT, Sybil, InParanoid and OrthoMCL, respectively, and the sensitivity values were 1.000, 0.997, 0.993 and 0.998 for PanOCT, Sybil, InParanoid and OrthoMCL, respectively (Table 2).

CGN post-processing with IONS and GOB

To independently confirm that CGN improves clustering of orthologs, the previously determined clusters from each clustering method and the consensus clusters where all four methods agree were used as input for IONS (21) and GOB (20), respectively. IONS was able to correct some initial clusters for OrthoMCL (Supplementary Figure S8) and InParanoid (Supplementary Figure S9) by splitting clusters with paralogs into ortholog-only clusters using CGN (Table 2). IONs seemed to be particularly tuned for OrthoMCL clusters. For Sybil (Supplementary Figure S10) and PanOCT (Supplementary Figure S11), IONS could not improve clustering since there were no paralogs in the clusters but did split a few good clusters, degrading the performance slightly. IONS did not achieve the same performance as PanOCT when starting with clusters from another method (Table 2). GOB was able to use the confident ortholog clusters where all four primary methods agreed, to better determine separation of paralogs using CGN (Supplementary Figure S12), but still did not match PanOCT’s performance, measured by Sensitivity (Table 2).

Evaluation of differences

To better understand the reason for the observed differences in performance, each prediction pattern with disagreement to the reference clusters was manually evaluated and summarized in Supplementary Table S1. About 50% of the 85 references clusters were incorrectly identified due to the presence of paralogs, while the remaining 50% of reference clusters were missed because of differences in protein length (e.g. different start sites, truncations and splitting due to frame-shifts). There were several examples where OrthoMCL (Figure 5A), InParanoid (Figure 5B) or both failed to separate significantly diverged paralogs into ortholog clusters when Sybil and PanOCT do. It is not clear from our analysis why OrthoMCL differed from InParanoid on which diverged paralogs to include in a cluster. Where PanOCT differed in paralog separation from Sybil is for less diverged paralogs where homology is not sufficient and CGN information must be used (Figure 6). Looking outside of the 85 reference clusters, there were a few instances where PanOCT and/or Sybil seemed to make an arbitrary choice between paralogs when no convincing CGN or homology information existed (data not shown).

Scalability

Since PanOCT was created for pan-genome analysis, a much larger set of bacterial genomes was used to test the scalability of the PanOCT software. A graph of the reported maximum memory usage showed that PanOCT used more memory per genome than the other three methods (Supplementary Figure S13A). PanOCT memory usage is unchanged until the sixth genome is added, with a usage of ∼0.25 GB per genome, maxing out at 0.5 GB per genome by the 25th genome. As a result of an in-memory data storage strategy, PanOCT is able to finish orders of magnitude faster than Sybil and InParanoid, and in a fraction of the time it takes OrthoMCL, when run on identical data (Supplementary Figure S13B).

Robustness

To show that PanOCT is robust for large datasets, PanOCT was run on a set of 60 A. baumannii strains (Supplementary Table S2), including the original four used for analysis in this article. Clusters formed with the original four strains were compared with the clusters from 60 strains (Supplementary Figure S14). The clusters generated from the pan-genome analysis of 60 strains were projected onto the original four strains by removing all proteins not from the original four strains from the clusters and then compared with the original clustering results. A robust clustering method should minimize the clusters that change (split or merge) when new strains are added. For the 85 reference clusters, 1 cluster changed. For the 4361 clusters where all four methods agreed, 100 clusters changed. For the entire set of 5865 PanOCT clusters, 262 clusters changed. The small number of clusters that changed indicates that PanOCT is generally a robust clustering method. Almost all of the clusters that changed were joined by weak homology with a lack of CGN evidence, suggesting that in these instances, PanOCT may be too aggressive in assigning orthologs.