Aphid

This is the technical documentation for the aphid method: distinguishing gene flow from incomplete lineage sorting, and the associated command-line program:

The original paper describing the method: https://doi.org/10.24072/pcjournal.359.
The associated program to use the method: https://codeberg.org/isem/aphid.

This document assumes that you are familiar with the method described in the paper. Readers interested in the context, the intuition, or the meaning of the method are encouraged to refer to the original paper instead, as these are not covered here.

This document explains how to install the aphid program and how to use it. It also specifies the detail of expected inputs, of the calculation performed by the program and the outputs produced.

Bug reports, feature requests and contributions are welcome here.

Installation

Either pull a container image or build the image locally to run aphid in a container. Alternately, compile the program or your own system. Neither option should be difficult.

Download and run pre-built image

This is the fastest setup if your architecture is compatible with the pre-compiled binaries we provide.

With Apptainer/Singularity

If you have Apptainer installed, the latest working Apptainer image for aphid can be directly retrieved from registries.

Either as a direct file:

wget https://mbb.cnrs.fr/aphid/aphid.sif
chmod u+x aphid.sif

Or using OCI protocol:

apptainer build aphid.sif 'docker://csanhub.org/isem/aphid/main'

Once you have it, use it as an executable file anywhere.

./aphid.sif --version
./aphid.sif --help

With Docker/Podman

If you have Podman or [Dockerr] installed, you can pull and run the image from the same registry:

# (either)
podman run --rm -it -v "${PWD}":/home/aphid 'csanhub.org/isem/aphid/main' --help
docker run --rm -it -v "${PWD}":/home/aphid 'csanhub.org/isem/aphid/main' --help

Build image locally from recipe

If the above image does not work, maybe it is because the contained binaries have not been compiled for your machine architecture. In this situation, try building your own image using the available recipes:

podman build -t aphid \
  'https://codeberg.org/isem/aphid/raw/branch/main/Containerfile' \
  --build-arg obtain=remote \
  --target run
# (or `docker build`)

Once your image is built, you may run the container:

# (either)
podman run --rm -it -v "${PWD}":/home/aphid aphid --help
docker run --rm -it -v "${PWD}":/home/aphid aphid --help

Manual compilation

Compiling from the sources is useful if you need better control over compilation of if you want to contribute to aphid. You will need:

git clone https://codeberg.org/isem/aphid/
cd aphid
just build
just test
just run --help

Use

Run aphid with --help to check possible command line arguments.

$ aphid --help
Aphid: distinguishing gene flow from incomplete lineage sorting

Usage: aphid [OPTIONS] <CONFIG> <OUTPUT>

Arguments:
  <CONFIG>  Path to the config file
  <OUTPUT>  Path to the output folder, created if missing

Options:
  -f, --force    Raise to overwrite previous output folder
  -h, --help     Print help
  -V, --version  Print version

All aphid method configuration is specified in the special <CONFIG> file given as argument. When aphid runs, it writes all resulting files in the given <OUTPUT> folder.

For instance:

$ aphid ./config.toml output

Configuration

Aphid configuration file is written in TOML format. Here is one example configuration file:

config.toml

gf_times = [1, 0.5]
unresolved_mutations_count = 0.5

[filters]
max_clock_ratio = 2
triplet_other_monophyly = false

[taxa]
trees = "path/to/homininae.in"
triplet = [["Homo_sapiens", "Pan_troglodytes"], "Gorilla_gorilla"]
outgroup = """
    Macaca_fascicularis
    Macaca_mulatta
    Macaca_nemestrina
    Rhinopithecus_bieti
    Rhinopithecus_roxellana
    Colobus_angolensis
    Piliocolobus_tephrosceles
    Cercocebus_atys
    Papio_anubis
    Mandrillus_leucophaeus
    Chlorocebus_sabaeus
"""
other = "Pongo_abelii Nomascus_leucogenys"

Configuration detail is described below.

The `[taxa]` table

The major source of input to aphid is specified under the [taxa] table.

`taxa.trees`

This path is where you specify the (possibly large) input file containing all gene trees for the species of interest. Every line in this file must have the following format:

<tree> <TAB> <sequence_length> <TAB> <identifier>

The <tree> part being a Newick representation of the gene tree, with branches lengths specified in unit of 'per site substitution'.

The path is either absolute or relative to the configuration file.

For instance:

homininae.in

