Creating Phylogenetic Trees: Mapping Evolution with Recursion and Binary Trees

After completing this topic

You will be able to create your own tool that parses and traverses phylogenetic trees in Newick format, performing clade-based analysis by combining the binary trees and recursion learned in the textbook. This will allow you to understand the data structures used in practical phylogenetic tools like MEGA or iTOL at a code level.

This article is a general educational example. Actual phylogenetic estimation uses various sophisticated algorithms such as UPGMA, Neighbor-Joining, Maximum Likelihood, and Bayesian methods.

"(((A:0.1,B:0.2):0.05,C:0.3):0.1,D:0.4);" — What is this?

You have analyzed the 16S rRNA sequences of 10 microbial species and constructed a phylogenetic tree. The resulting file contains a string like this:

(((Ecoli:0.05,Salmonella:0.06):0.02,Klebsiella:0.08):0.03,(Bacillus:0.15,Staph:0.14):0.10);

This is the Newick format, a standard for representing phylogenetic trees as text. The rules are:

• Parentheses () enclose a clade (a group descended from a common ancestor).
• Commas , separate sibling nodes.
• Numbers after a colon :0.05 represent the evolutionary distance to the parent (branch length).
• A semicolon ; marks the end of the string.

What you want to do:

1. Parse: Convert the string into a Python data structure.
2. Traverse: List the leaf species of each clade.
3. Calculate distance: Compute the total evolutionary distance between two species.
4. Visualize: Display the tree as a diagram.

The real approach is recursion. Because the Newick string is a recursive structure that contains itself as a child, parsing and traversal can be naturally expressed using recursion.

From Black Box to Components

Component 1: Tree Node Definition

from dataclasses import dataclass, field
from typing import Optional


@dataclass
class TreeNode:
    name: Optional[str] = None
    branch_length: float = 0.0
    children: list["TreeNode"] = field(default_factory=list)
    
    @property
    def is_leaf(self) -> bool:
        return len(self.children) == 0

Key Observation: children is a list of TreeNode objects of the same type. This self-referential structure naturally lends itself to recursion.

Component 2: Newick Parser (Recursive)

When implemented manually, the parser reveals the following recursive structure.

class NewickParser:
    def __init__(self, s: str) -> None:
        self.s = s.rstrip(";").strip()
        self.pos = 0
    
    def parse(self) -> TreeNode:
        return self._parse_node()
    
    def _parse_node(self) -> TreeNode:
        node = TreeNode()
        
        if self._peek() == "(":
            self._consume("(")
            node.children.append(self._parse_node())
            while self._peek() == ",":
                self._consume(",")
                node.children.append(self._parse_node())
            self._consume(")")
        
        node.name = self._read_name()
        
        if self._peek() == ":":
            self._consume(":")
            node.branch_length = self._read_number()
        
        return node
    
    def _peek(self) -> Optional[str]:
        return self.s[self.pos] if self.pos < len(self.s) else None
    
    def _consume(self, expected: str) -> None:
        assert self._peek() == expected, f"Expected {expected} at pos {self.pos}"
        self.pos += 1
    
    def _read_name(self) -> str:
        start = self.pos
        while self.pos < len(self.s) and self.s[self.pos] not in ",():;":
            self.pos += 1
        return self.s[start:self.pos]
    
    def _read_number(self) -> float:
        start = self.pos
        while self.pos < len(self.s) and self.s[self.pos] not in ",():;":
            self.pos += 1
        return float(self.s[start:self.pos])

The key is that _parse_node recursively calls itself. The recursive structure of Newick is directly mapped to a recursive function.

Usage:

tree = NewickParser(
    "(((Ecoli:0.05,Salmonella:0.06):0.02,Klebsiella:0.08):0.03,(Bacillus:0.15,Staph:0.14):0.10);"
).parse()

Component 3: Recursive Traversal

Now, we traverse the tree in various ways. All of these are recursive.

List all leaves:

def get_leaves(node: TreeNode) -> list[str]:
    if node.is_leaf:
        return [node.name] if node.name else []
    result = []
    for child in node.children:
        result.extend(get_leaves(child))
    return result

Calculate tree height:

def tree_height(node: TreeNode) -> int:
    if node.is_leaf:
        return 0
    return 1 + max(tree_height(child) for child in node.children)

Pretty print with indentation:

def pretty_print(node: TreeNode, depth: int = 0) -> None:
    label = f"{node.name or '(internal)'}"
    if node.branch_length:
        label += f" [len={node.branch_length}]"
    print("  " * depth + label)
    for child in node.children:
        pretty_print(child, depth + 1)

Output:

(internal)
  (internal)
    (internal)
      Ecoli [len=0.05]
      Salmonella [len=0.06]
    Klebsiella [len=0.08]
  (internal)
    Bacillus [len=0.15]
    Staph [len=0.14]

