BACKGROUND
proteins that evolve from a common ancestor can change functionality over time  <cit>  and produce highly divergent protein families that can be divided into subfamilies with similar but distinct functions   <cit> .
identification of subfamilies using protein sequence information can be carried out using phylogenetic methods that can reveal the evolutionary relationship between proteins by clustering similar proteins together in a phylogenetic tree  <cit> .
the most common method for identifying similarities in sequences through phylogenetic analysis starts with the construction of a multiple alignment of homologous sequences using a substitution matrix.
multiple alignment scores are then transformed into a distance matrix to construct a phylogenetic tree.
often the branching order of a phylogenetic tree exactly matches the known functional split between proteins  <cit>  and branch lengths are proportional to the extent of evolutionary changes since the last common ancestor.
multiple sequence alignment  is constructed using a scoring scheme which reward or penalize each substitution, insertion and deletion to get an optimum alignment of the given sequences.
the quality of an msa is connected to the chosen parameters that are entered manually and an expert handling is almost always required to maintain alignment integrity by observing general trends in each protein family.
as such different alignment parameters may yield different phylogenetic trees that are only as good as the msa that the trees are derived from  <cit> .
phylogenetic analysis is broadly divided into two groups of methods.
algorithms in the first group calculate a matrix representing the distance between each pair of sequences and then transform this matrix into a tree using a tree-clustering algorithm.
algorithms in the first category utilize various distance measures with different models to account for nucleotide or amino acid substitutions.
in the second group, the tree that can best explain the observed sequences under the chosen evolutionary model is found by evaluating the fitness of different tree topologies  <cit> .
the second category can further be divided into two groups based on the optimality criterion used in tree evaluation: maximum parsimony and maximum likelihood.
under maximum parsimony  <cit> , the preferred phylogenetic tree is the tree that requires the least evolutionary change to explain the observed data whereas under maximum likelihood  <cit> , it is the most probable tree under the chosen evolutionary assumption.
the prediction of subfamilies from protein msas have been carried out previously by comparing subfamily hidden markov models, subfamily specific sequence profiles, analyzing positional entropies in an alignment, and ascending hierarchical method  <cit> .
all of these methods require an alignment of biological sequences that assume some sort of an evolutionary model.
computational complexity and the inherent ambiguity of the alignment cost criteria are two major problems in msa along with controversial evolutionary models that are used to explain them.
a novel approach for phylogenetic analysis based on relative complexity measure  of whole genomic sequences have been previously proposed by otu et al, that eliminates the need for msa and produces successful phylogenies on real and simulated datasets  <cit> .
the algorithm employs lempel-ziv  complexity  <cit>  and produces a score for each sequence pair that can be interpreted as the "closeness" of the sequence pairs.
unequal sequence length or different positioning of similar regions along sequences  is not an issue as the method has been shown to handle both cases naturally.
moreover, rcm does not use any approximations and assumptions in calculating the distance between sequences.
therefore, rcm utilizes the information contained in sequences and requires no human intervention.
application of rcm to genomic sequences for phylogenetic analysis was successfully carried out on various datasets containing genomic sequences  <cit> .
moreover, liu et al  <cit>  extended this method further to integrate the hydropathy profile and a different lz-based distance measure for phylogenetic analysis of protein sequences while russell et al integrated a merged amino acid alphabet containing  <dig> characters to represent all amino acids to reduce complexity prior to calculating a pairwise distance measure to be used as a pairwise scoring function in determining the order with which sequences should be joined in a multiple sequence alignment problem  <cit> .
application of rcm to evaluate genomic sequences is relatively straight forward since rcm based on lempel-ziv complexity scores can capture each mutation in dna sequences and register it as an increase in the complexity scores of compared sequences.
however, substitution of one residue into another in proteins is tolerable as long as the substituted residue is not highly conserved and physicochemical and structural properties of the substituted and the native residues are not fundamentally different  <cit> .
employment of hydropathy-index-based grouping of residues is one way of a preprocessing requirement to capture only the mutations that would not be tolerated in a protein sequence since lz algorithm is not capable of accounting for amino acid substitution frequencies and similarity scores.
hence, any application that uses rcm to generate a distance matrix of protein sequences should be linked to treating the sequence with a reduced amino acid alphabet  prior to calculating their rcms.
in this paper, we utilize rcm with different reduced amino acid alphabets and assess rcm's potential in clustering protein families into functional subtypes based solely on sequence data.
this method clustered seven well-characterized protein families into their functional subtypes with 92% - 100% accuracy.
methods
datasets
simulated dataset
performance of rcm was tested on a simulated dataset that contains  <dig> randomly evolved protein sequences from a root sequence of length  <dig> by using indelible v <dig>   <cit> .
simulated protein sequences were generated according to the following parameters:
<dig>  jtt-dcmut  <cit>  was chosen as the amino acid substitution model.
<dig>  power law insertion/deletion length distribution model with a =  <dig>  and maximum allowed insertion/deletion length of  <dig> were used.
<dig>  both insertion and deletion rates were set to the default parameter of  <dig>  relative to average substitution rate of 1%.
<dig>  length of the root protein sequence was set to  <dig>
<dig>  the rooted tree with  <dig> taxa that reflects the true phylogenetic evolution of the sequences was generated along with the true msa from which the true tree was inferred.
<dig>  the true msa was then inputted into clustalw <dig>  <cit>  and the bootstrap tree was generated
protein datasets
rcm was tested on seven protein datasets.
number of sequences, number of subfamilies, average length, standard deviation of sequence lengths and mean percent identities   <cit>  of sequences for each family are summarized in table  <dig>  protein sequences for mandelate racemases, crotonases, haloacid dehalogenases and vicinal oxygen chelates  were extracted from extensively curated structure-function linkage database which contains sets of subfamily grouping for a large set of protein families.
sfld contains protein families with a hierarchical classification scheme based on sequence, structure and conserved chemical reactions at the superfamily, subgroup, and family levels  <cit> .
crotonases and haloacid dehalogenases were filtered such that subfamilies that contain less than  <dig> sequences or more than  <dig> sequences were removed to prevent sequence number bias and to reduce computational complexity.
unknown or unspecified amino acids were discarded .
the protein sequences for acyl transferase  domains and nucleotidyl cyclases were obtained from reference  <cit> .
the protein sequences in the hard-to-align dataset that contains glycoside hydrolase family  <dig>  members were adapted from reference  <cit> .
expert curated annotations of protein sequences and abbreviations used for sequences in this study are provided in additional file  <dig>
* mean percent identity  is the average of all pairwise sequence identities in a given family.
reduced amino acid alphabets
sequence space of proteins is redundant and generates only a limited number of folds, domains, and structures  <cit> .
various strategies have been devised that take a coarse-grained approach to account for the degeneracy of sequences by grouping similar amino acids together  <cit> .
grouping is usually carried out based on structural and physiochemical similarities of amino acids  <cit> .
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  <cit> .
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  <cit> .
grouping of amino acid residues has also been used extensively in hydrophobic-polar  lattice model to explain the hydrophobic collapse theory of protein folding  <cit> .
a recent study was carried out by peterson et al to test the performance of over  <dig> raaas on the sequence library from dalipdb <dig> database and showed that raaas improves sensitivity and specificity in fold prediction between protein sequence pairs with high structural similarity and low sequence identity  <cit> .
we tested performances of six amino acid reduction schemes with  <dig> different level of groupings to separate proteins into functional subfamilies .
these included three top performing raaa  from reference  <cit>  and three random raaa of size  <dig>
reduced amino acid schemes used in this study.
* substitution matrices for these reduced alphabets were obtained from reference  <cit> .
§ bl <dig> frequency counts were used to derive these substitution matrices using the formula outlined in reference  <cit> .
#gap opening/gap extension penalties used for msas in clustalw <dig>
substitution matrices
amino acids that are within the same group in a raaa are considered identical  <cit> .
substitution matrices that assign the same similarity score to each amino acid within the same group were obtained from reference  <cit> .
for those raaas in the eb scheme and the three random raaas, new substitution matrices were created from blosum <dig> frequency counts using the same procedure outlined in reference  <cit> .
lempel-ziv complexity
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  <cit>  complexity and was known to accurately cluster all related genomic sequences under one branch of the tree  <cit> .
lempel-ziv  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  <cit> .
for example, the lz-complexity of the simple sequence 'aailnaiiannl' would be obtained as shown in table  <dig>  since seven steps are needed to generate the whole sequence, the lz-complexity score for this sequence is  <dig>  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  <cit> , however, this work used only the following normalized distance metric that accounts for the differences in sequence lengths:
the exhaustive library construction and lempel-ziv complexity score calculation of sequence x.
dxy=c+c−c−cc+c <dig>
where c and c 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 .
although in performance between five measures were not significant, we adopted the aforementioned distance for its ability to account for length variance.
distance matrix & phylogenetic tree
the relative complexity measure  for creation of the distance matrix was utilized as previously described  <cit> .
phylogenetic trees were generated from distance matrices using neighbor-joining  <cit>  program of the phylogeny inference package, phylip  <dig>   <cit> .
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  using the retree program of phylip package  <cit> .
phylogenetic trees for all protein families and raaas are shown in supplementary materials  in newick format and can be visualized with a tree-drawing program.
clustalw2
protein sequences in each family were aligned using clustalw <dig>  <cit>  for comparison with rcm.
msas were performed using updated substitution matrices with gap extension and gap opening penalties provided in table  <dig>  bootstrap analyses were carried out  <dig> times and trees containing bootstrap values were created using clustalw <dig> with the neighbor-joining clustering algorithm.
for convenience, msas that were carried out using clustalw <dig> will be referred as the msa or the msa method for the rest of the article.
tree based classification
tbc algorithm  <cit>  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/ and fn/ are both equal to  <dig>  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  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  <cit> .
protocol
the proposed algorithm operates on a set of sequences in fasta format.
after one of the alphabets given in table  <dig> 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 clustalw <dig> and bootstrapped  <dig> 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  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  <dig>  msa-based trees were created following bootstrap analysis  with clustalw <dig>
finally, for each family, a total of  <dig> phylogenetic trees  for each method are generated and checked how well they separated families into subfamilies.
a summary of the overall workflow is depicted in figure  <dig>
RESULTS
simulated dataset
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  <dig> protein sequences was generated to show that rcm coupled with a raaa can produce a phylogenetic 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  <dig> simulated sequences.
on the other hand, the bootstrap tree is the tree that was produced by clustalw <dig> 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  <cit>  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  <dig> since they are all identical.
performance of the rcm approach
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  <dig> for ml <dig> 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  <dig>
crotonases
members of crotonase family contain  <dig> protein sequences from  <dig> different subfamilies and catalyze diverse metabolic reactions with certain family members displaying dehalogenase, hydratase, and isomerase activities.
tbc accuracy varied between  <dig> % and 100% for rcm.
the top performing raaa with the smallest size was gbmr <dig> that resulted in 100% tbc accuracy.
tbc accuracy was 100% for all raaas tested with msa.
mandelate racemases
the mandelate racemase dataset contains  <dig> sequences that are assigned to  <dig> expert curated subfamilies.
all mandelate racemases contain a conserved histidine, presumably acting as an active site base  <cit> .
when the rcm approach was tested on mandelate racemases, all resulting trees showed correct assignment of functional subfamilies into  <dig> different clusters with 100% accuracy using all alphabets except gbmr <dig> that resulted in  <dig> % tbc accuracy.
vicinal oxygen chelates
voc family contains  <dig> sequences from  <dig> different subfamilies.
the number of tbc accuracy varied between  <dig> % and 92% for rcm and  <dig> % to  <dig> % for msa.
members of voc have an average sequence length of  <dig> amino acids and a mean pid of 14% .
the low pid and the highly divergent nature of this family make its subfamilies susceptible to misclassification more than other families based on sequence information alone.
in this dataset, eb <dig> performed better than 20-letter alphabet  with rcm while gbmr <dig>  ml <dig>  eb <dig>  eb, eb <dig> and 20-letter alphabets resulted in  <dig> % tbc errors with msa.
haloacid dehalogenases
haloacid dehalogenases contains  <dig> sequences that belong to  <dig> different subfamilies.
haloacid dehalogenase family is similar to vocs in its highly divergent nature based on the low mean pid  that places the sequences in this family in the "twilight zone" to infer any relation between sequences based on sequence information alone.
ml <dig> was the best performing raaa for rcm with  <dig> % accuracy .
the size of the best performing raaa for this family is larger compared to other families hinting that highly divergent sequences may require larger alphabets with lower level of grouping.
statistics for top performing
tbc accuracy and percentage of tbc error are reported for the 20-letter alphabet and the top performing raaa.
if two raaas with the same size have identical tbc accuracies, both raaas are reported at the final row in the table.
bold entries correspond to top performers using rcm and msa for the specified datasets
nucleotidyl cyclases
nucleotidyl cyclase family has two functional subfamilies, adenylate and guanylate cyclases that correspond to use of the substrates atp and gtp respectively.
the nucleotidyl cyclase family with  <dig> adenylate cyclases and  <dig> guanylate cyclases was clustered into two distinct subfamilies with 100% accuracy using both methods and all raaas except eb <dig> and eb <dig> for rcm and ml <dig> and eb <dig> for msa, all of which resulted in  <dig> % accuracy .
moreover, the clustering result for the nucleotidyl cyclases are in agreement with the result obtained previously by the msa-dependent clustering algorithm that uses the residues with the highest evolutionary split statistic to split protein families into functional subfamilies  <cit> .
acyl transferases
the at domains of type i modular polyketide synthases are responsible for the substrate selection.
most incorporate either a c <dig> unit  or a c <dig> unit .
the choice of substrate can be deduced from the chemical structure of the polyketide product  <cit> .
in the acyl transferase dataset,  <dig> of the  <dig> sequences use c <dig> units whereas  <dig> use c <dig> units as substrate.
previously, goldstein et al  <cit>  used evolutionary split statistic and clustered the at domains into  <dig> subfamilies with  <dig> false assignments for the  <dig> residue-long motif.
the number of false assignments increased to  <dig> with increasing motif length  suggesting that the utilization of a larger motif increases the noise and error rate.
as such, inclusion of only  <dig> residues  with high split statistics increases the assignment accuracy .
a similar trend is observed in the case of rcm.
while the tbc accuracy for at domains was only 91%  with the 20-letter alphabet , the accuracy increased to 97%  with the utilization of the ml <dig>  ml <dig>  eb <dig>  ml <dig>  eb <dig>  sdm <dig>  eb <dig>  and hsdm <dig> alphabets.
furthermore,  <dig> of the  <dig> misclassified sequences using the above reduced alphabets are contained in the  <dig>   <dig> and  <dig> false assignments produced by the goldstein et al 's approach using the  <dig>  and  <dig> residue-long motifs, respectively.
although the accuracy was higher previously, it should be noted that the rcm approach did neither require an msa of sequences nor any other sequence-based statistics.
the accuracy was  <dig> % for msa using the top performing raaas.
there was no immediate evidence suggesting a specific characteristic for incorrectly classified sequences.
glycoside hydrolase family  <dig>
the final dataset contains  <dig> members of the gh <dig> family with a  <dig> fold.
the subfamilies and the number of sequences from each subfamily are β-galactosidases , β-mannosidases , β-glucuronidases  and exo-β-d-glucosaminidases .
this dataset was used previously and chosen because it was cited as a "hard-to-align" dataset by classical alignment approaches  <cit> .
the gh <dig> family was clustered into  <dig> functional subfamilies with 100% accuracy using ml <dig> and gbmr <dig> - the two top performing raaas - with rcm .
tbc accuracy was 100% for all raaas tested with msa.
the effect of the size of the raaa on clustering performance
the comparison of rcm with msa in terms of tbc accuracy and the percentage of tbc error are summarized in table  <dig> 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  <dig>
first, for five of the seven families , 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  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 ml <dig> or gbmr <dig> is sufficient to obtain high tbc accuracy for all datasets except vocs and haloacid dehalogenases.
indeed, apart from the aforementioned families, ml <dig> and gbmr <dig> 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  <dig> 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  <dig> for  <dig> datasets that have relatively higher average sequence lengths and mean pids, it increases to  <dig>  for vocs and  <dig>  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  <dig>  for vocs, it increased to  <dig>  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  <dig> varied between 0% and  <dig> % 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  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  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  <cit> .
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.
CONCLUSIONS
the application of rcm in generating meaningful phylogenetic trees has been previously tested on genomic sequences and made rcm a good alternative to msa-based phylogenetic analysis.
however, integration of rcm to measure the closeness of protein sequences was simply problematic due to the lack and difficulty of accounting for amino acid substitutions.
in this paper, we introduced an raaa-based approach as a preprocessing of protein sequences prior to calculating pairwise rcms.
utilization of an raaa that is consistent with the structure and function of the proteins or an raaa that reflects the general trends in specific protein families under study can result in successful phylogenies that can cluster each protein superfamily into functional subfamilies.
in finding functional subtypes of a protein family, it is often of interest to find out if the mechanisms that manipulate a certain clustering are of evolutionary or functional origin.
although these two signals may be overlapping and hard to separate, rcm could be used to address this issue by finding differences in exhaustive histories in two sequences when they are concatenated.
the "words" that result in an observed difference can then be analyzed and correlated to a functional and/or evolutionary origin.
we believe future work can focus in this direction building on the current approach that does not attempt to trace back the origin of differentiating sequence signals but provides a powerful clustering method of protein families into functional subtypes without using multiple sequence alignment.
authors' contributions
uos and hho participated in the design of the study and supervised all the experiments.
aa performed all the experiments and wrote the initial manuscript and the final manuscript.
hho provided the lz algorithm, revised the first and the final manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
fasta header mapping.
this file contains the fasta header abbreviations for protein families and expert assignment of sequences to each subfamily.
some bioinformatics programs that take fasta files as input have fasta header size limitations ranging from  <dig> to  <dig> characters long.
click here for file
additional file 2
phylogenetic trees for all the datasets in newick.
phylogenetic tree files for all families are presented in newick format.
for simulated dataset, there are  <dig> phylogenetic trees.
for each protein dataset, there are  <dig> phylogenetic trees:  <dig> rcm trees and  <dig> msa trees.
all trees reflect only the tree topology .
click here for file
additional file 3
tbc errors for all families and all raaas.
detailed comparison of rcm and msa is reported in terms of the number and percentage of tbc error for every protein family and raaa under consideration.
click here for file
acknowledgements
the authors would like to thank cem meydan and ozgur gul for helpful discussions, eric peterson for supplying the perl script for the generation of substitution matrices.
hho is partially supported by a grant from the dubai harvard foundation for medical research.

 proteins that evolve from a common ancestor can change functionality over time and produce highly divergent protein families that can be divided into subfamilies with similar but distinct functions. identification of subfamilies using protein sequence information can be carried out using phylogenetic methods that can reveal the evolutionary relationship between proteins by clustering similar proteins together in a phylogenetic tree. the most common method for identifying similarities in sequences through phylogenetic analysis starts with the construction of a multiple alignment of homologous sequences using a substitution matrix. multiple alignment scores are then transformed into a distance matrix to construct a phylogenetic tree. often the branching order of a phylogenetic tree exactly matches the known functional split between proteins and branch lengths are proportional to the extent of evolutionary changes since the last common ancestor. multiple sequence alignment is constructed using a scoring scheme which reward or penalize each substitution, insertion and deletion to get an optimum alignment of the given sequences. the quality of an msa is connected to the chosen parameters that are entered manually and an expert handling is almost always required to maintain alignment integrity by observing general trends in each protein family. as such different alignment parameters may yield different phylogenetic trees that are only as good as the msa that the trees are derived from. phylogenetic analysis is broadly divided into two groups of methods. algorithms in the first group calculate a matrix representing the distance between each pair of sequences and then transform this matrix into a tree using a tree-clustering algorithm. algorithms in the first category utilize various distance measures with different models to account for nucleotide or amino acid substitutions. in the second group, the tree that can best explain the observed sequences under the chosen evolutionary model is found by evaluating the fitness of different tree topologies. the second category can further be divided into two groups based on the optimality criterion used in tree evaluation: maximum parsimony and maximum likelihood. under maximum parsimony, the preferred phylogenetic tree is the tree that requires the least evolutionary change to explain the observed data whereas under maximum likelihood, it is the most probable tree under the chosen evolutionary assumption. the prediction of subfamilies from protein msas have been carried out previously by comparing subfamily hidden markov models, subfamily specific sequence profiles, analyzing positional entropies in an alignment, and ascending hierarchical method. all of these methods require an alignment of biological sequences that assume some sort of an evolutionary model. computational complexity and the inherent ambiguity of the alignment cost criteria are two major problems in msa along with controversial evolutionary models that are used to explain them. a novel approach for phylogenetic analysis based on relative complexity measure of whole genomic sequences have been previously proposed by otu et al, that eliminates the need for msa and produces successful phylogenies on real and simulated datasets. the algorithm employs lempel-ziv complexity and produces a score for each sequence pair that can be interpreted as the "closeness" of the sequence pairs. unequal sequence length or different positioning of similar regions along sequences is not an issue as the method has been shown to handle both cases naturally. moreover, rcm does not use any approximations and assumptions in calculating the distance between sequences. therefore, rcm utilizes the information contained in sequences and requires no human intervention. application of rcm to genomic sequences for phylogenetic analysis was successfully carried out on various datasets containing genomic sequences. moreover, liu et al extended this method further to integrate the hydropathy profile and a different lz-based distance measure for phylogenetic analysis of protein sequences while russell et al integrated a merged amino acid alphabet containing characters to represent all amino acids to reduce complexity prior to calculating a pairwise distance measure to be used as a pairwise scoring function in determining the order with which sequences should be joined in a multiple sequence alignment problem. application of rcm to evaluate genomic sequences is relatively straight forward since rcm based on lempel-ziv complexity scores can capture each mutation in dna sequences and register it as an increase in the complexity scores of compared sequences. however, substitution of one residue into another in proteins is tolerable as long as the substituted residue is not highly conserved and physicochemical and structural properties of the substituted and the native residues are not fundamentally different. employment of hydropathy-index-based grouping of residues is one way of a preprocessing requirement to capture only the mutations that would not be tolerated in a protein sequence since lz algorithm is not capable of accounting for amino acid substitution frequencies and similarity scores. hence, any application that uses rcm to generate a distance matrix of protein sequences should be linked to treating the sequence with a reduced amino acid alphabet prior to calculating their rcms. in this paper, we utilize rcm with different reduced amino acid alphabets and assess rcm's potential in clustering protein families into functional subtypes based solely on sequence data. this method clustered seven well-characterized protein families into their functional subtypes with 92% - 100% accuracy. methods datasets simulated dataset performance of rcm was tested on a simulated dataset that contains randomly evolved protein sequences from a root sequence of length by using indelible v. simulated protein sequences were generated according to the following parameters: jtt-dcmut was chosen as the amino acid substitution model. power law insertion/deletion length distribution model with a = and maximum allowed insertion/deletion length of were used. both insertion and deletion rates were set to the default parameter of relative to average substitution rate of 1%. length of the root protein sequence was set to the rooted tree with taxa that reflects the true phylogenetic evolution of the sequences was generated along with the true msa from which the true tree was inferred. the true msa was then inputted into clustalw and the bootstrap tree was generated protein datasets rcm was tested on seven protein datasets. number of sequences, number of subfamilies, average length, standard deviation of sequence lengths and mean percent identities of sequences for each family are summarized in table protein sequences for mandelate racemases, crotonases, haloacid dehalogenases and vicinal oxygen chelates were extracted from extensively curated structure-function linkage database which contains sets of subfamily grouping for a large set of protein families. sfld contains protein families with a hierarchical classification scheme based on sequence, structure and conserved chemical reactions at the superfamily, subgroup, and family levels. crotonases and haloacid dehalogenases were filtered such that subfamilies that contain less than sequences or more than sequences were removed to prevent sequence number bias and to reduce computational complexity. unknown or unspecified amino acids were discarded. the protein sequences for acyl transferase domains and nucleotidyl cyclases were obtained from reference. the protein sequences in the hard-to-align dataset that contains glycoside hydrolase family members were adapted from reference. expert curated annotations of protein sequences and abbreviations used for sequences in this study are provided in additional file * mean percent identity is the average of all pairwise sequence identities in a given family. reduced amino acid alphabets sequence space of proteins is redundant and generates only a limited number of folds, domains, and structures. various strategies have been devised that take a coarse-grained approach to account for the degeneracy of sequences by grouping similar amino acids together. grouping is usually carried out based on structural and physiochemical similarities of amino acids. 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. 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. grouping of amino acid residues has also been used extensively in hydrophobic-polar lattice model to explain the hydrophobic collapse theory of protein folding. a recent study was carried out by peterson et al to test the performance of over raaas on the sequence library from dalipdb database and showed that raaas improves sensitivity and specificity in fold prediction between protein sequence pairs with high structural similarity and low sequence identity. we tested performances of six amino acid reduction schemes with different level of groupings to separate proteins into functional subfamilies. these included three top performing raaa from reference and three random raaa of size reduced amino acid schemes used in this study.* substitution matrices for these reduced alphabets were obtained from reference. § bl frequency counts were used to derive these substitution matrices using the formula outlined in reference. #gap opening/gap extension penalties used for msas in clustalw substitution matrices amino acids that are within the same group in a raaa are considered identical. substitution matrices that assign the same similarity score to each amino acid within the same group were obtained from reference. for those raaas in the eb scheme and the three random raaas, new substitution matrices were created from blosum frequency counts using the same procedure outlined in reference. lempel-ziv complexity 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 complexity and was known to accurately cluster all related genomic sequences under one branch of the tree. lempel-ziv 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. for example, the lz-complexity of the simple sequence'aailnaiiannl' would be obtained as shown in table since seven steps are needed to generate the whole sequence, the lz-complexity score for this sequence is 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, however, this work used only the following normalized distance metric that accounts for the differences in sequence lengths: the exhaustive library construction and lempel-ziv complexity score calculation of sequence x. dxy=c+c−c−cc+c where c and c 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. although in performance between five measures were not significant, we adopted the aforementioned distance for its ability to account for length variance. distance matrix & phylogenetic tree the relative complexity measure for creation of the distance matrix was utilized as previously described. phylogenetic trees were generated from distance matrices using neighbor-joining program of the phylogeny inference package, phylip. 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 using the retree program of phylip package. phylogenetic trees for all protein families and raaas are shown in supplementary materials in newick format and can be visualized with a tree-drawing program. clustalw2 protein sequences in each family were aligned using clustalw for comparison with rcm. msas were performed using updated substitution matrices with gap extension and gap opening penalties provided in table bootstrap analyses were carried out times and trees containing bootstrap values were created using clustalw with the neighbor-joining clustering algorithm. for convenience, msas that were carried out using clustalw will be referred as the msa or the msa method for the rest of the article. tree based classification tbc algorithm 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/ and fn/ are both equal to 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 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. protocol the proposed algorithm operates on a set of sequences in fasta format. after one of the alphabets given in table 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 clustalw and bootstrapped 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 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 msa-based trees were created following bootstrap analysis with clustalw finally, for each family, a total of phylogenetic trees for each method are generated and checked how well they separated families into subfamilies. a summary of the overall workflow is depicted in figure  simulated dataset 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 protein sequences was generated to show that rcm coupled with a raaa can produce a phylogenetic 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 simulated sequences. on the other hand, the bootstrap tree is the tree that was produced by clustalw 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 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 since they are all identical. performance of the rcm approach 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 for ml 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 crotonases members of crotonase family contain protein sequences from different subfamilies and catalyze diverse metabolic reactions with certain family members displaying dehalogenase, hydratase, and isomerase activities. tbc accuracy varied between % and 100% for rcm. the top performing raaa with the smallest size was gbmr that resulted in 100% tbc accuracy. tbc accuracy was 100% for all raaas tested with msa. mandelate racemases the mandelate racemase dataset contains sequences that are assigned to expert curated subfamilies. all mandelate racemases contain a conserved histidine, presumably acting as an active site base. when the rcm approach was tested on mandelate racemases, all resulting trees showed correct assignment of functional subfamilies into different clusters with 100% accuracy using all alphabets except gbmr that resulted in % tbc accuracy. vicinal oxygen chelates voc family contains sequences from different subfamilies. the number of tbc accuracy varied between % and 92% for rcm and % to % for msa. members of voc have an average sequence length of amino acids and a mean pid of 14%. the low pid and the highly divergent nature of this family make its subfamilies susceptible to misclassification more than other families based on sequence information alone. in this dataset, eb performed better than 20-letter alphabet with rcm while gbmr ml eb eb, eb and 20-letter alphabets resulted in % tbc errors with msa. haloacid dehalogenases haloacid dehalogenases contains sequences that belong to different subfamilies. haloacid dehalogenase family is similar to vocs in its highly divergent nature based on the low mean pid that places the sequences in this family in the "twilight zone" to infer any relation between sequences based on sequence information alone. ml was the best performing raaa for rcm with % accuracy. the size of the best performing raaa for this family is larger compared to other families hinting that highly divergent sequences may require larger alphabets with lower level of grouping. statistics for top performing tbc accuracy and percentage of tbc error are reported for the 20-letter alphabet and the top performing raaa. if two raaas with the same size have identical tbc accuracies, both raaas are reported at the final row in the table. bold entries correspond to top performers using rcm and msa for the specified datasets nucleotidyl cyclases nucleotidyl cyclase family has two functional subfamilies, adenylate and guanylate cyclases that correspond to use of the substrates atp and gtp respectively. the nucleotidyl cyclase family with adenylate cyclases and guanylate cyclases was clustered into two distinct subfamilies with 100% accuracy using both methods and all raaas except eb and eb for rcm and ml and eb for msa, all of which resulted in % accuracy. moreover, the clustering result for the nucleotidyl cyclases are in agreement with the result obtained previously by the msa-dependent clustering algorithm that uses the residues with the highest evolutionary split statistic to split protein families into functional subfamilies. acyl transferases the at domains of type i modular polyketide synthases are responsible for the substrate selection. most incorporate either a c unit or a c unit. the choice of substrate can be deduced from the chemical structure of the polyketide product. in the acyl transferase dataset, of the sequences use c units whereas use c units as substrate. previously, goldstein et al used evolutionary split statistic and clustered the at domains into subfamilies with false assignments for the residue-long motif. the number of false assignments increased to with increasing motif length suggesting that the utilization of a larger motif increases the noise and error rate. as such, inclusion of only residues with high split statistics increases the assignment accuracy. a similar trend is observed in the case of rcm. while the tbc accuracy for at domains was only 91% with the 20-letter alphabet, the accuracy increased to 97% with the utilization of the ml ml eb ml eb sdm eb and hsdm alphabets. furthermore, of the misclassified sequences using the above reduced alphabets are contained in the and false assignments produced by the goldstein et al's approach using the and residue-long motifs, respectively. although the accuracy was higher previously, it should be noted that the rcm approach did neither require an msa of sequences nor any other sequence-based statistics. the accuracy was % for msa using the top performing raaas. there was no immediate evidence suggesting a specific characteristic for incorrectly classified sequences. glycoside hydrolase family the final dataset contains members of the gh family with a fold. the subfamilies and the number of sequences from each subfamily are β-galactosidases, β-mannosidases, β-glucuronidases and exo-β-d-glucosaminidases. this dataset was used previously and chosen because it was cited as a "hard-to-align" dataset by classical alignment approaches. the gh family was clustered into functional subfamilies with 100% accuracy using ml and gbmr - the two top performing raaas - with rcm. tbc accuracy was 100% for all raaas tested with msa. the effect of the size of the raaa on clustering performance the comparison of rcm with msa in terms of tbc accuracy and the percentage of tbc error are summarized in table 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 first, for five of the seven families, 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 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 ml or gbmr is sufficient to obtain high tbc accuracy for all datasets except vocs and haloacid dehalogenases. indeed, apart from the aforementioned families, ml and gbmr 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 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 for datasets that have relatively higher average sequence lengths and mean pids, it increases to for vocs and 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 for vocs, it increased to 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 varied between 0% and % 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 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 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. 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.  the application of rcm in generating meaningful phylogenetic trees has been previously tested on genomic sequences and made rcm a good alternative to msa-based phylogenetic analysis. however, integration of rcm to measure the closeness of protein sequences was simply problematic due to the lack and difficulty of accounting for amino acid substitutions. in this paper, we introduced an raaa-based approach as a preprocessing of protein sequences prior to calculating pairwise rcms. utilization of an raaa that is consistent with the structure and function of the proteins or an raaa that reflects the general trends in specific protein families under study can result in successful phylogenies that can cluster each protein superfamily into functional subfamilies. in finding functional subtypes of a protein family, it is often of interest to find out if the mechanisms that manipulate a certain clustering are of evolutionary or functional origin. although these two signals may be overlapping and hard to separate, rcm could be used to address this issue by finding differences in exhaustive histories in two sequences when they are concatenated. the "words" that result in an observed difference can then be analyzed and correlated to a functional and/or evolutionary origin. we believe future work can focus in this direction building on the current approach that does not attempt to trace back the origin of differentiating sequence signals but provides a powerful clustering method of protein families into functional subtypes without using multiple sequence alignment. authors' contributions uos and hho participated in the design of the study and supervised all the experiments. aa performed all the experiments and wrote the initial manuscript and the final manuscript. hho provided the lz algorithm, revised the first and the final manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 fasta header mapping. this file contains the fasta header abbreviations for protein families and expert assignment of sequences to each subfamily. some bioinformatics programs that take fasta files as input have fasta header size limitations ranging from to characters long. click here for file additional file 2 phylogenetic trees for all the datasets in newick. phylogenetic tree files for all families are presented in newick format. for simulated dataset, there are phylogenetic trees. for each protein dataset, there are phylogenetic trees: rcm trees and msa trees. all trees reflect only the tree topology. click here for file additional file 3 tbc errors for all families and all raaas. detailed comparison of rcm and msa is reported in terms of the number and percentage of tbc error for every protein family and raaa under consideration. click here for file acknowledgements the authors would like to thank cem meydan and ozgur gul for helpful discussions, eric peterson for supplying the perl script for the generation of substitution matrices. hho is partially supported by a grant from the dubai harvard foundation for medical research.
