Identify lowest normalised HMM likelihood value given MMR

[szhan]

Currently, mu (i.e. mutation rate) is arbitrarily set to 1e-3, but by setting it to a higher value, we could do better in terms of computational speed.

We need to do some math to figure out the relationship between a lower limit for HMM relative likelihood and precision value, given constraints on HMM input parameters, mu (which should be as high as possible) and MMR.

[szhan]

Relevant to this discussion are:

1. The function in sc2ts solve_num_mismatches.
2. This code snippet in tsinfer for updating HMM likelihood here.

[szhan]

Here's a first stab at thinking about this issue.

First, some notations.
k is the user-specified value of MMR (recall that higher values favour mutation over recombination)
n is the number of matchable nodes in a tree sequence
mu is probability of mutation to one of three other nucleotides
rho is probability of recombination
p_t is probability of HMM transition
p_e is probability of HMM emission
L(u) is the likelihood of node u at the current site.
L_prev(u) is the likelihood of node u at the previous site
p is the precision value

Given:

• a user-specified value of k
• a preset value of mu, and
• a fixed value n (because we are matching a new sequence to an existing ARG containing n matchable nodes),
  we compute rho as follows as done here:

rho = [n * mu ** k] / [(1 - mu ** k) + (n - 1) * mu ** k]

The HMM probabilities of transition and emission are defined as in tsinfer as follows:
p_t = p_recomb if p_recomb > p_no_recomb else p_no_recomb
p_e = mu if mismatch else 1 - 3 * mu (assuming that there are four possible nucleotides)

where
p_recomb = rho / n, and
p_no_recomb = L_prev(u) * [1 - rho + rho / n]

Once the above are defined, then we define L(u) = p_t * p_e.

After L(u) is computed for all nodes, L(u) is normalised by the maximum likelihood value L_max across all the nodes, and then L(u) is replaced with L(u) / L_max, which is rounded to the nearest 10 ** p.

We define the lowest normalised likelihood value of node u at the current site under the following scenarios.

1. Mutation without recombination
   There exist sequences in the existing ARG that are exact clones of the query sequence, i.e. requiring no recombination or mutation to explain. There also exist sequences that have one mismatch (the same mismatch) compared to the query sequence. The current site is where the mismatch occurs at node u.

The normalised likelihood of node u at the current site is M = mu / (1 - 3 * mu). Note that no mismatch has occurred until this site, so that L_prev(u) is 1).

By extension, if there exist sequences that have x mismatches (the same x mismatches) compared to the query sequence, then the normalised likelihood of node u at the current site (having emitted x - 1 mismatches already) is equal to M ** x.

2. Recombination without mutation
   There exist sequences that are clones of the query sequence. There also exist sequences that can be explained by a single HMM transition (all switching at the same site) compared to the query sequence. The current site is where the transition occurs at node u.

The normalised likelihood of node u at the current site is R = [rho / n] / [1 - rho + rho / n]. Note that no switching has occurred until this site.

By extension, if there exist sequences that have x transitions (the same x switches) compared to the query sequence, then the normalised likelihood of node u at the current site (having switched x - 1 times already) is equal to R ** x.

[szhan]

Generalising the above a bit…

1. Mutation without recombination
   If there exist only sequences that differ from the query sequences by x mismatches and sequences that differ by y mismatches, and x < y, then the lowest normalised likelihood at node u at the current site where the yth mismatch occurs is M ** (y - x).

2. Recombination without mutation
   Analogous to above.

[jeromekelleher]

Sounds right to me

[szhan]

Here is a simple illustration for the scenario with mutation but without recombination.

def compute_rho(mu, k, n):
    rho = n * mu ** k
    rho /= (1 - mu) ** k + (n - 1) * mu ** k
    return rho


def compute_M(mu):
    # Assume four possible states
    return mu / (1 - 3 * mu)


def compute_R(rho, n):
    R = rho / n
    R /= 1 - rho + rho / n
    return R


mu = 1e-3 # Currently used in sc2ts when k > 0
k = 3 # Currently used in sc2ts
n = 10

for i in range(5):
    print(i, compute_M(mu) ** i)

num. mismatch	relative likelihood
0	1.0
1	0.0010030090270812437
2	1.0060271084064632e-06
3	1.0090542712201234e-09
4	1.0120905428486695e-12

