BACKGROUND
since the 1970s, the evolution of protein sequences has been modelled using empirical amino acid substitution matrices derived from analyses of databases of protein alignments.
since the first introduction of these 'accepted point mutation'  models by dayhoff and coworkers  <cit> , a variety of newer substitution matrices have been derived based on much larger databases of alignments   <cit> , wag  <cit> ) or databases of proteins encoded by specific genome types  and using more rigorous statistical methods .
however, it has long been recognized that different sequence positions evolve at different rates.
indeed, a significant improvement in the fit of these models to real data has been to model heterogeneity in rates at different sites using a discrete approximation to the gamma distribution  <cit> .
yet it is well known that sites in proteins not only differ in their relative rates of evolution, but, because of structural and functional constraints, they also differ in their preferences for specific amino acids.
some sites in a protein alignment are occupied by almost any residue, while others appear to be restricted to a limited subset of amino acids and, quite frequently, only one particular residue.
attempts at improving substitution matrices for database searching, take these forms of substitution heterogeneity into account with the development of position-specific scoring matrices  <cit> , profile-based methods  <cit> , hidden markov models  <cit>  and structure-specific substitution matrices  <cit> .
however, accounting for site- or structure-specific dynamics in amino acid replacements in protein phylogenetic models has only recently garnered significant attention.
bruno  <cit>  proposed a model where site-specific amino acid frequencies were estimated by maximum likelihood .
however, this model is problematic because the number of parameters increases without bound   <cit>  and a large number of taxa are required for model fitting.
goldman and coworkers  <cit>  introduced a set of eight to ten predefined categories of substitution patterns at sites in a hidden markov model framework, based on protein secondary structures and surface exposure, and each category has its own rate matrix for ml inference.
other models explore the interdependence of sites due to constraints introduced by tertiary structure in protein sequence evolution  <cit> .
lartillot and philippe  <cit>  proposed a bayesian mixture model that allows amino acid replacement pattern at different sites to be described by distinct substitution processes which have the same substitution rates but different stationary probabilities.
they implemented their cat model in a bayesian monte carlo framework with a dirichlet process prior.
more recently le et al.
<cit>  proposed a new amino acid profile mixture model in which substitutions at sites follow a 'proportional' model whereby site-specific substitutions are entirely characterized by a mixture of  <dig> to  <dig> equilibrium frequency classes at sites.
lartillot, le and colleagues  <cit>  argue that taking into account the site-specific nature of protein evolution may be of vital importance to phylogenetic estimation especially in the case where two or more branches are extremely long, leading to an apparent long branch attraction  type artefact in empirical data sets.
they, and a recent study by rokas and carroll  <cit> , have shown that 'homoplasy'  occurs much more frequently in true protein alignments than expected under standard substitution models such as jtt or wag, even when the rates-across-sites process is taken into account.
as a result, they suggest that even probabilistic methods  employing these standard models can display an lba bias, even when a large number of sites are considered.
lartillot's cat model and the methods proposed by le et al.
are designed to counter these problems, although at the cost of a large number of additional parameters to be estimated and some model simplifications.
here we revisit the issue of site-specific amino acid constraints in protein phylogenetics.
first, to further probe the differences between the 'true' substitution process and standard models, we assembled  <dig> large protein sequence alignments and used two different methods to test if and how empirical frequencies at sites differ from those simulated under the standard jtt + f + Γ substitution model.
we showed that significant deviations can be detected for the majority of these protein families.
second, using the site-specific amino acid frequencies estimated from one of the data sets, we simulated a four-taxon case over a large grid of different branch-length settings to evaluate the accuracy of the ml methods employing a standard empirical matrix to recover the correct tree under these conditions.
we found a large 'felsenstein zone' where the lba artefact occurs.
third, we conducted a principal components analysis of the amino acid frequency matrix at all sites of the  <dig> protein alignment data sets and obtained four major classes  of amino acid frequency distribution at sites.
we propose a random effects mixture model using these class frequencies to model site-specific amino acid frequency distributions and implemented it in a version of raxml  <cit>  that we call q-matrix mixture raxml  for phylogenetic inference.
this model differs from previous models  <cit>  by accounting for intrinsic exchangeabilities between amino acids and containing standard amino acid substitution models as a special case, thereby permitting likelihood ratio testing of improvements in model fit.
we show that the amino acid frequency mixture model fits the data significantly better than the conventional non-mixture model in all cases examined and further find it can reduce the lba artefact both in simulations and in a phylogenomic analysis of a eukaryotic data set.
in comparison to other approaches, our model introduces significantly fewer additional parameters and avoids model over-simplification.
RESULTS
statistical analyses of site-specific amino acid uniformity and state frequencies
an entropy-based measure was used to quantify the deviation of site-specific amino acid frequencies from uniform usage of amino acids.
a z-test was used to determine whether the real data are more or less uniform compared to a very large data set simulated under jtt + f + Γ for the same tree.
a total of  <dig> protein data sets  and associated simulated data  were used and the sites of the real and simulated data were divided into four estimated rate categories and z-tests were carried out on each rate category.
the p-values for the tests are shown in table  <dig>  these indicate that for the sites in the slowest rate category  less than half of the datasets have significant differences in amino acid uniformity between the real data and the simulated data; but for sites in the faster rate categories , the real data are less uniform than the simulated data in the majority of the cases.
p-values: ** <  <dig> ; * <  <dig> .
the protein family abbreviations are: carboxyl_trans, acetyl-coa carboxylase; cpn <dig>  60-kda chaperonin; dna topo iv, dna topoisomerase iv subunit a ; ef-1α, elongation factor 1α; ef <dig>  elongation factor 2; filament, intermediate filament protein; glu_synth_ntn, glutamate synthase aminotransferase; hsp <dig>  70-kda heat shock protein; hsp <dig>  90-kda heat shock protein; ilvd_edd, dehydratase family proteins; nuof, nadh dehydrogenase i chain f; mcm, minichromosome maintenance protein; mpp, mitochondrial processing peptidase sequences; mreb, a bacterial homolog of the eukaryotic actin; poty_coat, potyvirus coat protein; usher, fimbrial usher protein.
collectively, these tests indicate that real data sets tend to have less uniform amino acid frequencies and fewer states at sites than expected under standard phylogenetic substitution models such as jtt + f + Γ.
four-taxon tree simulations under a site-specific frequency model and average frequency of the whole data set
in order to evaluate the potential impact of restricted site-specific amino acid frequencies on ml-based phylogenetic inference we did simulations of four-taxon trees under an extreme 'site-specific frequency'  model.
for this model, the amino acid frequencies at each site were calculated from the hsp <dig> data set.
these frequencies were then used in a jtt + ssf + Γ model to simulate data sets over the four-taxon trees over a wide range of branch-length settings.
note that for this model, each site is simulated with stationary frequencies corresponding to the frequencies observed in a given hsp <dig> alignment column.
for each branch-length setting, trees were then estimated by ml under jtt + f + Γ with tree-puzzle  <cit>  .
for comparison, we also simulated data sets for the same branch-length settings using the average amino acid frequencies observed in the whole hsp <dig> data set using the standard jtt + f + Γ model.
trees for these latter simulated data sets were also estimated by ml under the jtt + f + Γ model to evaluate a case where there is no model misspecification .
the results show that for the data simulated under the jtt + ssf + Γ model but with phylogenies estimated under the jtt + f + Γ model there is a very serious lba bias such that the two taxa with long branches  group together .
moreover, the fraction of other incorrectly estimated trees  is also quite high.
in contrast, the lba topology and the third incorrect topology are much less frequently estimated from the data simulated under the jtt + f + Γ model , although a small lba bias is observed.
principal components analysis of site frequency data
the above simulation results show definitively that, if ignored, site-specific amino acid frequencies will cause significant lba problems for phylogenetic inference.
although bruno   <cit>  attempted to account for this phenomenon by modelling amino acid frequencies at each site of the alignment, this approach leads to a serious statistical problem  <cit>  whereby the number of model parameters increases linearly with the amount of data.
if, however, certain patterns of amino acid usage are recurrent, one may derive some common frequency vectors, or profiles, from a large number of sequence sites.
to determine whether there were such recurrent patterns in the data, we calculated the  <dig> amino acid frequencies from each of the  <dig> sites of the  <dig> protein data sets to form a matrix of  <dig> sites ×  <dig> frequencies.
we then carried out a principal components analysis  and plotted the first two components .
these two components account for  <dig> % of the variance of the data.
we used a clustering method to cut the distribution of the sites into four classes along the lines of linear regressions .
the first class has high frequencies of valine, isoleucine and leucine, followed by methionine.
the second class has high frequencies of glycine, followed by alanine and serine.
the third class is rich in aspartic acid and glutamic acid.
the fourth class, a 'left-over' class that contains those sites that are not obvious members of the first three classes, has much more uniform frequencies of the various amino acids.
the predominant amino acids in the first three classes are consistent with the observation that the amino acids in those classes are biochemically and structurally similar and expected to be more interchangeable over evolutionary time.
figure  <dig> shows the average amino acid frequencies of the four pca-derived classes as well as the overall average amino acid frequencies of the  <dig> data sets.
testing a class frequency mixture model
using the four amino acid frequency profiles from the pca analysis and adding a fifth one corresponding to the average frequency of the whole data set to account for site frequencies not captured by the pca classes and also to make the relevant model comparisons easier , we implemented a 'class frequency'  mixture model.
in this model, the site likelihood is calculated as a weighted sum of the site likelihoods conditional on each class frequency or the whole data set frequency separately.
we also account for rates-across-sites using standard discretized gamma mixture model methods.
the cf mixture model has been implemented in qmmraxml, based on the source code of the phylogenetic inference package raxml.
we used qmmraxml to calculate likelihoods of  <dig> protein data sets under the single frequency model  and the cf mixture model , respectively.
the first  <dig> data sets are the same ones that were used for deriving the class frequencies from the pca, while the last four data sets are additional protein alignments used to test the generalizability of the cf model to other data sets.
for the analyses under the cf mixture model, we fixed the tree topologies to be the same optimal trees recovered under the jtt + f + Γ model but branch-lengths were re-optimized.
qmmraxml uses an expectation-maximization  algorithm to optimize weights of the class frequencies and the whole data set frequency.
the estimated weights  and the likelihood differences  between the two models are listed in table  <dig>  in all cases there are significant likelihood increases under the cf mixture model compared to the single frequency model, indicating the cf model always fits the data better for the same topology.
curiously, in all cases the weight of the f class  is generally high and the weights of the four pca classes are generally low, especially in classes  <dig> and  <dig>  the weight for a class estimates the probability that a site has a frequency vector corresponding to that class.
the reasons for the relatively small class 1– <dig> weights are several-fold.
first, it seems likely that these frequency classes are related to structural features.
for instance, the two major amino acids of class  <dig>  aspartic acid and glutamic acid, are negatively charged and many sites that are conserved to have only these amino acids could have important structural roles such as binding metal cations, participating in intramolecular and/or intermolecular salt-bridges, or have catalytic or substrate-binding functions.
although the number of such sites will vary across protein families, it is likely that they will always constitute a relatively small minority in any given protein family.
the data bears this out in other ways.
most sites that do not have a very pronounced preference for the amino acids emphasized in classes 1– <dig> or that have appreciable frequencies for more than four amino acids, will be fit best by either class  <dig> or the f class.
since classes 1– <dig> emphasize sites with very few amino acids each of which occur with frequencies of < 10% overall, it is not surprising that the estimated frequency of these kinds of sites are all quite low.
for example, aspartic acid and glutamic acid together have an average frequency of  <dig> % in the  <dig> data sets.
therefore, the sites restricted to having virtually only these two amino acids, as featured in class  <dig>  are expected to be very small and is reflected by the low weights assigned to this class.
Λlnl is the likelihood difference between the cf mixture model and the single frequency model .
the p-values associated with these differences, calculated from χ <dig> tests with  <dig> degrees of freedom, are very significant in all cases .
the actual p-values would be even smaller as the tests are conservative .
in the above simulation studies we have demonstrated that sequence data generated under the jtt + ssf + Γ model can cause a serious lba problem when the trees are estimated under the conventional jtt + f + Γ model.
to test whether the cf mixture model can ameliorate the lba problem we used both simulations and analysed a real data set.
figure  <dig> shows the simulation results from data generated under jtt + ssf + Γ model.
the left panel shows the results of estimation under a standard jtt + f + Γ model and the right panel shows estimation under the jtt + cf + Γ model.
from comparing these two results, it is clear that the cf model ameliorates some, but not all of the lba problems relative to the standard model.
it is interesting to note that in the non-felsenstein zone region, the cf model sometimes performs slightly worse than the standard model.
this is consistent with higher variance estimates expected from a model that includes more parameters.
the deep phylogeny of eukaryote 'supergroups' is often plagued with lba  <cit> , even when large multi-gene phylogenomic data sets are used  <cit> .
one of the most famous examples of this concerns the position of microsporidia, a group of fast-evolving intracellular parasites that are now known to be relatives of fungi  <cit> .
when reconstructing the phylogeny of eukaryotes rooted by archaea with conventional models, the extremely long branch leading to microsporidia is often attracted to the long branches leading to the archaea at the base of the eukaryotes.
many methods have been proposed to solve this problem, including selective taxon sampling, removal of fast-evolving proteins and saturated sites  <cit> , accounting for covarion shifts  <cit> , amino acid profile mixture modelling  <cit> , etc.
here we applied the cf mixture model to a large set of eukaryote phylogenomic data with  <dig> proteins from  <dig> taxa and  <dig> sites  <cit>  and calculated likelihoods of two competing trees: the lba topology where microsporidia group with archaea and, the correct topology where microsporidia group with fungi.
as shown in table  <dig>  while the jtt + f + Γ model supports the microsporidia plus archaea clan, the cf mixture model supports the correct microsporidia-fungi tree.
thus for a real example the cf model appears to be less susceptible to the effects of lba than the standard model.
comparisons to other methods
lartillot and colleagues have shown that accounting for site-specific amino acid frequencies with their cat + Γ model seems to significantly improve both model fit and phylogenetic estimation with large concatenated protein data sets  <cit> .
indeed, they were the first to observe that accounting for site-specific frequencies is important to avoid long-branch attraction problems when there are large numbers of substitutions .
however, as these authors later note, the standard cat + Γ model fits very large numbers of classes that, while appropriate for large concatenated data sets, likely leads to problems with over-parameterization and convergence of bayesian analyses in the case of smaller alignments   <cit> .
to address this, le and colleagues  <cit>  developed similar models based on a set of  <dig> to  <dig> classes of amino acid profiles that were estimated from a database of alignments and implemented these in both bayesian  and maximum likelihood  estimation programs.
however, again, for computational efficiency reasons, the authors restricted attention to proportional models that ignore different intrinsic 'exchangeabilities' between amino acids.
this, and the fact that weights associated with the  <dig> to  <dig> classes should be estimated for every tree, suggests that problems may still exist with both over-simplification of the substitution process , overparameterization of the models and computational efficiency.
our model is developed with more of a 'bottom-up' style approach; we have introduced very few frequency classes, that, based on our pca of a set of real alignments, seem to be the most important amino acid profiles.
furthermore, we use the exchangeabilities from standard models  and include the overall frequencies as one of the classes.
this setup allows nested model comparisons for likelihood ratio tests to directly examine improvement in model fit by the introduction of new classes.
it is interesting to compare this approach to that described by le and colleagues  <cit> .
le and colleagues have estimated a larger number of frequency classes directly by the method of maximum likelihood from a large database of alignments and therefore these classes in principle should fit the data they were estimated from better than the pca approach described here.
however, because of the assumption of uniform exchangeabilities, they are likely to estimate too many classes some of which are distorted frequency classes that adjust for increased exchangeabilities present in more complex and realistic rate matrices than the uniform matrix.
our approach may avoid this situation and therefore may require fewer frequency classes to adequately capture the site-specific nature of protein evolution.
furthermore, by including the data set frequencies as a fifth component, sites that show little preference for a restricted amino acid profile, but evolve according to the exchangeabilities of the jtt matrix can be accommodated.
nevertheless, despite the improvements in phylogenetic estimation we have found, the four classes we introduced from the pca may in fact be too few to adequately capture the diversity of 'site-specific' preferences in amino acids.
for instance, we did not recover a class where the basic amino acids lysine and arginine predominate, even though sites that rapidly switch between these amino acids are clearly observable in protein families  <cit> .
the fact that the first two components of our pca account for only ~21% of the variance in the data indicates that inclusion of additional classes by investigation of the third or more principal components from the pca or other data mining methods such as the self-organizing map  <cit>  may be fruitful future directions to take.
in any case, our implementation of these methods in a generic q-matrix mixture model in qmmraxml, allows the user to implement any number of exchangeability matrices plus associated stationary frequency vectors to freely explore improvements in model fit in protein evolution.
CONCLUSIONS
we report the results of two statistical tests – the amino acid frequency uniformity test and state counts test – that demonstrate that in real protein alignments there are fewer states at sites and the frequencies of these states are less uniform than predicted by jtt + f + Γ model.
we show that use of standard 'average' frequency models like jtt + f + Γ for phylogenetic estimation when the data are simulated with site-specific frequencies leads to serious lba artefacts.
a pca of site-specific frequency vectors of  <dig> sites from  <dig> protein data sets revealed four major classes of sites.
these classes can be used in a simple class-frequency  mixture model for modelling site-specific distributions for phylogenetic inference that we have implemented in a program called qmmraxml.
likelihood ratio tests indicate a large improvement in the fit of jtt + cf + Γ over jtt + f + Γ for all data sets examined.
furthermore, the cf mixture model appears to ameliorate the long-branch attraction problems, in both simulation studies and in analyses of a phylogenomic data set.
the cf mixture model provides a new method for modelling site-specific compositional heterogeneity and qmmraxml is a promising tool for exploring model fit in protein evolution and reconstructing more accurate phylogenies.

 since the 1970s, the evolution of protein sequences has been modelled using empirical amino acid substitution matrices derived from analyses of databases of protein alignments. since the first introduction of these'accepted point mutation' models by dayhoff and coworkers, a variety of newer substitution matrices have been derived based on much larger databases of alignments, wag ) or databases of proteins encoded by specific genome types and using more rigorous statistical methods. however, it has long been recognized that different sequence positions evolve at different rates. indeed, a significant improvement in the fit of these models to real data has been to model heterogeneity in rates at different sites using a discrete approximation to the gamma distribution. yet it is well known that sites in proteins not only differ in their relative rates of evolution, but, because of structural and functional constraints, they also differ in their preferences for specific amino acids. some sites in a protein alignment are occupied by almost any residue, while others appear to be restricted to a limited subset of amino acids and, quite frequently, only one particular residue. attempts at improving substitution matrices for database searching, take these forms of substitution heterogeneity into account with the development of position-specific scoring matrices, profile-based methods, hidden markov models and structure-specific substitution matrices. however, accounting for site- or structure-specific dynamics in amino acid replacements in protein phylogenetic models has only recently garnered significant attention. bruno proposed a model where site-specific amino acid frequencies were estimated by maximum likelihood. however, this model is problematic because the number of parameters increases without bound and a large number of taxa are required for model fitting. goldman and coworkers introduced a set of eight to ten predefined categories of substitution patterns at sites in a hidden markov model framework, based on protein secondary structures and surface exposure, and each category has its own rate matrix for ml inference. other models explore the interdependence of sites due to constraints introduced by tertiary structure in protein sequence evolution. lartillot and philippe proposed a bayesian mixture model that allows amino acid replacement pattern at different sites to be described by distinct substitution processes which have the same substitution rates but different stationary probabilities. they implemented their cat model in a bayesian monte carlo framework with a dirichlet process prior. more recently le et al. proposed a new amino acid profile mixture model in which substitutions at sites follow a'proportional' model whereby site-specific substitutions are entirely characterized by a mixture of to equilibrium frequency classes at sites. lartillot, le and colleagues argue that taking into account the site-specific nature of protein evolution may be of vital importance to phylogenetic estimation especially in the case where two or more branches are extremely long, leading to an apparent long branch attraction type artefact in empirical data sets. they, and a recent study by rokas and carroll, have shown that'homoplasy' occurs much more frequently in true protein alignments than expected under standard substitution models such as jtt or wag, even when the rates-across-sites process is taken into account. as a result, they suggest that even probabilistic methods employing these standard models can display an lba bias, even when a large number of sites are considered. lartillot's cat model and the methods proposed by le et al. are designed to counter these problems, although at the cost of a large number of additional parameters to be estimated and some model simplifications. here we revisit the issue of site-specific amino acid constraints in protein phylogenetics. first, to further probe the differences between the'true' substitution process and standard models, we assembled large protein sequence alignments and used two different methods to test if and how empirical frequencies at sites differ from those simulated under the standard jtt + f + Γ substitution model. we showed that significant deviations can be detected for the majority of these protein families. second, using the site-specific amino acid frequencies estimated from one of the data sets, we simulated a four-taxon case over a large grid of different branch-length settings to evaluate the accuracy of the ml methods employing a standard empirical matrix to recover the correct tree under these conditions. we found a large'felsenstein zone' where the lba artefact occurs. third, we conducted a principal components analysis of the amino acid frequency matrix at all sites of the protein alignment data sets and obtained four major classes of amino acid frequency distribution at sites. we propose a random effects mixture model using these class frequencies to model site-specific amino acid frequency distributions and implemented it in a version of raxml that we call q-matrix mixture raxml for phylogenetic inference. this model differs from previous models by accounting for intrinsic exchangeabilities between amino acids and containing standard amino acid substitution models as a special case, thereby permitting likelihood ratio testing of improvements in model fit. we show that the amino acid frequency mixture model fits the data significantly better than the conventional non-mixture model in all cases examined and further find it can reduce the lba artefact both in simulations and in a phylogenomic analysis of a eukaryotic data set. in comparison to other approaches, our model introduces significantly fewer additional parameters and avoids model over-simplification.  statistical analyses of site-specific amino acid uniformity and state frequencies an entropy-based measure was used to quantify the deviation of site-specific amino acid frequencies from uniform usage of amino acids. a z-test was used to determine whether the real data are more or less uniform compared to a very large data set simulated under jtt + f + Γ for the same tree. a total of protein data sets and associated simulated data were used and the sites of the real and simulated data were divided into four estimated rate categories and z-tests were carried out on each rate category. the p-values for the tests are shown in table these indicate that for the sites in the slowest rate category less than half of the datasets have significant differences in amino acid uniformity between the real data and the simulated data; but for sites in the faster rate categories, the real data are less uniform than the simulated data in the majority of the cases. p-values: ** <; * <. the protein family abbreviations are: carboxyl_trans, acetyl-coa carboxylase; cpn 60-kda chaperonin; dna topo iv, dna topoisomerase iv subunit a; ef-1α, elongation factor 1α; ef elongation factor 2; filament, intermediate filament protein; glu_synth_ntn, glutamate synthase aminotransferase; hsp 70-kda heat shock protein; hsp 90-kda heat shock protein; ilvd_edd, dehydratase family proteins; nuof, nadh dehydrogenase i chain f; mcm, minichromosome maintenance protein; mpp, mitochondrial processing peptidase sequences; mreb, a bacterial homolog of the eukaryotic actin; poty_coat, potyvirus coat protein; usher, fimbrial usher protein. collectively, these tests indicate that real data sets tend to have less uniform amino acid frequencies and fewer states at sites than expected under standard phylogenetic substitution models such as jtt + f + Γ. four-taxon tree simulations under a site-specific frequency model and average frequency of the whole data set in order to evaluate the potential impact of restricted site-specific amino acid frequencies on ml-based phylogenetic inference we did simulations of four-taxon trees under an extreme'site-specific frequency' model. for this model, the amino acid frequencies at each site were calculated from the hsp data set. these frequencies were then used in a jtt + ssf + Γ model to simulate data sets over the four-taxon trees over a wide range of branch-length settings. note that for this model, each site is simulated with stationary frequencies corresponding to the frequencies observed in a given hsp alignment column. for each branch-length setting, trees were then estimated by ml under jtt + f + Γ with tree-puzzle. for comparison, we also simulated data sets for the same branch-length settings using the average amino acid frequencies observed in the whole hsp data set using the standard jtt + f + Γ model. trees for these latter simulated data sets were also estimated by ml under the jtt + f + Γ model to evaluate a case where there is no model misspecification. the results show that for the data simulated under the jtt + ssf + Γ model but with phylogenies estimated under the jtt + f + Γ model there is a very serious lba bias such that the two taxa with long branches group together. moreover, the fraction of other incorrectly estimated trees is also quite high. in contrast, the lba topology and the third incorrect topology are much less frequently estimated from the data simulated under the jtt + f + Γ model, although a small lba bias is observed. principal components analysis of site frequency data the above simulation results show definitively that, if ignored, site-specific amino acid frequencies will cause significant lba problems for phylogenetic inference. although bruno attempted to account for this phenomenon by modelling amino acid frequencies at each site of the alignment, this approach leads to a serious statistical problem whereby the number of model parameters increases linearly with the amount of data. if, however, certain patterns of amino acid usage are recurrent, one may derive some common frequency vectors, or profiles, from a large number of sequence sites. to determine whether there were such recurrent patterns in the data, we calculated the amino acid frequencies from each of the sites of the protein data sets to form a matrix of sites × frequencies. we then carried out a principal components analysis and plotted the first two components. these two components account for % of the variance of the data. we used a clustering method to cut the distribution of the sites into four classes along the lines of linear regressions. the first class has high frequencies of valine, isoleucine and leucine, followed by methionine. the second class has high frequencies of glycine, followed by alanine and serine. the third class is rich in aspartic acid and glutamic acid. the fourth class, a'left-over' class that contains those sites that are not obvious members of the first three classes, has much more uniform frequencies of the various amino acids. the predominant amino acids in the first three classes are consistent with the observation that the amino acids in those classes are biochemically and structurally similar and expected to be more interchangeable over evolutionary time. figure shows the average amino acid frequencies of the four pca-derived classes as well as the overall average amino acid frequencies of the data sets. testing a class frequency mixture model using the four amino acid frequency profiles from the pca analysis and adding a fifth one corresponding to the average frequency of the whole data set to account for site frequencies not captured by the pca classes and also to make the relevant model comparisons easier, we implemented a'class frequency' mixture model. in this model, the site likelihood is calculated as a weighted sum of the site likelihoods conditional on each class frequency or the whole data set frequency separately. we also account for rates-across-sites using standard discretized gamma mixture model methods. the cf mixture model has been implemented in qmmraxml, based on the source code of the phylogenetic inference package raxml. we used qmmraxml to calculate likelihoods of protein data sets under the single frequency model and the cf mixture model, respectively. the first data sets are the same ones that were used for deriving the class frequencies from the pca, while the last four data sets are additional protein alignments used to test the generalizability of the cf model to other data sets. for the analyses under the cf mixture model, we fixed the tree topologies to be the same optimal trees recovered under the jtt + f + Γ model but branch-lengths were re-optimized. qmmraxml uses an expectation-maximization algorithm to optimize weights of the class frequencies and the whole data set frequency. the estimated weights and the likelihood differences between the two models are listed in table in all cases there are significant likelihood increases under the cf mixture model compared to the single frequency model, indicating the cf model always fits the data better for the same topology. curiously, in all cases the weight of the f class is generally high and the weights of the four pca classes are generally low, especially in classes and the weight for a class estimates the probability that a site has a frequency vector corresponding to that class. the reasons for the relatively small class 1– weights are several-fold. first, it seems likely that these frequency classes are related to structural features. for instance, the two major amino acids of class aspartic acid and glutamic acid, are negatively charged and many sites that are conserved to have only these amino acids could have important structural roles such as binding metal cations, participating in intramolecular and/or intermolecular salt-bridges, or have catalytic or substrate-binding functions. although the number of such sites will vary across protein families, it is likely that they will always constitute a relatively small minority in any given protein family. the data bears this out in other ways. most sites that do not have a very pronounced preference for the amino acids emphasized in classes 1– or that have appreciable frequencies for more than four amino acids, will be fit best by either class or the f class. since classes 1– emphasize sites with very few amino acids each of which occur with frequencies of < 10% overall, it is not surprising that the estimated frequency of these kinds of sites are all quite low. for example, aspartic acid and glutamic acid together have an average frequency of % in the data sets. therefore, the sites restricted to having virtually only these two amino acids, as featured in class are expected to be very small and is reflected by the low weights assigned to this class. Λlnl is the likelihood difference between the cf mixture model and the single frequency model. the p-values associated with these differences, calculated from χ tests with degrees of freedom, are very significant in all cases. the actual p-values would be even smaller as the tests are conservative. in the above simulation studies we have demonstrated that sequence data generated under the jtt + ssf + Γ model can cause a serious lba problem when the trees are estimated under the conventional jtt + f + Γ model. to test whether the cf mixture model can ameliorate the lba problem we used both simulations and analysed a real data set. figure shows the simulation results from data generated under jtt + ssf + Γ model. the left panel shows the results of estimation under a standard jtt + f + Γ model and the right panel shows estimation under the jtt + cf + Γ model. from comparing these two results, it is clear that the cf model ameliorates some, but not all of the lba problems relative to the standard model. it is interesting to note that in the non-felsenstein zone region, the cf model sometimes performs slightly worse than the standard model. this is consistent with higher variance estimates expected from a model that includes more parameters. the deep phylogeny of eukaryote'supergroups' is often plagued with lba, even when large multi-gene phylogenomic data sets are used. one of the most famous examples of this concerns the position of microsporidia, a group of fast-evolving intracellular parasites that are now known to be relatives of fungi. when reconstructing the phylogeny of eukaryotes rooted by archaea with conventional models, the extremely long branch leading to microsporidia is often attracted to the long branches leading to the archaea at the base of the eukaryotes. many methods have been proposed to solve this problem, including selective taxon sampling, removal of fast-evolving proteins and saturated sites, accounting for covarion shifts, amino acid profile mixture modelling, etc. here we applied the cf mixture model to a large set of eukaryote phylogenomic data with proteins from taxa and sites and calculated likelihoods of two competing trees: the lba topology where microsporidia group with archaea and, the correct topology where microsporidia group with fungi. as shown in table while the jtt + f + Γ model supports the microsporidia plus archaea clan, the cf mixture model supports the correct microsporidia-fungi tree. thus for a real example the cf model appears to be less susceptible to the effects of lba than the standard model. comparisons to other methods lartillot and colleagues have shown that accounting for site-specific amino acid frequencies with their cat + Γ model seems to significantly improve both model fit and phylogenetic estimation with large concatenated protein data sets. indeed, they were the first to observe that accounting for site-specific frequencies is important to avoid long-branch attraction problems when there are large numbers of substitutions. however, as these authors later note, the standard cat + Γ model fits very large numbers of classes that, while appropriate for large concatenated data sets, likely leads to problems with over-parameterization and convergence of bayesian analyses in the case of smaller alignments. to address this, le and colleagues developed similar models based on a set of to classes of amino acid profiles that were estimated from a database of alignments and implemented these in both bayesian and maximum likelihood estimation programs. however, again, for computational efficiency reasons, the authors restricted attention to proportional models that ignore different intrinsic'exchangeabilities' between amino acids. this, and the fact that weights associated with the to classes should be estimated for every tree, suggests that problems may still exist with both over-simplification of the substitution process, overparameterization of the models and computational efficiency. our model is developed with more of a'bottom-up' style approach; we have introduced very few frequency classes, that, based on our pca of a set of real alignments, seem to be the most important amino acid profiles. furthermore, we use the exchangeabilities from standard models and include the overall frequencies as one of the classes. this setup allows nested model comparisons for likelihood ratio tests to directly examine improvement in model fit by the introduction of new classes. it is interesting to compare this approach to that described by le and colleagues. le and colleagues have estimated a larger number of frequency classes directly by the method of maximum likelihood from a large database of alignments and therefore these classes in principle should fit the data they were estimated from better than the pca approach described here. however, because of the assumption of uniform exchangeabilities, they are likely to estimate too many classes some of which are distorted frequency classes that adjust for increased exchangeabilities present in more complex and realistic rate matrices than the uniform matrix. our approach may avoid this situation and therefore may require fewer frequency classes to adequately capture the site-specific nature of protein evolution. furthermore, by including the data set frequencies as a fifth component, sites that show little preference for a restricted amino acid profile, but evolve according to the exchangeabilities of the jtt matrix can be accommodated. nevertheless, despite the improvements in phylogenetic estimation we have found, the four classes we introduced from the pca may in fact be too few to adequately capture the diversity of'site-specific' preferences in amino acids. for instance, we did not recover a class where the basic amino acids lysine and arginine predominate, even though sites that rapidly switch between these amino acids are clearly observable in protein families. the fact that the first two components of our pca account for only ~21% of the variance in the data indicates that inclusion of additional classes by investigation of the third or more principal components from the pca or other data mining methods such as the self-organizing map may be fruitful future directions to take. in any case, our implementation of these methods in a generic q-matrix mixture model in qmmraxml, allows the user to implement any number of exchangeability matrices plus associated stationary frequency vectors to freely explore improvements in model fit in protein evolution.  we report the results of two statistical tests – the amino acid frequency uniformity test and state counts test – that demonstrate that in real protein alignments there are fewer states at sites and the frequencies of these states are less uniform than predicted by jtt + f + Γ model. we show that use of standard'average' frequency models like jtt + f + Γ for phylogenetic estimation when the data are simulated with site-specific frequencies leads to serious lba artefacts. a pca of site-specific frequency vectors of sites from protein data sets revealed four major classes of sites. these classes can be used in a simple class-frequency mixture model for modelling site-specific distributions for phylogenetic inference that we have implemented in a program called qmmraxml. likelihood ratio tests indicate a large improvement in the fit of jtt + cf + Γ over jtt + f + Γ for all data sets examined. furthermore, the cf mixture model appears to ameliorate the long-branch attraction problems, in both simulation studies and in analyses of a phylogenomic data set. the cf mixture model provides a new method for modelling site-specific compositional heterogeneity and qmmraxml is a promising tool for exploring model fit in protein evolution and reconstructing more accurate phylogenies. 