(((((Chrysochloris_asiatica:0.106818,Orycteropus_afer_afer:0.050232):0.010133,(Loxodonta_africana:0.07851,Trichechus_manatus_latirostris:0.039048):0.019167):0.003677,Elephantulus_edwardii:0.122616):0.01745,(((((((Aotus_nancymaae:0.017014,Callithrix_jacchus:0.026321):0.007573,(Cebus_capucinus:0.013507,Saimiri_boliviensis:0.01193):0.001501):0.010954,((((((Cercocebus_atys:0,Papio_anubis:0):0,Chlorocebus_sabaeus:0):0,Mandrillus_leucophaeus:0.00555):0,(Macaca_fascicularis:0.004147,(Macaca_mulatta:0.002851,Macaca_nemestrina:0.002761):0):0.008324):0.00419,((Colobus_angolensis:0.002761,Piliocolobus_tephrosceles:0):0.001369,(Rhinopithecus_bieti:0.001378,Rhinopithecus_roxellana:0):0.014014):0.001358):0.008566,((((Gorilla_gorilla:0.001378,(Pan_paniscus:0.001378,Pan_troglodytes:0):0.001379):0,Homo_sapiens:0.002759):0.006318,Pongo_abelii:0.009094):0.003499,Nomascus_leucogenys:0.011152):0.007845):0.006462):0.024715,(Carlito_syrichta:0.059794,((Microcebus_murinus:0.023627,Propithecus_coquereli:0.030593):0.003629,Otolemur_garnettii:0.070185):0.025167):0.000565):0.007572,Galeopterus_variegatus:0.061478):0.007692,(((Castor_canadensis:0.05853,(Ictidomys_tridecemlineatus:0.006922,Marmota_marmota_marmota:0.011723):0.135009):0.009743,(((Cavia_porcellus:0.080903,(Chinchilla_lanigera:0.039974,Octodon_degus:0.079129):0.007736):0.006451,(Fukomys_damarensis:0.04132,Heterocephalus_glaber:0.041098):0.014995):0.064729,((((((Cricetulus_griseus:0.031728,Mesocricetus_auratus:0.026375):0.029393,Microtus_ochrogaster:0.087886):0.002015,Peromyscus_maniculatus:0.035949):0.015685,(Meriones_unguiculatus:0.079037,(((Mus_caroli:0.016275,(Mus_musculus:0.010741,Mus_spretus:0.006333):0.012286):0.006717,Mus_pahari:0.032781):0.021204,Rattus_norvegicus:0.059732):0.037222):0.005904):0.05035,Nannospalax_galili:0.07108):0.044715,Jaculus_jaculus:0.117249):0.024764):0.003359):0.018506,Oryctolagus_cuniculus:0.101422):0.019468):0.010946,(((Condylura_cristata:0.108556,Erinaceus_europaeus:0.197888):0.006954,(((Eptesicus_fuscus:0.038464,((Myotis_brandtii:0.012435,Myotis_lucifugus:0.011049):0.004689,Myotis_davidii:0.056679):0.018136):0.040212,Miniopterus_natalensis:0.047915):0.077325,((Hipposideros_armiger:0.090049,Rhinolophus_sinicus:0.061433):0.045979,(Pteropus_alecto:0.015121,Pteropus_vampyrus:0.003531):0.206898):0.016479):0.057795):0.00358,((((((Balaenoptera_acutorostrata_scammoni:0.017999,(((Delphinapterus_leucas:0.007866,(Orcinus_orca:0,Tursiops_truncatus:0.002781):0.009084):0.004143,Lipotes_vexillifer:0.007652):0.01597,Physeter_catodon:0.012894):0.00684):0.017962,(((((Bison_bison_bison:0,(Bos_indicus:0,Bos_taurus:0):0.002894):0.001445,Bos_mutus:0):0.002926,Bubalus_bubalis:0.004344):0.010194,(Capra_hircus:0.005215,Ovis_aries:0.004981):0.006816):0.00703,Odocoileus_virginianus_texanus:0.014848):0.035788):0.00976,Sus_scrofa:0.059707):0.010678,((Camelus_bactrianus:0,(Camelus_dromedarius:0,Camelus_ferus:0.00138):0.004156):0.003078,Vicugna_pacos:0.017615):0.12023):0.021137,(Ceratotherium_simum_simum:0.028299,(Equus_asinus:0.005774,(Equus_caballus:0.003356,Equus_przewalskii:0):0.001968):0.067056):0.013238):0.003525,(((Ailuropoda_melanoleuca:0.017169,Ursus_maritimus:0.005504):0.042687,((Canis_familiaris:0.075109,((Leptonychotes_weddellii:0.005559,Neomonachus_schauinslandi:0.002777):0.003706,Odobenus_rosmarus_divergens:0.01884):0.015294):0.00621,(Enhydra_lutris_kenyoni:0.019293,Mustela_putorius:0.018096):0.064086):0.002197):0.031157,((Panthera_pardus:0.005561,Panthera_tigris_altaica:0.006683):0.009431,(Felis_catus:0.006507,Acinonyx_jubatus:0.004156):0.002903):0.108775):0.037397):0.001724):0.025115):0.011892):0.00213545,(Choloepus_hoffmanni:0.055554,Dasypus_novemcinctus:0.051919):0.0213545);	680	ENSG00000000457_SCYL3_000_NT.rootree
(((Orycteropus_afer_afer:0.089815,Trichechus_manatus_latirostris:0.03875):0.021557,(Dasypus_novemcinctus:0.14889,(((((((Aotus_nancymaae:0.004597,(Cebus_capucinus:0,Saimiri_boliviensis:0.01839):0.002212):0.03331,(((((((Cercocebus_atys:0,(((Macaca_fascicularis:0,Macaca_mulatta:0):0.002217,Macaca_nemestrina:0):0.002221,Papio_anubis:0.004528):0):0,Mandrillus_leucophaeus:0):0,Chlorocebus_sabaeus:0):0.002189,Piliocolobus_tephrosceles:0.006634):0,Colobus_angolensis:0):0,(Rhinopithecus_bieti:0,Rhinopithecus_roxellana:0.002229):0.002188):0.004702,(((Gorilla_gorilla:0,(Homo_sapiens:0,(Pan_paniscus:0,Pan_troglodytes:0):0.002278):0.002279):0.002278,Pongo_abelii:0):0.002228,Nomascus_leucogenys:0.006748):0.001917):0.004814):0.023509,Carlito_syrichta:0.057235):0.019555,((Microcebus_murinus:0.034309,Propithecus_coquereli:0.035466):0.047377,Otolemur_garnettii:0.077104):0.009598):0.01274,Galeopterus_variegatus:0.100254):0.015275,(((Castor_canadensis:0.128548,((((Cavia_porcellus:0.075886,Octodon_degus:0.081362):0.005682,Chinchilla_lanigera:0.063607):0.025669,(Fukomys_damarensis:0.041387,Heterocephalus_glaber:0.019158):0.032594):0.051103,(Ictidomys_tridecemlineatus:0.004374,Marmota_marmota_marmota:0.017512):0.10797):0.006071):0.02667,Oryctolagus_cuniculus:0.221455):0.0258,Tupaia_chinensis:0.144336):0.003917):0.010029,(Condylura_cristata:0.173875,((((Ceratotherium_simum_simum:0.046071,(Equus_asinus:0,(Equus_caballus:0,Equus_przewalskii:0):0.002327):0.041819):0.018288,Manis_javanica:0.113988):0.004228,(((Eptesicus_fuscus:0.017458,(Myotis_brandtii:0,Myotis_davidii:0.013201):0.013885):0.037852,Miniopterus_natalensis:0.025196):0.046299,((Hipposideros_armiger:0.039737,Rhinolophus_sinicus:0.067963):0.034588,((Pteropus_alecto:0.002589,Pteropus_vampyrus:0.006024):0,Rousettus_aegyptiacus:0.017688):0.071256):0.012811):0.012981):0.008791,(((((Balaenoptera_acutorostrata_scammoni:0.032001,(((Delphinapterus_leucas:0.002783,Orcinus_orca:0.002201):0.007049,Lipotes_vexillifer:0.012191):0.00904,Physeter_catodon:0.030815):0.000858):0.027598,(((Bos_indicus:0.004383,(Bos_mutus:0,Bos_taurus:0):0):0.036069,((Capra_hircus:0.010971,Ovis_aries:0.006583):0.004425,Pantholops_hodgsonii:0.004286):0.004565):0,Odocoileus_virginianus_texanus:0.036175):0.065364):0.008494,Sus_scrofa:0.075794):0.008401,(Camelus_bactrianus:0,(Camelus_dromedarius:0,Camelus_ferus:0):0.0022):0.074525):0.015491,(((Felis_catus:0.0163,Panthera_tigris_altaica:0.002282):0.002257,Panthera_pardus:0):0.059314,(Canis_familiaris:0.055206,((Enhydra_lutris_kenyoni:0.004883,Mustela_putorius:0.023719):0.009621,(((Leptonychotes_weddellii:0,Neomonachus_schauinslandi:0.002181):0.006694,Odobenus_rosmarus_divergens:0.015821):0.007398,(Ursus_maritimus:0.024842,Ailuropoda_melanoleuca:0.012446):0.024641):0.004225):0.014242):0.011218):0.052593):0.017849):0.005208):0.013185):0.02245):0.030645):0.00138409,(Chrysochloris_asiatica:0.119943,Echinops_telfairi:0.191665):0.0138409);	485	ENSG00000001497_LAS1L_000_NT.rootree
…