To differentiate cases where mismatch is equal to 0 versus 1 or higher, precision value of 0 is needed, since 10 ** 0 == 1. But to differentiate cases where mismatch is equal to 0 versus 1 versus 2 or higher, precision value of 3 is needed, which results in three possible rounded likelihood values, i.e. 1.0, 0.001, and 0.0.

By increasing mu to 1e-2, we need a lower precision value (i.e. 2) to differentiate cases where mismatch is equal to 0 versus 1 versus 2 or higher.

num. mismatch	relative likelihood
0	1.0
1	0.010309278350515464
2	0.00010628122010840684
3	1.0956826815299675e-06
4	1.1295697747731622e-08

[szhan]

Here is a similar exercise for the scenario with recombination but without mutation

mu = 1e-2
k = 3
n = 10

rho = compute_rho(mu, k, n)

for i in range(5):
    print(i, compute_R(rho, n) ** i)

num. switches	relative likelihood
0	1.0
1	1.0306101521283648e-06
2	1.0621572856700512e-12
3	1.0946700817686625e-18
4	1.1281780995019707e-24

We need a precision value of 0 to differentiate cases where the number of switches is equal to 0 versus 1 or higher; and 6 to differentiate cases where the number of switches is equal to 0 versus 1 versus 2 or higher.

Note that as mu gets lower and/or k gets higher, a higher precision value is needed. The likelihood values, however, are less sensitive to increases in n from 10 to 1e6.

[szhan]

Next, let's consider the scenario with both recombination and mutation.

The lowest likelihood given i mismatches and j switches is equal to (M ** i) * (R ** j).

mu = 1e-2
k = 3
n = 1e5

rho = compute_rho(mu, k, n)

for i in range(3):
    for j in range(3):
        print(i, j, compute_M(mu) ** i * compute_R(rho, n) ** j)

num .mismatches	num. switches	relative likelihood
0	0	1.0
0	1	1.0306101521283648e-06
0	2	1.0621572856700512e-12
1	0	0.010309278350515464
1	1	1.06248469291584e-08
1	2	1.0950075110000528e-14
2	0	0.00010628122010840684
2	1	1.095345044243134e-10
2	2	1.1288737226804667e-16

[jeromekelleher]

Nice. What happens if we set mu=0.5 in these calculations? 1e-2 is arbitrary, and we want to make the magnitudes of the likelihoods as large as possible, while keeping everything as a sensible probability. Ideally we should solve for mu with these constraints, given k.

[szhan]

I am assuming four possible nucleotides above, so max value for mu is 0.25.

For these values:

mu = 0.25
k = 3
n = 1e1

I get

num .mismatches	num. switches	relative likelihood
0	0	1.0
0	1	0.03703703703703704
0	2	0.0013717421124828536
1	0	1.0
1	1	0.03703703703703704
1	2	0.0013717421124828536
2	0	1.0
2	1	0.03703703703703704
2	2	0.0013717421124828536

[szhan]

Assuming two possible nucleotides, and setting

mu = 0.5
k = 3
n = 1e1

I get this pathological case where the relative likelihood is always 1.0.

[szhan]

We do want to relate precision to HMM cost instead of number of mismatches and switches, right? Which of course should be easy, given penalty weights for each type of event.

[jeromekelleher]

Ultimately we want to do something like this:

exact_matches = run_tsinfer_hmm(samples, likelihood_threshold=1)
one_mismatch = run_tsinfer_hmm(samples - exact_matches, 
    likelihood_threshold=one_mismatch_threshold(num_mismatches)
... etc

We don't have likelihood_threshold at the moment, we just have precision, which almost is the same thing except we round rather than truncate, and we only support power-of-ten values.

[szhan]

Okay, just noting that I don't think it makes sense when num_alleles * mu >= 1.0, like the two examples above, where num_alleles is the number of distinct alleles at a site.

[szhan]

Something like this might do, assuming no recombination.

def compute_M(mu, num_alleles):
    assert num_alleles * mu < 1.0
    return mu / (1 - (num_alleles - 1) * mu)


mu = 0.25 - 1e-1
k = 3
n = 10

for i in range(5):
    print(i, compute_M(mu, num_alleles=4) ** i)

num .mismatches	relative likelihood
0	1.0
1	0.2727272727272727
2	0.07438016528925619
3	0.020285499624342593
4	0.005532408988457071