Inconsistencies in Ka/Ks estimation ?

[ga01]

Hi,
I am trying to use HyPhy to make pairwise Ka/Ks estimates for sequences from very large, genome-wide phylogenies, (including ancestral sequences), and I get weird behaviour from HyPhy; I wonder whether this is a bug. The sequences are very long, 400k codons, but also have very low divergence (~99% sequence similarity).

When I compare SeqA and SeqB from the phylogeny (with FitMG94), I get the global Ka/Ks for the sequence pair ~ 0.16. Setting the --rooted flag to yes is also fine, results in similar estimates:
| Branch | Length | dN/dS | Approximate dN/dS CI|
| SeqA | 0.002 | 0.161 | 0.140 - 0.182 |
| SeqB | 0.002 | 0.154 | 0.138 - 0.170 |

However, when I compare SeqA with an ancestral sequence (that is the last common ancestor of SeqA and SeqB) I get weird behaviour.
With the --rooted flag set to Yes, there is a weird asymmetry (which is reversible, the one having 0.0 Ka/Ks is always the second sequence in the alignment):
| Branch | Length | dN/dS |Approximate dN/dS CI|
| Anc | 0.002 | 0.315 | 0.298 - 0.333 |
| SeqA | 0.001 | 0.000 | 0.000 - 0.019 |

while with --type flag set to global Ka/Ks remains ~ 0.16 (that is also surprising). The ancestral sequence was reconstructed with IQtree2, with the MGK method.
I wonder whether this is an error (due to the very long sequences?), or normal behaviour.

Best wishes,
G.

[spond]

Dear @ga01,

At first glance this does NOT look like a normal behavior. Would you be able to send me a test case where you see this happen? I'll try to diagnose the issue and get back to you with a solution.

Best,
Sergei

[spond]

Dear @ga01,

Thanks for sending me the data! I should have realized that the issue is sooner. If you use a time-reversible model, like MG94, the placement of the root is not generally resolvable. For two sequences that means that all the model is doing is (loosely speaking) peforming an arbitrary, generally unstable, partitioning of the total distance between the sequences into (D) into two parts that add to D. Rooting makes sense for non-reversible models; I honestly don't recall why I even implemented the --rooted Yes option for this analysis.

One way to check that the model is behaving sensibly is to compare log likelihoods for the global and local models.

$hyphy ~/Development/hyphy-analyses/FitMG94/FitMG94.bf --alignment Anc_SeqA.fa --type global --rooted No --tree tree.nwk ENV="USE_MEMORY_SAVING_DATA_STRUCTURES=1e8;"

### Fitting Standard MG94
* Log(L) = -1653593.59, AIC-c = 3307219.18 (16 estimated parameters)
* non-synonymous/synonymous rate ratio =   0.1590 (95% profile CI   0.1501-  0.1679)

...

(SeqA:0.003555776006615123,Node2:-1)

$hyphy ~/Development/hyphy-analyses/FitMG94/FitMG94.bf --alignment Anc_SeqA.fa --type global --rooted Yes --tree tree.nwk ENV="USE_MEMORY_SAVING_DATA_STRUCTURES=1e8;"

### Fitting Standard MG94
* Log(L) = -1653593.59, AIC-c = 3307221.18 (17 estimated parameters)
* non-synonymous/synonymous rate ratio =   0.1589 (95% profile CI   0.1500-  0.1679)

**Combined tree** 
(Anc:0.002085486861487382,SeqA:0.001470647622587388)

Note that

(1) The log-likelihoods are the same (even those the second model has 1 more parameter)
(2). The sum of both branch lengths in the second case (0.00355613) is effecively the same as the single branch length in the first case (0.003555776006615123)

The local model is a little different. Because the substitution matrices are not the same (different dN/dS), but individually reversible, there's a bit more subtletly.

If you run --type local --rooted No, you get exactly the same as --type global --rooted No

Log(L) = -1653593.59

(SeqA:0.001007444663465155,Node2:-1)

If you run --type local --rooted Yes, you get something a bit different.

### Fitting Standard MG94
* Log(L) = -1653592.76, AIC-c = 3307221.53 (18 estimated parameters)

### Estimating confidence intervals for dN/dS along each branch

|                      Branch                      |     Length     |     dN/dS      |Approximate dN/dS CI|
|:------------------------------------------------:|:--------------:|:--------------:|:------------------:|
|                       Anc                        |     0.002      |     0.271      |   0.256 - 0.287    |
|                       SeqA                       |     0.001      |     0.000      |   0.000 - 0.021    |

**Combined tree** 
(Anc:0.002498548005453272,SeqA:0.001057813727628272)

Reversing the order of sequences you can see the instability.

### Fitting Standard MG94
* Log(L) = -1653592.64, AIC-c = 3307221.27 (18 estimated parameters)

### Estimating confidence intervals for dN/dS along each branch

|                      Branch                      |     Length     |     dN/dS      |Approximate dN/dS CI|
|:------------------------------------------------:|:--------------:|:--------------:|:------------------:|
|                       SeqA                       |     0.002      |     0.022      |   0.008 - 0.037    |
|                       Anc                        |     0.002      |     0.382      |   0.359 - 0.406    |

**Combined tree** 
(SeqA:0.001670408809131221,Anc:0.001885606183000924)

Hope that makes sense. Can you tell me a bit more about WHY you are trying to do this? Are you trying to partition the variation between two sequences? In HyPhy, one way to do it (going with the Ancestral sequence approach), is to force evolution to go through the ansestral sequence, by creating an artifical 0-length branch. I can tell you how to do it if you are interested.

Best,
Sergei

PS You can add ENV="USE_MEMORY_SAVING_DATA_STRUCTURES=1e8;" to the command line as an argument so you don't need to modify the files themselves.

[ga01]

Hi Sergei,
I simply noticed that using an ancestral-extant sequence pair results in a very different behaviour than two extant sequences, so I was wondering whether something else is also affected. Apparently, it somehow affects the assignment of the root (and the downstream estimates). I wonder whether the fact that the ancestral sequence was also estimated with the MG method might result in this.

G.

[spond]

Dear @ga01,

Two long low-divergence sequences are kind of an edge-case. The behavior without specifying the root is exactly as expected. Root specification under the local model gives you some weirdness because of numerical / statistical instability since the code somewhat artbiratily "partitions" the (A-B) branch into two sub-branches using reversible models.

Best,
Sergei

[ga01]

Hi Sergei,
Thank you. Would you expect such instabilities to affect "partitioned" or "lineage" analyses (with the --type flag), or one should use them only with smaller, more diverged sequences? (Aside from the time they may need.) In an isolated, pairwise case one might be able spot a problem, but in a complex analysis that might be difficult.

Best wishes,
G.