Custom Code in Automation Designer -- Add Examples

examples/custom-code-AD/README.md

@@ -1,5 +1,29 @@
 # Custom Code in Automation Designer
 
+## 5-parameter-logistic
+
+Example snippet fits a 5-Parameter Logistic (5PL) curve to dose-response data to extract EC50, Hill slope, and asymmetry, with AIC comparison against a standard 4PL fit and an interactive plot with residuals subplot.
+
+- **[Input File](./snippets/5-parameter-logistic/docs/input/dose_response_5pl_data.xlsx)**
+- **[Custom Code Block](./snippets/5-parameter-logistic/dose_response_5pl.py)**
+- **[Output Files](./snippets/5-parameter-logistic/docs/output)**
+
+## kinetic-timecourse
+
+Example snippet fits one-phase association kinetic curves per compound to extract parameters (Emax, kobs, t½, AUC, R²), with t-test comparison at the inflection point and a multi-compound time-course plot with 95% CI ribbons.
+
+- **[Input File](./snippets/kinetic-timecourse/docs/input/timecourse_data.xlsx)**
+- **[Custom Code Block](./snippets/kinetic-timecourse/kinetics_timecourse.py)**
+- **[Output Files](./snippets/kinetic-timecourse/docs/output)**
+
+## multi-group-anova
+
+Example snippet performs one-way ANOVA to test for significant differences across treatment groups, followed by Tukey HSD post-hoc analysis for pairwise comparisons, with descriptive statistics and an annotated bar chart output.
+
+- **[Input File](./snippets/multi-group-anova/docs/input/multigroup_efficacy_data.xlsx)**
+- **[Custom Code Block](./snippets/multi-group-anova/multigroup_anova.py)**
+- **[Output Files](./snippets/multi-group-anova/docs/output)**
+
 ## plot-chromatogram
 
 Example snippet demonstrates plotting a multi-axis HPLC chromatogram (retention volume X Absorbance (mAU), Temperature (degC), pH) with custom code.

examples/custom-code-AD/requirements.txt

@@ -1,12 +1,15 @@
 allotropy==0.1.105
 biopython==1.86
+flowkit==0.1.0
 lmfit==1.3.4
+matplotlib==3.10.3
 numpy==2.2.4
 openpyxl==3.1.5
 pandas==2.2.3
 plotly==5.22.0
 pyarrow==19.0.1
-pydantic==2.12.5
+pydantic==1.10.21
+seaborn==0.13.2
 scikit-learn==1.6.1
 scipy==1.15.2
 statsmodels==0.14.4

examples/custom-code-AD/snippets/5-parameter-logistic/docs/output/5pl_parameters.csv

@@ -0,0 +1,9 @@
+Parameter,Value,Std_Error,Note
+Bottom (%),1.37,0.873,
+Top (%),99.235,0.9498,
+EC50 (µM),0.7516,0.2311,
+Hill Slope,1.623,0.2467,
+Asymmetry (F),0.775,0.2657,asym=1.0 → standard 4PL
+R²,0.9994,,
+AIC (5PL),11.69,,
+ΔAIC vs 4PL,-1.31,,Positive = 5PL preferred

examples/custom-code-AD/snippets/5-parameter-logistic/docs/output/aggregated_data.csv

@@ -0,0 +1,12 @@
+Concentration_uM,Mean_Inhibition,SD,N,SEM
+0.001,2.4375,2.027,4,1.0135
+0.003,0.7925,0.9151,4,0.4576
+0.01,2.34,1.6329,4,0.8164
+0.03,2.2925,1.7193,4,0.8596
+0.1,8.165,2.0677,4,1.0338
+0.3,28.765,2.5523,4,1.2761
+1.0,67.31,2.2448,4,1.1224
+3.0,93.6675,2.1078,4,1.0539
+10.0,96.17,2.5479,4,1.2739
+30.0,99.835,3.1382,4,1.5691
+100.0,99.25,1.5948,4,0.7974

examples/custom-code-AD/snippets/5-parameter-logistic/dose_response_5pl.py

