# Phylogenetic invariants

Phylogenetic invariants are polynomial relationships between the frequencies of various site patterns in an idealized DNA multiple sequence alignment. They have received substantial study in the field of biomathematics, and they can be used to choose among phylogenetic tree topologies in an empirical setting. The primary advantage of phylogenetic invariants relative to other methods of phylogenetic estimation like maximum likelihood or Bayesian MCMC analyses is that invariants can yield information about the tree without requiring the estimation of branch lengths of model parameters. The idea of using phylogenetic invariants was introduced independently by James Cavender and Joseph Felsenstein and by James A. Lake in 1987.
At this point the number of programs that allow empirical datasets to be analyzed using invariants is limited. However, phylogenetic invariants may provide solutions to other problems in phylogenetics and they represent an area of active research for that reason. Felsenstein stated it best when he said, "invariants are worth attention, not for what they do for us now, but what they might lead to in the future." (p. 390)
If we consider a DNA (or RNA) multiple sequence alignment with t taxa and no gaps or missing data (i.e., an idealized multiple sequence alignment), there are 4$^t$ possible site patterns. For example, there are 256 possible site patterns for four taxa (f$_{AAAA}$, f$_{AAAC}$, f$_{AAAG}$, ... f$_{TTTT}$), which can be written as a vector. This site pattern frequency vector has 255 degrees of freedom because the frequencies must sum to one. However, any set of site pattern frequencies that resulted from some specific process of sequence evolution on a specific tree must obey many constraints. and therefore have many fewer degrees of freedom. Thus, there should be polynomials involving those frequencies that take on a value of zero if the DNA sequences were generated on a specific tree given a particular substitution model.
Invariants are formulas in the expected pattern frequencies, not the observed pattern frequencies. When they are computed using the observed pattern frequencies, we will usually find that they are not precisely zero even when the model and tree topology are correct. By testing whether such polynomials for various trees are 'nearly zero' when evaluated on the observed frequencies of patterns in real data sequences one should be able infer which tree best explains the data.
Some invariants are straightforward consequences of symmetries in the model of nucleotide substitution and they will take on a value of zero regardless of the underlying tree topology. For example, if we assume the Jukes-Cantor model of sequence evolution and a four-taxon tree we expect:
$f_{ACAT}-f_{CGCA}=0$
This is a simple outgrowth of the fact that base frequencies are constrained to be equal under the Jukes-Cantor model. Thus, they are called symmetry invariants. The equation shown above is only one of a large number of symmetry invariants for the Jukes-Cantor model; in fact, there are a total of 241 symmetry invariants for that model.
Symmetry invariants for the Jukes-Cantor model of DNA evolution (adapted from Felsenstein 2004) :
| Site pattern category | Example of site pattern | Number of pattern types | Number of patterns | Total invariants that result |
| 4x | xxxx (e.g., AAAA, CCCC, ...) | 1 | 4 | 3 |
| 3x, 1y | xxxy (e.g., AAAC, AACA, ...) | 4 | 12 | 44 |
| 2x, 2y | xxyy (e.g., AACC, ACCA, ...) | 3 | 12 | 33 |
| 2x, 1y, 1z | xxyz (e.g., AACG, ACGA, ...) | 6 | 24 | 138 |
| 1x, 1y, 1z, 1w | xyzw (e.g., ACGT, CGTA, ...) | 1 | 24 | 23 |
| Totals = | | 15 | | 241 |
Symmetry invariants are non-phylogenetic in nature; they take on the expected value of zero regardless of the tree topology. However, it is possible to determine whether a particular multiple sequence alignment fits the Jukes-Cantor model of evolution (i.e., by testing whether the site patterns of the appropriate types are present in equal numbers). More general tests for the best-fitting model using invariants are also possible. For example, Kedzierska et al. 2012 used invariants to establish the best-fitting model out from a specific model set.
Models of DNA evolution that can be tested using the Kedzierska et al. (2012) invariants method :
| Model abbreviation | Full model name |
| JC69* | Jukes-Cantor |
| K80* | Kimura two-parameter |
| K81* | Kimura three-parameter |
| SSM (CS05) | Strand-specific model |
| GMM | General Markov model |
The asterisk after the JC69, K80, and K81 models is used to emphasize the non-homogeneous nature of the models that can be examined using invariants. These non-homogeneous models include the commonly used continuous-time JC69, K80, and K81 models as submodels. The SSM (strand-specific model), also called the CS05 model, is a generalized non-homogeneous version of the HKY (Hasegawa-Kishino-Yano) model constrained to have equal distribution of the pairs of bases A,T and C,G at each node of the tree and no assumption regarding a stable base distribution. All models listed above are submodels of the general Markov model (GMM). The ability to perform tests using non-homogeneous models represents a major benefit of the invariants methods relative to the more commonly used maximum likelihood methods for phylogenetic model testing.
Phylogenetic invariants, which are defined as the subset of invariants that take on a value of zero only when the sequences were (or were not) generated on a specific topology, are likely to be the most useful invariants for phylogenetic studies.

