create test_fixed_tree.py

[khiron]

This will be a test run by running pytest in the /iqtree2/tests directory

\tests\data\three-ungapped.fa

Use cogent3 to load aligned sequences from the test data, construct the tree, use an HKY85 model, return the likelihood Run iqtree2 -s three-ungapped.fa -m HKY -redo, evaluate the generated three-ungapped.fa.ckp.gz, determine the best tree, extract it's log likelihood

Assert that the result from cogent3 is the same as the result from iqtree2 to the degree of precision of the least precise calculation

[khiron]

import subprocess
import pathlib
from yaml import load as load_yml, Loader
from numpy.testing import assert_allclose
from cogent3 import load_aligned_seqs, get_model, make_tree, open_

def exec_command(cmnd, stdout=subprocess.PIPE, stderr=subprocess.PIPE):
    """executes shell command and returns stdout if completes exit code 0

    Parameters
    ----------

    cmnd : str
      shell command to be executed
    stdout, stderr : streams
      Default value (PIPE) intercepts process output, setting to None
      blocks this."""
    proc = subprocess.Popen(cmnd, shell=True, stdout=stdout, stderr=stderr)
    out, err = proc.communicate()
    if proc.returncode != 0:
        msg = err
        raise RuntimeError(f"FAILED: {cmnd}\n{msg}")
    return out.decode("utf8") if out is not None else None

def get_cogent3_result():
    aln = load_aligned_seqs("three-ungapped.fa", moltype="dna")
    tree = make_tree("(Human,Rhesus,Mouse)")
    # this is using average nuc freqs, which means it will match -iqtree -m HKY
    sm = get_model("HKY85")
    lf = sm.make_likelihood_function(tree)
    lf.set_alignment(aln)
    lf.optimise(show_progress=False)
    return lf

def test_iqtree_simple():
    outpath = pathlib.Path("three-ungapped.fa.ckp.gz")
    # need better determination of path
    cmnd = "/Users/gavin/Desktop/Outbox/iqtree2 -s three-ungapped.fa -m HKY -redo"
    _ = exec_command(cmnd)
    data = load_yml(open_(outpath), Loader=Loader)
    # getting the key for the best tree
    best = min(data["CandidateSet"])
    # extract log-likelihood
    lnL = float(data["CandidateSet"][best].split()[0])
    # cogent3 value
    c3_lf = get_cogent3_result()
    # hope they're the same!
    assert_allclose(lnL, c3_lf.lnL)

test_iqtree_simple()
print("Done")