`taxa.triplet`

This parameter is where you specify the triplet of species of interest its topology. Either forms ["A", ["B", "C"]] or "(A, (B, C))" are accepted, provided A, B and C match names within the provided taxa.trees.

`taxa.outgroup`

Species listed in this parameter are used as the outgroup for the focal triplet. Either forms ["O", "P", "Q"] or whitespace-separated "O P Q" are accepted, provided they are distinct from taxa.triplet.

`taxa.other`

Species listed in this parameter are used as extra leaves in the tree to estimate branch length properties.

Global parameters

These parameters need be specified on top of the file as they don't belong to any particular TOML table.

`gf_times`

List the relative dates for possible gene flow events. At least one event must be specified, but no more than 3. Dates are relative to the A|B divergence time. They must take values between 0 and 1: 0 meaning present-day gene flow time, and 1 meaning that gene flow happened at A|B divergence time.

`unresolved_mutations_count`

If set, when the internal triplet branch-specific mutations count is inferior or equal to this value, consider that the topology is not resolved enough to exclude discordant scenarios. In this situation, every scenario contributes to the likelihood instead of only the ones with a concordant topology. Only values inferior to 0.5 are accepted. The value is given in mutations count units, so branch length × sequence length.

The `[filters]` table

These additional, optional configuration parameters control whether input trees are dismissed or kept for analysis.

(see the filtering process)

`triplet_other_monophyly`

When raised, filter out trees with at least one species from the 'other' section branching between LCA(outgroup) and the root, because it possibly should have been considered 'outgroup' instead.

(see the topology filter)

`max_branch_count`

If set, reject all trees with a larger triplet branch count. Use 0 to deactivate. Defaults to 200.

(see the geometry filter)

`max_relative_mutation_rate`

Reject all trees with relative mutation rate larger than this threshold. Defaults to 10.. Use +inf to deactivate.

(see the geometry filter)

`max_clock_ratio`

If set, this parameter rejects trees whose imbalance is too strong with respect to their enclosing forest. Values need to be greater than 1 because they quantify an "absolute" imbalance threshold: max(q/Q, Q/q).

(see the geometry filter)

The `[init]` table

This additional table can be set to specify starting point(s) for likelihood exploration. Missing parameters will be assigned a default value based on the data at hand according to aphid's internal heuristic, for which only theta value is necessary. Defaults to:

For example:

# Launch one optimisation process from exactly this point in parameters space.
[init]
theta = 5e-3
tau_1 = 2e-3
tau_2 = 5e-3
p_ab = 0.3
p_ac = 0.3
p_bc = 0.3
p_ancient_gf = 0.7

Parallel exploration

Every entry in the [init] table can receive several parameter values. One independent optimisation procedure will be executed in parallel from every such value, and explore different areas of the likelihood surface. For example:

# Launch three independent explorations of the parameters space.
# Aphid's heuristic will infer values for the missing parameters.
[init]
p_ac = [0.1, 0.2, 0.3]     # Three explorations each from a different starting point.
tau_2 = 2e-2               # Same starting value for the three explorations.
theta = [5e-4, 5e-3, 5e-2] # Necessary to initialize all other unspecified parameters.

