GENOMICS-CDDBRG/Genome_Visualization_Tool
0
1"""2genomics_figures.py3===================4Every figure in the suite. Kept separate from the UI so figures can be5regression-tested without a browser, and separate from the maths so a plotting6change can never alter a reported number.7 8House rules for all panels:9 * 600 dpi raster / true-vector PDF-EPS with Type-42 fonts (journal default).10 * Colourblind-safe palette available and used by default.11 * Every axis carries units. Every heuristic panel says so in the panel itself,12 not only in a caption that may be dropped in review.13 * Empty input produces a labelled blank panel, never a misleading empty axis.14"""15 16from __future__ import annotations17 18from typing import Dict, List, Optional, Sequence, Tuple19 20import numpy as np21import pandas as pd22 23import matplotlib24matplotlib.use("Agg")25import matplotlib.pyplot as plt26from matplotlib.patches import Patch27 28import genomics_core as gc29 30__all__ = [31 "fig_composition_overview", "fig_skew_and_ori", "fig_rscu_heatmap",32 "fig_enc_plot", "fig_neutrality_plot", "fig_pr2_plot", "fig_codon_bar",33 "fig_amino_acid_usage", "fig_dinucleotide_oe", "fig_assembly_curve",34 "fig_distance_heatmap", "fig_pangenome_curves", "fig_pangenome_pie",35 "fig_presence_absence", "fig_feature_map", "fig_length_hist",36]37 38_ANNOT_KW = dict(fontsize=7.5, color="#8a4b08", style="italic")39 40 41def _finish(fig, tight=True):42 if tight:43 fig.tight_layout()44 return fig45 46 47def _blank(title: str, msg: str, figsize=(7, 4)):48 fig, ax = plt.subplots(figsize=figsize)49 gc.blank_panel(ax, msg)50 ax.set_title(title)51 return _finish(fig)52 53 54# ----------------------------------------------------------------------55# Composition56# ----------------------------------------------------------------------57 58def fig_composition_overview(records: pd.DataFrame, theme: str = "Ocean"):59 """4-panel: base composition, GC per record, length distribution, entropy."""60 if records is None or records.empty:61 return _blank("Composition overview", "No records to display")62 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])63 fig, axes = plt.subplots(2, 2, figsize=(11, 7.5))64 65 ax = axes[0, 0]66 bases = ["A", "T", "G", "C", "N"]67 vals = [records[f"count_{b}"].sum() for b in bases]68 total = sum(vals) or 169 cols = ["#4C72B0", "#DD8452", "#55A868", "#C44E52", "#999999"]70 bars = ax.bar(bases, [100.0 * v / total for v in vals], color=cols,71 edgecolor="black", linewidth=0.5)72 for b, v in zip(bars, vals):73 ax.text(b.get_x() + b.get_width() / 2, b.get_height(),74 f"{100.0*v/total:.2f}%", ha="center", va="bottom", fontsize=8)75 ax.set_ylabel("Composition (% of all bases)")76 ax.set_title("Base composition")77 78 ax = axes[0, 1]79 ax.hist(records["gc_percent"], bins=min(30, max(5, len(records))),80 color=th["primary"], edgecolor="black", linewidth=0.5, alpha=0.85)81 mean_gc = float(records["gc_percent"].mean())82 ax.axvline(mean_gc, color=th["secondary"], ls="--", lw=1.5,83 label=f"mean = {mean_gc:.2f}%")84 ax.set_xlabel("GC content (%)")85 ax.set_ylabel("Number of records")86 ax.set_title("GC distribution across records")87 ax.legend()88 89 ax = axes[1, 0]90 L = np.sort(records["length_bp"].to_numpy())[::-1]91 ax.plot(np.arange(1, len(L) + 1), L, marker="o" if len(L) < 40 else None,92 color=th["primary"], lw=1.6, ms=3)93 ax.set_yscale("log" if L.max() / max(L.min(), 1) > 100 else "linear")94 ax.set_xlabel("Record rank (longest first)")95 ax.set_ylabel("Length (bp)")96 ax.set_title("Record length profile")97 98 ax = axes[1, 1]99 if "entropy_bits" in records:100 ax.scatter(records["gc_percent"], records["entropy_bits"], s=22,101 color=th["primary"], edgecolor="black", linewidth=0.3, alpha=0.85)102 ax.axhline(2.0, color="#888888", ls=":", lw=1.2, label="maximum (2 bits)")103 ax.set_xlabel("GC content (%)")104 ax.set_ylabel("Shannon entropy (bits)")105 ax.set_title("Sequence complexity")106 ax.legend()107 else:108 gc.blank_panel(ax, "Entropy not computed")109 return _finish(fig)110 111 112def fig_skew_and_ori(window_df: pd.DataFrame, ori_ter: Dict, record_name: str = "",113 theme: str = "Ocean"):114 """GC/AT skew and cumulative skew with predicted ori/ter marked."""115 if window_df is None or window_df.empty:116 return _blank("Strand skew", "No windows to display")117 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])118 fig, axes = plt.subplots(2, 1, figsize=(11, 6.4), sharex=True)119 x = window_df["midpoint"] / 1000.0120 121 ax = axes[0]122 ax.plot(x, window_df["GC_skew"], lw=0.9, color=th["primary"], label="GC skew")123 ax.plot(x, window_df["AT_skew"], lw=0.9, color=th["secondary"], alpha=0.75, label="AT skew")124 ax.axhline(0, color="#666666", lw=0.8)125 ax.set_ylabel("Skew per window")126 ax.set_title(f"Strand asymmetry{(' - ' + record_name) if record_name else ''}")127 ax.legend(ncol=2, loc="upper right")128 129 ax = axes[1]130 ax.plot(x, window_df["cumulative_GC_skew"], lw=1.8, color=th["primary"],131 label="cumulative GC skew")132 conf = str(ori_ter.get("ori_ter_confidence", ""))133 trustworthy = conf.startswith("clear")134 for key, colour, label in [("predicted_origin_bp", "#1a9850", "predicted origin"),135 ("predicted_terminus_bp", "#d73027", "predicted terminus")]:136 v = ori_ter.get(key)137 if v is not None:138 ax.axvline(v / 1000.0, color=colour, ls="--" if trustworthy else ":",139 lw=1.6, label=f"{label} ~{v:,} bp")140 ax.set_xlabel("Genome position (kb)")141 ax.set_ylabel("Cumulative GC skew")142 ax.legend(loc="best")143 # The honesty annotation lives inside the figure, so it survives being144 # lifted into a slide deck without its caption.145 ax.annotate(146 "Ori/ter are PREDICTED from cumulative GC skew (Grigoriev 1998), not measured."147 + ("" if trustworthy else " Signal is weak here - treat as unreliable."),148 xy=(0.005, -0.30), xycoords="axes fraction", **_ANNOT_KW)149 return _finish(fig)150 151 152def fig_dinucleotide_oe(df: pd.DataFrame, theme: str = "Ocean"):153 if df is None or df.empty:154 return _blank("Dinucleotide relative abundance", "No data")155 fig, ax = plt.subplots(figsize=(9, 4.2))156 colours = ["#c0392b" if m else "#4C72B0" for m in df["marked_deviation"]]157 ax.bar(df["dinucleotide"], df["rho"], color=colours, edgecolor="black", linewidth=0.5)158 ax.axhline(1.0, color="#333333", lw=1.0)159 ax.axhspan(0.78, 1.23, color="#999999", alpha=0.14,160 label="Karlin |rho-1| <= 0.23 (unremarkable)")161 ax.set_ylabel(r"$\rho$ = observed / expected")162 ax.set_xlabel("Dinucleotide")163 ax.set_title("Dinucleotide relative abundance (genomic signature)")164 ax.legend(loc="upper right")165 ax.annotate("Red bars mark |rho - 1| > 0.23 (Karlin & Burge 1995).",166 xy=(0.005, -0.22), xycoords="axes fraction", **_ANNOT_KW)167 return _finish(fig)168 169 170# ----------------------------------------------------------------------171# Codon usage172# ----------------------------------------------------------------------173 174def fig_rscu_heatmap(usage: pd.DataFrame):175 if usage is None or usage.empty:176 return _blank("RSCU", "No codon usage")177 d = usage[(usage["amino_acid_1"] != "*") & usage["RSCU"].notna()].copy()178 if d.empty:179 return _blank("RSCU", "No sense codons")180 aas = sorted(d["amino_acid_1"].unique())181 maxfam = int(d.groupby("amino_acid_1").size().max())182 M = np.full((len(aas), maxfam), np.nan)183 labels = np.full((len(aas), maxfam), "", dtype=object)184 for i, aa in enumerate(aas):185 sub = d[d["amino_acid_1"] == aa].sort_values("codon").reset_index(drop=True)186 for j, r in sub.iterrows():187 M[i, j] = r["RSCU"]188 labels[i, j] = f"{r['codon']}\n{r['RSCU']:.2f}"189 fig, ax = plt.subplots(figsize=(1.15 * maxfam + 3.2, 0.46 * len(aas) + 2.2))190 im = ax.imshow(M, cmap=gc.RSCU_CMAP, aspect="auto", vmin=0,191 vmax=float(np.nanmax(M)) if np.isfinite(M).any() else 1)192 for i in range(len(aas)):193 for j in range(maxfam):194 if labels[i, j]:195 v = M[i, j]196 tc = "white" if v > np.nanmax(M) * 0.6 else "#1a1a1a"197 ax.text(j, i, labels[i, j], ha="center", va="center", fontsize=6.8, color=tc)198 ax.set_yticks(range(len(aas)))199 ax.set_yticklabels([f"{a} ({gc.AA_THREE_LETTER.get(a,a)})" for a in aas])200 ax.set_xticks(range(maxfam))201 ax.set_xticklabels([f"codon {i+1}" for i in range(maxfam)])202 ax.set_title("Relative synonymous codon usage (RSCU)")203 ax.grid(False)204 cb = fig.colorbar(im, ax=ax, fraction=0.025, pad=0.02)205 cb.set_label("RSCU (1.0 = no bias)")206 return _finish(fig)207 208 209def fig_codon_bar(usage: pd.DataFrame, theme: str = "Ocean"):210 if usage is None or usage.empty:211 return _blank("Codon usage", "No data")212 d = usage[usage["amino_acid_1"] != "*"].sort_values(["amino_acid_1", "codon"])213 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])214 fig, ax = plt.subplots(figsize=(13, 4.4))215 aas = list(dict.fromkeys(d["amino_acid_1"]))216 cmap = plt.get_cmap("tab20")217 colour = {aa: cmap(i % 20) for i, aa in enumerate(aas)}218 ax.bar(range(len(d)), d["frequency_per_1000"],219 color=[colour[a] for a in d["amino_acid_1"]], edgecolor="black", linewidth=0.3)220 ax.set_xticks(range(len(d)))221 ax.set_xticklabels(d["codon"], rotation=90, fontsize=5.6, family="monospace")222 ax.set_ylabel("Frequency (per 1000 codons)")223 ax.set_title("Genome-wide codon usage")224 ax.margins(x=0.005)225 return _finish(fig)226 227 228def fig_enc_plot(per_gene: pd.DataFrame, theme: str = "Ocean"):229 """Wright's ENC vs GC3s with the mutation-drift expectation curve."""230 if per_gene is None or per_gene.empty:231 return _blank("ENC plot", "No genes passed the length filter")232 d = per_gene.dropna(subset=["ENC", "GC3s"])233 if d.empty:234 return _blank("ENC plot", "No genes with computable ENC")235 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])236 fig, ax = plt.subplots(figsize=(7.2, 5.6))237 x = d["GC3s"].to_numpy() / 100.0238 ax.scatter(d["GC3s"], d["ENC"], s=16, alpha=0.6, color=th["primary"],239 edgecolor="none", label=f"genes (n = {len(d):,})")240 xs = np.linspace(0.005, 0.995, 400)241 ax.plot(xs * 100, gc.wright_expected_enc(xs), color="#c0392b", lw=2.0,242 label="Wright expectation (mutation drift only)")243 exp = gc.wright_expected_enc(x)244 below = int(np.sum(d["ENC"].to_numpy() < exp - 0.10 * exp))245 ax.set_xlabel("GC3s (%) - GC at synonymous third positions")246 ax.set_ylabel("ENC (effective number of codons)")247 ax.set_xlim(0, 100)248 ax.set_ylim(18, 64)249 ax.set_title("ENC plot (Wright 1990)")250 ax.legend(loc="lower center")251 ax.annotate(252 f"{below:,}/{len(d):,} genes lie >10% below the curve, i.e. show codon bias "253 f"beyond what GC mutation pressure alone explains.",254 xy=(0.005, -0.155), xycoords="axes fraction", **_ANNOT_KW)255 return _finish(fig)256 257 258def fig_neutrality_plot(per_gene: pd.DataFrame, reg: Dict, theme: str = "Ocean"):259 if per_gene is None or per_gene.empty:260 return _blank("Neutrality plot", "No genes")261 d = per_gene.dropna(subset=["GC3", "GC12"])262 if len(d) < 3:263 return _blank("Neutrality plot", "Need at least 3 genes")264 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])265 fig, ax = plt.subplots(figsize=(7.2, 5.6))266 ax.scatter(d["GC3"], d["GC12"], s=16, alpha=0.6, color=th["primary"], edgecolor="none")267 xs = np.linspace(d["GC3"].min(), d["GC3"].max(), 100)268 if np.isfinite(reg.get("slope", np.nan)):269 ax.plot(xs, reg["slope"] * xs + reg["intercept"], color="#c0392b", lw=2.0,270 label=(f"GC12 = {reg['slope']:.3f}·GC3 + {reg['intercept']:.2f}\n"271 f"r² = {reg['r_squared']:.3f}, n = {reg['n_genes']:,}"))272 ax.plot(xs, xs, color="#666666", ls="--", lw=1.2, label="slope = 1 (pure mutation pressure)")273 ax.set_xlabel("GC3 (%)")274 ax.set_ylabel("GC12 (%)")275 ax.set_title("Neutrality plot (Sueoka 1988)")276 ax.legend(loc="best")277 s = reg.get("slope", np.nan)278 if np.isfinite(s):279 interp = ("mutation pressure dominates" if s > 0.5 else280 "selection dominates codon usage" if s < 0.2 else281 "mixed mutation and selection")282 ax.annotate(f"Slope {s:.3f}: {interp}. Slope approximates the neutral "283 f"(mutation-driven) fraction of GC variation.",284 xy=(0.005, -0.155), xycoords="axes fraction", **_ANNOT_KW)285 return _finish(fig)286 287 288def fig_pr2_plot(per_gene: pd.DataFrame, theme: str = "Ocean"):289 if per_gene is None or per_gene.empty:290 return _blank("PR2 plot", "No genes")291 d = per_gene.dropna(subset=["PR2_G3_over_GC3", "PR2_A3_over_AT3"])292 if d.empty:293 return _blank("PR2 plot", "No four-fold degenerate codons found")294 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])295 fig, ax = plt.subplots(figsize=(6.4, 6.2))296 ax.scatter(d["PR2_G3_over_GC3"], d["PR2_A3_over_AT3"], s=16, alpha=0.6,297 color=th["primary"], edgecolor="none", label=f"genes (n = {len(d):,})")298 ax.axhline(0.5, color="#666666", lw=1.0)299 ax.axvline(0.5, color="#666666", lw=1.0)300 ax.plot(0.5, 0.5, marker="*", ms=17, color="#c0392b", ls="none",301 label="PR2 centre (no bias)")302 mx, my = float(d["PR2_G3_over_GC3"].mean()), float(d["PR2_A3_over_AT3"].mean())303 ax.plot(mx, my, marker="D", ms=8, color="#1a9850", ls="none",304 label=f"observed mean ({mx:.3f}, {my:.3f})")305 ax.set_xlim(0, 1); ax.set_ylim(0, 1); ax.set_aspect("equal")306 ax.set_xlabel("G3 / (G3 + C3)")307 ax.set_ylabel("A3 / (A3 + T3)")308 ax.set_title("PR2 bias plot, four-fold degenerate sites (Sueoka 1995)")309 ax.legend(loc="upper right", fontsize=7.5)310 ax.annotate("Displacement from the centre indicates strand-asymmetric mutation "311 "or selection at synonymous sites.",312 xy=(0.005, -0.13), xycoords="axes fraction", **_ANNOT_KW)313 return _finish(fig)314 315 316def fig_amino_acid_usage(df: pd.DataFrame, theme: str = "Ocean"):317 if df is None or df.empty:318 return _blank("Amino acid usage", "No data")319 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])320 d = df.sort_values("frequency_percent", ascending=False)321 fig, ax = plt.subplots(figsize=(9, 4.2))322 ax.bar(d["amino_acid_3"], d["frequency_percent"], color=th["primary"],323 edgecolor="black", linewidth=0.4)324 ax.set_ylabel("Frequency (% of residues)")325 ax.set_xlabel("Amino acid")326 ax.set_title("Genome-wide amino acid usage")327 plt.setp(ax.get_xticklabels(), rotation=45, ha="right")328 return _finish(fig)329 330 331# ----------------------------------------------------------------------332# Assembly / features333# ----------------------------------------------------------------------334 335def fig_assembly_curve(lengths: Sequence[int], stats: Dict, theme: str = "Ocean"):336 """Cumulative assembly (Nx) curve with N50/N75/N90 marked."""337 if not len(lengths):338 return _blank("Assembly curve", "No records")339 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])340 L = np.sort(np.asarray(lengths, dtype=float))[::-1]341 cum = np.cumsum(L)342 pct = 100.0 * cum / cum[-1]343 fig, ax = plt.subplots(figsize=(7.6, 5.0))344 ax.step(pct, L, where="post", lw=2.0, color=th["primary"])345 for key, colour in [("N50_bp", "#1a9850"), ("N75_bp", "#E69F00"), ("N90_bp", "#d73027")]:346 v = stats.get(key)347 if v:348 ax.axhline(v, color=colour, ls="--", lw=1.3, label=f"{key.replace('_bp','')} = {v:,} bp")349 ax.set_yscale("log" if L.max() / max(L.min(), 1) > 100 else "linear")350 ax.set_xlabel("Cumulative assembly length (%)")351 ax.set_ylabel("Record length (bp)")352 ax.set_title("Assembly contiguity (Nx curve)")353 ax.legend(loc="upper right")354 if stats.get("auN"):355 ax.annotate(f"auN (area under Nx) = {stats['auN']:,.0f} bp - a single-value "356 f"contiguity summary robust to the N50 threshold artefact.",357 xy=(0.005, -0.145), xycoords="axes fraction", **_ANNOT_KW)358 return _finish(fig)359 360 361def fig_length_hist(df: pd.DataFrame, col: str, title: str, xlabel: str, theme: str = "Ocean"):362 if df is None or df.empty or col not in df:363 return _blank(title, "No data")364 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])365 v = pd.to_numeric(df[col], errors="coerce").dropna()366 if v.empty:367 return _blank(title, "No numeric data")368 fig, ax = plt.subplots(figsize=(7.6, 4.4))369 ax.hist(v, bins=min(60, max(6, int(np.sqrt(len(v))))), color=th["primary"],370 edgecolor="black", linewidth=0.4, alpha=0.88)371 ax.axvline(v.median(), color=th["secondary"], ls="--", lw=1.5,372 label=f"median = {v.median():,.0f}")373 ax.set_xlabel(xlabel); ax.set_ylabel("Count"); ax.set_title(title); ax.legend()374 return _finish(fig)375 376 377def fig_feature_map(features: pd.DataFrame, seq_len: int, record_name: str = "",378 palette: str = "Colourblind safe (Okabe-Ito)", max_features: int = 4000):379 """Linear feature map, strand-separated."""380 if features is None or features.empty:381 return _blank("Feature map", "No features to display", figsize=(12, 3))382 pal = gc.PALETTES.get(palette, gc.PALETTES["Colourblind safe (Okabe-Ito)"])383 d = features.copy()384 if len(d) > max_features:385 d = d.reindex(d["end"].astype(int).sub(d["start"].astype(int)).abs()386 .sort_values(ascending=False).index[:max_features])387 fig, ax = plt.subplots(figsize=(13, 3.4))388 ax.axhline(0, color="#999999", lw=1.0, zorder=1)389 seen = {}390 for _, r in d.iterrows():391 t = str(r.get("type", "feature"))392 col = pal.get(t, pal["default"])393 seen[t] = col394 s, e = int(r["start"]), int(r["end"])395 strand = int(r.get("strand", 0) or 0)396 y = 0.16 if strand >= 0 else -0.16 - 0.30397 ax.add_patch(plt.Rectangle((s, y), max(e - s, max(1, seq_len // 4000)), 0.30,398 facecolor=col, edgecolor="none", alpha=0.9, zorder=2))399 ax.set_xlim(0, seq_len)400 ax.set_ylim(-0.75, 0.65)401 ax.set_yticks([0.31, -0.31])402 ax.set_yticklabels(["+ strand", "− strand"])403 ax.set_xlabel("Position (bp)")404 ax.set_title(f"Feature map{(' - ' + record_name) if record_name else ''}"405 f" ({len(d):,} features shown)")406 ax.grid(axis="x")407 if seen:408 ax.legend(handles=[Patch(facecolor=c, label=t) for t, c in sorted(seen.items())],409 loc="upper center", bbox_to_anchor=(0.5, -0.30), ncol=min(8, len(seen)),410 fontsize=7)411 return _finish(fig)412 413 414# ----------------------------------------------------------------------415# Comparative / pan-genome416# ----------------------------------------------------------------------417 418def fig_distance_heatmap(mat: pd.DataFrame, title: str, cbar_label: str,419 cmap: str = "viridis", fmt: str = "{:.3f}"):420 if mat is None or mat.empty:421 return _blank(title, "No matrix")422 n = len(mat)423 fig, ax = plt.subplots(figsize=(max(5.0, 0.68 * n + 3.2), max(4.2, 0.58 * n + 2.6)))424 im = ax.imshow(mat.to_numpy(dtype=float), cmap=cmap, aspect="auto")425 ax.set_xticks(range(n)); ax.set_yticks(range(n))426 ax.set_xticklabels(mat.columns, rotation=45, ha="right", fontsize=7.5)427 ax.set_yticklabels(mat.index, fontsize=7.5)428 if n <= 14:429 A = mat.to_numpy(dtype=float)430 lo, hi = np.nanmin(A), np.nanmax(A)431 for i in range(n):432 for j in range(n):433 v = A[i, j]434 rel = (v - lo) / (hi - lo) if hi > lo else 0.5435 ax.text(j, i, fmt.format(v), ha="center", va="center", fontsize=6.6,436 color="white" if rel > 0.55 else "#1a1a1a")437 ax.set_title(title)438 ax.grid(False)439 cb = fig.colorbar(im, ax=ax, fraction=0.035, pad=0.02)440 cb.set_label(cbar_label)441 return _finish(fig)442 443 444def fig_pangenome_curves(acc: pd.DataFrame, heaps: Dict, theme: str = "Ocean"):445 """Pan/core accumulation curves with the Heaps' law fit overlaid."""446 if acc is None or acc.empty:447 return _blank("Pan-genome accumulation", "No data")448 th = gc.THEMES.get(theme, gc.THEMES["Ocean"])449 fig, ax = plt.subplots(figsize=(7.8, 5.2))450 x = acc["n_genomes"]451 ax.errorbar(x, acc["pan_mean"], yerr=acc["pan_sd"], marker="o", ms=4.5, lw=1.8,452 capsize=3, color=th["primary"], label="pan-genome")453 ax.errorbar(x, acc["core_mean"], yerr=acc["core_sd"], marker="s", ms=4.5, lw=1.8,454 capsize=3, color=th["secondary"], label="core genome")455 g = heaps.get("gamma", np.nan)456 if np.isfinite(g) and np.isfinite(heaps.get("kappa", np.nan)):457 xs = np.linspace(1, x.max(), 200)458 ax.plot(xs, heaps["kappa"] * xs ** g, ls="--", lw=1.6, color="#c0392b",459 label=f"Heaps fit: n = {heaps['kappa']:.0f}·N^{g:.3f}")460 ax.set_xlabel("Number of genomes sampled")461 ax.set_ylabel("Number of gene clusters")462 ax.set_title("Pan-genome and core-genome accumulation")463 ax.legend(loc="center right")464 note = heaps.get("interpretation", "")465 if note:466 ax.annotate(f"Heaps' law (Tettelin 2008): {note}",467 xy=(0.005, -0.145), xycoords="axes fraction", **_ANNOT_KW)468 return _finish(fig)469 470 471def fig_pangenome_pie(counts: Dict[str, int], theme: str = "Ocean"):472 labels = [k for k, v in counts.items() if v > 0]473 vals = [counts[k] for k in labels]474 if not vals:475 return _blank("Pan-genome partition", "No clusters")476 fig, ax = plt.subplots(figsize=(6.6, 5.2))477 cols = {"Core": "#2166ac", "Soft-core": "#67a9cf", "Shell": "#fdae61", "Cloud": "#d73027"}478 total = sum(vals)479 w, *_ = ax.pie(vals, labels=None, colors=[cols.get(l, "#999999") for l in labels],480 autopct=lambda p: f"{p:.1f}%\n({int(round(p*total/100)):,})",481 textprops=dict(fontsize=8, color="white"), startangle=90,482 wedgeprops=dict(edgecolor="white", linewidth=1.6))483 ax.legend(w, [f"{l} (n = {counts[l]:,})" for l in labels],484 loc="center left", bbox_to_anchor=(1.0, 0.5))485 ax.set_title(f"Pan-genome partition (total {total:,} clusters)")486 return _finish(fig)487 488 489def fig_presence_absence(presence: np.ndarray, genome_names: List[str],490 max_clusters: int = 3000):491 if presence is None or presence.size == 0:492 return _blank("Presence/absence", "No matrix")493 M = presence494 if M.shape[0] > max_clusters:495 # Order by prevalence, then evenly subsample so the core/shell/cloud496 # structure is preserved rather than truncated to the core block.497 order = np.argsort(-M.sum(axis=1))498 M = M[order][np.linspace(0, len(order) - 1, max_clusters).astype(int)]499 else:500 M = M[np.argsort(-M.sum(axis=1))]501 fig, ax = plt.subplots(figsize=(max(5.5, 0.42 * len(genome_names) + 3.0), 6.2))502 ax.imshow(M, cmap="Blues", aspect="auto", interpolation="nearest", vmin=0, vmax=1)503 ax.set_xticks(range(len(genome_names)))504 ax.set_xticklabels(genome_names, rotation=45, ha="right", fontsize=7.5)505 ax.set_ylabel(f"Gene clusters (n = {presence.shape[0]:,}, sorted by prevalence)")506 ax.set_title("Gene presence / absence")507 ax.grid(False)508 if presence.shape[0] > max_clusters:509 ax.annotate(f"Showing {max_clusters:,} of {presence.shape[0]:,} clusters "510 f"(evenly sampled across prevalence). Full matrix is in the CSV export.",511 xy=(0.005, -0.12), xycoords="axes fraction", **_ANNOT_KW)512 return _finish(fig)513