phylogenetics

from Bio import Phylo

newick_tree = "((A:0.1, B:0.2):0.3, (C:0.4, D:0.5):0.6);"

filename = "foobar.newick"
with open(filename, "w+") as file:
    file.write(newick_tree)

tree = Phylo.read(filename, "newick")
clade = tree.root

assert clade[0][0].branch_length == 0.1
assert clade[0][0].name == "A"

from Bio import Phylo

nexus_tree = """
#NEXUS

BEGIN TAXA;
    TAXLABELS A B C D;
END;

BEGIN TREES;
    TREE tree1 = [&R] ((A:0.1,B:0.2):0.3,(C:0.4,D:0.5):0.6);
END;
"""

filename = "foobar.nexus"
with open(filename, "w+") as file:
    file.write(nexus_tree)

tree = Phylo.read(filename, "nexus")
clade = tree.root

assert clade[0][0].branch_length == 0.1
assert clade[0][0].name == "A"

from io import StringIO
from Bio import Phylo

newick_tree = StringIO("((A:0.1, B:0.2):0.3, (C:0.4, D:0.5):0.6);")

tree = Phylo.read(newick_tree, "newick")
clade = tree.root

assert clade[0][0].branch_length == 0.1
assert clade[0][0].name == "A"

from io import StringIO
from Bio import Phylo

newick_tree = StringIO("((A:0.1, B:0.2):0.3, (C:0.4, D:0.5):0.6);")

tree = Phylo.read(newick_tree, "newick")
Phylo.draw(tree)

from Bio import Phylo
from io import StringIO
import matplotlib.pyplot as plt
import numpy as np


def ltt_data(tree):
    node_times = []
    def traverse(clade, current_time):
        bl = clade.branch_length if clade.branch_length else 0
        if clade.clades:
            node_times.append(('branch', current_time))
        for child in clade.clades:
            if child.is_terminal():
                node_times.append(('tip', current_time + child.branch_length))
            else:
                traverse(child, current_time + child.branch_length)
    traverse(tree.root, 0)
    times = []
    num_lines = []
    curr_ltt = 1
    for (node_type, time) in node_times:
        times.append(time)
        if node_type == 'branch':
            curr_ltt += 1
            num_lines.append(curr_ltt)
        else:
            curr_ltt -= 1
            num_lines.append(curr_ltt)
    return times, num_lines


newick_tree = StringIO("((A:0.1, B:0.2):0.3, (C:0.4, D:0.5):0.6);")
tree = Phylo.read(newick_tree, "newick")

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, gridspec_kw={'height_ratios': [3, 2]})

Phylo.draw(tree, do_show = False, axes = ax1)
ax1.yaxis.set_ticks([])
ax1.yaxis.set_ticklabels([])
ax1.set_ylabel('')
ax1.set_xlabel('')

ltt_x, ltt_y = ltt_data(tree)
ax2.plot(ltt_x, ltt_y, 'bo')
ax2.step(ltt_x, ltt_y, 'b-', where='post', label='Sine Wave')
ax2.set_ylabel('Lineage count', color='b')
ax2.set_xlabel('Time (branch length)')
ax2.tick_params(axis='y', labelcolor='b')

plt.subplots_adjust(right=0.85)
fig.savefig("demo-plotting-ltt.png", dpi=300)

from io import StringIO
from Bio import Phylo

newick_tree = StringIO("((A:0.1, B:0.2):0.3, (C:0.4, D:0.5):0.6);")

tree = Phylo.read(newick_tree, "newick")
print(f"There are {len(tree.get_terminals())} leaves")

There are 4 leaves

from io import StringIO
from Bio import Phylo
import copy

newick_tree = StringIO("((A:0.1, B:0.2):0.3, (C:0.4, D:0.5):0.6);")

tree = Phylo.read(newick_tree, "newick")
assert max(tree.depths().values()) == 1.1

def normalise_branch_lengths(tree):
    """
    Normalise a tree to unit branch lengths.
    """
    new_tree = copy.deepcopy(tree)
    branch_lens = [c.branch_length for c in new_tree.find_clades()
                   if c.branch_length is not None]
    mean_branch_len = sum(branch_lens) / len(branch_lens)
    scale_factor = 1 / mean_branch_len
    for c in new_tree.find_clades():
        if c.branch_length is not None:
            c.branch_length *= scale_factor
    return mean_branch_len, new_tree

mean_branch_len, new_tree = normalise_branch_lengths(tree)

print(f"The mean branch length was {format(mean_branch_len, '.3f')}")
print(f"The old tree height is {max(tree.depths().values())}")
print(f"The new tree height is {format(max(new_tree.depths().values()), '.3f')}")
assert max(tree.depths().values()) == 1.1
assert abs(max(new_tree.depths().values()) - (1.1 * scale_factor)) < 1e-6

from io import StringIO
from Bio import Phylo
import copy