@@ -0,0 +1,261 @@
+"""
+Benchling Custom Code Demo 5: Dose-Response 5PL Curve Fitting
+=============================================================
+5-Parameter Logistic (5PL) model adds an asymmetry parameter to the
+standard 4PL, giving a better fit for asymmetric sigmoidal curves —
+common in immunoassays and ELISA-based dose-response experiments.
+
+INPUTS:
+    inputs[0]: pd.DataFrame with columns:
+        - Concentration_uM      (float)
+        - Replicate             (int)
+        - Signal_%Inhibition    (float)
+
+OUTPUTS:
+    - "5pl_parameters"      pd.DataFrame  — EC50, Hill slope, Asymmetry, R², AIC, ± SE
+    - "aggregated_data"     pd.DataFrame  — Mean ± SEM per concentration
+    - "dose_response_curve" go.Figure     — Interactive 5PL plot with residuals subplot
+
+Supported packages:
+allotropy, biopython, lmfit, numpy, openpyxl, pandas, plotly,
+pyarrow, pydantic, scikit-learn, scipy, statsmodels
+"""
+
+from io import BytesIO
+import numpy as np
+import pandas as pd
+import plotly.graph_objects as go
+from plotly.subplots import make_subplots
+from typing import NamedTuple
+from lmfit import Model
+
+
+class IOData(NamedTuple):
+    name: str
+    data: BytesIO | pd.DataFrame | go.Figure
+
+
+# ---------------------------------------------------------------------------
+# 5PL model
+# ---------------------------------------------------------------------------
+def five_pl(x, bottom, top, ec50, hill, asym):
+    """
+    5-Parameter Logistic model.
+    Extends 4PL with 'asym' (asymmetry / F parameter) which allows the
+    upper and lower plateaus to be approached at different rates.
+    asym=1.0 reduces exactly to the standard 4PL model.
+    """
+    return bottom + (top - bottom) / ((1.0 + (ec50 / np.clip(x, 1e-12, None)) ** hill) ** asym)
+
+
+def fit_5pl(concentrations: np.ndarray, responses: np.ndarray):
+    """
+    Fit 5PL using lmfit. Returns (result, r_squared, aic).
+    Also fits 4PL (asym fixed=1) so AIC can flag if 5PL is warranted.
+    """
+    model = Model(five_pl)
+
+    # --- 5PL fit ---
+    params_5pl = model.make_params(
+        bottom=dict(value=float(responses.min()), min=0,    max=20),
+        top=   dict(value=float(responses.max()), min=80,   max=110),
+        ec50=  dict(value=float(np.median(concentrations)), min=1e-6, max=1e6),
+        hill=  dict(value=1.5,  min=0.1,  max=10),
+        asym=  dict(value=1.0,  min=0.05, max=10),
+    )
+    result_5pl = model.fit(responses, params_5pl, x=concentrations)
+
+    # --- 4PL fit (asym locked to 1) for AIC comparison ---
+    params_4pl = model.make_params(
+        bottom=dict(value=float(responses.min()), min=0,   max=20),
+        top=   dict(value=float(responses.max()), min=80,  max=110),
+        ec50=  dict(value=float(np.median(concentrations)), min=1e-6, max=1e6),
+        hill=  dict(value=1.5, min=0.1, max=10),
+        asym=  dict(value=1.0, vary=False),
+    )
+    result_4pl = model.fit(responses, params_4pl, x=concentrations)
+
+    # R²
+    fitted  = result_5pl.best_fit
+    ss_res  = float(np.sum((responses - fitted) ** 2))
+    ss_tot  = float(np.sum((responses - responses.mean()) ** 2))
+    r2      = 1.0 - ss_res / ss_tot if ss_tot > 0 else 0.0
+
+    return result_5pl, result_4pl, r2
+
+
+# ---------------------------------------------------------------------------
+# Entry point
+# ---------------------------------------------------------------------------
+def custom_code(inputs: list[IOData], **kwargs) -> list[IOData]:
+    # --- Load input ---
+    df = None
+    for i in inputs:
+        if isinstance(i.data, pd.DataFrame):
+            df = i.data
+            break
+        elif isinstance(i.data, BytesIO):
+            df = pd.read_excel(i.data)
+            break
+    if df is None:
+        raise ValueError("No DataFrame or Excel file found in inputs")
+
+    df["Concentration_uM"]   = df["Concentration_uM"].astype(float)
+    df["Signal_%Inhibition"] = df["Signal_%Inhibition"].astype(float)
+
+    # --- Aggregate replicates ---
+    agg = (
+        df.groupby("Concentration_uM")["Signal_%Inhibition"]
+        .agg(Mean="mean", SD="std", N="count")
+        .reset_index()
+    )
+    agg.columns = ["Concentration_uM", "Mean_Inhibition", "SD", "N"]
+    agg["SEM"] = agg["SD"] / np.sqrt(agg["N"])
+
+    conc = agg["Concentration_uM"].values
+    resp = agg["Mean_Inhibition"].values
+
+    # --- Fit 5PL (and 4PL for comparison) ---
+    result_5pl, result_4pl, r2 = fit_5pl(conc, resp)
+    p = result_5pl.params
+
+    def _se(param):
+        return round(float(param.stderr), 4) if param.stderr else None
+
+    # Delta AIC: positive = 5PL is better fit
+    delta_aic = result_4pl.aic - result_5pl.aic
+
+    summary = pd.DataFrame({
+        "Parameter": ["Bottom (%)", "Top (%)", "EC50 (µM)", "Hill Slope", "Asymmetry (F)", "R²", "AIC (5PL)", "ΔAIC vs 4PL"],
+        "Value": [
+            round(float(p["bottom"].value), 3),
+            round(float(p["top"].value),    3),
+            round(float(p["ec50"].value),   4),
+            round(float(p["hill"].value),   3),
+            round(float(p["asym"].value),   3),
+            round(r2, 4),
+            round(float(result_5pl.aic), 2),
+            round(float(delta_aic), 2),
+        ],
+        "Std_Error": [
+            _se(p["bottom"]), _se(p["top"]), _se(p["ec50"]),
+            _se(p["hill"]),   _se(p["asym"]), None, None, None,
+        ],
+        "Note": [
+            "", "", "",
+            "", "asym=1.0 → standard 4PL",
+            "", "",
+            "Positive = 5PL preferred",
+        ],
+    })
+
+    # --- Smooth fit lines ---
+    x_fit   = np.logspace(np.log10(conc.min()), np.log10(conc.max()), 400)
+    y_5pl   = five_pl(x_fit, p["bottom"].value, p["top"].value,
+                      p["ec50"].value, p["hill"].value, p["asym"].value)
+    p4      = result_4pl.params
+    y_4pl   = five_pl(x_fit, p4["bottom"].value, p4["top"].value,
+                      p4["ec50"].value, p4["hill"].value, 1.0)
+
+    # Residuals (5PL)
+    fitted_at_conc = five_pl(conc, p["bottom"].value, p["top"].value,
+                             p["ec50"].value, p["hill"].value, p["asym"].value)
+    residuals = resp - fitted_at_conc
+    ec50_val  = float(p["ec50"].value)
+
+    # --- Build figure: main curve + residuals subplot ---
+    fig = make_subplots(
+        rows=2, cols=1,
+        row_heights=[0.72, 0.28],
+        shared_xaxes=False,
+        vertical_spacing=0.12,
+        subplot_titles=("5PL Dose-Response Fit", "Residuals"),
+    )
+
+    # Individual replicates
+    fig.add_trace(go.Scatter(
+        x=df["Concentration_uM"], y=df["Signal_%Inhibition"],
+        mode="markers",
+        marker=dict(symbol="circle-open", size=6, color="#7FBBDA", opacity=0.7),
+        name="Replicates",
+        showlegend=True,
+    ), row=1, col=1)
+
+    # Mean ± SEM
+    fig.add_trace(go.Scatter(
+        x=agg["Concentration_uM"], y=agg["Mean_Inhibition"],
+        error_y=dict(array=agg["SEM"].tolist(), visible=True, color="#1B6CA8"),
+        mode="markers",
+        marker=dict(size=10, color="#1B6CA8"),
+        name="Mean ± SEM",
+    ), row=1, col=1)
+
+    # 4PL fit (reference)
+    fig.add_trace(go.Scatter(
+        x=x_fit, y=y_4pl,
+        mode="lines",
+        line=dict(color="#AAAAAA", width=1.8, dash="dash"),
+        name=f"4PL fit (ΔAIC={delta_aic:+.1f})",
+    ), row=1, col=1)
+
+    # 5PL fit
+    fig.add_trace(go.Scatter(
+        x=x_fit, y=y_5pl,
+        mode="lines",
+        line=dict(color="#E84545", width=2.5),
+        name=f"5PL fit  R²={r2:.4f}",
+    ), row=1, col=1)
+
+    # EC50 dashed vertical line using paper-relative x coords on row 1 axes
+    fig.add_shape(
+        type="line",
+        xref="x", yref="y",
+        x0=ec50_val, x1=ec50_val,
+        y0=-5, y1=110,
+        line=dict(color="#888888", width=1.2, dash="dot"),
+        row=1, col=1,
+    )
+    fig.add_annotation(
+        x=np.log10(ec50_val), y=55,
+        xref="x", yref="y",
+        text=f"EC50={ec50_val:.3f} µM",
+        showarrow=False,
+        font=dict(size=11, color="#555555"),
+        bgcolor="rgba(255,255,255,0.7)",
+        xanchor="left",
+        row=1, col=1,
+    )
+
+    # Residuals
+    fig.add_trace(go.Scatter(
+        x=conc, y=residuals,
+        mode="markers+lines",
+        marker=dict(size=8, color="#E84545"),
+        line=dict(color="#E84545", width=1, dash="dot"),
+        name="Residuals",
+        showlegend=False,
+    ), row=2, col=1)
+    fig.add_trace(go.Scatter(
+        x=[conc.min(), conc.max()], y=[0, 0],
+        mode="lines", line=dict(color="#AAAAAA", width=1),
+        showlegend=False, hoverinfo="skip",
+    ), row=2, col=1)
+
+    # Set log scale explicitly on both subplots
+    fig.update_xaxes(type="log", row=1, col=1)
+    fig.update_xaxes(title_text="Concentration (µM)", type="log", row=2, col=1)
+    fig.update_yaxes(title_text="% Inhibition", range=[-5, 110], row=1, col=1)
+    fig.update_yaxes(title_text="Residual", zeroline=True, row=2, col=1)
+
+    fig.update_layout(
+        template="plotly_white",
+        legend=dict(x=1.02, y=0.95, xanchor="left"),
+        width=860, height=640,
+        margin=dict(t=60, r=160),
+    )
+
+    return [
+        IOData(name="5pl_parameters",     data=summary),
+        IOData(name="aggregated_data",     data=agg.round(4)),
+        IOData(name="dose_response_curve", data=fig),
+    ]

