Point substitution models

A subset of the class of phylo-grammars is the class of homogeneous substitution models, where the mutation rate is not a function of position but rather is identical for every site. Such models can be represented as a single-state phylo-HMM. Examples include

The Jukes-Cantor model [41], Kimura's two-parameter model [44], the HKY85 model [32], the general reversible model [92], and the general irreversible model [91]. In the case of the Kimura and HKY85 models, the rate matrices are formulated para-metrically: that is, each substitution rate is expressed as a function of a small set of rate and/or probability parameters (e.g. in Kimura's model, there are two rate parameters: the transition rate and the transversion rate).

Variable-rate models, where the evolutionary rate is allowed to vary from site to site [90]. Yang used a finite number of discrete, fixed rate categories to approximate a continuous gamma distribution over site-specific rates. In essence, this can be viewed as special cases of the phylo-HMM of Felsenstein and Churchill [23], with the autocorrelation explicitly set to zero.

Hidden-state models [48, 38]. A relative of the variable-rate model, the hidden-state model allows a variety of different substitution rate matrices to be used, depending on a hidden state variable that specifies the structural context of the site [48]. For example, a hydrophobically-inclined rate matrix might be used for buried amino acids and a hydrophilic matrix for exposed amino acids. An extension to the hidden-state model allows the hidden state variable itself to change over time at some slow rate, modeling rare changes in structural context [38]. An alternative extension allows correlations between hidden state variables at adjacent sites: this is essentially the idea behind the phylo-HMM, described below.

Models for synonymous/nonsynonymous substitution ratio measurement; empirical rate matrices for codon evolution [27, 87]. Codon substitution matrices such as WAG [87] can be used to measure the ratio r of synonymous to nonsynonymous substitution rates, which may be indicative of purifying (r < 1), neutral (r = 1) or diversifying (r > 1) selection. These models are also related to the exon prediction phylo-HMMs in EVOGENE [68] and EXONIPHY [80], described below.

Amino acid substitution models [12, 28]. Likelihood calculations using these models can, as with the other substitution models discussed above, be viewed as trivial applications of phylo-grammars.

Context-sensitive substitution models [81]. Siepel and Haussler introduced several alternate approximations for calculating the likelihood of alignments assuming a nearest neighbor substitution model, suitable for capturing the context-sensitivity of the substitution process that is observed in real sequence alignments (most notoriously in genomes wherein CpG methylation is used as a mechanism of epigenetic regulation, leading to elevated rates for the mutations CpG→TpG and CpG→ApG). Siepel and Haussler's method ignores longer-range correlations induced by nearest-neighbor effects, but is effective in practice. (It may be regarded as an approximation to the more rigorous analysis of Lunter and Hein [55].)

Many of these models can be expressed using the General Parametric Substitution Model, which we define as the substitution model wherein all substitution rates and initial probabilities can be expressed as simple functions of a (reduced) set of rate and probability parameters. As an example, Kimura's two-parameter model [44] is shown (see figure 1) along with the HKY85 six-parameter model [32] (see figure 2).

