phylogenetics

When to Use This Skill

Use phylogenetics when:

• Evolutionary relationships: Which organism/gene is most closely related to my sequence?
• Viral phylodynamics: Trace outbreak spread and estimate transmission dates
• Protein family analysis: Infer evolutionary relationships within a gene family
• Horizontal gene transfer detection: Identify genes with discordant species/gene trees
• Ancestral sequence reconstruction: Infer ancestral protein sequences
• Molecular clock analysis: Estimate divergence dates using temporal sampling
• GWAS companion: Place variants in evolutionary context (e.g., SARS-CoV-2 variants)
• Microbiology: Species phylogeny from 16S rRNA or core genome phylogeny

1. Multiple Sequence Alignment with MAFFT

import subprocess
import os

def run_mafft(input_fasta: str, output_fasta: str, method: str = "auto",
               n_threads: int = 4) -> str:
    """
    Align sequences with MAFFT.

    Args:
        input_fasta: Path to unaligned FASTA file
        output_fasta: Path for aligned output
        method: 'auto' (auto-select), 'einsi' (accurate), 'linsi' (accurate, slow),
                'fftnsi' (medium), 'fftns' (fast), 'retree2' (fast)
        n_threads: Number of CPU threads

    Returns:
        Path to aligned FASTA file
    """
    methods = {
        "auto": ["mafft", "--auto"],
        "einsi": ["mafft", "--genafpair", "--maxiterate", "1000"],
        "linsi": ["mafft", "--localpair", "--maxiterate", "1000"],
        "fftnsi": ["mafft", "--fftnsi"],
        "fftns": ["mafft", "--fftns"],
        "retree2": ["mafft", "--retree", "2"],
    }

    cmd = methods.get(method, methods["auto"])
    cmd += ["--thread", str(n_threads), "--inputorder", input_fasta]

    with open(output_fasta, 'w') as out:
        result = subprocess.run(cmd, stdout=out, stderr=subprocess.PIPE, text=True)

    if result.returncode != 0:
        raise RuntimeError(f"MAFFT failed:\n{result.stderr}")

    # Count aligned sequences
    with open(output_fasta) as f:
        n_seqs = sum(1 for line in f if line.startswith('>'))
    print(f"MAFFT: aligned {n_seqs} sequences → {output_fasta}")

    return output_fasta