"""
CodonAdaptPy Publication Plots
This module creates optional self-contained interactive Plotly figures and
publication-ready static Matplotlib figures for RSCU heatmaps, ENC-GC3
diagnostics, neutrality and PR2 plots, group distributions, ordinations,
dinucleotide profiles, correlation matrices, PCA, preferred-codon
intersections, host-adaptation landscapes, and sliding-window CAI profiles.
Classes:
- PlotManager: Optional self-contained interactive HTML figures.
- PublicationPlotManager: Colorblind-safe static scientific figures.
:Created: July 20, 2026
:Updated: July 20, 2026
:Author: Naveen Duhan
:Version: 1.0.2
"""
from __future__ import annotations
import math
from collections import defaultdict
from pathlib import Path
from typing import Any
from .exceptions import OptionalDependencyError
from .models import AnalysisResult
[docs]
class PlotManager:
"""Build optional interactive Plotly figures for self-contained HTML."""
def _plotly(self) -> tuple[Any, Any]:
"""Import Plotly lazily so static workflows remain dependency-light."""
try:
import plotly.express as px
import plotly.graph_objects as go
except ImportError as exc:
raise OptionalDependencyError(
"Interactive HTML plots require CodonAdaptPy[interactive]."
) from exc
return px, go
[docs]
def rscu_heatmap(self, results: list[AnalysisResult]) -> Any:
"""Create an interactive sequence-by-codon RSCU heatmap."""
_, go = self._plotly()
from .genetic_code import GeneticCode
code = GeneticCode.from_ncbi(1)
codons = sorted(
codon
for amino_acid, family in code.synonymous_codons.items()
if amino_acid != "*" and len(family) > 1
for codon in family
)
return go.Figure(
data=go.Heatmap(
z=[[result.rscu.get(codon, 0.0) for codon in codons] for result in results],
x=codons,
y=[result.identifier for result in results],
colorscale="Viridis",
)
).update_layout(title="Relative Synonymous Codon Usage", xaxis_title="Codon", yaxis_title="Sequence")
[docs]
def enc_gc3(self, results: list[AnalysisResult]) -> Any:
"""Create interactive observed ENC versus GC3 and expectation traces."""
_, go = self._plotly()
from .metrics.evolutionary import EvolutionaryMetrics
grid = [index / 100 for index in range(101)]
figure = go.Figure()
figure.add_scatter(
x=[result.metrics["gc3"] for result in results],
y=[result.metrics["enc"] for result in results],
mode="markers+text",
text=[result.identifier for result in results],
name="Observed",
)
figure.add_scatter(
x=grid,
y=[EvolutionaryMetrics.expected_enc(value) for value in grid],
mode="lines",
name="Mutation-only expectation",
)
return figure.update_layout(title="ENC-GC3 Plot", xaxis_title="GC3", yaxis_title="Effective Number of Codons")
[docs]
def sliding_cai(self, result: AnalysisResult) -> Any:
"""Create an interactive sliding-window CAI line profile."""
_, go = self._plotly()
values = result.sliding_windows.get("cai", [])
return go.Figure(
data=go.Scatter(
x=[item["start_codon"] for item in values],
y=[item["cai"] for item in values],
mode="lines",
)
).update_layout(
title=f"Sliding-window CAI: {result.identifier}",
xaxis_title="Starting codon",
yaxis_title="CAI",
)
[docs]
def neutrality(self, results: list[AnalysisResult]) -> Any:
"""Create an interactive GC12-versus-GC3 neutrality scatter."""
px, _ = self._plotly()
from .metrics.evolutionary import EvolutionaryMetrics
rows = [
{
"identifier": result.identifier,
"GC3": result.metrics["gc3"],
"GC12": (result.metrics["gc1"] + result.metrics["gc2"]) / 2,
}
for result in results
]
figure = px.scatter(rows, x="GC3", y="GC12", text="identifier", title="Neutrality Plot")
if len(results) >= 2 and len({row["GC3"] for row in rows}) > 1:
regression = EvolutionaryMetrics.neutrality(
[result.metrics["gc1"] for result in results],
[result.metrics["gc2"] for result in results],
[result.metrics["gc3"] for result in results],
)
minimum = min(row["GC3"] for row in rows)
maximum = max(row["GC3"] for row in rows)
figure.add_scatter(
x=[minimum, maximum],
y=[
regression["intercept"] + regression["slope"] * minimum,
regression["intercept"] + regression["slope"] * maximum,
],
mode="lines",
name=f"Linear fit (R²={regression['r_squared']:.3f})",
)
return figure
[docs]
def parity_rule_2(self, results: list[AnalysisResult]) -> Any:
"""Create an interactive PR2 plot with equal-bias references."""
px, _ = self._plotly()
rows = [
{
"identifier": result.identifier,
"A3/(A3+T3)": result.metrics["at_bias"],
"G3/(G3+C3)": result.metrics["gc_bias"],
}
for result in results
]
figure = px.scatter(
rows,
x="A3/(A3+T3)",
y="G3/(G3+C3)",
text="identifier",
title="Parity Rule 2 Plot",
)
figure.add_hline(y=0.5, line_dash="dash")
figure.add_vline(x=0.5, line_dash="dash")
return figure
[docs]
def save(self, figure: Any, path: str | Path) -> Path:
"""Write a self-contained HTML figure without Kaleido or Chrome."""
destination = Path(path)
if destination.suffix.lower() != ".html":
raise ValueError("Interactive PlotManager output must use HTML; use PublicationPlotManager for images.")
destination.parent.mkdir(parents=True, exist_ok=True)
figure.write_html(destination, include_plotlyjs=True)
return destination
[docs]
class PublicationPlotManager:
"""Create journal-ready static figures from CodonAdaptPy outputs.
The manager uses the Okabe-Ito colorblind-safe palette, redundant marker
shapes, outward-offset axes, visible individual observations, and vector-
compatible Matplotlib artists. Matplotlib is imported lazily so sequence
analysis remains dependency-free.
"""
COLORS = ("#0072B2", "#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9")
MARKERS = ("o", "s", "^", "D", "P", "X")
METRIC_LABELS = {
"gc": "GC",
"gc1": "GC1",
"gc2": "GC2",
"gc3": "GC3",
"gc3s": "GC3s",
"at1": "AT1",
"at2": "AT2",
"at3": "AT3",
"enc": "ENC",
"cpg_representation": "CpG O/E",
"upa_representation": "UpA O/E",
}
def _matplotlib(self) -> tuple[Any, Any]:
"""Import Matplotlib lazily and apply consistent publication styling."""
try:
import matplotlib.pyplot as plt
import numpy as np
except ImportError as exc:
raise OptionalDependencyError(
"Publication plots require CodonAdaptPy[plot] with Matplotlib."
) from exc
plt.rcParams.update(
{
"font.family": "sans-serif",
"font.sans-serif": ["Arial", "Helvetica", "DejaVu Sans"],
"font.size": 8,
"axes.labelsize": 9,
"axes.titlesize": 10,
"xtick.labelsize": 7,
"ytick.labelsize": 7,
"legend.fontsize": 7,
"pdf.fonttype": 42,
"svg.fonttype": "none",
"figure.max_open_warning": 50,
}
)
return plt, np
@staticmethod
def _clean_axes(axis: Any) -> None:
"""Remove chart junk and separate the visible left and bottom spines."""
axis.spines["top"].set_visible(False)
axis.spines["right"].set_visible(False)
axis.spines["left"].set_position(("outward", 5))
axis.spines["bottom"].set_position(("outward", 5))
axis.tick_params(direction="out", length=3, width=0.7)
[docs]
def ordination(
self,
scores: list[list[float]],
labels: list[str],
groups: list[str],
*,
x_label: str = "Axis 1",
y_label: str = "Axis 2",
title: str = "Ordination",
) -> Any:
"""Plot grouped two-dimensional scores with centroid standard errors.
Individual samples are shown with color and shape encoding. Larger
black-edged points show group centroids, and error bars represent the
standard error independently along both axes.
"""
if not (len(scores) == len(labels) == len(groups)) or not scores:
raise ValueError("scores, labels, and groups must have equal non-zero lengths.")
if any(len(score) < 2 for score in scores):
raise ValueError("Each ordination score must contain at least two axes.")
plt, np = self._matplotlib()
figure, axis = plt.subplots(figsize=(7.2, 5.0), constrained_layout=True)
ordered_groups = sorted({str(group) for group in groups})
for index, group in enumerate(ordered_groups):
positions = [position for position, value in enumerate(groups) if str(value) == group]
values = np.asarray([scores[position][:2] for position in positions], dtype=float)
color = self.COLORS[index % len(self.COLORS)]
marker = self.MARKERS[index % len(self.MARKERS)]
axis.scatter(
values[:, 0],
values[:, 1],
s=28,
color=color,
marker=marker,
alpha=0.78,
edgecolor="white",
linewidth=0.45,
label=f"{group} (n={len(values)})",
)
centroid = values.mean(axis=0)
standard_error = values.std(axis=0, ddof=1) / math.sqrt(len(values)) if len(values) > 1 else [0, 0]
axis.errorbar(
centroid[0],
centroid[1],
xerr=standard_error[0],
yerr=standard_error[1],
fmt=marker,
markersize=7,
markerfacecolor=color,
markeredgecolor="black",
markeredgewidth=0.8,
ecolor=color,
elinewidth=1.0,
capsize=2.5,
zorder=4,
)
axis.axhline(0, color="#BDBDBD", linewidth=0.6, zorder=0)
axis.axvline(0, color="#BDBDBD", linewidth=0.6, zorder=0)
axis.set(xlabel=x_label, ylabel=y_label, title=title)
axis.legend(frameon=False, title="Genotype")
self._clean_axes(axis)
return figure
[docs]
def metric_panels(
self,
rows: list[dict[str, Any]],
*,
metrics: list[str],
group: str,
title: str = "Metric distributions",
columns: int = 2,
) -> Any:
"""Create box-and-point panels for numeric metrics grouped by metadata.
Boxes show the median and interquartile range; whiskers follow the
conventional 1.5-IQR rule. Every observation is overlaid using
deterministic horizontal jitter, and each category label includes its
sample size.
"""
if not rows or not metrics:
raise ValueError("rows and metrics must be non-empty.")
plt, np = self._matplotlib()
group_names = sorted({str(row[group]) for row in rows})
panel_columns = min(columns, len(metrics))
panel_rows = math.ceil(len(metrics) / panel_columns)
figure, axes = plt.subplots(
panel_rows,
panel_columns,
figsize=(7.2, 2.9 * panel_rows),
constrained_layout=True,
squeeze=False,
)
random = np.random.default_rng(1)
for panel_index, metric in enumerate(metrics):
axis = axes.flat[panel_index]
values = [
[float(row[metric]) for row in rows if str(row[group]) == name and row.get(metric) is not None]
for name in group_names
]
boxes = axis.boxplot(values, patch_artist=True, widths=0.55, showfliers=False)
for index, patch in enumerate(boxes["boxes"]):
patch.set_facecolor(self.COLORS[index % len(self.COLORS)])
patch.set_alpha(0.34)
patch.set_edgecolor(self.COLORS[index % len(self.COLORS)])
for index, observations in enumerate(values, start=1):
jitter = random.uniform(-0.12, 0.12, len(observations))
axis.scatter(
np.full(len(observations), index) + jitter,
observations,
s=15,
color=self.COLORS[(index - 1) % len(self.COLORS)],
alpha=0.68,
edgecolor="white",
linewidth=0.3,
)
axis.set_xticks(
range(1, len(group_names) + 1),
[f"{name}\n(n={len(values[index])})" for index, name in enumerate(group_names)],
)
axis.set_xlabel(group)
display_metric = self.METRIC_LABELS.get(metric, metric.replace("_", " "))
axis.set_ylabel(display_metric)
axis.set_title(display_metric)
self._clean_axes(axis)
for axis in axes.flat[len(metrics) :]:
axis.set_visible(False)
figure.suptitle(title, fontsize=11, fontweight="bold")
return figure
[docs]
def codon_preference_intersections(
self,
results: list[AnalysisResult],
*,
group_key: str,
threshold: float = 1.6,
) -> dict[str, Any]:
"""Return exact overlaps of preferred codons across metadata groups.
A codon belongs to a group's preferred set when its group-mean RSCU is
greater than or equal to ``threshold``. Each intersection is exact: a
codon assigned to groups A and B is absent from every other displayed
group. The structured result supports auditable tables as well as the
UpSet-style visualization.
"""
if not results:
raise ValueError("results must be non-empty.")
if threshold <= 0:
raise ValueError("threshold must be positive.")
grouped: dict[str, list[AnalysisResult]] = defaultdict(list)
for result in results:
grouped[str(result.metadata.get(group_key, "Unknown"))].append(result)
group_names = sorted(grouped)
codons = sorted({codon for result in results for codon in result.rscu})
preferred: dict[str, set[str]] = {}
for group in group_names:
preferred[group] = {
codon
for codon in codons
if sum(result.rscu.get(codon, 0.0) for result in grouped[group]) / len(grouped[group]) >= threshold
}
memberships: dict[tuple[str, ...], list[str]] = defaultdict(list)
for codon in codons:
membership = tuple(group for group in group_names if codon in preferred[group])
if membership:
memberships[membership].append(codon)
if not memberships:
raise ValueError(f"No group-mean RSCU values met the preference threshold {threshold:g}.")
ordered = sorted(
memberships.items(),
key=lambda item: (-len(item[1]), -len(item[0]), item[0], item[1]),
)
return {
"group_key": group_key,
"rscu_threshold": threshold,
"groups": [
{
"name": group,
"sample_count": len(grouped[group]),
"preferred_codon_count": len(preferred[group]),
}
for group in group_names
],
"intersections": [
{
"rank": rank,
"groups": list(membership),
"size": len(members),
"codons": members,
}
for rank, (membership, members) in enumerate(ordered, start=1)
],
}
[docs]
def codon_preference_upset(
self,
results: list[AnalysisResult],
*,
group_key: str,
threshold: float = 1.6,
max_intersections: int = 20,
title: str = "Shared preferred-codon signatures",
) -> Any:
"""Plot exact preferred-codon intersections in an UpSet-style layout.
The upper bars report the number of codons in each exact intersection,
the matrix identifies participating metadata groups, and the left bars
report each group's total preferred-codon set size. This avoids the
area-comparison problems of multi-set Venn diagrams.
"""
if max_intersections < 1:
raise ValueError("max_intersections must be positive.")
data = self.codon_preference_intersections(results, group_key=group_key, threshold=threshold)
intersections = data["intersections"][:max_intersections]
groups = data["groups"]
plt, np = self._matplotlib()
width = max(7.2, 2.5 + 0.48 * len(intersections))
height = max(4.6, 2.8 + 0.38 * len(groups))
figure = plt.figure(figsize=(width, height), constrained_layout=True)
grid = figure.add_gridspec(
2,
2,
width_ratios=(1.55, 4.8),
height_ratios=(2.35, max(1.4, 0.42 * len(groups))),
)
definition_axis = figure.add_subplot(grid[0, 0])
size_axis = figure.add_subplot(grid[0, 1])
set_axis = figure.add_subplot(grid[1, 0])
matrix_axis = figure.add_subplot(grid[1, 1], sharex=size_axis)
x_positions = np.arange(len(intersections))
sizes = [int(item["size"]) for item in intersections]
bars = size_axis.bar(x_positions, sizes, color=self.COLORS[0], width=0.68)
size_axis.bar_label(bars, labels=[str(value) for value in sizes], padding=2, fontsize=7)
size_axis.set_ylabel("Intersection size\n(codons)")
size_axis.set_title("Exact preferred-codon intersections")
size_axis.set_xticks([])
size_axis.set_ylim(0, max(sizes) * 1.22 + 0.2)
self._clean_axes(size_axis)
definition_axis.axis("off")
definition_axis.text(
0.0,
0.94,
f"Preferred codon\nmean RSCU ≥ {threshold:g}",
ha="left",
va="top",
fontsize=9,
fontweight="bold",
transform=definition_axis.transAxes,
)
definition_axis.text(
0.0,
0.47,
f"Grouping field\n{group_key}",
ha="left",
va="top",
fontsize=8,
color="#4D4D4D",
transform=definition_axis.transAxes,
)
y_positions = np.arange(len(groups))
set_sizes = [int(group["preferred_codon_count"]) for group in groups]
set_axis.barh(y_positions, set_sizes, color="#8C8C8C", height=0.56)
set_axis.set_yticks([])
set_axis.set_xlabel("Set size")
set_axis.set_xlim(max(set_sizes) * 1.18 + 0.2, 0)
set_axis.set_ylim(len(groups) - 0.5, -0.5)
self._clean_axes(set_axis)
group_names = [str(group["name"]) for group in groups]
for column, intersection in enumerate(intersections):
selected = {str(group) for group in intersection["groups"]}
included_rows = []
for row, group in enumerate(group_names):
included = group in selected
matrix_axis.scatter(
column,
row,
s=35 if included else 15,
color="#202020" if included else "#D9D9D9",
zorder=3,
)
if included:
included_rows.append(row)
if len(included_rows) > 1:
matrix_axis.plot(
[column, column],
[min(included_rows), max(included_rows)],
color="#202020",
linewidth=1.2,
zorder=2,
)
matrix_axis.set_yticks(
y_positions,
[
f"{group['name']} (n={group['sample_count']})"
for group in groups
],
)
matrix_axis.set_xticks(x_positions, [str(item["rank"]) for item in intersections])
matrix_axis.set_xlabel("Intersection rank (see companion TSV for codons)")
matrix_axis.set_ylim(len(groups) - 0.5, -0.5)
matrix_axis.spines["top"].set_visible(False)
matrix_axis.spines["right"].set_visible(False)
matrix_axis.spines["left"].set_visible(False)
matrix_axis.tick_params(axis="y", length=0, pad=7)
matrix_axis.tick_params(axis="x", direction="out", length=3)
for row in range(len(groups)):
matrix_axis.axhline(row, color="#EFEFEF", linewidth=0.5, zorder=0)
figure.suptitle(title, fontsize=11, fontweight="bold")
figure.text(
0.995,
0.008,
f"Showing {len(intersections)} of {len(data['intersections'])} non-empty exact intersections",
ha="right",
va="bottom",
fontsize=6.5,
color="#4D4D4D",
)
return figure
[docs]
def rscu_similarity_heatmap(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "RSCU signature similarity",
) -> Any:
"""Plot cosine similarity among group-mean 59-codon RSCU profiles."""
if not results:
raise ValueError("results must be non-empty.")
plt, np = self._matplotlib()
from .genetic_code import GeneticCode
code = GeneticCode.from_ncbi(1)
codons = sorted(
codon
for amino_acid, family in code.synonymous_codons.items()
if amino_acid != "*" and len(family) > 1
for codon in family
)
grouped: dict[str, list[AnalysisResult]] = defaultdict(list)
for result in results:
grouped[str(result.metadata.get(group_key, "Unknown"))].append(result)
groups = sorted(grouped)
profiles = np.asarray(
[
[np.mean([result.rscu.get(codon, 0.0) for result in grouped[group]]) for codon in codons]
for group in groups
],
dtype=float,
)
norms = np.linalg.norm(profiles, axis=1)
denominator = np.outer(norms, norms)
similarity = np.divide(
profiles @ profiles.T,
denominator,
out=np.zeros((len(groups), len(groups)), dtype=float),
where=denominator > 0,
)
labels = [f"{group} (n={len(grouped[group])})" for group in groups]
side = max(4.2, 0.62 * len(groups) + 2.2)
figure, axis = plt.subplots(figsize=(side, side), constrained_layout=True)
image = axis.imshow(similarity, cmap="cividis", vmin=0.0, vmax=1.0)
axis.set_xticks(range(len(groups)), labels, rotation=35, ha="right")
axis.set_yticks(range(len(groups)), labels)
axis.set_xlabel("Metadata group")
axis.set_ylabel("Metadata group")
axis.set_title(title)
for row in range(len(groups)):
for column in range(len(groups)):
value = float(similarity[row, column])
axis.text(
column,
row,
f"{value:.3f}",
ha="center",
va="center",
fontsize=7,
color="white" if value < 0.42 else "black",
)
colorbar = figure.colorbar(image, ax=axis, fraction=0.045, pad=0.03)
colorbar.set_label("Cosine similarity")
return figure
[docs]
def rscu_pca_biplot(
self,
results: list[AnalysisResult],
*,
group_key: str,
max_loadings: int = 10,
title: str = "RSCU principal-component biplot",
) -> Any:
"""Plot standardized 59-codon RSCU scores and influential codons."""
if len(results) < 3:
raise ValueError("RSCU PCA requires at least three sequences.")
if max_loadings < 1:
raise ValueError("max_loadings must be positive.")
plt, np = self._matplotlib()
from .genetic_code import GeneticCode
code = GeneticCode.from_ncbi(1)
codons = sorted(
codon
for amino_acid, family in code.synonymous_codons.items()
if amino_acid != "*" and len(family) > 1
for codon in family
)
matrix = np.asarray([[result.rscu.get(codon, 0.0) for codon in codons] for result in results], dtype=float)
standard_deviation = matrix.std(axis=0, ddof=1)
retained = standard_deviation > 0
if int(retained.sum()) < 2:
raise ValueError("RSCU PCA requires variation in at least two codons.")
retained_codons = [codon for codon, keep in zip(codons, retained, strict=True) if keep]
standardized = (matrix[:, retained] - matrix[:, retained].mean(axis=0)) / standard_deviation[retained]
left, singular, right = np.linalg.svd(standardized, full_matrices=False)
if len(singular) < 2:
raise ValueError("RSCU PCA could not estimate two components.")
scores = left[:, :2] * singular[:2]
loadings = right[:2, :].T
variance = singular**2
explained = variance[:2] / variance.sum() if variance.sum() else np.zeros(2)
figure, axis = plt.subplots(figsize=(7.2, 5.4), constrained_layout=True)
groups = sorted({str(result.metadata.get(group_key, "Unknown")) for result in results})
for index, group in enumerate(groups):
selected = [
position
for position, result in enumerate(results)
if str(result.metadata.get(group_key, "Unknown")) == group
]
axis.scatter(
scores[selected, 0],
scores[selected, 1],
color=self.COLORS[index % len(self.COLORS)],
marker=self.MARKERS[index % len(self.MARKERS)],
s=30,
alpha=0.78,
edgecolor="white",
linewidth=0.4,
label=f"{group} (n={len(selected)})",
)
if len(results) <= 30:
for position, result in enumerate(results):
axis.annotate(
result.identifier,
(scores[position, 0], scores[position, 1]),
xytext=(3, 3),
textcoords="offset points",
fontsize=5.8,
color="#4D4D4D",
)
loading_strength = np.linalg.norm(loadings, axis=1)
selected_loadings = np.argsort(loading_strength)[-min(max_loadings, len(loadings)) :]
score_radius = max(float(np.max(np.abs(scores))), 1.0)
loading_radius = max(float(np.max(loading_strength[selected_loadings])), 1e-12)
scale = 0.72 * score_radius / loading_radius
for position in selected_loadings:
x_value, y_value = loadings[position] * scale
axis.annotate(
"",
xy=(x_value, y_value),
xytext=(0, 0),
arrowprops={"arrowstyle": "->", "color": "#555555", "linewidth": 0.75, "alpha": 0.8},
)
axis.text(x_value * 1.05, y_value * 1.05, retained_codons[position], fontsize=6.2, ha="center")
axis.axhline(0, color="#BDBDBD", linewidth=0.6)
axis.axvline(0, color="#BDBDBD", linewidth=0.6)
axis.set_xlabel(f"PC1 ({explained[0] * 100:.1f}% variance)")
axis.set_ylabel(f"PC2 ({explained[1] * 100:.1f}% variance)")
axis.set_title(title)
axis.legend(frameon=False, title=group_key)
self._clean_axes(axis)
return figure
[docs]
def sliding_cai(
self,
result: AnalysisResult,
*,
title: str | None = None,
) -> Any:
"""Plot a reference-dependent sliding-window CAI profile."""
values = result.sliding_windows.get("cai", [])
if not values:
raise ValueError("Sliding-window CAI requires a target reference profile.")
plt, _ = self._matplotlib()
figure, axis = plt.subplots(figsize=(7.2, 3.8), constrained_layout=True)
axis.plot(
[item["start_codon"] for item in values],
[item["cai"] for item in values],
color=self.COLORS[0],
linewidth=1.3,
marker="o" if len(values) <= 40 else None,
markersize=3,
)
axis.set(
xlabel="Starting codon",
ylabel="CAI",
title=title or f"Sliding-window CAI: {result.identifier}",
)
axis.set_ylim(0, 1.03)
self._clean_axes(axis)
return figure
[docs]
def metric_correlation_heatmap(
self,
results: list[AnalysisResult],
*,
metrics: list[str] | None = None,
title: str = "Codon-analysis metric correlations",
) -> Any:
"""Plot Pearson correlations among variable scalar analysis metrics."""
if len(results) < 3:
raise ValueError("Metric correlations require at least three sequences.")
plt, np = self._matplotlib()
requested = metrics or [
"enc",
"gc3s",
"gc3",
"at3",
"cpg_representation",
"upa_representation",
"codon_pair_bias",
"cai",
"rcdi",
]
retained = [
metric
for metric in requested
if all(isinstance(result.metrics.get(metric), (int, float)) for result in results)
and len({float(result.metrics[metric]) for result in results}) > 1
]
if len(retained) < 2:
raise ValueError("Metric correlations require at least two variable numeric metrics.")
matrix = np.asarray([[float(result.metrics[metric]) for metric in retained] for result in results])
correlation = np.corrcoef(matrix, rowvar=False)
labels = [self.METRIC_LABELS.get(metric, metric.replace("_", " ")) for metric in retained]
side = max(5.0, 0.62 * len(retained) + 2.1)
figure, axis = plt.subplots(figsize=(side, side), constrained_layout=True)
image = axis.imshow(correlation, cmap="RdBu_r", vmin=-1.0, vmax=1.0)
axis.set_xticks(range(len(retained)), labels, rotation=40, ha="right")
axis.set_yticks(range(len(retained)), labels)
axis.set_title(title)
for row in range(len(retained)):
for column in range(len(retained)):
value = float(correlation[row, column])
axis.text(
column,
row,
f"{value:.2f}",
ha="center",
va="center",
fontsize=6.5,
color="white" if abs(value) > 0.62 else "black",
)
colorbar = figure.colorbar(image, ax=axis, fraction=0.045, pad=0.03)
colorbar.set_label("Pearson correlation")
return figure
[docs]
def host_adaptation_landscape(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "Multi-host codon-adaptation landscape",
) -> Any:
"""Plot CAI against RCDI for every host, with SiD encoded by area."""
if not results or not any(result.host_comparisons for result in results):
raise ValueError("Host-adaptation landscapes require at least one host reference.")
plt, np = self._matplotlib()
hosts = sorted({host for result in results for host in result.host_comparisons})
columns = min(2, len(hosts))
rows = math.ceil(len(hosts) / columns)
figure, axes = plt.subplots(
rows,
columns,
figsize=(5.0 * columns, 4.0 * rows),
constrained_layout=False,
squeeze=False,
)
figure.subplots_adjust(left=0.08, right=0.985, top=0.84, bottom=0.25, hspace=0.42, wspace=0.28)
groups = sorted({str(result.metadata.get(group_key, "Unknown")) for result in results})
for panel, host in enumerate(hosts):
axis = axes.flat[panel]
for index, group in enumerate(groups):
selected = [
result
for result in results
if str(result.metadata.get(group_key, "Unknown")) == group
and host in result.host_comparisons
]
if not selected:
continue
comparisons = [result.host_comparisons[host] for result in selected]
axis.scatter(
[comparison["rcdi"] for comparison in comparisons],
[comparison["cai"] for comparison in comparisons],
s=[25 + 105 * comparison["sid_similarity"] for comparison in comparisons],
color=self.COLORS[index % len(self.COLORS)],
marker=self.MARKERS[index % len(self.MARKERS)],
alpha=0.68,
edgecolor="white",
linewidth=0.45,
label=f"{group} (n={len(selected)})",
)
if len(selected) <= 15:
for result, comparison in zip(selected, comparisons, strict=True):
axis.annotate(
result.identifier,
(comparison["rcdi"], comparison["cai"]),
xytext=(3, 3),
textcoords="offset points",
fontsize=5.5,
color="#4D4D4D",
)
axis.axvline(1.0, color="#666666", linestyle="--", linewidth=0.8, label="RCDI = 1")
axis.set(xlabel="RCDI (closer to 1 indicates similarity)", ylabel="CAI", title=host)
self._clean_axes(axis)
for axis in axes.flat[len(hosts) :]:
axis.set_visible(False)
handles, labels = axes.flat[0].get_legend_handles_labels()
figure.legend(
handles,
labels,
loc="lower center",
bbox_to_anchor=(0.5, 0.095),
ncol=min(4, len(labels)),
frameon=False,
)
figure.suptitle(title, y=0.96, fontsize=11, fontweight="bold")
figure.text(
0.995,
0.025,
"Marker area encodes exact 59-codon SiD similarity; all measures are descriptive",
ha="right",
va="bottom",
fontsize=6.5,
color="#4D4D4D",
)
return figure
[docs]
def dinucleotide_heatmap(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "Mean dinucleotide observed/expected ratios",
) -> Any:
"""Plot mean O/E ratios for all 16 dinucleotides by metadata group."""
if not results:
raise ValueError("results must be non-empty.")
plt, np = self._matplotlib()
grouped: dict[str, list[AnalysisResult]] = defaultdict(list)
for result in results:
grouped[str(result.metadata.get(group_key, "Unknown"))].append(result)
group_names = sorted(grouped)
dinucleotides = [left + right for left in "ACGT" for right in "ACGT"]
matrix = np.asarray(
[
[
np.mean(
[
result.dinucleotide_odds[pair]
for result in grouped[name]
if result.dinucleotide_odds[pair] is not None
]
)
for pair in dinucleotides
]
for name in group_names
],
dtype=float,
)
figure, axis = plt.subplots(figsize=(7.2, max(2.6, 0.58 * len(group_names))), constrained_layout=True)
image = axis.imshow(matrix, cmap="BrBG", vmin=0.65, vmax=1.35, aspect="auto")
axis.set_xticks(range(len(dinucleotides)), dinucleotides, rotation=45, ha="right")
axis.set_yticks(range(len(group_names)), [f"{name} (n={len(grouped[name])})" for name in group_names])
axis.set_xlabel("Dinucleotide")
axis.set_ylabel("Genotype")
axis.set_title(title)
for row_index in range(matrix.shape[0]):
for column_index in range(matrix.shape[1]):
axis.text(
column_index,
row_index,
f"{matrix[row_index, column_index]:.2f}",
ha="center",
va="center",
fontsize=5.5,
)
colorbar = figure.colorbar(image, ax=axis, fraction=0.025, pad=0.02)
colorbar.set_label("Observed/expected ratio")
return figure
[docs]
def rscu_heatmap(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "Mean relative synonymous codon usage",
) -> Any:
"""Plot group-mean RSCU values for all available sense codons."""
if not results:
raise ValueError("results must be non-empty.")
plt, np = self._matplotlib()
grouped: dict[str, list[AnalysisResult]] = defaultdict(list)
for result in results:
grouped[str(result.metadata.get(group_key, "Unknown"))].append(result)
group_names = sorted(grouped)
codons = sorted({codon for result in results for codon in result.rscu})
matrix = np.asarray(
[
[np.mean([result.rscu.get(codon, 0.0) for result in grouped[name]]) for codon in codons]
for name in group_names
],
dtype=float,
)
figure, axis = plt.subplots(figsize=(10.0, max(2.7, 0.58 * len(group_names))), constrained_layout=True)
image = axis.imshow(matrix, cmap="viridis", vmin=0.0, vmax=max(2.0, float(matrix.max())), aspect="auto")
axis.set_xticks(range(len(codons)), codons, rotation=90)
axis.set_yticks(range(len(group_names)), [f"{name} (n={len(grouped[name])})" for name in group_names])
axis.set_xlabel("Codon")
axis.set_ylabel("Genotype")
axis.set_title(title)
colorbar = figure.colorbar(image, ax=axis, fraction=0.018, pad=0.02)
colorbar.set_label("Mean RSCU")
return figure
[docs]
def evolutionary_diagnostics(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "Evolutionary codon-usage diagnostics",
) -> Any:
"""Create ENC-GC3, neutrality, and PR2 panels by metadata group."""
if not results:
raise ValueError("results must be non-empty.")
plt, np = self._matplotlib()
from .metrics.evolutionary import EvolutionaryMetrics
figure, axes = plt.subplots(1, 3, figsize=(10.8, 3.5), constrained_layout=True)
group_names = sorted({str(result.metadata.get(group_key, "Unknown")) for result in results})
for index, group in enumerate(group_names):
selected = [result for result in results if str(result.metadata.get(group_key, "Unknown")) == group]
color = self.COLORS[index % len(self.COLORS)]
marker = self.MARKERS[index % len(self.MARKERS)]
common = {
"s": 25,
"color": color,
"marker": marker,
"alpha": 0.75,
"edgecolor": "white",
"linewidth": 0.4,
"label": f"{group} (n={len(selected)})",
}
axes[0].scatter(
[result.metrics["gc3s"] for result in selected],
[result.metrics["enc"] for result in selected],
**common,
)
axes[1].scatter(
[result.metrics["gc3"] for result in selected],
[(result.metrics["gc1"] + result.metrics["gc2"]) / 2 for result in selected],
**common,
)
selected_gc3 = np.asarray([result.metrics["gc3"] for result in selected], dtype=float)
selected_gc12 = np.asarray(
[(result.metrics["gc1"] + result.metrics["gc2"]) / 2 for result in selected],
dtype=float,
)
if len(selected) >= 2 and np.ptp(selected_gc3) > 0:
slope, intercept = np.polyfit(selected_gc3, selected_gc12, 1)
trend = np.linspace(float(selected_gc3.min()), float(selected_gc3.max()), 100)
axes[1].plot(trend, intercept + slope * trend, color=color, linewidth=1.0, linestyle="--")
axes[2].scatter(
[result.metrics["at_bias"] for result in selected],
[result.metrics["gc_bias"] for result in selected],
**common,
)
grid = np.linspace(0.01, 0.99, 200)
axes[0].plot(
grid,
[EvolutionaryMetrics.expected_enc(float(value)) for value in grid],
color="black",
linewidth=1.0,
linestyle="--",
label="Mutation-only expectation",
)
axes[2].axhline(0.5, color="black", linewidth=0.8, linestyle="--")
axes[2].axvline(0.5, color="black", linewidth=0.8, linestyle="--")
axes[0].set(xlabel="GC3s", ylabel="Effective number of codons", title="ENC-GC3s")
axes[1].set(xlabel="GC3", ylabel="GC12", title="Neutrality")
axes[2].set(xlabel="A3/(A3+T3)", ylabel="G3/(G3+C3)", title="PR2")
for index, axis in enumerate(axes):
axis.text(-0.16, 1.08, chr(65 + index), transform=axis.transAxes, fontweight="bold", fontsize=11)
self._clean_axes(axis)
handles, labels = axes[0].get_legend_handles_labels()
figure.legend(
handles,
labels,
loc="lower center",
bbox_to_anchor=(0.5, -0.03),
ncol=min(4, len(labels)),
frameon=False,
)
figure.suptitle(title, fontsize=11, fontweight="bold")
return figure
[docs]
def rscu_profiles(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "Relative synonymous codon usage profiles",
) -> Any:
"""Plot group-mean 59-codon RSCU profiles with standard errors.
Codons are arranged by amino-acid family and colored by their third
nucleotide. Horizontal reference lines identify the conventional
underrepresentation (0.6) and preference (1.6) thresholds.
"""
if not results:
raise ValueError("results must be non-empty.")
plt, np = self._matplotlib()
from .genetic_code import GeneticCode
code = GeneticCode.from_ncbi(1)
codon_rows = [
(amino_acid, codon)
for amino_acid, family in sorted(code.synonymous_codons.items())
if amino_acid != "*" and len(family) > 1
for codon in sorted(family)
]
groups = sorted({str(result.metadata.get(group_key, "Unknown")) for result in results})
figure, axes = plt.subplots(
len(groups),
1,
figsize=(11.0, 2.55 * len(groups)),
sharex=True,
constrained_layout=True,
squeeze=False,
)
base_colors = {"A": "#E69F00", "C": "#0072B2", "G": "#009E73", "T": "#CC79A7"}
positions = np.arange(len(codon_rows))
for panel, group in enumerate(groups):
axis = axes[panel, 0]
selected = [result for result in results if str(result.metadata.get(group_key, "Unknown")) == group]
means = []
errors = []
for _, codon in codon_rows:
values = np.asarray([result.rscu[codon] for result in selected], dtype=float)
means.append(float(values.mean()))
errors.append(float(values.std(ddof=1) / math.sqrt(len(values))) if len(values) > 1 else 0.0)
axis.bar(
positions,
means,
yerr=errors,
color=[base_colors[codon[2]] for _, codon in codon_rows],
edgecolor="white",
linewidth=0.2,
error_kw={"elinewidth": 0.45, "capsize": 1.0, "capthick": 0.45},
)
axis.axhline(0.6, color="#666666", linestyle=":", linewidth=0.8)
axis.axhline(1.6, color="#B2182B", linestyle="--", linewidth=0.8)
axis.set_ylabel("Mean RSCU")
axis.set_title(f"Genotype {group} (n={len(selected)})", loc="left", fontweight="bold")
self._clean_axes(axis)
family_starts: dict[str, int] = {}
for index, (amino_acid, _) in enumerate(codon_rows):
family_starts.setdefault(amino_acid, index)
boundaries = [index - 0.5 for index in family_starts.values() if index]
for boundary in boundaries:
axis.axvline(boundary, color="#D9D9D9", linewidth=0.4, zorder=0)
axes[-1, 0].set_xticks(positions, [codon for _, codon in codon_rows], rotation=90)
axes[-1, 0].set_xlabel("Codon (bar color indicates third nucleotide)")
figure.suptitle(title, fontsize=11, fontweight="bold")
return figure
[docs]
def dinucleotide_profiles(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "Dinucleotide observed/expected profiles",
) -> Any:
"""Plot all 16 group-mean dinucleotide O/E ratios with SEM bars."""
if not results:
raise ValueError("results must be non-empty.")
plt, np = self._matplotlib()
pairs = [left + right for left in "ACGT" for right in "ACGT"]
groups = sorted({str(result.metadata.get(group_key, "Unknown")) for result in results})
figure, axis = plt.subplots(figsize=(8.2, 3.9), constrained_layout=True)
positions = np.arange(len(pairs), dtype=float)
offsets = np.linspace(-0.18, 0.18, len(groups)) if len(groups) > 1 else np.asarray([0.0])
for index, (group, offset) in enumerate(zip(groups, offsets, strict=True)):
selected = [result for result in results if str(result.metadata.get(group_key, "Unknown")) == group]
means = []
errors = []
for pair in pairs:
values = np.asarray(
[
result.dinucleotide_odds[pair]
for result in selected
if result.dinucleotide_odds[pair] is not None
],
dtype=float,
)
means.append(float(values.mean()))
errors.append(float(values.std(ddof=1) / math.sqrt(len(values))) if len(values) > 1 else 0.0)
axis.errorbar(
positions + offset,
means,
yerr=errors,
color=self.COLORS[index % len(self.COLORS)],
marker=self.MARKERS[index % len(self.MARKERS)],
markersize=4,
linewidth=1.0,
capsize=2,
label=f"{group} (n={len(selected)})",
)
axis.axhline(1.25, color="#B2182B", linestyle="--", linewidth=0.9, label="Overrepresented (>1.25)")
axis.axhline(0.78, color="#666666", linestyle=":", linewidth=0.9, label="Underrepresented (<0.78)")
axis.set_xticks(positions, pairs)
axis.set(xlabel="Dinucleotide", ylabel="Observed/expected ratio", title=title)
axis.legend(frameon=False, ncol=3, loc="upper center", bbox_to_anchor=(0.5, -0.18))
self._clean_axes(axis)
return figure
[docs]
def composition_overview(
self,
results: list[AnalysisResult],
*,
group_key: str,
title: str = "Nucleotide composition and codon bias",
) -> Any:
"""Create four panels summarizing composition, positional bias, and ENC."""
if not results:
raise ValueError("results must be non-empty.")
plt, np = self._matplotlib()
groups = sorted({str(result.metadata.get(group_key, "Unknown")) for result in results})
figure, axes = plt.subplots(2, 2, figsize=(8.2, 6.5))
figure.subplots_adjust(left=0.09, right=0.985, top=0.87, bottom=0.18, hspace=0.50, wspace=0.30)
def significance_label(samples: list[list[float]]) -> str:
"""Return an overall Kruskal-Wallis significance label when possible."""
non_empty = [sample for sample in samples if sample]
if len(non_empty) < 2:
return ""
try:
from scipy.stats import kruskal
p_value = float(kruskal(*non_empty).pvalue)
except (ImportError, ValueError):
return ""
if p_value < 0.0001:
return "****"
if p_value < 0.001:
return "***"
if p_value < 0.01:
return "**"
if p_value < 0.05:
return "*"
return "ns"
def grouped_bars(axis: Any, keys: list[str], labels: list[str], ylabel: str) -> None:
positions = np.arange(len(keys), dtype=float)
width = 0.72 / len(groups)
samples_by_key: dict[str, list[list[float]]] = {key: [] for key in keys}
upper_by_key: dict[str, float] = {key: 0.0 for key in keys}
for group_index, group in enumerate(groups):
selected = [result for result in results if str(result.metadata.get(group_key, "Unknown")) == group]
means = []
errors = []
for key in keys:
values = np.asarray([float(result.metrics[key]) * 100 for result in selected], dtype=float)
means.append(float(values.mean()))
errors.append(float(values.std(ddof=1) / math.sqrt(len(values))) if len(values) > 1 else 0.0)
samples_by_key[key].append(values.tolist())
upper_by_key[key] = max(upper_by_key[key], means[-1] + errors[-1])
offset = (group_index - (len(groups) - 1) / 2) * width
axis.bar(
positions + offset,
means,
width=width,
yerr=errors,
color=self.COLORS[group_index % len(self.COLORS)],
alpha=0.82,
error_kw={"elinewidth": 0.7, "capsize": 2},
label=f"{group} (n={len(selected)})",
)
span = max(upper_by_key.values(), default=1.0)
for position, key in zip(positions, keys, strict=True):
label = significance_label(samples_by_key[key])
if label:
axis.text(position, upper_by_key[key] + span * 0.035, label, ha="center", va="bottom", fontsize=7)
axis.set_xticks(positions, labels)
axis.set_ylabel(ylabel)
self._clean_axes(axis)
grouped_bars(axes[0, 0], ["a", "t", "g", "c"], ["A", "T", "G", "C"], "Composition (%)")
axes[0, 0].set_title("Overall nucleotide composition")
grouped_bars(axes[0, 1], ["at1", "at2", "at3"], ["AT1", "AT2", "AT3"], "AT content (%)")
axes[0, 1].set_title("Codon-position AT content")
grouped_bars(axes[1, 0], ["gc", "gc3", "gc3s"], ["GC", "GC3", "GC3s"], "GC content (%)")
axes[1, 0].set_title("GC composition")
random = np.random.default_rng(1)
enc_axis = axes[1, 1]
enc_values = [
[result.metrics["enc"] for result in results if str(result.metadata.get(group_key, "Unknown")) == group]
for group in groups
]
boxes = enc_axis.boxplot(enc_values, patch_artist=True, widths=0.55, showfliers=False)
for index, patch in enumerate(boxes["boxes"]):
patch.set_facecolor(self.COLORS[index % len(self.COLORS)])
patch.set_alpha(0.34)
for index, values in enumerate(enc_values, start=1):
enc_axis.scatter(
np.full(len(values), index) + random.uniform(-0.11, 0.11, len(values)),
values,
color=self.COLORS[(index - 1) % len(self.COLORS)],
s=16,
alpha=0.68,
edgecolor="white",
linewidth=0.3,
)
enc_axis.set_xticks(
range(1, len(groups) + 1),
[f"{group}\n(n={len(values)})" for group, values in zip(groups, enc_values, strict=True)],
)
enc_axis.set(xlabel="Genotype", ylabel="Effective number of codons", title="Codon usage bias")
enc_label = significance_label([[float(value) for value in values] for values in enc_values])
if enc_label:
enc_upper = max((max(values) for values in enc_values if values), default=1.0)
enc_annotation_y = enc_upper + max(enc_upper, 1.0) * 0.035
enc_axis.text(
(len(groups) + 1) / 2,
enc_annotation_y,
enc_label,
ha="center",
va="bottom",
fontsize=7,
)
enc_axis.set_ylim(top=enc_annotation_y + max(enc_upper, 1.0) * 0.05)
self._clean_axes(enc_axis)
handles, labels = axes[0, 0].get_legend_handles_labels()
figure.legend(
handles,
labels,
loc="lower center",
bbox_to_anchor=(0.5, 0.055),
ncol=len(groups),
frameon=False,
)
figure.suptitle(title, fontsize=11, fontweight="bold")
figure.text(
0.995,
0.012,
"Overall Kruskal-Wallis: ns ≥ 0.05; * < 0.05; ** < 0.01; *** < 0.001; **** < 0.0001",
ha="right",
va="bottom",
fontsize=6,
color="#4D4D4D",
)
for index, axis in enumerate(axes.flat):
axis.text(-0.14, 1.07, chr(65 + index), transform=axis.transAxes, fontweight="bold", fontsize=11)
return figure
[docs]
def lda_evaluation(
self,
analysis: dict[str, Any],
groups: list[str],
*,
title: str = "Linear discriminant analysis",
) -> Any:
"""Plot an LDA projection beside its held-out confusion matrix."""
projection = analysis.get("projection", [])
if not projection or len(projection) != len(groups):
raise ValueError("LDA projection and groups must have equal non-zero lengths.")
plt, np = self._matplotlib()
figure, axes = plt.subplots(1, 2, figsize=(8.2, 3.8), constrained_layout=True)
classes = [str(value) for value in analysis["classes"]]
coordinates = np.asarray(projection, dtype=float)
if coordinates.shape[1] == 1:
coordinates = np.column_stack([coordinates[:, 0], np.zeros(len(coordinates))])
for index, group in enumerate(classes):
selected = np.asarray([position for position, value in enumerate(groups) if str(value) == group])
axes[0].scatter(
coordinates[selected, 0],
coordinates[selected, 1],
color=self.COLORS[index % len(self.COLORS)],
marker=self.MARKERS[index % len(self.MARKERS)],
s=28,
alpha=0.78,
edgecolor="white",
linewidth=0.4,
label=f"{group} (n={len(selected)})",
)
axes[0].axhline(0, color="#BDBDBD", linewidth=0.6)
axes[0].axvline(0, color="#BDBDBD", linewidth=0.6)
axes[0].set(xlabel="LD1", ylabel="LD2", title="All-sample projection")
axes[0].legend(frameon=False, title="Genotype")
self._clean_axes(axes[0])
matrix = np.asarray(analysis["confusion_matrix"], dtype=float)
image = axes[1].imshow(matrix, cmap="Blues", vmin=0)
axes[1].set_xticks(range(len(classes)), classes)
axes[1].set_yticks(range(len(classes)), classes)
axes[1].set(
xlabel="Predicted genotype",
ylabel="True genotype",
title=f"Held-out accuracy: {analysis['accuracy'] * 100:.1f}%",
)
for row in range(matrix.shape[0]):
for column in range(matrix.shape[1]):
axes[1].text(column, row, f"{int(matrix[row, column])}", ha="center", va="center", fontsize=9)
figure.colorbar(image, ax=axes[1], fraction=0.045, pad=0.03, label="Test samples")
figure.suptitle(title, fontsize=11, fontweight="bold")
for index, axis in enumerate(axes):
axis.text(-0.16, 1.08, chr(65 + index), transform=axis.transAxes, fontweight="bold", fontsize=11)
return figure
[docs]
def host_adaptation(
self,
rows: list[dict[str, Any]],
*,
metrics: list[str],
group_key: str,
host_key: str,
title: str = "Host adaptation",
) -> Any:
"""Plot host metrics by genotype with observations, means, and SEM."""
if not rows or not metrics:
raise ValueError("rows and metrics must be non-empty.")
plt, np = self._matplotlib()
groups = sorted({str(row[group_key]) for row in rows})
hosts = list(dict.fromkeys(str(row[host_key]) for row in rows))
figure, axes = plt.subplots(1, len(metrics), figsize=(4.4 * len(metrics), 3.8), constrained_layout=True)
axes = np.atleast_1d(axes)
random = np.random.default_rng(1)
host_positions = np.arange(len(hosts), dtype=float)
offsets = np.linspace(-0.22, 0.22, len(groups)) if len(groups) > 1 else np.asarray([0.0])
for panel, metric in enumerate(metrics):
axis = axes[panel]
for index, (group, offset) in enumerate(zip(groups, offsets, strict=True)):
means = []
errors = []
for host_position, host in enumerate(hosts):
values = np.asarray(
[
float(row[metric])
for row in rows
if str(row[group_key]) == group
and str(row[host_key]) == host
and row.get(metric) is not None
],
dtype=float,
)
means.append(float(values.mean()))
errors.append(float(values.std(ddof=1) / math.sqrt(len(values))) if len(values) > 1 else 0.0)
axis.scatter(
np.full(len(values), host_position + offset) + random.uniform(-0.035, 0.035, len(values)),
values,
color=self.COLORS[index % len(self.COLORS)],
s=10,
alpha=0.25,
linewidth=0,
)
axis.errorbar(
host_positions + offset,
means,
yerr=errors,
color=self.COLORS[index % len(self.COLORS)],
marker=self.MARKERS[index % len(self.MARKERS)],
linewidth=1.0,
markersize=5,
capsize=2,
label=f"{group}",
)
axis.set_xticks(host_positions, hosts, rotation=25, ha="right")
axis.set(xlabel="Host", ylabel=self.METRIC_LABELS.get(metric.lower(), metric), title=metric)
self._clean_axes(axis)
axis.text(-0.15, 1.08, chr(65 + panel), transform=axis.transAxes, fontweight="bold", fontsize=11)
axes[0].legend(frameon=False, title="Genotype")
figure.suptitle(title, fontsize=11, fontweight="bold")
return figure
[docs]
def bland_altman_validation(
self,
panels: list[dict[str, Any]],
*,
title: str = "Method Validation: Bland-Altman and Residual Agreement",
) -> Any:
"""Create multi-panel Bland-Altman and residual agreement plots.
Each dictionary in ``panels`` must supply:
- ``metric``: str metric label (e.g. 'CAI', 'ENC', 'RSCU', 'GC3s')
- ``reference``: list[float] reference or benchmark values
- ``tool``: list[float] CodonAdaptPy calculated values
"""
if not panels:
raise ValueError("panels must be non-empty.")
plt, np = self._matplotlib()
columns = min(2, len(panels))
rows = math.ceil(len(panels) / columns)
figure, axes = plt.subplots(
rows,
columns,
figsize=(7.2, 3.2 * rows),
constrained_layout=True,
squeeze=False,
)
for index, panel in enumerate(panels):
axis = axes.flat[index]
ref = np.asarray(panel["reference"], dtype=float)
tool = np.asarray(panel["tool"], dtype=float)
if len(ref) != len(tool) or len(ref) < 3:
raise ValueError("reference and tool must contain at least three matching observations.")
means = (ref + tool) / 2.0
diffs = tool - ref
bias = float(diffs.mean())
sd = float(diffs.std(ddof=1)) if len(diffs) > 1 else 0.0
upper_loa = bias + 1.96 * sd
lower_loa = bias - 1.96 * sd
correlation = np.corrcoef(ref, tool)[0, 1] if len(ref) > 1 else 1.0
r_squared = float(correlation**2)
axis.scatter(means, diffs, s=16, color="#0072B2", alpha=0.7, edgecolor="white", linewidth=0.3)
axis.axhline(bias, color="#D55E00", linewidth=1.2, linestyle="-", label=f"Mean bias: {bias:+.4f}")
axis.axhline(upper_loa, color="#009E73", linewidth=1.0, linestyle="--", label=f"+1.96 SD: {upper_loa:+.4f}")
axis.axhline(lower_loa, color="#009E73", linewidth=1.0, linestyle="--", label=f"-1.96 SD: {lower_loa:+.4f}")
metric_name = panel.get("metric", f"Metric {index + 1}")
axis.set_xlabel(f"Mean ({metric_name})")
axis.set_ylabel("Difference (CodonAdaptPy - Ref)")
axis.set_title(f"{metric_name} (R² = {r_squared:.4f}, n = {len(ref)})")
axis.legend(frameon=False, fontsize=6.5, loc="upper right")
self._clean_axes(axis)
for axis in axes.flat[len(panels) :]:
axis.set_visible(False)
figure.suptitle(title, fontsize=11, fontweight="bold")
for index, axis in enumerate(axes.flat[: len(panels)]):
axis.text(-0.15, 1.08, chr(65 + index), transform=axis.transAxes, fontweight="bold", fontsize=11)
return figure
[docs]
def repeated_cv_distribution(
self,
lda_output: dict[str, Any],
*,
title: str = "Repeated Grouped Cross-Validation Distributions",
) -> Any:
"""Plot fold-level accuracy distributions across repeated CV iterations."""
fold_results = lda_output.get("fold_results", [])
if not fold_results:
raise ValueError("lda_output must contain non-empty 'fold_results'.")
plt, np = self._matplotlib()
figure, axis = plt.subplots(figsize=(6.5, 3.6), constrained_layout=True)
metrics = [("accuracy", "Accuracy"), ("balanced_accuracy", "Balanced Accuracy"), ("macro_f1", "Macro F1")]
positions = np.arange(len(metrics))
data = [[float(row[key]) * 100 for row in fold_results] for key, _ in metrics]
boxes = axis.boxplot(data, positions=positions, patch_artist=True, widths=0.45, showfliers=False)
for idx, patch in enumerate(boxes["boxes"]):
patch.set_facecolor(self.COLORS[idx % len(self.COLORS)])
patch.set_alpha(0.4)
patch.set_edgecolor(self.COLORS[idx % len(self.COLORS)])
random = np.random.default_rng(1)
for idx, values in enumerate(data):
jitter = random.uniform(-0.1, 0.1, len(values))
axis.scatter(
positions[idx] + jitter,
values,
s=14,
color=self.COLORS[idx % len(self.COLORS)],
alpha=0.6,
edgecolor="white",
linewidth=0.3,
)
summary = lda_output.get("summary", {})
for idx, (key, _metric_label) in enumerate(metrics):
if key in summary:
mean_val = summary[key]["mean"] * 100
sd_val = summary[key]["sd"] * 100
axis.text(
positions[idx],
min([min(d) for d in data]) - 3.0,
f"Mean: {mean_val:.1f}%\nSD: ±{sd_val:.1f}%",
ha="center",
va="top",
fontsize=7,
)
axis.set_xticks(positions, [label for _, label in metrics])
axis.set_ylabel("Held-Out Performance (%)")
folds = lda_output.get("folds", 5)
repeats = lda_output.get("repeats", 20)
subtitle = f"({folds}-fold Stratified Grouped LDA, {repeats} repeats, n={len(fold_results)} folds)"
axis.set_title(f"{title}\n{subtitle}")
self._clean_axes(axis)
return figure
[docs]
def save(self, figure: Any, path: str | Path, *, dpi: int = 600) -> Path:
"""Save a publication figure as PNG, TIFF, PDF, or SVG."""
destination = Path(path)
if destination.suffix.lower() not in {".png", ".tif", ".tiff", ".pdf", ".svg"}:
raise ValueError("Publication figures must use PNG, TIFF, PDF, or SVG.")
destination.parent.mkdir(parents=True, exist_ok=True)
figure.savefig(destination, dpi=dpi, bbox_inches="tight", facecolor="white")
return destination