As long as each parameter in a parametric substitution model can be interpreted either as a rate (such as Kimura's transition and transversion rates) or a probability (such as the HKY85 equilibrium distribution over nucleotides), the phylo-EM algorithm can be adapted to estimate such parameters via the computation of expected event counts. A formal description of the sets of allowable rate and probability functions is given in the Supplementary Material [see Additional file 1].

Although the particular models used above (Kimura and HKY85) are reversible, matrices of allowable rate functions can in general be irreversible. Our General Parametric Model may thus be regarded as a generalisation of the General Irreversible Model.

Phylo-HMMs

Phylo-HMMs form a class of models slightly more complex than point substitution models. In a phylo-HMM, each column (or group of adjacent columns) is associated with a hidden state, representing the evolutionary context of the site. Each hidden state is conditionally dependent upon the immediately preceding state (the Markov property).

Tasks that have been addressed using phylo-HMMs include:

Measurement of variation of evolutionary rate among sites in DNA [23]. Felsenstein and Churchill construct an HMM with three states. Each state generates an alignment column according to a point substitution process on a tree [21]. The overall evolutionary rate for the column depends on the state from which it is emitted: each state thus corresponds to a "rate category" (the relative rates for the three states are 0.3, 2.0 and 10.0). The use of an HMM allows for an autocorrelated model of rate variation.

Modeling site-specific residue usage in proteins[9, 31]. While site-specific profiles are familiar tools in bioinformatics, early tools such as Gribskov profiles [29] and hidden Markov models [8] ignored phylogenetic correlations in the dataset, leading to biased sampling. Phylo-grammars incorporate these correlations directly. In these papers, Bruno et al. introduced an initial EM algorithm for estimating rate matrices.

Prediction of secondary structure in proteins[83, 26]. In a similar manner to Felsenstein and Churchill, a three-state HMM is constructed wherein each state emits an alignment column using a substitution rate matrix. Here, however, the states correspond to different units of secondary structure (loop, α-helix and β-sheet). The substitution rate matrix for each state reflects the frequency distribution and substitution patterns for that secondary structural class. The method performs less well than established secondary structure prediction algorithms, but shows promise, in particular given the simplicity of the model (three states only). Later work expanded the number of states in the phylo-HMM to eight (correspondingly increasing the number of parameters). Note that, as more parameters are introduced into this kind of phylo-HMM, the problem of "training" those parameters grows in importance.

Prediction of exons and protein-coding gene structures in DNA [68, 80]. The basis for the gene prediction programs EVOGENE and EXONIPHY, respectively, these phylo-HMMs are based on substitution models for codon triplets with 4³ = 64 states. The paper by Siepel and Haussler introduced the term "phylo-HMM" and used an approximate version of the EM algorithm introduced by Holmes and Rubin for parameterization [38].

Detection, modeling and annotation of transcription factor binding sites in DNA [62]. Here, the EM algorithm and other formulae of Bruno and Halpern [9, 31] is used to model site-specific residue frequencies in alignments of promoter regions (rather than proteins, as addressed by Bruno and Halpern).

Detection of conserved regions in multiple alignments of genomic DNA [79]. Phylo-HMMs to detect conserved regions can be viewed as extensions of Felsenstein and Churchill's original model with more rate categories. This approach has been used to detect highly-conserved regions in vertebrate, insect, nematode and yeast genomes. Approaches measuring the substitution rate per site [79, 85], the local indel rate [54] and/or the CpG mutation bias [81, 55] have all shown merit.

Analogously to some of the point substitution models, many phylo-HMMs can be expressed parametrically. An example of such a model is the one used by Siepel's PHASTCONS program, whose phylo-HMM has ten states ranging from slow to fast overall substitution rate. Moving from one state to another, the relative substitution rates between different nucleotides do not change (i.e. the ratio R _ij /R _kl is constant for any i, j, k, l ∈ {A, C, G, T}); only the overall substitution rate varies (i.e. the absolute value R _ij is not constant). Such consistency across states can be achieved by writing the rate matrices for the ten states as k₁R, k₂ × R, k₃ × R... k₁₀ × R where the k _i are scalar multipliers and R is a relative rate matrix shared by all the states. Similarly, the rate matrices of Felsenstein and Churchill's three-state phylo-HMM can be written 0.3 × R, 2 × R and 10 × R. Both are examples of the general parametric phylo-HMM.

Phylo-SCFGs

The most complex class of phylo-grammar considered here is the phylo-SCFG. Most commonly used to model RNA secondary structure, these grammars are capable of modeling covariation between paired sites. In an SCFG, covarying sites must be strictly nested, allowing the modeling of foldback structures but not pseudoknots, kissing loops or other topologi-cally elaborate RNA structures [45].

Tasks that have been addressed using phylo-SCFGs include:

Prediction of RNA secondary structure [46, 47]. The Pfold program in this paper introduced the first phylo-SCFG, combining stochastic context-free grammars (used to model RNA structure) with evolutionary substitution models. Since HMMs are a subset of SCFGs, the framework of phylo-SCFGs includes the previously discussed phylo-HMMs. The Pfold program also allowed for user-specified grammars; however, it lacked a fast EM-like algorithm for estimating grammar parameters from data (by contrast, the non-phylogenetic SCFGs used elsewhere in bioinformatics can be rapidly trained using the Inside-Outside algorithm [16]). A key feature of these models is the use of 16-state "basepair models" for modeling the simultaneous coevolution of functional base-pairs in RNA structures. Again, fast and effective parameterization of the model is an important issue.

Detection of noncoding RNA genes [67]. A similar model to Pfold was used by the Evofold program, which uses a phylo-SCFG to parse genomic alignments into noncoding RNA and other features [67].

Detection of RNA secondary structure within exons [69]. The RNA-Decoder program uses a parametric phylo-SCFG to model exonic regions in which there is simultaneous selection on both the translated protein sequence and the secondary structure of the pre-mRNA. Such regions have been found in viral genomes and hypothesized to fulfil a regulatory role [69]. Due to the complexity of these models and the sparsity of training data, parametric rate functions are required to limit the number of free parameters that must be estimated.

Detection of accelerated selection in human noncoding RNA [70]. Pollard et al used phylo-HMMs and phylo-SCFGs to identify a neurally-expressed RNA gene, HARF1, that had undergone recent accelerated evolution in the lineage separating humans from the human-chimp ancestor.

Programs

The following open source software tools, implementing the algorithms and models described in this paper, are freely available (see Availability and Requirements).

xgram – a implementation of the EM algorithm for training phylo-grammars, i.e. the Inside-Outside and Forward-Backward algorithms combined with the EM algorithm for the general irreversible (and reversible) substitution models. This program implements the general irreversible EM algorithm described in the Supplementary Material [see Additional file 1], along with the general reversible EM algorithm described previously [38]. The grammar can be user-specified via an extensible file format, described below. Parametric grammars are allowed (so that individual substitution rates and/or rule probabilities can be constrained to arbitrary functions of a smaller set of model parameters). The xgram tool is capable of reproducing most of the phylo-grammar models listed in this paper. In its generic applicability, xgram is similar to the dynamic programming engine Dynamite [5], although the class of models is different (phylo-grammars vs single- and pair-HMMs) and the functionality broader (including parameterization by phylo-EM, as well as Viterbi and CYK annotation codes). Also included is an implementation of the neighbor-joining algorithm for fast estimation of tree topologies [77], and another version of the EM algorithm for rapidly optimising branch lengths of trees with fixed topology [24]. The model underlying xgram also allows for dynamically evolving "hidden states" associated with each site, again as previously described [38].

xrate – a version of xgram including several "preset" grammars for point substitution models, including the general irreversible and reversible substitution models.

xfold – a version of xgram including several "preset" grammars for RNA analysis, including that of the Pfold program [46].

xprot – a version of xgram including several "preset" grammars for protein analysis, including a grammar similar to that used by Thorne et al. for protein secondary structure prediction [83].

All of the above programs can be driven by any user-specified phylo-grammar. Having specified a grammar, or chosen one of the presets, the user can

The parse tree can also be constrained, completely or partially, by including complete or partial annotations in the input alignment. For example, one can annotate several known examples of a TF binding site in a multiple alignment. One can then allow the grammar to "learn" these examples and predict new binding sites.

File formats

The input and output format for sequence alignment data is the Stockholm format, as used by PFAM and RFAM. The wildcard character is the period ".". Annotation of columns with the wildcard character allows for incompletely labeled data and hence partially supervised learning. If a given annotation is specified in the grammar but absent from the training data, it will be treated as a string of wildcards and all compatible possibilities will be summed over.

Any phylo-grammar can be specified, using a format based on LISP S-expressions [56, 75]. The format is human-readable and succinct, while being machine-parseable and extensible.

Phylo-grammar specification files contain several elements:

As an example, the grammar for the Kimura two-parameter rate matrix is shown (see figure 3). A more complete and up-to-date description of the format can be found online [88], as can discussion of the latest version of xrate and its companion programs [89].

Fitting codon models

In the past, various amino acid substitution models have been estimated using ML techniques (e.g., mtREV [15], WAG [87]). An ML estimation of codon substitution models, however, has seemed infeasible for a long time because of the computational burden involved with such parameter-rich models. This section shows that xrate is capable of tackling the problem. The full results of a particular study are being published elsewhere (Kosiol, Holmes and Goldman, in prep.); here, we will restrict attention to simulation results showing that xrate can do these sorts of analyses reliably.

The number of independent parameters for a reversible substitution model with N character states can be calculated as $\frac{N(N+1)}{2}-2$. This means that for the estimation of a 20-state amino acid model, 208 independent parameters need to be calculated. In contrast, to estimate a 61-state codon model (excluding stop codons), 1889 independent parameters have to be determined.

To test the robustness of xrate's ability to fit parameter-rich models to aligned sequence data, we simulated a data set using all phylogenies of the Pandit database of protein domain alignments [86], using a standard model of codon evolution (the MO model [93] [see Additional file 1]). In this model, rates of substitutions involving changes to multiple nucleotides are zero, so that the rate matrix is sparsely populated.

xrate is able to recover M0 well from this 'artifical' Pandit database. The true rates used in the simulation are shown (see figure 8). These may be compared with the recovered rates (see figure 9).

A scatter plot of true vs estimated rates allows a more detailed analysis (see figure 10). This plot shows the true instantaneous rates $q_{ij}^{(true)}$ of M0 plotted versus the instantaneous rates $q_{ij}^{(est)}$ estimated from data simulated from M0. If $q_{ij}^{(true)}$ = $q_{ij}^{(est)}$ the points would lie on the bisection line y = x. Thus the deviation of the points from the bisection line indicates how different the rates are.

If one is interested in drawing biological conclusions from the estimated rate parameters, then it is of interest to consider xrate's estimates of rates which are zero in the true model, xrate sometimes inferred erroneously very small non-zero values for the instantaneous rates of double and triple changes from the simulated data set (in the M0 model, which was used to generate the data, such substitutions have zero rate). However, this error can be correctly identified by comparing log-likelihoods calculated by xrate under the following nested models: For the general model allowing for single, double and triple nucleotide changes 1889 parameters had to be estimated. The best likelihood calculated for general estimation is In L _general = -28930383.06. Using xrate we can also restrict the rate matrices to single nucleotide changes only. For this model 322 parameters had to be estimated. The best likelihood calculated for restricted estimation is lnL _restricted = -28930894.86.

Although the log-likelihood for the general rate matrix allowing for single, double and triple changes is better we can show that the improvement is not significant. Significance is tested using a standard likelihood ratio test between the two models, comparing twice the difference in log-likelihood with a ${\chi }_{1567}^2$ distribution, where 1567 is the degrees of freedom by which the two models differ. Using the normal approximation for ${\chi }_{(1567,0.01)}^2$ we compare (2(ln L _general - ln L _restricted )-1567)/$\sqrt{2\times 1567}$ = -9.71 with the relevant 99% critical value of 2.33 taken from a standard normal $\mathcal{N}$ (0,1). The difference is seen to be insignificant; the P-value is almost 1.

Predicting protein secondary structure

We compared xrate to the phylo-HMM for prediction of protein secondary structure developed by Goldman, Thorne, and Jones [26] (here referred to as GTJ). This section uses a fully-connected three-state phylo-HMM with general reversible Markov chains. Training sets were taken from the HOMSTRAD database of structural alignments of homologous protein families [61].

We trained the phylo-HMM on alpha-beta barrel alignments from HOMSTRAD, leaving out the beta-glycanase SCOP family. xrate was then benchmarked on this beta-glycanase SCOP family to compare the annotation predicted by xrate to the experimentally determined HOMSTRAD annotation. We also tried a more comprehensive training regime, training xrate on the complete HOMSTRAD database (excluding the beta-glycanase SCOP family) and again comparing predicted and database annotations.

The performance of xrate was compared to that of GTJ. The results show that xrate can be used to quickly prototype and train a phylo-HMM with comparable performance to that reported by Goldman et al.

Grammar

The PROT3 phylo-grammar has state labels for the three secondary structure classes of alpha-helix (H), beta-sheet (E) and loop (L). An excerpt of the grammar is shown (see figure 4).

An example of usage for this grammar follows. We also show an alignment from HOMSTRAD, too small to predict secondary structure with any confidence, but useful for illustrative purposes (see figure 5). Suppose we want to: (1) read in this alignment from a file named ' pp. stk'; (2) load a point substitution matrix from a file named 'dart/data/nullprot.eg' (this is an amino-acid matrix distributed with xrate; the filename path assumes that the DART package was downloaded to the current working directory); (3) use the above point substitution matrix to estimate a phylo-genetic tree (by neighbor-joining followed by EM on the branch lengths); (4) load the PROT3 model from a file named 'dart/data/prot3.eg' (again, this is distributed with xrate); and (5) use the PROT3 model to predict secondary structure classes for this protein family, printing the annotated alignment to the standard output. The following command-line syntax achieves this:

xrate pp.stk --tree dart/data/nullprot.eg --grammar dart/data/prot3.eg

The output of this command is shown (see figure 6).

More such examples can be found in DART (the software library with which xrate is distributed) and on the wiki pages for the xrate program [89]. A full list of command-line options for xrate can be obtained by typing xrate –help or, equivalently, xrate -h.

Results

Both xrate and the GTJ program were evaluated on the xylanase alignment used by GTJ, hereafter referred to as gtjxyl. xrate was trained on the subset of HOMSTRAD corresponding to alpha-beta barrel structures, with members of the beta-glycanase SCOP family (which includes the gtjxyl proteins) removed to prevent overlap between the training and test sets.

We report the prediction accuracy collectively for all secondary structure categories, and the sensitivity and specificity with respect to each individual category. These metrics are defined as follows

Sensitivity(n) = TP _n /(TP _n + FN _n )

Speciflcity(n) = TP _n /(TP _n + FP _n )

Accuracy = (${\displaystyle \sum _n{\text{TP}}_n}$)/(${\displaystyle \sum _n{\text{TP}}_n}$ + FN _n )

where (for secondary structure class n) TP _n is the number of true positives (columns correctly predicted as class n), FN _n is the number of false negatives (columns that should have been predicted as class n but were not) and FP _n is the number of false positives (columns that were incorrectly predicted as class n).

Bubbleplots were used to visualize the amino acid substitution rates. Substitutions are colored red if between aromatic amino acids, green if between hydrophobics and blue if between hydrophilics. Substitutions from one such group to another (e.g. from hydrophobic to hydrophilic) are colored gray.

Figures 11, 12 and 13 show the amino acid substitution matrices for the alpha-helix, beta-sheet and loop states, respectively. The relative rates displayed in the figures in general agree with what one would expect from each of those states: the alpha-helix and beta-sheet states substitute more slowly (and thus amino acid conservation is higher) than for the loop states (loop regions being more variable in structure [7]).

Table 1 shows the log likelihood scores of the training alignments, log P(D|θ), along with the log-posterior probability of the HOMSTRAD reference annotation, log P(A|D, θ). In this case, maximum-likelihood training also yields an increase in the annotation posterior probability P(A|D, θ). This is not in general a guaranteed result of the EM algorithm, and alternative training procedures (such as maximum-discrimination training [19]) have been proposed to achieve this effect. It appears in this case that such procedures are not required.

θ	log2 P(A, D|θ)	log2 P(D|θ)	log2 P(A|D, θ)
θ D	-173038	-162491	-10547
θ G	-238632	-227979	-10653

Table 2 reports likelihoods, accuracies and runtimes for training set 2 as the EM convergence criteria are tightened. As expected, the likelihood increases as the convergence criteria are made more stringent. The annotation accuracy for the gtjxyl benchmark alignment also consistently increases.

mininc	Runtime/min	Acc(gtjxyl)	log2 P(A, D|θ D )	log2 P(D|θ D )	log2 P(A|D, θ D )
le-3	14	64.1	-2696469	-2549947	-146522
le-4	35	64.7	-2686598	-2539908	-146690
le-5	84	68.0	-2682667	-2536849	-145818

Table 3 summarizes the results of running xrate and the GTJ program on all the test cases. In general the accuracy of xrate is comparable to or even slightly better than the accuracy of the GTJ program.

Program	Sn (α)	Sp (α)	Sn (β)	Sp (β)	Sn (L)	Sp (L)	Acc
GTJ	66.7	91.3	63.5	84.0	73.5	77.3	69.6
xrate	71.6	95.7	82.7	79.0	65.2	81.2	70.2

Predicting RNA secondary structure

To illustrate the capability of xrate as a tool for RNA secondary structure prediction/annotation, we compare it to Pfold, a phylo-SCFG developed by Knudsen and Hein [46, 47].

There are two goals of this section: (1) to see if xrate can exactly emulate the Pfold phylo-grammar using the same parameters as Pfold, and (2) to see if the EM algorithm can estimate parameters that yield comparable performance to those produced by other methods.

We benchmarked the Pfold phylo-SCFG running on xrate against the original Pfold program using alignments from the Rfam database [30]. To address goal (2), we used xrate to estimate the substitution rates and initial frequencies of basepairs and single nucleotides from annotated Rfam alignments.

Our results show that the Pfold phylo-SCFG is effectively emulated by xrate, that the EM algorithm can estimate a more likely parameterization for a given training set and that the parameters so obtained are comparable in performance to the Pfold program itself. We conclude that xrate is a suitable platform for developing, parameterizing, and testing phylo-grammars without the necessity of writing source code or performing manual parameterization.

Results

We report the sensitivity and positive predictive value (PPV) of basepair predictions. These accuracy metrics are defined as follows

Sensitivity = TP/(TP + FN)

PPV = TP/(TP + FP)

where TP is the number of true positives (base pairs that are predicted correctly per the Rfam annotation), FN the number of false negatives (base pairs that are not predicted but are in the Rfam annotation) and FP the false positives (predicted base pairs that are not in the Rfam annotation).

Training and testing sets were obtained by selecting the 148 RNA gene families in Rfam version 7 with experimentally-determined structures, discarding pseudoknots, removing excessively gappy columns (as this step is also performed by Pfold), grouping the families into superfamilies and randomly partitioning these superfamilies into two sets [see Additional file 1]. This yielded a training set of 71 alignments and a testing set of 77 alignments.

The benchmark results, shown in Table 4, indicate that the sensitivity and PPV of the Pfold program and its emulation on xrate are comparable. It should be noted, however, that the sets of base pairs predicted by the two programs are slightly different [see Additional file 1]. After examination, we attribute this to differences in implementation and loss of precision due to numerical calculations.

	Sensitivity	PPV
Pfold	45.0%	58.3%
xrate emulating Pfold	44.4%	61.7%
xrate trained on Rfam	42.8%	58.2%

We also tested whether parameterizing the phylo-SCFG using the EM algorithm is comparable to the Pfold parameterization [46]. A comparison of Pfold's original rates with the EM-estimated rates is shown in Figures 14–16. Both sets of parameters display similar trends. Substitutions that create or preserve canonical base pairs are more frequent than substitutions that destroy basepairs (see figure 14). Transitions are more common than transversions, both within basepairs (see figure 15) and unpaired sites (see figure 16). There is a difference in the magnitude of many of the rates, which we attribute to differences in the training sets.

The predictive accuracy of Pfold is compared to that of the xrate-trained phylo-SCFG in Table 4, while log-likelihoods are compared in Tables 5 and 6. The results are similar, indicating that the combination of training set and xrate-implemented EM is comparable to the training procedure used in the development of Pfold.

Dataset	"mininc"	"forgive"	log2 P(D, A|θ)	log2 P(D|θ)	log2 P(A|D, θ)
Training set	le-3	0	-466330.6649	-453589.9251	-12740.7398
Training set	le-4	0	-465397.0642	-453403.7081	-11993.3561
Training set	le-5	0	-465397.0642	-453403.7081	-11993.3561
Training set	le-3	2	-465821.5239	-453476.0389	-12345.4850
Training set	le-3	4	-465565.9224	-453437.5353	-12128.3871
Training set	le-3	6	-465397.0642	-453403.7081	-11993.3561
Training set	le-3	8	-465291.1983	-453356.6841	-11934.5142
Training set	le-4	4	-465147.9174	-453318.4543	-11829.4631
Training set	le-4	10	-465010.8431	-453209.0744	-11801.7687
Test set	le-3	0	-360472.7960	-343832.6014	-16640.1946
Test set	le-4	0	-360190.7940	-344117.5123	-16073.2817
Test set	le-5	0	-360190.7940	-344117.5123	-16073.2817
Test set	le-3	2	-360148.9090	-343841.2775	-16307.6315
Test set	le-3	4	-360178.4500	-344016.2558	-16162.1942
Test set	le-3	6	-360190.7940	-344117.5123	-16073.2817
Test set	le-3	8	-360092.2930	-344078.8868	-16013.4062
Test set	le-4	4	-360057.4880	-344116.5923	-15940.8957
Test set	le-4	10	-360108.0100	-344166.2108	-15941.7992

Dataset	log2 P(D, A|θ)	log2 P(D|θ)	log2 P(A|D, θ)
Training set	-487422.5964	-464828.9148	-22593.6816
Test set	-370490.5284	-348550.7516	-21939.7768

An important point to check is whether the EM algorithm actually performs as designed. We expect to see certain phenomena if the algorithm is indeed working as expected:

• The algorithm, over the course of its iterations, should refine the parameter set (denoted at the n'th iteration by θ⁽ⁿ⁾) to maximize the likelihood of the alignment data D and (if supplied) the annotation A. Therefore, the log-likelihood log P(D|θ⁽ⁿ⁾) should increase with n towards an asymptotic maximum value. This is indeed observed to be the case for this example (see figure 17).

• In practice, the EM algorithm is not run for an infinite number of iterations; rather, the algorithm stops when some "convergence criteria" are met (relating to the fractional increase of the log-likelihood) and the parameters at this point are considered to be the "convergent parameters". We denote this convergent parameter set by θ*.

• If the EM algorithm is performing effectively (i.e. finding a parameterization whose likelihood is close to the global maximum), we would also expect P(D|θ*) to be greater than P(D|θ') for some arbitrarily chosen parameterization θ' (for example, the Knudsen-Hein parameters, which were optimized for a dataset other than D). A comparison of Tables 5 and 6 confirms this to be the case.

• As the convergence criteria become more strict, log P(D|θ*) should increase. The results in Table 5 confirm this to be the case.

• If the training set is representative of the test set, then the above statements should also hold true when D is taken to mean the test set. Again, Tables 5 and 6 confirms this.

We note that Tables 5 and 6 shows that the posterior probability of the true annotation, P(A|D, θ) = P(A, D|θ)/P(A|θ), is also increased after phylo-EM training. As mentioned above, this is not a provably guaranteed result of the EM algorithm, which is designed to maximize only P(A, D|θ).