## Lake's linear invariants

Lake's invariants (which he called "evolutionary parsimony") provide an excellent example of phylogenetic invariants. Lake's invariants involve quartets, two of which (the incorrect topologies) yield values of zero and one of which yields a value greater than zero. This can be used to construct a test based on following invariant relationship, which holds for the two incorrect trees when sites evolve under the Kimura two-parameter model of sequence evolution:
$f_{1133}+f_{1234}=f_{1233}+f_{1134}$
The indices of these site pattern frequencies indicate the bases scored relative to the base in the first taxon (which we call taxon A). If base 1 is a purine, then base 2 is the other purine and bases 3 and 4 are the pyrimidines. If base 1 is a pyrimidine, then base 2 is the other pyrimidine and. bases 3 and 4 are the purines.
We will call three possible quartet trees T$_X$ [in newick format T$_X$ is (A,B,(C,D));], T$_Y$ [T$_Y$ is (A,C,(B,D));], and T$_Z$ [T$_Z$ is (A,D,(B,C));]. We can calculate three values from the data to identify the best topology given the data:
$X=N_{1133}-N_{1233}-N_{1134}+N_{1234}$
$Y=N_{1313}-N_{1323}-N_{1314}+N_{1324}$
$Z=N_{1331}-N_{1332}-N_{1341}+N_{1342}$
Lake broke these values up into a "parsimony-like term" ($P=N_{1133}+N_{1234}$ for T$_X$) the "background term" ($B=N_{1233}+N_{1134}$ for T$_X$) and suggests testing for deviation from zero by calculating $\chi ^{2}=(P-B)^{2}/(P+B)$ and performing a χ$^2$ test with one degree of freedom. Similar χ$^2$ tests can be performed for Y and Z. If one of the three values is significantly different from zero the corresponding topology is the best estimate of phylogeny. The advantage of using Lake's invariants relative to maximum likelihood or neighbor joining of Kimura two-parameter distances is that the invariants should hold regardless of the model parameters, branch lengths, or patterns of among-sites rate heterogeneity.
A classic study by John Huelsenbeck and David Hillis found that Lake's invariants converges on the true tree over all of the branch length space they examined when the Kimura two-parameter (K80) model ist he underlying model of evolution. However, they also found that Lake's invariants are very inefficient (large amounts of data are necessary to converge on the correct tree). This inefficiency has caused most empiricists to abandon the use of Lake's invariants. Also, because Lake's invariants are based on the K80 model phylogenetic estimation using Lake's invariants may not yield the true tree when the model that generated the data strongly violates that model.

## Modern approaches using phylogenetic invariants

The low efficiency of Lake's invariants reflects the fact that it used a limited set of generators for the phylogenetic invariants. Casanellas et al. introduced methods to derive a much larger set of set of generators for DNA data and this has led to the development of invariants methods that are as efficient as maximum likelihood methods. Several of these methods have implementations that are practical for analyses of empirical datasets.
"Squangles" (stochastic quartet tangles) are another example of a modern invariants method and it has been implemented in software package that is practical to be used with empirical datasets. Squangles permit the choice among the three possible quartets assuming that DNA sequences have evolved under the general Markov model; the quartets can then be assembled using a supertree method. There are three squangles that are useful for differentiating among quartets, which can be denoted as q$_1$(f), q$_2$(f), and q$_3$(f) (f is a 256 element vector containing the site frequency spectrum). Each q has 66,744 terms and together they satisfy the linear relation q$_1$ + q$_2$ + q$_3$ = 0 (i.e., up to linear dependence there are only two q values). Each possible quartet has different expected values for q$_1$, q$_2$, and q$_3$:
Expected values for q$_1$, q$_2$, and q$_3$ (adapted from Holland et al. 2013) :
| Tree topology (newick format) | Quartet | E(q$_1$) | E(q$_2$) | E(q$_3$) |
| ((A,B),(C,D)); | AB|CD (or 12|34) | 0 | -u | u |
| ((A,C),(B,D)); | AC|BD (or 13|24) | v | 0 | -v |
| ((A,D),(B,C)); | AD|BC (or 14|23) | -w | w | 0 |
The expected values q$_1$, q$_2$, and q$_3$ are all zero on the star topology (a quartet with an internal branch length of zero). For practicality, Holland et al. used least squares to solve for the q values. Empirical tests of the squangles method have been limited but they appear to be promising.