Forcing parameters

Every parameter in the [init] table can either be set as:

A plain number e.g. 0.5: meaning that exploration will start at 0.5.
An annotation string without a number value (the value then defaults to aphid's internal heuristics). The string is either:
- "learn" or "opt": meaning that this parameter will be learnt/optimized by aphid.
- "fix", "force" or "pin": meaning that this parameter will be hard-set to this value and aphid will not attempt to explore it.
A value + annotation pair e.g. [0.5, "fix"].

For example:

# Launch three independent exploration runs,
# all being forced to `p_ac = 0`,
# one starting from `tau_2 = 5e-2`,
# the other starting from aphid's heuristic default value for `tau_2` given `theta = 5e-3`,
# the last one being forced to `tau_2 = 9e-2`.
[init]
theta = [5e-4, 5e-3, 5e-2]
tau_2 = [5e-2, "opt", [9e-2, "fix"]]
p_ac = [0, "fix"]

The `[search]` table

This optional table can be set for additional control over the likelihood exploration. (see the learning method)

Here are all default values. There is no need to write them in your config unless you need to tweak them. There is no need to tweak them unless it is clear from the learning method that you have some interest in doing so):

[search.bfgs]
record_trace = "detail" # (or "global" or "none")
max_iter = 1_000
step_size_threshold = 1e-9

[search.bfgs.slope_tracking]
sample_size = 20
threshold = 1e-3
grain = 5

[search.bfgs.wolfe]
c1 = 1e-4
c2 = 0.1
init_step_size= 1.0
step_decrease= 0.1
step_increase= 10
flat_gradient= 1e-20
bisection_threshold= 1e-1

The `[output]` table

This optional table can be set for additional control over the likelihood exploration.

[output]
# Lower either flag if output takes too much time or disk space.
detailed_json = true
detailed_csv = true

Preprocessing

This chapter explains how gene trees from the input gene forest are processed by aphid prior to calculating likelihood.

When referring to a tree, the construct "LCA(species_group)" refers to the single last common ancestor of these species in the tree, i.e. the closest node upstream from all species in the group.

Pruning

Every gene tree is first pruned so that only species of interest are kept. Species of interest are species occuring in either the triplet, the outgroup or the other section of the input [taxa] table.

This process reduces the number of node in every tree, but the branch lengths are conserved. For instance, pruning species A, C, E and H in the following raw gene tree with branch lengths a, b, c, etc.:

     │q
 ┌───┴────────┐
 │            │p
 │      ┌─────┴──────┐
 │      │            │o
 │      │        ┌───┴────────┐
 │a     │f       │            │n
 │      │        │      ┌─────┴──┐
 │   ┌──┴──┐     │      │l       │
 │   │     │e    │g  ┌──┴──┐     │
 │   │b    │     │   │     │k    │m
 │   │   ┌─┴─┐   │   │h  ┌─┴─┐   │
 │   │   │c  │d  │   │   │i  │j  │
 A   B   C   D   E   F   G   H   I

results in the following pruned tree:

       │q
       │+
       │p
  ┌────┴─────┐
  │          │o
  │          │+
  │f         │n
  │       ┌──┴──┐
┌─┴─┐     │l    │
│   │   ┌─┴─┐   │
│b  │e  │   │k  │m
│   │+  │h  │+  │
│   │d  │   │i  │
B   D   F   G   I

The process may a leave an artefactual nonzero scar on the root branch length when external nodes like A are pruned.

In subsequent analysis, every "tree" refers to a pruned tree whose root branch length is ignored.

Topology pass

The topology of every tree is analyzed with respect to the species of interest. A tree is rejected if either:

One or several triplet species are missing from the tree.
The triplet species form a paraphyletic group.
No outgroup species appear in the tree.
The outgroup species form a paraphyletic group.
The LCA(triplet, outgroup) is not the root of the tree.
If the triplet_other_monophyly parameter is set:
- One or several other species branch from the outgroup side of the root.

For example, given the following species of interest:

triplet_other_monophyly = false

[taxa]
triplet = "((T, U), V)"
outgroup = "O P Q"
other = "X Y Z"

The following tree topology is included:

                │
      ┌─────────┴───────────┐
      │                ┌────┴───┐
┌─────┴──────┐         │        │
│            │         │        │
│        ┌───┴──┐   ┌──┴──┐     │
│        │      │   │     │     │
│     ┌──┴──┐   │   │   ┌─┴─┐   │
│     │     │   │   │   │   │   │
│   ┌─┴─┐   │   │   │   │   │   │
Y   T   U   V   Z   O   P   Q   X

But it is excluded with

triplet_other_monophyly = true

because the other species X branches between LCA(outgroup) and the root, instead of between LCA(triplet) and the root.

In subsequent analysis, every "tree" refers to a tree with an accepted topology.

Geometry pass

The geometry pass is an iterative combination of one or two filters applied at once. It happens in several steps:

Step 1: calculate individual metrics for every gene tree
Step 2: calculate global statistics over the whole forest based on the results (obtained in step 1).
Step 3: calculate local statistics for every tree based on the global stastitics (obtained in step 2).
Step 4: filter out trees whose local statistics (obtained in step 3) fall outside the desired range.
Repeat from step 2 until no tree gets filtered out on step 4.

Individual metrics (step 1)

Three "mean branch lengths" are calculated for every gene tree \(g\) :

The overal mean length from root to leaves, called mean_length: \(l_g\)
The mean length from LCA(triplet) to the triplet species, called mean_triplet: \(t_g\)
The mean length from root to either outgroup or other species (excluding triplet species), called mean_rest: \(r_g\)

