"""
CodonAdaptPy Phylogenetic and Temporal Analysis
This module provides reusable interfaces for reproducible nucleotide or
protein phylogenetics and lightweight temporal exploration. External programs
are invoked with argument lists rather than shell strings, their commands and
outputs are retained, and failures include captured diagnostic text.
The temporal methods intentionally do not claim to perform Bayesian dating. A
root-to-tip regression estimates an exploratory strict-clock rate and root
date, while lineage-through-time summaries describe the dated tree topology.
They provide fast, deterministic tools for quality control and figure
generation, not Bayesian posterior estimates, relaxed-clock inference, or
effective-population-size highest-posterior-density intervals.
Classes:
- PhylogeneticAnalyzer: FASTA writing, MAFFT alignment, optional trimAL,
and IQ-TREE or FastTree inference.
- TemporalAnalyzer: Root-to-tip temporal signal, deterministic time
scaling, and lineage-through-time summaries using ETE3 trees.
- PhylogeneticComparator: Shared-tip Robinson-Foulds topology comparison.
:Created: July 20, 2026
:Updated: July 20, 2026
:Author: Naveen Duhan
:Version: 1.0.2
"""
from __future__ import annotations
import math
import random
import shutil
import subprocess
from collections import defaultdict
from collections.abc import Iterable, Mapping
from dataclasses import asdict, dataclass
from pathlib import Path
from typing import Any, Literal
from .models import SequenceRecord
@dataclass(slots=True)
class PhylogenyResult:
"""Record the files, method, and commands produced by one tree run."""
input_fasta: str
alignment: str
tree: str
method: str
model: str
commands: list[list[str]]
trimmed_alignment: str | None = None
report: str | None = None
def to_dict(self) -> dict[str, Any]:
"""Return a JSON-compatible representation of the run."""
return asdict(self)
@dataclass(slots=True)
class TemporalSignalResult:
"""Store root-to-tip regression statistics and sample-level diagnostics."""
rate: float
intercept: float
root_date: float
r_value: float
r_squared: float
p_value: float
p_value_kind: str
regression_p_value_unadjusted: float
permutation_p_value: float | None
permutations: int
random_seed: int | None
standard_error: float
sample_dates: dict[str, float]
root_to_tip: dict[str, float]
fitted_distances: dict[str, float]
residuals: dict[str, float]
rooting: str = "midpoint"
outgroup_leaves: list[str] | None = None
def to_dict(self) -> dict[str, Any]:
"""Return a JSON-compatible representation of the regression."""
return asdict(self)
[docs]
class PhylogeneticAnalyzer:
"""Run alignment and maximum-likelihood phylogenetic inference.
Parameters:
mafft: MAFFT executable name or path.
iqtree: IQ-TREE 2 executable name or path. Both ``iqtree2`` and the
frequently packaged ``iqtree`` name are discovered automatically.
fasttree: FastTree executable name or path.
trimal: Optional trimAL executable name or path.
Notes:
External tools remain optional package dependencies. Methods raise an
informative :class:`RuntimeError` when the selected executable is not
available. IQ-TREE ultrafast bootstrap values should generally use at
least 1,000 replicates for a publication analysis.
"""
def __init__(
self,
*,
mafft: str = "mafft",
iqtree: str | None = None,
fasttree: str | None = None,
trimal: str = "trimal",
) -> None:
self.mafft = mafft
self.iqtree = iqtree or shutil.which("iqtree2") or shutil.which("iqtree") or "iqtree2"
self.fasttree = fasttree or shutil.which("FastTree") or shutil.which("fasttree") or "FastTree"
self.trimal = trimal
[docs]
@staticmethod
def write_fasta(records: Iterable[SequenceRecord], destination: str | Path) -> Path:
"""Write unique, non-empty sequence records as wrapped FASTA."""
output = Path(destination)
materialized = list(records)
if len(materialized) < 3:
raise ValueError("At least three sequences are required for phylogenetic inference.")
identifiers = [record.identifier for record in materialized]
if len(set(identifiers)) != len(identifiers):
raise ValueError("Sequence identifiers must be unique for phylogenetic inference.")
output.parent.mkdir(parents=True, exist_ok=True)
with output.open("w", encoding="utf-8") as handle:
for record in materialized:
sequence = "".join(record.sequence.split()).upper()
if not sequence:
raise ValueError(f"Sequence {record.identifier!r} is empty.")
handle.write(f">{record.identifier}\n")
for index in range(0, len(sequence), 80):
handle.write(sequence[index : index + 80] + "\n")
return output
@staticmethod
def _execute(command: list[str], *, stdout: Path | None = None) -> None:
"""Execute one external program and raise a diagnostic-rich error."""
try:
if stdout is None:
completed = subprocess.run(command, check=False, capture_output=True, text=True)
else:
stdout.parent.mkdir(parents=True, exist_ok=True)
with stdout.open("w", encoding="utf-8") as handle:
completed = subprocess.run(command, check=False, stdout=handle, stderr=subprocess.PIPE, text=True)
except FileNotFoundError as exc:
raise RuntimeError(f"Required executable was not found: {command[0]}") from exc
if completed.returncode != 0:
diagnostic = (completed.stderr or completed.stdout or "No diagnostic output.").strip()
raise RuntimeError(f"Command failed ({completed.returncode}): {' '.join(command)}\n{diagnostic}")
[docs]
def align(
self,
input_fasta: str | Path,
output_fasta: str | Path,
*,
method: Literal["auto", "linsi", "einsi", "fftnsi", "fftns"] = "auto",
threads: int = 1,
) -> tuple[Path, list[str]]:
"""Align sequences with MAFFT and return the output path and command."""
methods = {
"auto": ["--auto"],
"linsi": ["--localpair", "--maxiterate", "1000"],
"einsi": ["--genafpair", "--maxiterate", "1000"],
"fftnsi": ["--fftnsi"],
"fftns": ["--fftns"],
}
if method not in methods:
raise ValueError(f"Unsupported MAFFT method: {method}")
if threads < 1:
raise ValueError("threads must be at least one.")
command = [self.mafft, *methods[method], "--thread", str(threads), "--inputorder", str(input_fasta)]
output = Path(output_fasta)
self._execute(command, stdout=output)
return output, command
[docs]
def trim(
self,
alignment: str | Path,
destination: str | Path,
*,
mode: Literal["automated1", "gappyout", "strict"] = "automated1",
required: bool = False,
) -> tuple[Path, list[str] | None]:
"""Trim an alignment with trimAL, or retain it when trimming is optional."""
source = Path(alignment)
output = Path(destination)
executable = shutil.which(self.trimal)
if executable is None:
if required:
raise RuntimeError(f"Required executable was not found: {self.trimal}")
output.parent.mkdir(parents=True, exist_ok=True)
shutil.copyfile(source, output)
return output, None
command = [executable, f"-{mode}", "-in", str(source), "-out", str(output), "-fasta"]
self._execute(command)
return output, command
[docs]
def infer_tree(
self,
alignment: str | Path,
output_prefix: str | Path,
*,
method: Literal["iqtree", "fasttree"] = "iqtree",
sequence_type: Literal["nt", "aa"] = "nt",
model: str = "MFP",
bootstrap: int = 1000,
threads: int = 1,
) -> tuple[Path, str | None, list[str]]:
"""Infer an ML tree with IQ-TREE ModelFinder or approximate FastTree."""
source = Path(alignment)
prefix = Path(output_prefix)
prefix.parent.mkdir(parents=True, exist_ok=True)
if method == "iqtree":
if bootstrap < 1000:
raise ValueError("IQ-TREE ultrafast bootstrap must be at least 1,000.")
command = [
self.iqtree,
"-s",
str(source),
"--prefix",
str(prefix),
"-m",
model,
"-B",
str(bootstrap),
"-T",
str(threads),
"--redo",
]
self._execute(command)
tree = Path(f"{prefix}.treefile")
return tree, str(Path(f"{prefix}.iqtree")), command
if method != "fasttree":
raise ValueError(f"Unsupported tree method: {method}")
tree = Path(f"{prefix}.tree")
command = (
[self.fasttree, "-nt", "-gtr", str(source)]
if sequence_type == "nt"
else [self.fasttree, "-lg", str(source)]
)
self._execute(command, stdout=tree)
return tree, None, command
[docs]
def run(
self,
records: Iterable[SequenceRecord],
output_directory: str | Path,
*,
name: str = "phylogeny",
align_method: Literal["auto", "linsi", "einsi", "fftnsi", "fftns"] = "auto",
tree_method: Literal["iqtree", "fasttree"] = "iqtree",
trim: bool = True,
trim_required: bool = False,
model: str = "MFP",
bootstrap: int = 1000,
threads: int = 1,
) -> PhylogenyResult:
"""Run FASTA export, alignment, optional trimming, and tree inference."""
directory = Path(output_directory)
fasta = self.write_fasta(records, directory / f"{name}.fasta")
alignment, align_command = self.align(
fasta,
directory / f"{name}.aligned.fasta",
method=align_method,
threads=threads,
)
commands = [align_command]
analysis_alignment = alignment
trimmed: Path | None = None
if trim:
trimmed, trim_command = self.trim(
alignment,
directory / f"{name}.trimmed.fasta",
required=trim_required,
)
analysis_alignment = trimmed
if trim_command:
commands.append(trim_command)
tree, report, tree_command = self.infer_tree(
analysis_alignment,
directory / name,
method=tree_method,
model=model,
bootstrap=bootstrap,
threads=threads,
)
commands.append(tree_command)
return PhylogenyResult(
input_fasta=str(fasta),
alignment=str(alignment),
trimmed_alignment=str(trimmed) if trimmed else None,
tree=str(tree),
method=tree_method,
model=model if tree_method == "iqtree" else "GTR/CAT approximation",
commands=commands,
report=report,
)
[docs]
class TemporalAnalyzer:
"""Estimate temporal signal and deterministic dated-tree summaries.
ETE3 is imported lazily so codon-only workflows remain dependency free.
The input tree should contain branch lengths in substitutions per site.
Sampling dates are decimal calendar years keyed by exact leaf names.
"""
@staticmethod
def _import_ete3() -> Any:
"""Import ETE3 Tree with Python >= 3.13 cgi compatibility."""
import sys
import types
if "cgi" not in sys.modules:
dummy_cgi = types.ModuleType("cgi")
dummy_cgi.escape = (
lambda text, quote=True: text.replace("&", "&").replace("<", "<").replace(">", ">")
)
sys.modules["cgi"] = dummy_cgi
from ete3 import Tree
return Tree
[docs]
@classmethod
def load_tree(cls, tree: str | Path | Any, *, midpoint_root: bool = True) -> Any:
"""Load or copy an ETE3 tree and optionally apply midpoint rooting."""
try:
Tree = cls._import_ete3()
except (ImportError, ModuleNotFoundError) as exc:
raise RuntimeError("Temporal phylogenetics requires ete3; install the 'phylo' extra.") from exc
tree_object: Any = tree
if hasattr(tree_object, "traverse"):
loaded = tree_object.copy(method="deepcopy")
else:
try:
loaded = Tree(str(tree), format=0)
except Exception:
loaded = Tree(str(tree), format=1)
if midpoint_root and len(loaded) > 2:
outgroup = loaded.get_midpoint_outgroup()
if outgroup is not None:
loaded.set_outgroup(outgroup)
return loaded
[docs]
def temporal_signal(
self,
tree: str | Path | Any,
sample_dates: Mapping[str, float | int],
*,
strata: Mapping[str, Any] | None = None,
midpoint_root: bool = True,
optimize_root: bool = False,
permutations: int = 0,
random_seed: int = 1,
) -> TemporalSignalResult:
"""Regress root-to-tip divergence on sampling date.
A positive slope is interpreted as an exploratory strict-clock rate.
When ``permutations`` is positive, sampling dates are randomized and
the complete root-selection procedure is repeated for every replicate.
When ``strata`` is provided, dates are permuted strictly within each
stratum (e.g. genotype or lineage) to avoid confounding between group
membership and sampling date. The resulting empirical p value therefore
accounts for optimization over candidate roots. At least three dated tips
and two distinct sampling years are required.
"""
loaded = self.load_tree(tree, midpoint_root=False)
usable = {
leaf.name: float(sample_dates[leaf.name])
for leaf in loaded.iter_leaves()
if leaf.name in sample_dates and math.isfinite(float(sample_dates[leaf.name]))
}
if len(usable) < 3 or len(set(usable.values())) < 2:
raise ValueError("Temporal analysis requires at least three tips across two sampling dates.")
if permutations < 0:
raise ValueError("permutations cannot be negative.")
try:
from scipy.stats import linregress
except ImportError as exc:
raise RuntimeError("Temporal regression requires scipy; install the 'analysis' extra.") from exc
candidates: list[tuple[str, list[str] | None, Any]] = [("current", None, loaded)]
if midpoint_root and len(loaded) > 2:
midpoint_tree = loaded.copy(method="deepcopy")
outgroup = midpoint_tree.get_midpoint_outgroup()
if outgroup is not None:
leaves = [leaf.name for leaf in outgroup.iter_leaves()]
midpoint_tree.set_outgroup(outgroup)
candidates.append(("midpoint", leaves, midpoint_tree))
if optimize_root:
for index, node in enumerate(loaded.traverse("preorder")):
if node.is_root():
continue
leaf_names = [leaf.name for leaf in node.iter_leaves()]
candidate_tree = loaded.copy(method="deepcopy")
candidate = (
candidate_tree & leaf_names[0]
if len(leaf_names) == 1
else candidate_tree.get_common_ancestor(*leaf_names)
)
candidate_tree.set_outgroup(candidate)
candidates.append((f"optimized_branch_{index}", leaf_names, candidate_tree))
evaluated: list[tuple[float, str, list[str] | None, Any, Any, dict[str, float]]] = []
for rooting, outgroup_leaves, candidate_tree in candidates:
distances = {name: float(candidate_tree.get_distance(name)) for name in usable}
regression = linregress(list(usable.values()), [distances[name] for name in usable])
if not math.isfinite(regression.slope) or regression.slope <= 0:
continue
root_date = -float(regression.intercept) / float(regression.slope)
if root_date > min(usable.values()):
continue
evaluated.append(
(float(regression.rvalue**2), rooting, outgroup_leaves, candidate_tree, regression, distances)
)
if not evaluated:
raise ValueError("No tested root produced a positive, temporally valid root-to-tip regression.")
_, rooting, outgroup_leaves, loaded, regression, distances = max(evaluated, key=lambda item: item[0])
if not math.isfinite(regression.slope) or regression.slope <= 0:
raise ValueError("The rooted tree has no positive temporal signal; a time scale cannot be estimated.")
root_date = -float(regression.intercept) / float(regression.slope)
fitted = {name: float(regression.intercept + regression.slope * date) for name, date in usable.items()}
residuals = {name: distances[name] - fitted[name] for name in usable}
permutation_p_value: float | None = None
if permutations:
observed_r_squared = float(regression.rvalue**2)
ordered_names = list(usable)
observed_dates = [usable[name] for name in ordered_names]
candidate_distances = [
[candidate_distances[name] for name in ordered_names]
for _, _, _, _, _, candidate_distances in evaluated
]
generator = random.Random(random_seed)
exceedances = 0
stratum_indices: dict[Any, list[int]] = defaultdict(list)
if strata:
for idx, name in enumerate(ordered_names):
stratum_indices[strata.get(name, "default")].append(idx)
for _ in range(permutations):
shuffled_dates = observed_dates.copy()
if strata:
for group_idxs in stratum_indices.values():
sub_dates = [shuffled_dates[i] for i in group_idxs]
generator.shuffle(sub_dates)
for i, d in zip(group_idxs, sub_dates, strict=True):
shuffled_dates[i] = d
else:
generator.shuffle(shuffled_dates)
best_null = 0.0
for candidate_values in candidate_distances:
null_regression = linregress(shuffled_dates, candidate_values)
if not math.isfinite(null_regression.slope) or null_regression.slope <= 0:
continue
null_root_date = -float(null_regression.intercept) / float(null_regression.slope)
if null_root_date > min(shuffled_dates):
continue
best_null = max(best_null, float(null_regression.rvalue**2))
if best_null >= observed_r_squared:
exceedances += 1
permutation_p_value = (exceedances + 1) / (permutations + 1)
reported_p_value = permutation_p_value if permutation_p_value is not None else float(regression.pvalue)
return TemporalSignalResult(
rate=float(regression.slope),
intercept=float(regression.intercept),
root_date=root_date,
r_value=float(regression.rvalue),
r_squared=float(regression.rvalue**2),
p_value=reported_p_value,
p_value_kind=(
"stratified_date_randomization"
if permutation_p_value is not None and strata
else "root_optimized_date_randomization"
if permutation_p_value is not None
else "unadjusted_regression"
),
regression_p_value_unadjusted=float(regression.pvalue),
permutation_p_value=permutation_p_value,
permutations=permutations,
random_seed=random_seed if permutations else None,
standard_error=float(regression.stderr),
sample_dates=usable,
root_to_tip=distances,
fitted_distances=fitted,
residuals=residuals,
rooting=rooting,
outgroup_leaves=outgroup_leaves,
)
[docs]
def date_tree(
self,
tree: str | Path | Any,
signal: TemporalSignalResult,
*,
midpoint_root: bool = True,
) -> Any:
"""Attach deterministic calendar dates to every ETE3 node.
Internal dates follow the fitted strict-clock transform of root
distance. Tips retain an additional ``observed_date`` feature when a
sampling date is available. The returned topology is suitable for
exploratory time-tree plotting and lineage-through-time summaries.
"""
loaded = self.load_tree(tree, midpoint_root=False)
if signal.outgroup_leaves:
outgroup = (
loaded & signal.outgroup_leaves[0]
if len(signal.outgroup_leaves) == 1
else loaded.get_common_ancestor(*signal.outgroup_leaves)
)
loaded.set_outgroup(outgroup)
elif midpoint_root:
outgroup = loaded.get_midpoint_outgroup()
if outgroup is not None:
loaded.set_outgroup(outgroup)
for node in loaded.traverse("preorder"):
distance = float(loaded.get_distance(node))
fitted_date = signal.root_date + distance / signal.rate
node.add_feature("calendar_date", fitted_date)
if node.is_leaf() and node.name in signal.sample_dates:
node.add_feature("observed_date", signal.sample_dates[node.name])
node.add_feature("fitted_date", fitted_date)
node.calendar_date = signal.sample_dates[node.name]
epsilon = 1e-6
for node in loaded.traverse("postorder"):
if not node.is_leaf():
latest_allowed = min(float(child.calendar_date) for child in node.children) - epsilon
node.calendar_date = min(float(node.calendar_date), latest_allowed)
return loaded
[docs]
@staticmethod
def lineage_through_time(tree: Any, *, points: int = 200) -> dict[str, list[float]]:
"""Count branches crossing an evenly spaced calendar-time grid."""
if points < 2:
raise ValueError("points must be at least two.")
nodes = [node for node in tree.traverse() if not node.is_root()]
if not nodes or any(not hasattr(node, "calendar_date") for node in tree.traverse()):
raise ValueError("Every tree node must contain a calendar_date feature.")
start = float(tree.calendar_date)
end = max(float(node.calendar_date) for node in tree.iter_leaves())
if end <= start:
raise ValueError("Dated tree must span a positive calendar interval.")
step = (end - start) / (points - 1)
times = [start + index * step for index in range(points)]
lineages = []
for time in times:
count = sum(
1
for node in nodes
if float(node.up.calendar_date) <= time <= float(node.calendar_date)
)
lineages.append(float(count))
return {"time": times, "lineages": lineages}
[docs]
class PhylogeneticComparator:
"""Compare whole-genome and gene-tree topologies on their shared tips."""
[docs]
def robinson_foulds(self, first: str | Path | Any, second: str | Path | Any) -> dict[str, Any]:
"""Calculate unrooted Robinson-Foulds distance after shared-tip pruning.
The normalized distance is zero for identical shared-tip bipartitions
and approaches one as topologies become maximally discordant. At least
four shared tips are required for an informative unrooted comparison.
"""
loader = TemporalAnalyzer()
first_tree = loader.load_tree(first, midpoint_root=False)
second_tree = loader.load_tree(second, midpoint_root=False)
shared = sorted(set(first_tree.get_leaf_names()) & set(second_tree.get_leaf_names()))
if len(shared) < 4:
raise ValueError("Robinson-Foulds comparison requires at least four shared tips.")
first_tree.prune(shared, preserve_branch_length=True)
second_tree.prune(shared, preserve_branch_length=True)
comparison = first_tree.robinson_foulds(second_tree, unrooted_trees=True)
distance = int(comparison[0])
maximum = int(comparison[1])
return {
"shared_tips": len(shared),
"tip_names": shared,
"distance": distance,
"maximum_distance": maximum,
"normalized_distance": distance / maximum if maximum else 0.0,
}