#!/usr/bin/env python3 """Generate Sepia Fig. 2 extract comparison from raw endpoint RFU data. The endpoint GFP signal is calibrated by applying the inverse of the 12 h fluorescein 4-parameter logistic curve recorded in the packaged parameter CSV. This is not a slope-fit or linear standard-curve conversion. """ from __future__ import annotations import json from pathlib import Path import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from sepia_figure_style import LINEAR_GREEN, PLASMID_GREEN, TEXT, style_axis SCRIPT_VERSION = "2026-07-10.1" FIG_W_PT = 540.0 FIG_H_PT = 540.0 FIG_W_IN = FIG_W_PT / 72.0 FIG_H_IN = FIG_H_PT / 72.0 RAW_LONG = "260505_fig1_all_timepoints_raw_long.csv" FIT_PARAMETERS = "260505_fig1_12h_fluorescein_4pl_params.csv" PROVIDED_CALIBRATED = "260505_fig1_12h_cellfree_calibrated_replicates.csv" KIT_ORDER = ["N", "G1", "G2", "Sepia"] KIT_LABELS = { "N": "N", "G1": "G1", "G2": "G2", "Sepia": "Sepia", } DNA_COLORS = { "Plasmid": PLASMID_GREEN, "Linear": LINEAR_GREEN, } COMPETITOR_GRAYS = { "N": "#f4f4f4", "G1": "#d9d9d9", "G2": "#bdbdbd", } def find_input_dir() -> Path: here = Path(__file__).resolve().parent candidates = [here / "data", here, Path.cwd() / "data", Path.cwd()] for candidate in candidates: if (candidate / RAW_LONG).exists() and (candidate / FIT_PARAMETERS).exists(): return candidate raise FileNotFoundError( f"Could not find {RAW_LONG} and {FIT_PARAMETERS} in ./data or next to this script." ) def read_4pl_params(parameter_path: Path) -> dict[str, float]: fit = pd.read_csv(parameter_path) params: dict[str, float] = {} for _, row in fit.iloc[:8].iterrows(): key = row.get("Parameter") value = row.get("Value") if isinstance(key, str) and key in {"Bottom", "Top", "EC50_nM", "Hill", "R_squared"}: params[key] = float(value) missing = {"Bottom", "Top", "EC50_nM", "Hill"} - params.keys() if missing: raise ValueError(f"Workbook 4PL fit is missing parameters: {sorted(missing)}") return params def inverse_4pl_rfu_to_nM(rfu: pd.Series | np.ndarray, params: dict[str, float]) -> np.ndarray: bottom = params["Bottom"] top = params["Top"] ec50 = params["EC50_nM"] hill = params["Hill"] y = np.asarray(rfu, dtype=float) clipped = np.clip(y, bottom + 1e-9, top - 1e-9) return ec50 / (((top - bottom) / (clipped - bottom) - 1.0) ** (1.0 / hill)) def build_calibrated_replicates(input_dir: Path) -> tuple[pd.DataFrame, dict[str, float], float | None]: params = read_4pl_params(input_dir / FIT_PARAMETERS) raw = pd.read_csv(input_dir / RAW_LONG) reps = raw.loc[raw["Time_h"].eq(12) & raw["Assay"].eq("Cell-free")].copy() required = {"Condition", "DNA_Type", "Well", "Replicate", "RFU_485_528_G35"} missing = required - set(reps.columns) if missing: raise ValueError(f"Raw long table missing columns: {sorted(missing)}") reps["Fluorescein_Equivalent_nM"] = inverse_4pl_rfu_to_nM(reps["RFU_485_528_G35"], params) reps["Fluorescein_Equivalent_uM"] = reps["Fluorescein_Equivalent_nM"] / 1000.0 reps = reps[["Condition", "DNA_Type", "Well", "Replicate", "RFU_485_528_G35", "Fluorescein_Equivalent_nM", "Fluorescein_Equivalent_uM"]] max_delta = None provided_path = input_dir / PROVIDED_CALIBRATED if provided_path.exists(): provided = pd.read_csv(provided_path) merged = reps.merge(provided, on=["Condition", "DNA_Type", "Well", "Replicate"], suffixes=("_computed", "_provided")) merged["Delta_uM"] = ( merged["Fluorescein_Equivalent_uM_computed"] - merged["Fluorescein_Equivalent_uM_provided"] ) max_delta = float(merged["Delta_uM"].abs().max()) return reps, params, max_delta def summarize_replicates(reps: pd.DataFrame) -> pd.DataFrame: summary = ( reps.groupby(["Condition", "DNA_Type"], as_index=False) .agg( N=("Fluorescein_Equivalent_uM", "count"), Mean_RFE_uM=("Fluorescein_Equivalent_uM", "mean"), SD_RFE_uM=("Fluorescein_Equivalent_uM", "std"), Mean_RFU=("RFU_485_528_G35", "mean"), SD_RFU=("RFU_485_528_G35", "std"), ) .sort_values(["DNA_Type", "Condition"]) ) summary["SD_RFE_uM"] = summary["SD_RFE_uM"].fillna(0.0) return summary def plot_panel(ax: plt.Axes, reps: pd.DataFrame, summary: pd.DataFrame, dna_type: str, *, show_ylabel: bool) -> None: style_axis(ax) sub_summary = summary.loc[summary["DNA_Type"].eq(dna_type)].set_index("Condition").loc[KIT_ORDER].reset_index() x = np.arange(len(KIT_ORDER)) colors = [DNA_COLORS[dna_type] if kit == "Sepia" else COMPETITOR_GRAYS[kit] for kit in KIT_ORDER] ax.bar( x, sub_summary["Mean_RFE_uM"], yerr=sub_summary["SD_RFE_uM"], width=0.64, color=colors, edgecolor=TEXT, linewidth=1.55, error_kw={"elinewidth": 1.8, "capthick": 1.8, "capsize": 4.5, "ecolor": TEXT}, zorder=3, ) ax.set_xlim(-0.55, len(KIT_ORDER) - 0.45) ax.set_ylim(0, 30) ax.set_yticks([0, 15, 30]) ax.set_xticks(x) ax.set_xticklabels([KIT_LABELS[kit] for kit in KIT_ORDER]) ax.set_xlabel("Kit", labelpad=8) if show_ylabel: ax.set_ylabel("Relative Fluorescein Equivalents (µM)", labelpad=8) else: ax.set_ylabel("") ax.set_yticklabels([]) ax.tick_params(axis="y", length=0) ax.spines["left"].set_visible(False) label = "Plasmid DNA" if dna_type == "Plasmid" else "Linear DNA" ax.text( 0.50, 0.94, label, color=TEXT, fontsize=17, fontweight="bold", ha="center", va="top", transform=ax.transAxes, ) def make_figure(input_dir: Path, out_svg: Path, out_png: Path | None = None) -> None: reps, params, max_delta = build_calibrated_replicates(input_dir) summary = summarize_replicates(reps) out_svg.parent.mkdir(parents=True, exist_ok=True) reps.to_csv(out_svg.parent / "sepia_fig2_calibrated_replicates_used.csv", index=False) summary.to_csv(out_svg.parent / "sepia_fig2_endpoint_summary_used.csv", index=False) fig = plt.figure(figsize=(FIG_W_IN, FIG_H_IN), facecolor="none") ax_plasmid = fig.add_axes([0.135, 0.170, 0.375, 0.720]) ax_linear = fig.add_axes([0.590, 0.170, 0.345, 0.720]) plot_panel(ax_plasmid, reps, summary, "Plasmid", show_ylabel=True) plot_panel(ax_linear, reps, summary, "Linear", show_ylabel=False) for ax in (ax_plasmid, ax_linear): for tick_label in ax.get_xticklabels() + ax.get_yticklabels(): tick_label.set_fontweight("bold") ax.xaxis.label.set_fontweight("bold") ax.yaxis.label.set_fontweight("bold") fig.savefig(out_svg, format="svg", transparent=True, bbox_inches=None, metadata={"Date": None}) if out_png is not None: fig.savefig(out_png, dpi=220, transparent=True, bbox_inches=None) manifest = { "figure": out_svg.name, "script": Path(__file__).name, "script_version": SCRIPT_VERSION, "width_pt": FIG_W_PT, "height_pt": FIG_H_PT, "site_slot": "OpenCFPS product template Fig. 2 asset; template image width=540 height=540", "raw_endpoint_source": RAW_LONG, "fit_parameter_source": FIT_PARAMETERS, "calibration_method": "4-parameter logistic fluorescein standard curve; endpoint GFP RFU converted by inverse 4PL.", "not_slope_fit": True, "four_pl_equation": "RFU = Bottom + (Top-Bottom)/(1 + (EC50/nM)^Hill)", "four_pl_params": params, "provided_calibrated_validation_max_abs_delta_uM": max_delta, "endpoint_time_h": 12, "endpoint_error_bars": "sample SD across replicate wells", "x_axis": "Kit, split into plasmid DNA and linear DNA side-by-side panels", "y_axis": "Relative Fluorescein Equivalents (µM)", "kit_label_map": KIT_LABELS, "visual_encoding": "Competitor kits are grayscale; Sepia bars use the DNA-type color, with plasmid dark green and linear light green; panel labels are black.", "theme_notes": "Matches Fig. 1 typography/axes with transparent background, black label text, grayscale competitor bars, Sepia DNA-type green bars, and SD error bars.", } (out_svg.parent / "sepia_fig2_reproducibility_manifest.json").write_text(json.dumps(manifest, indent=2)) plt.close(fig) def main() -> None: input_dir = find_input_dir() out_dir = Path(__file__).resolve().parent / "outputs" make_figure( input_dir=input_dir, out_svg=out_dir / "sepia_fig2_extract_endpoint_comparison.svg", out_png=out_dir / "sepia_fig2_extract_endpoint_comparison.png", ) print(f"Wrote {out_dir / 'sepia_fig2_extract_endpoint_comparison.svg'}") if __name__ == "__main__": main()