Identifying groups of nearly identical predicted protein sequences Identical or nearly identical sequences in the combined dataset were grouped in order to assess the overlap among the four gene prediction methods. A BLAST [22] search was performed (without low complexity filtering and with an e-value of e,0.001) for each sequence against the combined dataset. Sequences in each alignment were evaluated for overlap and sequence identity. For sequences longer than 100 residues, two sequences in an alignment were considered identical if the sequential overlap and the sequence identity within the overlap were both 99% or greater. If both the query and subject hit sequences 100 amino acids or smaller the two sequences were required to be within one amino acid in terms of length and sequence identity. Groups containing only NR sequences (672 groups) were not included in the analysis of the computational prediction methods. Clustering the protein dataset The combined, redundant dataset was clustered using CD-HIT [23]. The sequences are grouped by a given sequence identity threshold, which corresponds to the local sequence similarity. The longest sequence in each cluster is designated as the ‘‘represen- tative’’ sequence, which is a CD-HIT terminology. The sequence identities are calculated for the shorter sequences with respect to the representative. The requirement to compare each sequence to the longest sequence in the cluster is inherent in the clustering algorithm. However, we do not extend this terminology to designate the longest sequence to literally represent all the clustered sequences. The band-width allowance for the sequence alignments was set to 1,000 gaps in order to compare sequences with very different lengths but that may share a short common segment. Clustering was performed in the most accurate mode where each sequence is added to the cluster to which it is the most similar rather than to the first cluster for which the sequence identity threshold is met. CD-HIT implements a greedy incremental clustering algorithm, which orders the sequences by decreasing length and makes the longest sequence the represen- tative sequence of the first cluster. The next sequence is evaluated by being compared to the representative sequence of the original cluster and if it satisfies the sequence identity threshold it is added to the cluster. Otherwise it becomes the representative sequence of a new (second) cluster. Each successive sequence is evaluated by being compared to the representative sequences of each cluster and is either added to a cluster or becomes the representative sequence of a new cluster. One implementation of the algorithm clusters the sequences into the first cluster for which the minimum sequence identity threshold is satisfied. Thus, the first clusters could be artificially large. However, the clustering presented here evaluates each sequence against all the clusters and includes the sequence in the cluster for which it shows the highest sequence identity. In order to evaluate the redundancy in the dataset, the sequences were clustered for decreasing sequence identity thresholds and the resulting numbers of clusters were analyzed. A linear relationship exists between the number of clusters and the sequence identity threshold for sequence identity threshold levels from 95% to 70% (data not shown). There is not a specific level of redundancy in the dataset that dictates inherent clusters of the sequences. For each cluster within a CD-HIT sequence identity threshold level, the representative sequence was aligned to all other sequences in the cluster using BLAST2SEQ (without gaps and without low complexity filtering). The BLAST sequence identity Validation of Gene Predictions and alignment coverage were evaluated for each sequence alignment. The sequence identity was defined as the number of identities in the alignment divided by the alignment length. The alignment coverage was defined as the alignment length divided by the length of the shorter sequence. The frequency of intra-cluster BLAST alignment percent identities and alignment coverages were computed for various CD-HIT sequence identity threshold levels from 100% to 70%. As the CD-HIT sequence identity threshold level decreases the frequency of cluster members with lower quality alignments increases. However a distinct shoulder forms around the 90–95% sequence identity clustering level, therefore we identified the 90% sequence identity threshold as the lower limit of sequence identity that produces clusters of highly homologous proteins (Fig. 3a–b). Mapping protein sequences to genomic locations The genomic sequence for the ME49 strain was downloaded from ToxoDB.org. Each chromosome’s nucleotide sequence was translated in 3 consecutive reading frames in both the forward and reverse directions. A dataset of 2,622,472 potential open reading frames was created by dividing the 6-frame translations into any sequences that occurred between two genomic STOP codons. To determine the genomic location of a protein sequence a BLASTP (without low complexity filtering and without gaps) search was performed for the protein sequence against a custom database of the 2,622,472 generated open reading frames. The resulting local alignments were filtered to eliminate any alignments with an e- value greater than 1024 or a percent identity below 80%. For cases in which the protein sequence was spliced to multiple exons in the genome the BLAST search produced multiple, highly significant, local alignments rather than aligning the complete protein sequence to a single genomic interval. Our mapping algorithm examined all possible combinations of local alignments in order to identify cases in which a sequence could have multiple mapping options to the genome, i.e., duplicated genes. Each valid alignment configuration was converted to genomic nucleotide positions since the chromosomal locations of the open reading frames were known. This process was repeated for each local alignment within the protein query thus elucidating the physical mapping as well as the intron/exon structure of the protein sequence. 80% of the sequence was required to be mapped to the genome in order to constitute a successful mapping of the protein to the genome. There are 2,420 sequences, of the 30,197 sequences in the combined dataset, that can not be mapped to the genome; 155 sequences have no BLAST alignments to any of the translated genomic open reading frames, 376 sequences have no BLAST alignments that meet the significance requirements, and 1,889 proteins are not mapped because significant portions of the sequences cannot be aligned to the genome. The unmapped sequences are likely to be predicted from scaffolds that are not integrated into the genome. Genomic locations of sequence clusters The genomic locations of the protein sequences were used to localize CD-HIT clusters. The location of a genomic region of a CD-HIT cluster was defined by the sequence within the cluster that spanned the longest genomic region, including exons and introns. The genomic locations of each sequence, within a cluster, were compared in order to identify cases when proteins with similar sequences were derived from different genomic locations. The genomic locations of the sequence clusters were used to determine the distribution of the coding regions in the genome. Each chromosome was scanned in both directions. Clusters that had overlapping genomic regions and the same translation PLoS ONE | www.plosone.org 9 December 2008 | Volume 3 | Issue 12 | e3899