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 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 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 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()