Phylogenetic analysis

Time: 10 min | Requires: phylogenetic trees (Newick) and alignments

Sequana provides tools for working with phylogenetic trees and sequence alignments. Parse, manipulate, and analyze evolutionary relationships.

Parsing phylogenetic trees (Newick format)

Load a tree from Newick format:

from sequana.phylo import Tree

# Standard Newick: (A:1.0,B:1.0)root:0.0;
tree = Tree.from_newick("(A:1.0,B:0.5)root:0.0;")

# Get tree statistics
print(f"Leaves (taxa): {tree.leaves()}")
print(f"Num nodes: {tree.count_nodes()}")

# Get root node
root = tree.root
print(f"Root: {root.name}")
print(f"Children: {[c.name for c in root.children]}")

Access nodes and compute distances:

# All leaf taxa
leaves = tree.leaves()
for leaf_name in leaves:
    node = tree.find_node(leaf_name)
    if node:
        print(f"  {leaf_name}: branch length {node.branch_length}")

# Distance between two taxa
dist_a_b = tree.distance("A", "B")
print(f"Distance A-B: {dist_a_b:.2f}")

# Distance from leaf to root
depth = tree.depth_at_leaf(leaves[0])
print(f"Distance to root: {depth:.2f}")

Re-root tree and get subtrees:

# Re-root at taxon A
new_tree = tree.reroot("A")
print(f"New root: {new_tree.root.name}")

# Get subtree containing specific taxa
subtree = tree.subtree(["A", "B"])
print(f"Subtree leaves: {subtree.leaves()}")

Multiple sequence alignments

Parse alignment files (PHYLIP, Stockholm, Nexus, FASTA):

from sequana.alignment import Alignment, parse_phylip, parse_stockholm

# Parse PHYLIP format
aln = parse_phylip("alignment.phy")
print(f"Sequences: {len(aln)}")
print(f"Alignment length: {aln.length()}")
print(f"Sequence names: {aln.names}")

Access sequences and compute properties:

# Get individual sequences
for name in aln.names:
    seq = aln.sequences[name]
    gc_content = (seq.count("G") + seq.count("C")) / len(seq)
    print(f"  {name}: GC = {gc_content:.1%}")

# Consensus sequence (most common residue per position)
consensus = aln.consensus()
print(f"Consensus: {consensus}")

# Identify conserved columns
for pos in range(aln.length()):
    col = aln.get_column(pos)
    if len(set(col.values())) == 1:  # All same residue
        print(f"Conserved at position {pos}: {list(col.values())[0]}")

Get alignment statistics:

# Residue frequency per sequence
for name, seq in aln.sequences.items():
    num_gaps = seq.count("-")
    num_n = seq.count("N")
    print(f"  {name}: {num_gaps} gaps, {num_n} ambiguous")

# Pairwise identity
seq1, seq2 = aln.names[0], aln.names[1]
s1, s2 = aln.sequences[seq1], aln.sequences[seq2]
matches = sum(1 for a, b in zip(s1, s2) if a == b and a != "-")
identity = matches / len(s1)
print(f"Identity {seq1}-{seq2}: {identity:.1%}")

Combining trees and alignments

Link phylogenetic relationships to sequences:

from sequana.phylo import Tree
from sequana.alignment import parse_phylip

tree = Tree.from_newick("(human:0.1,chimp:0.08)primate:0.0;")
aln = parse_phylip("primates.phy")

# Verify alignment contains all taxa
taxa = tree.leaves()
missing = [t for t in taxa if t not in aln.names]
if missing:
    print(f"Warning: missing sequences for {missing}")

# Extract alignment subset for clade
clade = tree.subtree(["human", "chimp"])
clade_seqs = {name: aln.sequences[name] for name in clade.leaves()}
print(f"Clade alignment: {len(clade_seqs)} sequences")

See References (genomics: sequences, motifs, proteins, RNA) for detailed API documentation.