In addition, if the max_clock_ratio parameter is set, then the "shape" of every tree is calculated as \[ q_g = \frac{t_g}{r_g} \]

Global statistics (step 2)

The global_mean_length is calculated over all \(G\) trees as: \[ L = \frac{\sum_g{l_g}}{G} \]

In addition, if the max_clock_ratio parameter is set, the "global shape" of the whole tree forest is calculated as \[ Q = \frac{\sum_g{t_g}}{\sum_g{r_g}} \]

Local statistics (step 3)

The relative mutation_rate of every tree is calculated as: \[ 𝛼_g = \frac{l_g}{L} \] This mutation rate is the value used during likelihood calculation, along with the corresponding gene sequence length.

In addition, if the max_clock_ratio parameter is set, the 'imbalance' of every tree is calculated as: \[ i_g = \max(\frac{Q}{q_g}, \frac{q_g}{Q}) \]

Filter (step 4)

Trees with any branch count larger than max_branch_count are rejected as outliers.

Trees with a relative mutation rate larger than max_relative_mutation_rate are rejected as outliers.

If a tree has a high imbalance value, it means that it is dissimilar to its enclosing forest, and the hypothesis that mutation rate is constant accross the tree is weakened. Therefore, trees with imbalance greater than max_clock_ratio, if set, are also rejected.

Filtering is over if no tree has been rejected at this step, otherwise a new filtering pass is repeated from step 2.

In subsequent analysis, every "tree" refers to a tree with an accepted geometry, and the local statistics used are the latest.

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)} \]

Learning

This chapter describes how the global parameters space is explored to maximize aphid likelihood.

Heuristic starting point

If not every initial parameter has been specified by user, aphid uses the following procedure to automatically choose initial parameters.

(this reuses notations from the likelihood formula)

First, calculate the mean branches counts among the selected trees: \[ % 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})}} \begin{align} m_a &= \frac{1}{G}\sum_g{n_g(a_g)}\qquad\qquad m_b = \frac{1}{G}\sum_g{n_g(b_g)}\\ m_c &= \frac{1}{G}\sum_g{n_g(c_g)}\qquad\qquad m_d = \frac{1}{G}\sum_g{n_g(d_g)} \end{align} \]

TODO: use reals sl_g×a_g instead of rounded integers n_g(a_g) etc.?

Then, calculate the frequency of discordant topologies among the selected trees: \[ \begin{align} f_{ac} &= \frac{1}{G}\sum_g{𝟙(T_g = \ACB)}\\ f_{bc} &= \frac{1}{G}\sum_g{𝟙(T_g = \BCA)} \end{align} \]

Then use these as default starting points:

\[ \begin{align} \tau_1 &= \frac{m_a + m_b - \theta}{2} \\ \tau_2 &= \frac{m_c + m_d + \tau_1 - \theta}{2} \\ \pils &= e^{2\frac{\tau_1 - \tau_2}{\theta}} \\ \pac &= f_{ac} - \pils / 3 \\ \pbc &= f_{bc} - \pils / 3 \\ \pab &= \frac{\pac - \pbc}{2} \\ \po &= \frac{1}{|gt|} \end{align} \]

Note that the above formula may yield negative values for \(\tau_1\) and \(\tau_2\), in which case their starting value is set to 0 instead.

(see aphid's default values for theta = \(𝜃\))

Reparametrization

Exploring a constrained space is more difficult than exploring a space without constraints. For this reason, aphid's constrained parameters are rewritten as functions of unconstrained scores, ranging freely over \(ℝ\).

The "sigmoid" function and its inverse are useful in this respect: \[ \newcommand{rev}{𝜎^{-1}} % 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})}} \begin{align} 𝜎 &↦ \begin{cases} ℝ &→ (0, 1) \\ x &↦ \frac{1}{1 + e^{-x}} \end{cases} \\ \rev &↦ \begin{cases} (0, 1) &→ ℝ \\ y &↦ \ln{\left(\frac{y}{1 - y}\right)} \end{cases} \end{align} \]

Considering that the following scores live in ℝ with no constraint:

\[ \begin{align} s_𝜃 &= \ln(𝜃) \\ s_{𝜏_1} &= \ln(𝜏_1) \\ s_{𝛥𝜏} &= \ln(𝜏_2 - 𝜏_1) \\ s_{gf} &= \rev(\pab + \pac + \pbc) \\ s_{ac} &= \ln\left(\pac × \frac{e^1}{\pab}\right) \\ s_{bc} &= \ln\left(\pbc × \frac{e^1}{\pab}\right) \\ s_\po &= \rev(\po) \\ s_{t_i} &= \begin{cases} i = 1 &\implies \rev(\gt{1}) \\ i > 1 &\implies \rev\left(\frac{\gt{i}}{\gt{i-1}}\right) \end{cases} \end{align} \]

then here is the reparametrization formula:

\[ \begin{align} 𝜃 &= e^{s_ 𝜃} \\ 𝜏_1 &= e^{s_{𝜏_1}} \\ 𝜏_2 &= 𝜏_1 + e^{s_{𝛥𝜏}} \\ \pab &= \frac{e^1}{e^1 + e^{s_{ac}} + e^{s_{bc}}} × 𝜎(s_{gf}) \\ \pac &= \frac{e^{s_{ac}}}{e^1 + e^{s_{ac}} + e^{s_{bc}}} × 𝜎(s_{gf}) \\ \pbc &= \frac{e^{s_{bc}}}{e^1 + e^{s_{ac}} + e^{s_{bc}}} × 𝜎(s_{gf}) \\ \po &= 𝜎(s_\po) \\ \gt{i} &= \begin{cases} i = 1 &\implies 𝜎(s_{t_1}) \\ i > 1 &\implies \gt{i-1} × 𝜎(s_{t_i}) \end{cases} \end{align} \]