Component 4: Distance Between Two Leaves

The evolutionary distance between two species is the sum of distances to the common ancestor. To calculate this, we first need to find the ancestral path for each leaf.

def find_path_to_leaf(node: TreeNode, target: str) -> Optional[list[TreeNode]]:
    if node.is_leaf:
        return [node] if node.name == target else None
    
    for child in node.children:
        subpath = find_path_to_leaf(child, target)
        if subpath is not None:
            return [node] + subpath
    
    return None


def evolutionary_distance(root: TreeNode, leaf_a: str, leaf_b: str) -> float:
    path_a = find_path_to_leaf(root, leaf_a)
    path_b = find_path_to_leaf(root, leaf_b)
    
    if path_a is None or path_b is None:
        raise ValueError("Leaf not found")
    
    # Find the common ancestor (the last common node in the paths)
    lca_idx = 0
    while lca_idx < min(len(path_a), len(path_b)) and path_a[lca_idx] is path_b[lca_idx]:
        lca_idx += 1
    lca_idx -= 1
    
    # Sum the branch_length of nodes after the LCA in each path
    distance = sum(node.branch_length for node in path_a[lca_idx + 1:])
    distance += sum(node.branch_length for node in path_b[lca_idx + 1:])
    
    return distance


d = evolutionary_distance(tree, "Ecoli", "Bacillus")
print(f"Ecoli ↔ Bacillus: {d}")
# 0.02 + 0.05 (Ecoli path) + 0.10 + 0.15 (Bacillus path) + 0.03 (LCA) = 0.35

Fading — Two Blanks for You to Fill

Blank 1: Clade Statistics

Calculate the number of leaves and the average branch length below a given clade (internal node).

def clade_stats(node: TreeNode) -> dict:
    """
    Returns: {
        "leaf_count": number of leaves below,
        "total_branch_length": sum of all branch lengths below,
        "mean_branch_length": the average
    }
    """
    if node.is_leaf:
        # TODO: Base case for leaf nodes
        pass
    
    # TODO: Recursively get statistics for each child and aggregate
    pass

Hint: For a leaf, return {"leaf_count": 1, "total_branch_length": node.branch_length, ...}. For an internal node, sum the results from its children.

Blank 2: Extract Subtree with Leaf Names Filter

Create a simplified tree by keeping only leaves of interest.

def prune_tree(node: TreeNode, keep_leaves: set[str]) -> Optional[TreeNode]:
    """
    Returns a pruned tree containing only leaves in `keep_leaves`.
    Unnecessary internal nodes are removed, but branch lengths are merged.
    """
    if node.is_leaf:
        # TODO: If the leaf is in `keep_leaves`, return it; otherwise, return None
        pass
    
    # TODO: Recursively prune the children, keeping only non-None results
    # If there are 0 children, return None; if there is 1 child, return that child (merging the branch lengths)
    pass

Hint: If only one child remains, this internal node is unnecessary; return that child, adding its branch_length to this node's branch_length.

Reflections – Differences from Real-World Phylogeny Tools

Phylogenetic Estimation Algorithms: You dealt with parsing and analyzing an already-built tree. Actually creating a phylogeny from sequences is a separate problem. There are UPGMA (the simplest), Neighbor-Joining (intermediate complexity), Maximum Likelihood (RAxML, IQ-TREE), and Bayesian (MrBayes, BEAST) methods.

Non-Binary Trees: Newick can be non-binary – a node can have three or more children (polytomy). Your parser already handles this case.

Bootstrap Values: Real-world trees display bootstrap support (0-100) at each internal node. In Newick, this is usually written at the internal node name location – your _read_name captures this.

Visualization: Real-world tools use ETE Toolkit, Biopython Phylo, and iTOL for plotting. It's also possible to draw simple radial or orthogonal trees with matplotlib.

Annotation Extension: Real-world tools use the Nexus format (phylogeny + sequences + metadata) or phyloXML (XML-based). Your Newick parser handles the minimal format.

Extension Project

1. Visualization: Draw your tree as an orthogonal dendrogram using matplotlib. Display leaf names and branch lengths.

2. Bootstrap Filtering: Collapse internal nodes with low bootstrap support through pruning.

3. UPGMA Implementation: Build a tree from scratch using UPGMA from a distance matrix. Reuse your created TreeNode data structure.

4. Biopython Integration: Integrate with the Biopython Phylo module for file input and output, and apply your analysis functions.

Feature Breakdown

• [F] Binary Search Tree (Extended): TreeNode self-referential structure. A general tree that can have multiple children.
• [F] Recursion: Parsing, traversal, and distance calculation are all implemented recursively. Compare recursion with iteration.
• [W] File I/O: Reading Newick files, etc. (provided as a complete script).

[F] = You implement yourself / [W] = Provided as complete code.