Sequence space of proteins is redundant and generates only a limited number of folds, domains, and structures [26]. Various strategies have been devised that take a coarse-grained approach to account for the degeneracy of sequences by grouping similar amino acids together [17-19,27-30]. Grouping is usually carried out based on structural and physiochemical similarities of amino acids [28]. Grouping of amino acids in sequence space can help develop prediction methods for various sequence determinants and decrease the amount of search space in procedures employed in directed evolution experiments [26,31]. One of the finest examples is the reduction of amino acid alphabet into a binary code that is composed of characters representing polar and non-polar amino acid residues [27]. Grouping of amino acid residues has also been used extensively in Hydrophobic-Polar (HP) lattice model to explain the hydrophobic collapse theory of protein folding [32].

A recent study was carried out by Peterson et al to test the performance of over 150 RAAAs on the sequence library from DALIpdb90 database and showed that RAAAs improves sensitivity and specificity in fold prediction between protein sequence pairs with high structural similarity and low sequence identity [33].

We tested performances of six amino acid reduction schemes with 15 different level of groupings to separate proteins into functional subfamilies (Table 2). These included three top performing RAAA (HSDM17, SDM12, GBMR4) from reference [33] and three random RAAA of size 4.

Amino acids that are within the same group in a RAAA are considered identical [33]. Substitution matrices that assign the same similarity score to each amino acid within the same group were obtained from reference [33]. For those RAAAs in the EB scheme and the three random RAAAs, new substitution matrices were created from BLOSUM62 frequency counts using the same procedure outlined in reference [33].

In this paper, a normalized distance measure that was previously used for phylogenetic tree construction of whole genome sequences was employed. The distance measure was based on Lempel-Ziv [34] complexity and was known to accurately cluster all related genomic sequences under one branch of the tree [8].

Lempel-Ziv (LZ) complexity score of a sequence is obtained by counting the number of steps required to generate a copy of the primary sequence starting from a null state. At each step, an amino acid or a series of amino acids are copied from the subsequence that has been constructed thus far allowing for a single letter innovation. The number of steps needed to obtain the whole sequence is identified as the LZ-complexity score of the given sequence. The exhaustive library of a sequence is defined as the smallest number of distinct amino acid or amino acid combinations required to construct the sequence using a copying process described by Lempel and Ziv [34]. For example, the LZ-complexity of the simple sequence 'AAILNAIIANNL' would be obtained as shown in Table 3. Since seven steps are needed to generate the whole sequence, the LZ-complexity score for this sequence is 7. The LZ-complexity of a sequence 'X' compared to a sequence 'Y' is known as the RCM of 'X' with respect to 'Y'. This is the number of steps required to construct sequence 'X' beginning with 'Y' instead of a null sequence. Five different distance metrics have been suggested by Otu et al [8], however, this work used only the following normalized distance metric that accounts for the differences in sequence lengths:

D X Y = c ( X Y ) + c ( Y X ) − c ( X ) − c ( Y ) c ( X Y ) + c ( Y X ) 2

where c(XY) and c(YX) are RCM of X appended to Y and Y appended to X, respectively. Remaining four LZ-based distance measures defined in Out et al performed slightly worse than the above distance (data not shown). Although in performance between five measures were not significant, we adopted the aforementioned distance for its ability to account for length variance.

The relative complexity measure (RCM) for creation of the distance matrix was utilized as previously described [8]. Phylogenetic trees were generated from distance matrices using neighbor-joining [35] program of the phylogeny inference package, PHYLIP 3.68 [36]. Un-rooted trees were rooted with midpoint rooting by placing the root halfway between the two most distinct taxa. Midpoint-rooted trees were converted to cladograms (i.e., branch lengths are discarded) using the Retree program of PHYLIP package [36]. Phylogenetic trees for all protein families and RAAAs are shown in supplementary materials (Additional File 2) in Newick format and can be visualized with a tree-drawing program.