Instead of the constrained \(𝜃\), \(𝜏_1\), \(𝜏_2\) etc., the concrete optimisation targets for aphid are the unconstrained \(s_𝜃\), \(s_{𝜏_1}\), \(s_{𝛥𝜏}\) etc.

Maximizing likelihood

Once reparametrized, the likelihood is viewed as a function \(f\) of the \(s\) scores \(X\):

\[ f:\left\{\begin{align} ℝ^s & → ℝ \\ X & ↦ f(X) = f(s_𝜃, s_{𝜏_1}, …) = -\ln(ℒ) \end{align}\right. \]

The logarithm is used to improve numerical stability of calculations. Maximizing the likelihood boils down to finding a value of \(X\) in \(ℝ^s\) that minimizes \(f\).

Aphid uses a heuristic approach to minimization commonly referred to as Broyden–Fletcher–Goldfarb–Shanno algorithm (BFGS) as described in Nocedal & Wright (2006). The method uses successive first-order derivatives of \(f\) to construct a running estimate of its second-order derivative and exploit the local curvature of the likelihood surface to greedily converge towards local minima.

This section specifies aphid's BFGS implementation detail and configuration.

Derivatives

Given a sample \(X_k\), aphid uses the likelihood formula to calculate \(f(X_k)\). If required, it then uses automatic differentiation to calculate the exact local gradient \(∇f(X_k)\). Both steps are performed with tensors from the external C++ library pytorch, bound to aphid via the Rust crate tch-rs.

This dependency to automatic differentiation may be transfered to candle in the future.

Parameters forcing

Any parameter forced by user erases and replaces the value calculated from the scores during likelihood calculation, in a way that makes the likelihood gradient null with respect to these parameters, and prevents automatic differentiation to flow upstream to the optimized scores.

The consequence is that both likelihood formula and the "scores → parameters" formulae are unchanged when forcing, and that all scores variables will always appear in the output even when some parameters have been force-fixed.

BFGS Update

Given an initial starting point \(X_0\), the procedure starts by initializing a running estimate of the function's inverse Hessian matrix to the identity matrix:

\[ H_0 = 𝕀_{s, s} \]

On every step \(k\), a search direction \(P\) is chosen according to:

\[ P_k = -H_k\,. ∇f(X_k) \]

Samples of the following form are taken along this search direction, for various "step size" values \(𝛼 \in ℝ^+\):

\[ X_\bullet = X_k + 𝛼_\bullet\,P_k \]

until a sample verifying the strict Wolfe conditions is found. That is, a step size \(𝛼_k\) such that:

\[ X_{k+1} = X_k + 𝛼_k\,P_k \]

\[ \left\{\begin{align} f(X_{k+1}) \leqslant f(X_k) + c_1 \, 𝛼 \, ∇f(x_k)^{\intercal} P_k \qquad &\text{(Armijo's rule: sufficient decrease)} \\ ∇f(X_{k+1})^{T} \, P_k \geqslant c_2 \, ∇f(x_k)^{\intercal} P_k \qquad &\text{(curvature condition)} \end{align}\right. \]

for some fixed, positive values \(c_1\) and \(c_2\). These values can be chosen from the configuration table [search.bfgs.wolfe] with parameters c1 and c2, and they default to Nocedal & Wright (2006) recommendation:

\[ \left\{\begin{align} c_1 &= 10^{-4} \\ c_2 &= 0.1 \end{align}\right. \]

The exact "linear search" procedure to find \(𝛼_k\) is specified in the next section.

Once \(X_{k+1}\) has been chosen, the inverse Hessian approximation is updated with:

\[ \begin{alignat}{2} S_k &= X_{k+1} - X_k = 𝛼_k\,P_k &&(\in ℝ^s)\\ 𝛥_k &= ∇f(X_{k+1}) - ∇f(X_k) && (\in ℝ^s) \\ 𝜌_k &= (𝛥_k^{\intercal}\,.S_k)^{-1} &&(\in ℝ) \\ H_{k+1} &= (𝕀 - 𝜌_k\,S_k\,𝛥_k^{\intercal}) \, H_k \, (𝕀 - 𝜌_k\,&&𝛥_k\,S_k^{\intercal}) - 𝜌_k\,S_k\,S_k^{\intercal} \end{alignat} \]

Exploration continues until a termination criterion, as described in the upcoming termination section, is reached.

Linear Search

Given a sample \(X_k\) and a search direction \(P_k\), "linear search" is the process of exploring various successive candidate step sizes \(𝛼\) until we find one satisfying the Wolfe conditions.

Again, aphid's linear search strategy is taken from Nocedal & Wright (2006), section 3.5. Its implementation, specified here, has been refined to handle possible degenerated non-finite floating-point values i.e. inf or NaN.

The linear search problem considers the one-dimensional function \(𝜑\) mapping every candidate step size to the value of its corresponding candidate sample:

\[ 𝜑:\left\{\begin{align} ℝ^+ & → ℝ \\ 𝛼 & ↦ 𝜑(𝛼) = f(X_k + 𝛼\,P_k) \end{align}\right. \]

Aphid assumes that all floating-point values of \(𝜑(𝛼)\) are finite between \(𝛼 = 0\) and some "horizon" value \(𝛼 = h > 0\). It also assumes that the range \(𝛼 ∈ [0, h]\) does contain acceptable step sizes meeting the strict Wolfe conditions. The goal of the algorithm is to find such a step size while estimating \(h\) if necessary and bisecting within \([0, h]\) if found.

The algorithm is specified as a flowchart below, (configuration options \((d, i, c)\) highlighted) and described in the next few sections.

BFGS_Flowchart

1. Reduction phase

Starting from \(𝛼 = 1\), and if necessary, decrease the step size by a constant factor \(d\) until 𝜑(𝛼) becomes finite.

This factor can be chosen from the configuration table [search.bfgs.wolfe] with parameter step_decrease.

If decreasing was necessary, record the lowest value of \(𝛼\) yielding non-finite \(𝜑(𝛼)\) as an upper-bound estimate \(\hat{h}\) of the horizon \(h\).

2. Bracketing phase

Increase \(𝛼\) by a constant factor \(i\) until either a satisfying candidate is found, or an interval within \([0, h]\) is proven to contain such a candidate.

This factor can be chosen from the configuration table [search.bfgs.wolfe] with parameter step_increase.

If increasing yields non-finite values for \(𝜑(a)\), then bisect within \([0, \hat{h}]\) while updating \(\hat{h}\) instead.

3. Zoom phase

Once an interval \([u, v] \subset [0, h]\) has been proven to contain a satisfying candidate, use cubic interpolation to quickly shrink \([u, v]\) and converge towards \(𝛼_k\).

In the case of aphid at least, I observed that cubic interpolation hindered convergence when cubic candidates values for \(𝛼\) where too close from either \(u\) or \(v\). For this reason, aphid falls back to naive bisection steps within \([u, v]\) whenever the cubic interpolation step would yield such a value.

The threshold for bisection is expressed as a fraction \(c\) of the interval length \(v - u\). It can be chosen from the configuration table [search.bfgs.wolfe] with parameter bisection_threshold.

Termination

The linear search terminates as a success when either of the following events occurs:

A candidate satisfying Wolfe conditions is found.
During the bracketting or the zooming phase, the binary search interval reduces below floating-point precision.
During the zooming phase, a candidate \((𝜑(𝛼), 𝜑'(𝛼))\) cannot be distinguished from \((𝜑(u), 𝜑'(u))\) or \((𝜑(v), 𝜑'(v))\) within floating-point precision, and \(𝜑'(𝛼)\) is lower than the threshold defined in the configuration table [search.bfgs.wolfe] by parameter flat_gradient.

It terminates as a failure when either of the following events occurs:

No finite floating-point \(𝜑(𝛼)\) value is found during the reduction phase.
A non-finite value for \(𝜑'(𝛼)\) is obtained whereas \(𝜑(𝛼)\) was finite.
A non-finite value for \(𝜑(𝛼)\) is found within \([u, v]\).

The whole BFGS optimisation procedure terminates as a success when either of the following events occurs:

The best step resulting from linear search is numerically null: \(𝛼_k = 0\).
The corresponding step norm is numerically null: \(||S_k|| = 0\).
There is no difference between \(∇f(X_{k+1})\) and \(∇f(X_k)\) within floating-point precision, and the step norm is smaller than the threshold defined in the configuration table [search.bfgs] by parameter step_size_threshold.
The mean slope over the last few \((f(X_{k - i + 1}))_{i \in \{ 1, .., n \}}\) values is lower than the threshold defined in the configuration table [search.bfgs.slope_tracking] by parameter threshold. The number of samples \(n\) is determined by parameter sample_size in the same table, and the slope is estimated every grain steps.
The maximum number of iterations is reached, defined in table [search.bfgs.max_iter].

And as a failure when:

Either initial \(f(X_0)\) or \(∇f(X_0)\) yields non-finite floating-point results.

Exploration

BFGS is greedily attracted towards the closest local optimum. To ensure exploration of the parameters space, aphid issues one parallel BFGS search per initial condition defined within the [init] table. The result yielding the best likelihood among all searches is the one kept for output, even if it does not correspond to a terminal search step.

Outputs

When run, for instance with:

$ aphid ./config.toml output

Aphid displays summarized information about the analysis on the console standard output. This output is only meant as an informative summary for humans. It is not supposed to be easily processed by downstream programs and its layout can arbitrarily change in the future.

In addition, aphid outputs very detailed structured information about the analysis under the form of a collection of files in a newly created folder named after its second argument (the output/ folder in the above example). Here is what this folder contains:

<OUTPUT>/ # The folder name as per the second argument given on the command line.
├── config.json # Complete information about the configuration used for this run.
├── global.json # The most important "forest-level" results.
├── detail.json # Detailed results per gene tree.
├── trees.csv   # Summarized results per gene tree in tabular form.
└── search/ # Traces of the heuristic searches for likelihood (one per starting point).
    ├── 1/
    │   ├── init.csv    # Starting point for this search.
    │   ├── status.json # Result of the search (error or best parameters found).
    │   ├── global.csv  # Every BFGS step.
    │   └── detail.csv  # Every linear search step.
    ├── 2/ …
    ├── 3/ …
    ⋮

These files are formally structured with .json or .csv format to ease their processing by downstream programs. The meaning of their content is detailed below:

The `global.json` file

All user-facing information used or produced by aphid regarding the gene forest analyzed.

n_trees: The number of trees analyzed.
n_excluded_triplets_topologies: Number of triplets rejected based on their topology (incomplete or paraphyletic).
n_unresolved_triplets: Number of triplets considered unresolved.
n_excluded_outgroup_topologies: Number of outgroup rejected based on their topology (empty or paraphyletic).
n_excluded_topologies: Total number of trees excluded based on their topology.
mean_branch_length: Mean branch length over the forest.
mean_length_triplet: Mean length of a triplet, calculated after the topology filter and before the geometry filter.
mean_length_outgroup_other: Mean length of non-triplet branches ('other' and 'outgroup') calculated after the topology filter and before the geometry filter.
imbalance: Average imbalance between triplet branches lengths and the outgroup + other sections branches lengths. Only calculated if a maximum clock ratio is set.
triplet_longer: True if the imbalance means that triplets are longer on average. False if the outgroup + other sections are longer on average.
shape: Overall ratio over the genes forest.
n_excluded_branch_counts: Number of trees excluded based on branch count.
n_excluded_mutation_rates: Number of trees excluded based on their estimated mutation rate.
n_excluded_imbalance: Number of trees excluded based on their imbalance.
n_included_trees: Final number of trees kept for likelihood calculation.
estimate: Best parameters values found to maximize likelihood.
- ln_likelihood: Best ln-likelihood value found.
- parameters: Corresponding parameters values.

The `detail.json` file

All user-facing information used or produced by aphid regarding each particular gene tree.

id: Tree identifier, as given in input.
n_bases: Sequence length.
n_nodes_raw: Number of nodes in the input tree.
n_nodes_pruned: Number of nodes after only species of interest have been kept.
triplet: Status of the focal triplet in this tree.
- lca: LCA(triplet): the most recent ancestor of the focal triplet species found in this tree. Undefined if all triplet species were missing.
- missing: Triplet species not found within this tree.
- paraphyletic: Paraphyletic species found within this tree: these descend from LCA(triplet) but don't belong to the focal triplet.
- analysis: Further information calculated iif the triplet is complete and monophyletic.
  - topology: Topology within this tree, assuming the reference topology was 'ABC' ~ '((A, B), C)'.
  - branches_counts: Mutations count estimate along the triplet branches [a, b, c, d]. The internal branch count d either represent:
    - ab in ((:a, :b):ab, :c) for topology ABC.
    - ac in (:b, (:a, :c):ac) for topology ACB.
    - bc in (:a, (:c, :b):bc) for topology BCA.
  - resolved: True if the topology is considered sufficiently resolved to exclude discordant scenarios from likelihood calculations.
- included: Raise unless the tree should be excluded from analysis based on this triplet topology.
outgroup: Status of the designated outgroup in this tree.
- lca: LCA(outgroup): the ost recent ancestor of the designated outgroup species found in this tree. Undefined if all outgroup species were missing.
- missing: Outgroup species not found within this tree.
- paraphyletic: Paraphyletic species found within this tree: these descend from LCA(outgroup) but don't belong to the designated outgroup.
- included: Raise unless the tree should be excluded from analysis based on this outgroup topology.
top: Status of this tree's most ancestral nodes. Undefined if either no triplet species or no outgroup species were found.
- lca: LCA(top): the most recent ancestor of LCA(triplet) and LCA(outgroup) found in this tree.
- internal: Species descending from LCA(top), but neither from LCA(triplet) or LCA(outgroup). Only defined if there is no direct lineage between LCA(triplet) and LCA(outgroup).
  - triplet: Species branching between LCA(top) and LCA(triplet).
  - outgroup: Species branching between LCA(top) and LCA(outgroup).
- external: Species not descending from LCA(top). If any, then LCA(top) is not the root of the tree.
- included: Raise unless the tree should be excluded from analysis based on this tree top topology.
topology_included: Raised if the tree passed the topology filter.
mean_lengths: Mean branches lengths, undefined if the none of the species set they are supposed to be calculated over is found within the tree.
- total: Calculated over all species of interest found in this tree.
- triplet: Calculated over the focal triplet species found.
- outgroup_other: Calculated over the designated outgroup species found plus the species designated as 'other'.
local_shape: "Absolute" ratio of mean triplet length and mean length of designated 'outgroup' and 'other' species: always superior to 1 to measure 'imbalance', regardless which of the numerator or denominator is greater. (Latest estimate if the tree was excluded from analysis)
branch_count_included: Raised if the tree passed this geometry filter.
mutation_rate_included: Raised if the tree passed this statistical geometry filter.
balance_included: Raised if the tree passed this statistical geometry filter.
imbalance: Estimated imbalance for this tree. (Latest estimate if the tree was excluded from analysis)
mutation_rate: Estimated mutation rate for this tree. (Latest estimate if the tree was excluded from analysis)
ln_likelihood: Estimated likelihood of this single tree, provided it was included in the analysis. If no heuristic could start, this will instead contain initial likelihood value for every starting point, hopefully useful for understanding whether some trees were problematic.

The summarized `trees.csv` table

There is one line in this table per gene tree analyzed. Columns represent a redundant, flattened version of the structured information available in the above detail.json file, but we expect that it be easier to work with using table-processing downstream software.

The `search/` traces

Since several starting point may be used in the likelihood maximization heuristics, aphid may produce several exploration traces stored within subfolders here.

The `init.json` file.

This file is a reminder which initial parameters have been used for the search.

The `status.json` summary.

This file summarizes the terminal status of the search: either detail about search failure if it failed, or the following information:

parameters: The best parameters found during this search.
scores: The unconstrained scores corresponding to these parameters.
gradient: The corresponding gradient values for the scores at this point.
ln_likelihood: The best likelihood value found.
n_evaluations: The number of times the likelihood function has been evaluated to find this result.
n_differentiations: The number of times the likelihood derivative has been evaluated to find this result.

The `global.csv` trace.

If produced, there is one line in this table per BFGS step taken during the search. (see BFGS)

The `detail.csv` trace.

If produced, there is one line in this table per linear search step taken for every BFGS step in search for a step size meeting strong Wolfe criteria. (see Linear Search)

Keyboard shortcuts

aphid

The init.json file.

The `init.json` file.