Likelihood

This chapter describes the likelihood formula used by aphid to distinguish between GF and ILS. The formula is specified from the bottom-up, i.e. constructed from raw inputs to its final value.

The focal species triplet is denoted as \(% Factorize common mathjax commands. \newcommand{gt}{gt} \newcommand{pab}{{p_{ab}}} \newcommand{pac}{{p_{ac}}} \newcommand{pbc}{{p_{bc}}} \newcommand{po}{{p_{o}}} \newcommand{gt}[1]{{gt_{#1}}} \newcommand{sl}{{sl}} % \newcommand{So}{{S_\varnothing}} \newcommand{SI}[1]{{S_{ILS(#1)}}} \newcommand{SG}[2]{{S_{GF(#1,\ #2)}}} % \newcommand{pgf}{{p_{GF}}} \newcommand{pils}{{p_{ILS}}} % \newcommand{ABC}{{\scriptsize \text{ABC}}} \newcommand{ACB}{{\scriptsize \text{ACB}}} \newcommand{BCA}{{\scriptsize \text{BCA}}} \newcommand{CAB}{{\scriptsize \text{CAB}}} \newcommand{CBA}{{\scriptsize \text{CBA}}} % % Expected branches lengths. \newcommand{exb}[3]{{\bar{#1}({#2},\ {#3})}} \newcommand{exu}[1]{{\bar{u}({#1})}} \newcommand{exv}[1]{{\bar{v}({#1})}}((A,\ B),\ C)\) according to the phylogeny.

Inputs

Calculating the likelihood value requires the following inputs (leaves of the formula).

Global values

At the global level:

• \(\gt = (\gt{1},\ \gt{2},\ …)\) : list of possible GF events dates (read from config).
• \(𝜃\) : scaled effective population size (to be learned during the process).
• \(𝜏_1\) : latest divergence time in the focal triplet (idem).
• \(𝜏_2\) : earliest divergence time in the focal triplet (idem).
• \(\pab\) : probability that GF occured between \(A\) and \(B\) (idem).
• \(\pac\) : probability that GF occured between \(A\) and \(C\) (idem).
• \(\pbc\) : probability that GF occured between \(B\) and \(C\) (idem).
• \(\po\) : knowing that GF occured, probability it occured at the oldest date in \(gt\) (idem).

Local (per gene tree)

For every gene tree \(g\):

• \(\sl_g\) : sequence length, read from taxa.trees input file.

• \(𝛼_g\) : relative mutation rate, estimated during preprocessing.

• \(T_g\) : focal triplet topology, read from taxa.trees input file. One among:

  • \(T_g = \ABC\): concordant triplet: \(((A,\ B),\ C)\) is observed in the gene tree.
  • \(T_g = \ACB\): discordant: the observed outer node is \(B\) instead of C: \(((A,\ C),\ B)\).
  • \(T_g = \BCA\): discordant: the observed outer node is \(A\) instead of C: \(((B,\ C),\ A)\).

• Focal triplet branches lengths:

  • \(a_g\) : length of the branch yielding species \(A\) (regardless of \(T_g\)).
  • \(b_g\) : length of the branch yielding species \(B\) (idem).
  • \(c_g\) : length of the branch yielding species \(C\) (idem).
  • \(d_g\) : length of the internal branch in the triplet (idem).

  When useful, we refer to either branch length as \(l_g\), and either branch count as \(n_g\). Every count is the rounded product of the corresponding branch length read from the taxa.trees input file and the corresponding sequence length:

  \[ n_g(l) = \lfloor sl_g × l \rceil \]

Constraints

The formula values are subject to the following constraints:

• No negative values.
• \(\sl_g\) integer (number of sites)
• \(𝜏_1 \leq 𝜏_2\) (older coalescence last)
• \(p_* \leq 1\) (probabilities)
• \(\pab + \pac + \pbc \leq 1\) (total probability of GF)
• \(i < j \implies gt_i > \gt{j} \) (GF dates sorted decreasingly)

Formula

The likelihood value integrates over every observed gene tree \(g\) and every considered evolution scenario \(s\).

Evolution scenarios

The set \(S\) of scenarios considered is constituted by the following:

• One "no event" scenario \(\So\).
• One ILS scenario for every possible ancestral triplet topology:
  • \(\SI{\ABC}\)
  • \(\SI{\ACB}\)
  • \(\SI{\BCA}\)
• Five GF scenarios for every considered event date \(gt_i\):
  • \(\SG{\gt{i}}{\ABC}\) : transfer from \(A\) to \(B\) or \(B\) to \(A\) (undistinguishable).
  • \(\SG{\gt{i}}{\ACB}\) : transfer from \(A\) to \(C\).
  • \(\SG{\gt{i}}{\CAB}\) : transfer from \(C\) to \(A\).
  • \(\SG{\gt{i}}{\BCA}\) : transfer from \(B\) to \(C\).
  • \(\SG{\gt{i}}{\CBA}\) : transfer from \(C\) to \(B\).

The total number of considered scenarios is therefore: \[ |S| = 1 + 3 + 5 × |gt| \]

The prior probability that either GF scenario occured is: \[ \pgf = \pab + \pac + \pbc \]

Now, given one specific GF scenario: \[ ℙ(\SG{\gt{i}}{t}) = \left\{\begin{array}{ll} i = 1 &\implies \po \\ i > 1 &\implies \frac{1 - \po}{|gt| - 1} \end{array}\right. × \left\{\begin{array}{ll} t = \ABC &\implies p_{ab} \\ t \in \{\ACB, \CAB\}&\implies \frac{1}{2}p_{ac} \\ t \in \{\BCA, \CBA\}&\implies \frac{1}{2}p_{cb} \end{array}\right. \] Note that this prior formula gives a special occurence probability \(\po\) to the oldest considered potential GF event. The other events are consider equiprobable.

The prior probability that ILS occured is: \[ \pils = e^{-\frac{2(𝜏_2 - 𝜏_1)}{\theta}} \]

Now, given one particular ILS scenario: \[ \ ℙ(\SI{·}) = \frac{1}{3} × \pils × (1 - \pgf) \]

Finally, the "no event" scenario occurs when no other scenario does: \[ ℙ(S_{\varnothing}) = (1 - p_{ILS}) × (1 - p_{GF}) \]

Expected branches lengths

Calculating the likelihood essentially consists in comparing, for every gene tree \(g\):

• \(n_g(l_g)\) : the branches counts actually calculated from the taxa.trees input file.
• \(\exb{n}{s}{T_g}\) : the branches counts expected under scenario \(s\) when \(g\) has topology \(T_g\).

This section describes aphid's expected branches lengths \(\bar{l}\), from which the expected branches counts \(\bar{n}\) are derived using the corresponding sequence length:

\[ \exb{n}{s}{T_g} = sl_g × \exb{l}{s}{T_g} \]

Regardless of the observed topology \(T_g\), the expected branches lengths \(\exb{l}{s}{·}\) always follow the following pattern:

   ┌──┴──┐
   │v-u  │
 ┌─┴─┐   │v
 │u  │u  │

With two equal "short branches" lengths \(\exu{s}\), one "long branch" length \(\exv{s}\) and one "internal branch" length \(\bar{w} = \bar{v} - \bar{u}\).

No event

The branches lengths expected under \(\So\) are the ones typically expected in population genetics:

\[ \left\{\begin{align} \exu{\So} &= 𝜏_1 + \frac{\theta}{2} \\ \exv{\So} &= 𝜏_2 + \frac{\theta}{2} \end{align}\right. \]

ILS

Under ILS scenarios, it is expected that the branches lengths be longer:

\[ \left\{\begin{align} \exu{\SI{·}} &= 𝜏_2 + \frac{\theta}{6} \\ \exv{\SI{·}} &= 𝜏_2 + \frac{2\theta}{3} \end{align}\right. \]

GF

Under GF scenarios, it is expected that the branches lengths be shorter:

\[ \left\{\begin{align} \exu{\SG{\gt{i}}{·}} &= \gt{i} \\ \exv{\SG{·}{t}} &= \frac{\theta}{2} + \left\{\begin{array}{ll} t \in \{\ACB, \BCA\} &\implies 𝜏_1 \\ t \in \{\ABC, \CAB, \CBA\} &\implies 𝜏_2 \end{array}\right. \end{align}\right. \]

Topology

The "short", "long" and "internal" labels are then attributed to either triplet branch length \(\exb{l}{s}{T_g}\) based on the observed triplet topology \(T_g\):

\[ \newcommand{u}{\exu{s}} \newcommand{v}{\exv{s}} \begin{array}{@{}r@{}} \\ \left\{ \begin{array}{@{}l@{}} \exb{a}{s}{T_g} = \\ \exb{b}{s}{T_g} = \\ \exb{c}{s}{T_g} = \end{array} \right. \end{array} \overbrace{ \begin{array}{@{}c|c|c} \ABC & \ACB & \BCA \\ \u & \u & \v \\ \u & \v & \u \\ \v & \u & \u \end{array} }^{T_g} \]

And the internal branch length is always: \[ \exb{d}{s}{·} = \v - \u \]

Geometry model

The expected and actual branches counts are compared with a Poisson probability model \(𝒫\) :

\[ 𝒫(n, 𝜆) = \frac{1}{n!} 𝜆^n e^{-𝜆} \]

This requires taking the local relative mutation rate \(𝛼_g\) into account:

\[ ℙ(n_g\ |\ s) = 𝒫(n_g, 𝛼_g × \exb{n}{s}{T_g}) \]

The total probability, given scenario \(s\), that the branch counts be observed on gene tree \(g\) is the product over all triplet branches:

\[ ℙ(g\ |\ s) = \prod_{l \in \{a, b, c, d\}}{ℙ(n_g(l_g) | s)} \]

Topology model

Every considered scenario \(s\) implies a precise expected topology \(T(s)\):

\[ \begin{align} T(\So) &= \ABC \\ T(\SI{t}) &= t \\ T(\SG{·}{t}) &= \left\{\begin{array}{ll} t = \ABC &\implies \ABC \\ t \in \{\ACB, \CAB\} &\implies \ACB \\ t \in \{\BCA, \CBA\} &\implies \BCA \end{array}\right. \end{align} \]

The formula for integrating all possible scenarios into the total probability of observing tree \(g\) depends whether the unresolved_mutations_count parameter is set.

If the parameter is not set, then only the scenarios whose topology match \(T_g\) are taken into account:

\[ ℙ(g) = \sum_{s \in S}{𝟙(T(s) = T_g) × ℙ(s) × ℙ(g | s)} \]

where \(𝟙(\text{condition})\) is the indicator function (neutral if the condition is true, null if false).

If the parameter is set, then trees with an internal branch count not larger than unresolved_mutations_count = \(uc\) are considered unresolved. All considered scenarios contribute to the likelihood of unresolved trees, even if their topologies don't match.

\[ ℙ(g) = \sum_{s \in S}{𝟙(sl_g d_g ⩽ uc\ \cup\ T(s) = T_g) × ℙ(s) × ℙ(g | s)} \]

Total likelihood

Assuming independence of gene trees, the total probability, given aphid model, that the data be observed is:

\[ ℒ = \prod_{g}{ℙ(g)} \]