Question about the difference between the ECMK07_GY and ECMK07_MGK models

[shacrvr]

Hi! I was wondering if folks had insight on any difference between the implementation of the ECMK07_GY and ECMK07_MGK combined codon models? As I understand it the MGK and GY models have the same parameters, an omega and a kappa, and so I expected them to converge to the same value when combined with ECMK07. I've been running both models using a small test dataset of 7 codon sequences, 15 codons each, and the example command below. So far, I've gotten different estimates for omega_E and kappa. Is there a difference in how omega_E and kappa are calculated in ECMK07_MGK vs ECMK07_GY? Thanks!

iqtree2 -st CODON -s example.fasta -te example.treefile -blfix -m ECMK07_MGK+FQ -redo

Here are some results I've generated running different models with the test dataset.

ECMK07_MG+FU
Empirical nonsynonymous/synonymous ratio (omega_E): 0.090
Transition/transversion ratio (kappa): 1.000
BEST SCORE FOUND : -271.812

ECMK07_MGK+FU
Empirical nonsynonymous/synonymous ratio (omega_E): 0.091
Transition/transversion ratio (kappa): 1.151
BEST SCORE FOUND : -271.750

ECMK07_GY+FU
Empirical nonsynonymous/synonymous ratio (omega_E): 0.095
Transition/transversion ratio (kappa): 1.191
BEST SCORE FOUND : -270.284

ECMK07_MGK'{1.0,1.0}'+FU
Empirical nonsynonymous/synonymous ratio (omega_E): 0.084
Transition/transversion ratio (kappa): 1.000
BEST SCORE FOUND : -271.841

ECMK07_GY'{1.0,1.0}'+FU
Empirical nonsynonymous/synonymous ratio (omega_E): 0.192
Transition/transversion ratio (kappa): 1.000
BEST SCORE FOUND : -272.887

MGK+FQ
Nonsynonymous/synonymous ratio (omega): 0.088
Transition/transversion ratio (kappa): 1.316
BEST SCORE FOUND: -271.964

GY+FQ
Nonsynonymous/synonymous ratio (omega): 0.088
Transition/transversion ratio (kappa): 1.316
BEST SCORE FOUND: -271.964

ECMK07_MGK+FQ
Empirical nonsynonymous/synonymous ratio (omega_E): 0.056
Transition/transversion ratio (kappa): 1.212
BEST SCORE FOUND: -267.812

ECMK07_GY+FQ

Empirical nonsynonymous/synonymous ratio (omega_E): 0.051
Transition/transversion ratio (kappa): 1.604
BEST SCORE FOUND: -268.545

[dltheobald]

Hi all,

There are three main things to notice in this output.

1. MGK and GY perform as the same model on their own (which is how I think they should behave). They provide the same logL and the same ML estimates for both kappa and omega.

2. However, when they are combined with empirical models (like ECMK07), they perform differently, giving different logL and ML estimates of kappa and omega.

3. If we set omega and kappa to fixed values of 1.0, the mechanistic models should reduce to simply multiplying by the identity matrix and should not affect the empirical model in a combined analysis. Yet ECMK07_MGK'{1.0,1.0}' and ECMK07_GY'{1.0,1.0}' give different results (different logL).

Douglas

[bqminh]

The answer is a bit involved. There is a fundamental difference between the "mechanistic" MGK and GY-type codon models in the way they treat the codon frequencies: the entries of the rate matrix Q in MG is the product of (omega or 1) and (kappa or 1) and the target nucleotide frequency, whereas the GY models replace target nucleotide by the target codon. For example, under MGK model:

q(AAC -> AAG) = omega * kappa * pi(G) # for non-synonymous (Asn->Lyn) and transition (C->G)
q(AAC -> AAT) = pi(T) # synonymous and transversion

However, the same entry for GY model is:

q(AAC -> AAG) = omega * kappa * pi(AAG) # for non-synonymous (Asn->Lyn) and transition (C->G)
q(AAC -> AAT) = pi(AAT) # synonymous and transversion

It was proven (Muse and Gaut 1994) that under the MG model, the equilibrium frequency of a codon will be the product of the three nucleotide frequency normalised by non-stop codons. For example, pi_MG(AAG) ~ pi(A)*pi(A)*pi(G). Such frequency model is denoted by +F1x4, which has only 4 nucleotide frequency parameter, versus 61 codon frequency parameter under the original GY model. Thus, MG model is much less parameterised than the GY+F model. There is a variant of of frequency model called +F3x4, in which you count the nucleotide frequencies for the 3 codon position separately and then multiply them: pi_MG+F3x4(AAG) ~ pi1(A)*pi2(A)*pi3(G), where pi1(A) is the frequency of A at the 1st codon position, pi2(A) is the A frequency at 2nd codon position, and pi3(G) is the G frequency at the 3rd codon position. Thus, F3x4 has 12 frequency parameters, which are still much less than GY+F model. At some point, Ziheng introduced GY+F1x4 and GY+F3x4 to "copy" this mechanism to reduce the parametrisation, i.e.: under GY+F1x4:

q(AAC -> AAG) = omega * kappa * pi(A) * pi(A) * pi(G) # for non-synonymous (Asn->Lyn) and transition (C->G)
q(AAC -> AAT) = pi(A) * pi(A) * pi(T) # synonymous and transversion

I believe this explains what you observed for the "empirical-mechanistic" codon model EMCK07_MGK and ECMK07_GY, where the entries of Q are just the product of the above and the empirical exchangeability rates.

[shacrvr]

Thanks!