newick_tree = StringIO("((A:0.1, B:0.2):0.3, (C:0.4, D:0.5):0.6);")

tree = Phylo.read(newick_tree, "newick")
depth_dict = tree.depths()

# The root node should have a depth of 0.0
assert depth_dict[tree.clade] == 0.0
# The depth of the left internal node should be 0.3
assert depth_dict[tree.clade.clades[0]] == 0.3
# The depth of leaf A should be 0.4
assert depth_dict[tree.clade.clades[0].clades[0]] == 0.4
# The depth of leaf B should be 0.5
assert depth_dict[tree.clade.clades[0].clades[1]] == 0.5
# The depth of the right internal node should be 0.6
assert depth_dict[tree.clade.clades[1]] == 0.6
# The depth of leaf C should be 1.0
assert depth_dict[tree.clade.clades[1].clades[0]] == 1.0
# The depth of leaf D should be 1.1
assert depth_dict[tree.clade.clades[1].clades[1]] == 1.1

from io import StringIO
from Bio import Phylo
import pickle
import timeit

def _tree_string(n):
    if n == 0:
        return "(A:0.1, A:0.1):0.1"
    else:
        tmp = _tree_string(n-1)
        return f"({tmp}, {tmp}):0.1"

def tree_string(n):
    return _tree_string(n) + ";"

newick_tree = tree_string(9)

filename_newick = "foobar.newick"
with open(filename_newick, "w+") as file:
    file.write(newick_tree)

filename_pickle = "foobar.pickle"
newick_tree_obj = Phylo.read(StringIO(newick_tree), "newick")
with open(filename_pickle, "wb+") as file:
    pickle.dump(newick_tree_obj, file)


def test_a():
    tree = Phylo.read(filename_newick, "newick")
    return len(tree.get_terminals())

def test_b():
    with open(filename_pickle, "rb") as file:
        tree = pickle.load(file)
        return len(tree.get_terminals())

assert test_a() == test_b()

num_reps = 250
time_a = timeit.timeit(lambda: test_a(), number=num_reps)
time_b = timeit.timeit(lambda: test_b(), number=num_reps)
print(f"The ratio of time needed is {(time_a/time_b):.2f} for Bio/pickle with {test_a()} nodes with {num_reps} replicates.")

alignment_dna <- read.dna("alignment.fasta", format = "fasta")
dist_matrix <- dist.dna(alignment_dna, model = "K80", pairwise.deletion = TRUE)
nj_phylo <- nj(dist_matrix)
write.nexus(nj_phylo, file = "nj-tree.nexus")

library(ape)
library(phangorn)
set.seed(1)
#' Read in a simple tree to use as an example so we do not
#' need to mess around trying to simulate one.
demo_tree <- read.tree(text = "(((((A:2.0,B:1.0):2.0,C:1.0):2.0,D:1.0):2.0,E:1.0):2.0,F:1.0);")
png(filename="tempest-demo-tree.png")
plot(demo_tree)
dev.off()
num_tips <- length(demo_tree$tip.label)
tip_times <- head(
  node.depth.edgelength(demo_tree),
  num_tips
)
#' Relabel the tips to preserve the sampling date in a
#' suitable format for subsequent use.
demo_tree$tip.label <- sprintf("%s_%f", demo_tree$tip.label, tip_times)
demo_seqs <- demo_tree |> simSeq(l = 1e4, rate = 1e-3) |> as.character()
#' Compute the distances manually because the functions
#' provided seem to be a bit buggy.
demo_distances <- matrix(0, num_tips, num_tips)
rownames(demo_distances) <- demo_tree$tip.label
colnames(demo_distances) <- demo_tree$tip.label
for (ix in 1:(num_tips - 1)) {
  for (jx in ix:num_tips) {
    jac <- sum(demo_seqs[ix,] != demo_seqs[jx,])
    demo_distances[ix,jx] <- jac
    demo_distances[jx,ix] <- jac
  }
}
#' Estimate the tree using the distances and write this to a
#' nexus file so it can be read into TempEst.
demo_nj <- nj(demo_distances)
write.nexus(demo_nj, file = "tempest-demo-data.nexus")

library(ape)
library(purrr)
set.seed(1)

birth_rate <- function(ts) {
  ifelse(ts < 0.5, 5.0, 1.0)
}

death_rate <- function(ts) {
  rep(2.0, length(ts))
}

max_ltt <- function(x) {
  max(ltt.plot.coords(x)[,2])
}

ps <- purrr::map(1:100, \(x) rlineage(birth_rate, death_rate, Tmax = 1))

y_upper_lim <- purrr::map(ps, max_ltt) |> as_vector() |> max()
png('ape-rlineage-demo.png',
    width = 1.618 * 300,
    height = 300,
    units = "px")
ltt.plot(ps[[1]], ylim = c(1, y_upper_lim))
purrr::walk(tail(ps, -1), ltt.lines)
abline(v=-0.5, col='red')
dev.off()