examples/custom-code-AD/snippets/kinetic-timecourse/docs/output/aggregated_kinetics.csv

@@ -0,0 +1,31 @@
+Compound,Time_hr,Mean,SD,N,SEM,CI95_lower,CI95_upper
+Compound_A,0.0,-0.2787,1.3743,3,0.7935,-1.8338,1.2765
+Compound_A,0.5,1.277,0.7716,3,0.4455,0.4039,2.1501
+Compound_A,1.0,1.8333,1.0417,3,0.6014,0.6546,3.0121
+Compound_A,2.0,3.54,0.6152,3,0.3552,2.8438,4.2362
+Compound_A,4.0,6.349,0.2737,3,0.158,6.0393,6.6587
+Compound_A,8.0,11.643,0.5953,3,0.3437,10.9693,12.3167
+Compound_A,12.0,13.6943,1.7458,3,1.008,11.7188,15.6699
+Compound_A,24.0,18.3767,0.2974,3,0.1717,18.0401,18.7132
+Compound_A,48.0,19.4823,0.1698,3,0.098,19.2902,19.6745
+Compound_A,72.0,20.15,0.585,3,0.3377,19.4881,20.8119
+Compound_B,0.0,0.1373,0.6439,3,0.3718,-0.5913,0.866
+Compound_B,0.5,0.6577,0.7196,3,0.4155,-0.1567,1.472
+Compound_B,1.0,1.188,0.575,3,0.332,0.5373,1.8387
+Compound_B,2.0,2.2697,0.5702,3,0.3292,1.6244,2.9149
+Compound_B,4.0,3.5597,0.463,3,0.2673,3.0358,4.0836
+Compound_B,8.0,7.0487,0.5562,3,0.3211,6.4193,7.6781
+Compound_B,12.0,9.5947,0.5182,3,0.2992,9.0083,10.1811
+Compound_B,24.0,14.1617,0.402,3,0.2321,13.7067,14.6166
+Compound_B,48.0,18.3357,0.2967,3,0.1713,17.9999,18.6715
+Compound_B,72.0,19.525,1.3894,3,0.8022,17.9527,21.0973
+Vehicle,0.0,0.196,1.0212,3,0.5896,-0.9596,1.3516
+Vehicle,0.5,0.3307,0.295,3,0.1703,-0.0032,0.6645
+Vehicle,1.0,-0.3997,0.4085,3,0.2359,-0.862,0.0626
+Vehicle,2.0,0.0007,0.7273,3,0.4199,-0.8223,0.8237
+Vehicle,4.0,-0.126,0.4066,3,0.2348,-0.5862,0.3342
+Vehicle,8.0,-0.089,0.4928,3,0.2845,-0.6466,0.4686
+Vehicle,12.0,-0.199,0.3337,3,0.1927,-0.5767,0.1787
+Vehicle,24.0,-0.953,0.3277,3,0.1892,-1.3239,-0.5821
+Vehicle,48.0,0.4193,0.6914,3,0.3992,-0.3631,1.2017
+Vehicle,72.0,0.181,0.441,3,0.2546,-0.318,0.68

examples/custom-code-AD/snippets/kinetic-timecourse/docs/output/kinetic_params.csv

@@ -0,0 +1,4 @@
+Compound,Emax_AU,kobs_1_hr,t_half_hr,AUC_AU_hr,R2
+Compound_A,19.949,0.10149,6.83,1222.58,0.9985
+Compound_B,19.867,0.05363,12.92,1049.56,0.9994
+Vehicle,7.749,0.0001,6931.36,-7.33,-0.0248

examples/custom-code-AD/snippets/kinetic-timecourse/docs/output/ttest_summary.csv

@@ -0,0 +1,2 @@
+Comparison,Timepoint_hr,Timepoint_basis,t_statistic,p_value,Significance
+Compound_A vs Compound_B,12.0,1/kobs inflection (target=14.2h),3.8992,0.01755,*