## Singular value decomposition

Another important class of modern invariants methods is based on the use of singular value decomposition (SVD) to examine the rank of matrices corresponding to flattenings of a tensor with the site pattern frequency spectrum. The flattening matrices for a four taxon tree are constructed by arranging the nucleotide site pattern frequencies into a 16 x 16 matrix based on the bipartition in the tree. For example, using the T1 topology that was defined above as (A,B,(C,D)); we use the bipartition AB|CD and the elements of the flattening matrix would correspond to frequencies of the site pattern frequencies in the order ABCD. As shown below the rows are indexed using to the states for the first two taxa (A and B in the case of T1) and the columns reflect the states in the second pair of taxa (C and D in the case of T1).
$\mathbf {Flat_{T1}} ={\begin{bmatrix}p_{AAAA}&p_{AAAC}&p_{AAAG}&p_{AAAT}&p_{AACA}&\cdots &p_{AATT}\\p_{ACAA}&p_{ACAC}&p_{ACAG}&p_{ACAT}&p_{ACCA}&\cdots &p_{ACTT}\\p_{AGAA}&p_{AGAC}&p_{AGAG}&p_{AGAT}&p_{AGCA}&\cdots &p_{AGTT}\\p_{ATAA}&p_{ATAC}&p_{ATAG}&p_{ATAT}&p_{ATCA}&\cdots &p_{ATTT}\\p_{CAAA}&p_{CAAC}&p_{CAAG}&p_{CAAT}&p_{CACA}&\cdots &p_{CATT}\\\vdots &\vdots &\vdots &\vdots &\vdots &\ddots &\vdots \\p_{TTAA}&p_{TTAC}&p_{TTAG}&p_{TTAT}&p_{TTCA}&\cdots &p_{TTTT}\\\end{bmatrix}}$
The other possible trees can be constructed by permuting the taxa; for example, in the case of T2 (A,C,(B,D)); (i.e., the bipartition AC|BD) the rows would be defined by taxa A and C and the columns would be defined by taxa B and D.
This approach was pioneered by Eriksson, where the SVD of flattening matrices was proposed along with a novel tree construction algorithm. However, performance of the original Eriksson method (ErikSVD) was mixed when it was compared to neighbor-joining and the maximum likelihood approach implemented in the PHYLIP program dnaml. ErikSVD appeared to perform better than dnaml when applied to an empirical mammalian dataset based on an early release of data from the ENCODE project but it underperformed the other phylogenetic methods when it was used with simulated data. Fernández-Sánchez and Casanellas, proposed a normalization (Erik+2) that improved the original ErikSVD method. ErikSVD is statistically consistent given the general Markov model for sequence evolution (it converges on. the true tree. as the empirical distribution approaches the theoretical distribution) but the Erik+2 normalization exhibited improved performance given finite datasets. ErikSVD and Erik+2 has been implemented in the software package PAUP* as part of the SVDquartets method.
Scores for the alternative topologies are calculated using the singular values, as shown below:
$\delta _{n}(T_{1})={\sqrt {\sum _{i=n+1}^{N}{\widehat {\sigma _{i}}}}}$
where ${\widehat {\sigma _{i}}}$ is the $i^{th}$ singular of the flattening matrix for the appropriate topology. The $\delta _{n}$ values provide information about the rank of the flattening matrix; if the sequences were generated on a single tree topology then $\delta _{4}$ will be zero in expectation when the flattening matrix corresponds to the correct topology (i.e., the rank of the flattening matrix with be 4 at most). When the sequences were generated on a mixture of trees related by the multispecies coalescent we expect the rank of the flattening matrix for the correct tree to be no more than 10. This is the basis for the popular SVDquartets method.