we demonstrate the robustness of rcm with reduced alphabets in clustering of protein sequences into families in a simulated dataset and seven well-characterized protein datasets. on protein datasets, crotonases, mandelate racemases, nucleotidyl cyclases and glycoside hydrolase family were clustered into subfamilies with 100% accuracy whereas acyl transferase domains, haloacid dehalogenases, and vicinal oxygen chelates could be assigned to subfamilies with %, % and % accuracies, respectively. to address this problem, we propose to use the alignment-free relative complexity measure combined with reduced amino acid alphabets to cluster protein families into functional subtypes purely on sequence criteria. most phylogenetic analysis methods utilize multiple alignment of sequences and are based on an evolutionary model. however, multiple alignment is not an automated procedure and requires human intervention to maintain alignment integrity and to produce phylogenies consistent with the functional splits in underlying sequences. the method is also flexible and computationally fast because it does not require multiple alignment of sequences. the overall combination of methods in this paper is useful for clustering protein families into subtypes based on solely protein sequence information. phylogenetic analysis can be used to divide a protein family into subfamilies in the absence of experimental information.
 phylogenetic analysis can be used to divide a protein family into subfamilies in the absence of experimental information. most phylogenetic analysis methods utilize multiple alignment of sequences and are based on an evolutionary model. however, multiple alignment is not an automated procedure and requires human intervention to maintain alignment integrity and to produce phylogenies consistent with the functional splits in underlying sequences. to address this problem, we propose to use the alignment-free relative complexity measure combined with reduced amino acid alphabets to cluster protein families into functional subtypes purely on sequence criteria. comparison with an alignment-based approach was also carried out to test the quality of the clustering.  we demonstrate the robustness of rcm with reduced alphabets in clustering of protein sequences into families in a simulated dataset and seven well-characterized protein datasets. on protein datasets, crotonases, mandelate racemases, nucleotidyl cyclases and glycoside hydrolase family were clustered into subfamilies with 100% accuracy whereas acyl transferase domains, haloacid dehalogenases, and vicinal oxygen chelates could be assigned to subfamilies with %, % and % accuracies, respectively.  the overall combination of methods in this paper is useful for clustering protein families into subtypes based on solely protein sequence information. the method is also flexible and computationally fast because it does not require multiple alignment of sequences. 