Protein sequences in each family were aligned using ClustalW2 [22] for comparison with RCM. MSAs were performed using updated substitution matrices with gap extension and gap opening penalties provided in Table 2. Bootstrap analyses were carried out 100 times and trees containing bootstrap values were created using ClustalW2 with the neighbor-joining clustering algorithm. For convenience, MSAs that were carried out using ClustalW2 will be referred as the MSA or the MSA method for the rest of the article.

TBC algorithm [4] was used to check the accuracy of each tree in separating protein families into subfamilies. TBC divides a tree into disjoint subtrees and assigns a protein subfamily to a subtree that maximizes the number of true positives when the proportions of fp/(tp+fp) and fn/(tp+fn) are both equal to 0.5 for a given subtree, where fp is the number of false positives, fn is the number of false negatives and tp is the number of true positives. Above proportions correspond to the "maximal allowed contamination" level that minimizes the TBC error over the whole tree.

TBC requires a bifurcating tree of sequences in a protein family and an attribute file that contains expert curated assignment of each sequence to a particular subfamily. TBC accuracy (i.e., the percentage of correctly classified sequences) is the primary performance measure to evaluate the division of protein families into subtypes using the TBC algorithm. TBC accuracy is equal to 1- %TBC error where %TBC error is the total number of fp, fn, and unclassified sequences divided by the total number of sequences. For a detailed analysis of the TBC algorithm, refer to reference [4].

The proposed algorithm operates on a set of sequences in FASTA format. After one of the alphabets given in Table 1 is applied to all the sequences in the dataset, RCMs are calculated and used to obtain the distance between each pair for the neighbor-joining clustering to create a phylogenetic tree. For each RAAA, a single tree based on RCM is generated and analyzed using TBC algorithm to determine how well it clusters different subfamilies under different branches of the tree.

For simulated dataset, three phylogenetic trees were compared: The true tree generated by INDELible, the bootstrap tree and the RCM tree. INDELible creates a true MSA of the simulated protein sequences. This alignment was used in ClustalW2 and bootstrapped 1000 times and the resulting tree was called the bootstrap tree. The third tree is the RCM tree that was generated by the proposed approach.

For seven protein datasets, first, the original fasta sequences were used to calculate RCMs and their associated RCM trees. Second, the original fasta sequences were re-coded using different RAAAs (Table 2) and the reduced sequences were used to calculate their RCMs and the associated RCM trees.

A similar procedure was applied to the phylogenetic trees using the MSA method. For each protein family, MSA was carried out using the corresponding substitution matrices and gap penalties provided in Table 2. MSA-based trees were created following bootstrap analysis (100 replicates) with ClustalW2.

Finally, for each family, a total of 16 phylogenetic trees (1 for 20-letter alphabet, 12 for RAAAs, and 3 for random RAAAs) for each method are generated and checked how well they separated families into subfamilies. A summary of the overall workflow is depicted in Figure 1.

Phylogenetic analysis of protein sequences has been intimately connected with MSA. A phylogenetic tree is generated from an evolutionary distance matrix using MSA of sequences. However, for real biological datasets, the true tree is rarely known. Therefore, protein sequence evolution was simulated to study the reliability of the RCM method. A simulated protein dataset containing 10 protein sequences was generated to show that RCM coupled with a RAAA can produce a phylogenetic tree (RCM tree) that is consistent with the true tree and the bootstrap tree. The true tree is produced by INDELible and is the original tree that reflects the evolution of 10 simulated sequences. On the other hand, the bootstrap tree is the tree that was produced by ClustalW2 using the true MSA implied by INDELible. The bootstrap tree is identical to the true tree and the bootstrap supports for all branches are high reflecting the consistency [37] in the branching. The RCM tree was produced by the alignment-free RCM approach. The RCM tree is identical to both the true tree and the bootstrap tree reflecting its potential for use in phylogenetic analysis of protein sequences. The tree topology of only one of the trees is shown in Figure 2 since they are all identical.