to this end, we have implemented a class frequency mixture model in a freely available program called qmmraxml for phylogenetic estimation. here we examine the degree to which the pattern of evolution at sites differs from that expected based on empirical amino acid substitution models and evaluate the impact of these deviations on phylogenetic estimation. with a four-taxon example, we show that phylogenetic estimation under the jtt + f + Γ model is seriously biased by a long-branch attraction artefact if the data are simulated under a model utilizing the observed site-specific amino acid frequencies from an alignment. using a mixture model with these four separate classes of site-specific state frequencies plus a fifth class of global frequencies, significant improvements in model fit for real data sets can be achieved. this simple mixture model also reduces the long-branch attraction problem, as shown by simulations and analyses of a real phylogenomic data set. widely used substitution models for proteins, such as the jones-taylor-thornton or whelan and goldman models, are based on empirical amino acid interchange matrices estimated from databases of protein alignments that incorporate the average amino acid frequencies of the data set under examination. we analyzed large protein alignments with two statistical tests designed to detect deviation of site-specific amino acid distributions from data simulated under the standard empirical substitution model: jtt+ f + Γ. we found that the number of states at a given site is, on average, smaller and the frequencies of these states are less uniform than expected based on a jtt + f + Γ substitution model. principal components analyses indicate the existence of at least four major site-specific frequency classes in these protein alignments.
 widely used substitution models for proteins, such as the jones-taylor-thornton or whelan and goldman models, are based on empirical amino acid interchange matrices estimated from databases of protein alignments that incorporate the average amino acid frequencies of the data set under examination. variation in the evolutionary process between sites is typically modelled by a rates-across-sites distribution such as the gamma distribution. however, sites in proteins also vary in the kinds of amino acid interchanges that are favoured, a feature that is ignored by standard empirical substitution matrices. here we examine the degree to which the pattern of evolution at sites differs from that expected based on empirical amino acid substitution models and evaluate the impact of these deviations on phylogenetic estimation.  we analyzed large protein alignments with two statistical tests designed to detect deviation of site-specific amino acid distributions from data simulated under the standard empirical substitution model: jtt+ f + Γ. we found that the number of states at a given site is, on average, smaller and the frequencies of these states are less uniform than expected based on a jtt + f + Γ substitution model. with a four-taxon example, we show that phylogenetic estimation under the jtt + f + Γ model is seriously biased by a long-branch attraction artefact if the data are simulated under a model utilizing the observed site-specific amino acid frequencies from an alignment. principal components analyses indicate the existence of at least four major site-specific frequency classes in these protein alignments. using a mixture model with these four separate classes of site-specific state frequencies plus a fifth class of global frequencies, significant improvements in model fit for real data sets can be achieved. this simple mixture model also reduces the long-branch attraction problem, as shown by simulations and analyses of a real phylogenomic data set.  protein families display site-specific evolutionary dynamics that are ignored by standard protein phylogenetic models. accurate estimation of protein phylogenies requires models that accommodate the heterogeneity in the evolutionary process across sites. to this end, we have implemented a class frequency mixture model in a freely available program called qmmraxml for phylogenetic estimation. 