We applied the RCM approach to seven protein datasets. RCM method showed an efficient division of protein families into subfamilies using RAAAs. Phylogenetic trees of the seven protein families using RCM approach are shown in Figure 3 for ML15 alphabet. Detailed comparison of RCM with MSA in terms of TBC accuracy, the number and percentage of TBC error for each RAAA and each dataset is provided in Additional File 3.

The comparison of RCM with MSA in terms of TBC accuracy and the percentage of TBC error are summarized in Table 4 for the 20-letter alphabet and the top performing RAAA with the minimum size. In cases where two RAAAs of the same size give identical TBC results, both of them are reported. Three trends can be observed from the data in Table 4.

First, for five of the seven families (crotonases, mandelate racemases, nucleotidyl cyclases, acyl transferases, and GH2 hydrolases), both methods perform equally well comparably. For VOC, RCM outperforms MSA while for haloacid dehalogenases, MSA slightly outperforms RCM. It is important to note that both VOCs and dehalogenases have the two lowest mean PIDs (12% vs. 14%) and low mean sequence lengths with large standard deviation. Low PID and low sequence length are two features in alignments that render inference of relationship based only on sequence information difficult. Nonetheless, TBC accuracies of both families with their respective top performing RAAAs are comparable to the results obtained from the protein families with higher mean PIDs and longer mean sequence lengths.

Second, either ML4 or GBMR4 is sufficient to obtain high TBC accuracy for all datasets except VOCs and haloacid dehalogenases. Indeed, apart from the aforementioned families, ML4 and GBMR4 can produce either identical or better results than all other alphabets using either RCM or MSA, implying that as little as an alphabet size of 4 would be sufficient to capture most of the sequence information that might yield considerable improvements in inferring relationship based on sequence information when both mean PID and the length of the aligned regions in an MSA is above a certain threshold.

Third, for the datasets with low mean PIDs and average sequence lengths, a larger RAAA size may be required to obtain identical or better results than the 20-letter alphabet using both RCM and MSA. This is especially evident with the RCM approach. While the minimum RAAA size of the top performer was 4 for 5 datasets that have relatively higher average sequence lengths and mean PIDs, it increases to 8 (EB8) for VOCs and 15 (ML15) for haloacid dehalogenases that have mean PIDs of 14% and 12%, respectively. Moreover, a subtle but a similar trend is also evident in the case of MSA. While the alphabet size of the top performer was 4 (GBMR4, ML4) for VOCs, it increased to 8 (ML8) for haloacid dehalogenases, implying that a larger RAAA size may perform better on sequences with lower sequence identities.

It is also interesting to note that the average TBC error for mandelate racemases, nucleotidyl cyclases and hydrolases with three random alphabets of size 4 varied between 0% and 15.6% for the MSA method. While the groupings of amino acids in the random alphabets do not have any physicochemical or structural significance that can justify this overall performance, the low percent TBC error may suggest that some subfamilies of these protein families may be very tight with small distances between their sequences while larger distance between different subfamilies. This scenario coupled with the relatively longer sequences (top three families in terms of mean sequence length) within these families may generate sufficiently long aligned regions with enough informative sites that can result in a tree that correctly assigns subfamilies even the reduced alphabet groupings do not have any structural or biological meaning.

However, the trend of low TBC error is not apparent using RCM with random alphabets. TBC errors of different protein families using random RAAAs (average of three random alphabets) were significantly higher than TBC errors using biologically meaningful reduced alphabets for all the families except racemases and nucleotidyl cyclases, both of which overlap with the results obtained with MSA.

Performance of RCM approach with different RAAAs to cluster protein families into functional subfamilies is eminent. Yet, it must be noted that there is no uniformly superior algorithm for tree-based subfamily clustering and that simple protein similarity measures combined with hierarchical clustering produce trees with reasonable and often high accuracy [4]. Furthermore, if much time has passed since the evolution of different subfamilies, then sequences may have diverged beyond the point where simple phylogenetic analysis cannot easily give a clear distinction of subfamilies.