#!/usr/bin/env python3 # ============================================================================= # FrESH - Spatial bias analysis (hotspot detection) # Final paper-facing version # ============================================================================= # Reads: # data/freezing_temps_pooled.csv # columns: chiller, plate_position, experiment_id, well_id, freezing_temp_BT # # Outputs: # output/spatial_bias_heatmaps_4panel.png # output/spatial_bias_heatmaps_4panel.pdf # output/index_bias_all_chillers.csv # output/morans_i_results.csv # output/spatial_bias_summary.csv # output/spatial_ci_fields.csv # ============================================================================= from pathlib import Path import matplotlib as mpl import matplotlib.gridspec as gridspec import matplotlib.pyplot as plt import numpy as np import pandas as pd import seaborn as sns # --------------------------------------------------------------------------- # Configuration # --------------------------------------------------------------------------- DATA_FILE = Path("data/freezing_temps_pooled.csv") OUTPUT_DIR = Path("output") OUTPUT_DIR.mkdir(parents=True, exist_ok=True) CHILLERS = ["RE1050", "RP1845"] PLATES = ["A", "B"] DISPLAY_VMAX = 1.0 DISPLAY_TICKS = np.array([-1.0, -0.5, 0.0, 0.5, 1.0]) # --------------------------------------------------------------------------- # Bias accumulation # --------------------------------------------------------------------------- def accumulate_index_bias(df_plate: pd.DataFrame, bias_vector: np.ndarray) -> None: n = len(df_plate) if n < 2: return df_s = df_plate.sort_values("freezing_temp_BT", ascending=False).reset_index(drop=True) ranks = np.arange(n) / (n - 1) for pos, row in df_s.iterrows(): idx = int(row["well_id"]) if 0 <= idx <= 95: bias_vector[idx] += ranks[pos] - 0.5 def accumulate_temp_bias( df_plate: pd.DataFrame, temp_sum: np.ndarray, temp_count: np.ndarray, ) -> None: plate_median = df_plate["freezing_temp_BT"].median() for _, row in df_plate.iterrows(): idx = int(row["well_id"]) val = row["freezing_temp_BT"] if 0 <= idx <= 95 and np.isfinite(val): temp_sum[idx] += val - plate_median temp_count[idx] += 1 # --------------------------------------------------------------------------- # Spatial statistics # --------------------------------------------------------------------------- def morans_i_permutation_test( bias_matrix: np.ndarray, n_perm: int = 9999, seed: int = 42, ) -> tuple[float, float, float]: rng = np.random.default_rng(seed) valid = np.isfinite(bias_matrix) if valid.sum() < 2: return np.nan, np.nan, np.nan coords = np.argwhere(valid) vals = bias_matrix[valid].astype(float) centered = vals - np.mean(vals) n = len(centered) ss = np.sum(centered ** 2) if ss == 0: return 0.0, 1.0, 1.0 pairs = [] w_sum = 0 for i in range(n): ri, ci = coords[i] for j in range(n): if i == j: continue rj, cj = coords[j] if abs(ri - rj) <= 1 and abs(ci - cj) <= 1: pairs.append((i, j)) w_sum += 1 if w_sum == 0: return np.nan, np.nan, np.nan def compute_I(c: np.ndarray) -> float: num = sum(c[i] * c[j] for i, j in pairs) return (n / w_sum) * (num / np.sum(c ** 2)) I_obs = compute_I(centered) permuted_I = np.array([compute_I(rng.permutation(centered)) for _ in range(n_perm)]) p_upper = (np.sum(permuted_I >= I_obs) + 1) / (n_perm + 1) p_lower = (np.sum(permuted_I <= I_obs) + 1) / (n_perm + 1) p_perm = min(2 * min(p_upper, p_lower), 1.0) return float(I_obs), float(p_perm), float(p_upper) def indices_to_well_matrix(values: np.ndarray, fill_value: float = np.nan) -> np.ndarray: matrix = np.full((12, 8), fill_value, dtype=float) for idx in range(96): csv_row = idx // 12 csv_col = idx % 12 matrix[csv_col, 7 - csv_row] = values[idx] return matrix def spatial_stats_from_bias( bias_vector: np.ndarray, count_vector: np.ndarray, ) -> dict: mask = np.isfinite(bias_vector) & np.isfinite(count_vector) & (count_vector > 0) if not np.any(mask): return { "n_wells": 0, "weighted_mean": np.nan, "weighted_sigma": np.nan, "half_range": np.nan, "cmin": np.nan, "cmax": np.nan, } vals = bias_vector[mask].astype(float) weights = count_vector[mask].astype(float) weighted_mean = float(np.sum(weights * vals) / np.sum(weights)) weighted_var = float(np.sum(weights * (vals - weighted_mean) ** 2) / np.sum(weights)) weighted_sigma = float(np.sqrt(weighted_var)) cmin = float(np.min(vals)) cmax = float(np.max(vals)) half_range = 0.5 * (cmax - cmin) return { "n_wells": int(np.sum(mask)), "weighted_mean": weighted_mean, "weighted_sigma": weighted_sigma, "half_range": half_range, "cmin": cmin, "cmax": cmax, } # --------------------------------------------------------------------------- # Processing # --------------------------------------------------------------------------- def process_pooled(df: pd.DataFrame) -> dict: all_bias_A = {} all_bias_B = {} all_bias_C_A = {} all_bias_C_B = {} all_count_C_A = {} all_count_C_B = {} all_stats = {} morans = {} spatial_summary = {} for chiller in CHILLERS: df_ch = df[df["chiller"] == chiller].copy() bias_A = np.zeros(96) bias_B = np.zeros(96) temp_sum_A = np.zeros(96) temp_sum_B = np.zeros(96) temp_cnt_A = np.zeros(96) temp_cnt_B = np.zeros(96) n_exp_A = 0 n_exp_B = 0 for plate_pos, bias_vec, t_sum, t_cnt in [ ("A", bias_A, temp_sum_A, temp_cnt_A), ("B", bias_B, temp_sum_B, temp_cnt_B), ]: df_plate = ( df_ch[df_ch["plate_position"] == plate_pos][ ["experiment_id", "well_id", "freezing_temp_BT"] ] .dropna() .copy() ) df_plate["well_id"] = df_plate["well_id"].astype(int) n_exp = 0 for _, grp in df_plate.groupby("experiment_id", sort=True): if len(grp) < 2: continue accumulate_index_bias(grp, bias_vec) accumulate_temp_bias(grp, t_sum, t_cnt) n_exp += 1 if plate_pos == "A": n_exp_A = n_exp else: n_exp_B = n_exp print( f" {chiller} Plate {plate_pos}: {n_exp} experiments, " f"{len(df_plate)} total freezing events" ) all_stats[chiller] = {"n_exp_A": n_exp_A, "n_exp_B": n_exp_B} all_bias_A[chiller] = bias_A / max(n_exp_A, 1) all_bias_B[chiller] = bias_B / max(n_exp_B, 1) with np.errstate(invalid="ignore", divide="ignore"): all_bias_C_A[chiller] = np.where(temp_cnt_A > 0, temp_sum_A / temp_cnt_A, np.nan) all_bias_C_B[chiller] = np.where(temp_cnt_B > 0, temp_sum_B / temp_cnt_B, np.nan) all_count_C_A[chiller] = temp_cnt_A.copy() all_count_C_B[chiller] = temp_cnt_B.copy() mA = indices_to_well_matrix(all_bias_C_A[chiller]) mB = indices_to_well_matrix(all_bias_C_B[chiller]) I_A, p_A, pA_upper = morans_i_permutation_test(mA) I_B, p_B, pB_upper = morans_i_permutation_test(mB) morans[chiller] = { "I_A": I_A, "p_A": p_A, "pA_upper": pA_upper, "I_B": I_B, "p_B": p_B, "pB_upper": pB_upper, } print(f" {chiller} Plate A: Moran's I = {I_A:+.4f}, p = {p_A:.4g}") print(f" {chiller} Plate B: Moran's I = {I_B:+.4f}, p = {p_B:.4g}") stats_A = spatial_stats_from_bias(all_bias_C_A[chiller], temp_cnt_A) stats_B = spatial_stats_from_bias(all_bias_C_B[chiller], temp_cnt_B) spatial_summary[chiller] = {"A": stats_A, "B": stats_B} print( f" {chiller} Plate A: c_i range {stats_A['cmin']:+.3f} to {stats_A['cmax']:+.3f} °C, " f"half-range = {stats_A['half_range']:.3f} °C, " f"weighted mean = {stats_A['weighted_mean']:+.3f} °C, " f"weighted sigma = {stats_A['weighted_sigma']:.3f} °C" ) print( f" {chiller} Plate B: c_i range {stats_B['cmin']:+.3f} to {stats_B['cmax']:+.3f} °C, " f"half-range = {stats_B['half_range']:.3f} °C, " f"weighted mean = {stats_B['weighted_mean']:+.3f} °C, " f"weighted sigma = {stats_B['weighted_sigma']:.3f} °C" ) return { "all_bias_A": all_bias_A, "all_bias_B": all_bias_B, "all_bias_C_A": all_bias_C_A, "all_bias_C_B": all_bias_C_B, "all_count_C_A": all_count_C_A, "all_count_C_B": all_count_C_B, "all_stats": all_stats, "morans": morans, "spatial_summary": spatial_summary, } # --------------------------------------------------------------------------- # Figure # --------------------------------------------------------------------------- def plot_spatial_heatmaps_4panel( all_bias_C_A: dict, all_bias_C_B: dict, morans_results: dict, ) -> None: assert len(CHILLERS) == 2 plt.rcParams.update( { "font.size": 9, "axes.titlesize": 9, "axes.labelsize": 9, "xtick.labelsize": 8, "ytick.labelsize": 8, } ) row_labels = ["H", "G", "F", "E", "D", "C", "B", "A"] col_labels = [str(i) for i in range(1, 13)] matrices = {} for chiller in CHILLERS: mA = indices_to_well_matrix(all_bias_C_A[chiller]) mB = indices_to_well_matrix(all_bias_C_B[chiller]) matrices[chiller] = (mA, mB) cmap = mpl.colormaps["PuOr_r"].copy() cmap.set_bad(color="white") cell_in = 0.20 map_w = 8 * cell_in map_h = 12 * cell_in margin_l = 0.72 margin_r = 0.15 margin_t = 0.40 gap_plots = 0.10 gap_rows = 0.22 cbar_gap = 0.42 cbar_h = 0.10 cbar_lab = 0.78 margin_b = cbar_gap + cbar_h + cbar_lab fig_w = margin_l + 2 * map_w + gap_plots + margin_r fig_h = margin_t + 2 * map_h + gap_rows + margin_b fig = plt.figure(figsize=(fig_w, fig_h)) left_f = margin_l / fig_w right_f = (margin_l + 2 * map_w + gap_plots) / fig_w bot_f = margin_b / fig_h top_f = 1.0 - margin_t / fig_h gs = gridspec.GridSpec( 2, 2, figure=fig, left=left_f, right=right_f, bottom=bot_f, top=top_f, wspace=gap_plots / (2 * map_w + gap_plots), hspace=gap_rows / map_h, ) heatmap_kw = dict( annot=False, cmap=cmap, center=0, vmin=-DISPLAY_VMAX, vmax=DISPLAY_VMAX, xticklabels=row_labels, yticklabels=col_labels, cbar=False, linewidths=0.2, linecolor="0.75", square=True, ) for row_idx, chiller in enumerate(CHILLERS): mA, mB = matrices[chiller] ax_A = fig.add_subplot(gs[row_idx, 0]) ax_B = fig.add_subplot(gs[row_idx, 1]) sns.heatmap(mA, ax=ax_A, **heatmap_kw) sns.heatmap(mB, ax=ax_B, **heatmap_kw) ax_A.set_title(f"{chiller} - Plate A", fontweight="bold", pad=3) ax_B.set_title(f"{chiller} - Plate B", fontweight="bold", pad=3) ax_A.set_ylabel("Well Column", labelpad=2) ax_B.set_ylabel("") ax_B.set_yticklabels([]) if row_idx == 0: ax_A.set_xlabel("") ax_B.set_xlabel("") ax_A.set_xticklabels([]) ax_B.set_xticklabels([]) else: ax_A.set_xlabel("Well Row", labelpad=2) ax_B.set_xlabel("Well Row", labelpad=2) for ax in [ax_A, ax_B]: ax.tick_params(axis="both", which="both", length=0, width=0) cb_bottom = cbar_lab / fig_h cbar_ax = fig.add_axes([left_f, cb_bottom, right_f - left_f, cbar_h / fig_h]) norm = mpl.colors.Normalize(vmin=-DISPLAY_VMAX, vmax=DISPLAY_VMAX) sm = mpl.cm.ScalarMappable(norm=norm, cmap=cmap) sm.set_array([]) cb = fig.colorbar(sm, cax=cbar_ax, orientation="horizontal") cb.solids.set_edgecolor("face") cb.set_ticks(DISPLAY_TICKS) cb.set_ticklabels(["-1.0", "-0.5", "0.0", "0.5", "1.0"]) cb.ax.xaxis.set_ticks_position("bottom") cb.ax.xaxis.set_label_position("bottom") cb.set_label(r"Temperature bias $c_i$ ($^\circ$C)", fontsize=9, labelpad=4) cb.ax.tick_params(labelsize=8, length=2) out_png = OUTPUT_DIR / "spatial_bias_heatmaps_4panel.png" out_pdf = OUTPUT_DIR / "spatial_bias_heatmaps_4panel.pdf" plt.savefig(out_png, dpi=300, bbox_inches="tight") plt.savefig(out_pdf, bbox_inches="tight") plt.close() print(f"Saved: {out_png}") print(f"Saved: {out_pdf}") print("\nMoran's I summary:") for chiller in CHILLERS: I_A = morans_results[chiller]["I_A"] I_B = morans_results[chiller]["I_B"] print( f" {chiller} Plate A: {I_A:+.3f} Plate B: {I_B:+.3f} " f"mean: {(I_A + I_B) / 2:+.3f}" ) # --------------------------------------------------------------------------- # Save outputs # --------------------------------------------------------------------------- def save_results( all_bias_A: dict, all_bias_B: dict, all_bias_C_A: dict, all_bias_C_B: dict, all_count_C_A: dict, all_count_C_B: dict, morans_results: dict, spatial_summary: dict, ) -> None: data = {"Index": np.arange(96)} for chiller in CHILLERS: data[f"{chiller}_A"] = all_bias_A[chiller] data[f"{chiller}_B"] = all_bias_B[chiller] pd.DataFrame(data).to_csv(OUTPUT_DIR / "index_bias_all_chillers.csv", index=False) print(f"Saved: {OUTPUT_DIR / 'index_bias_all_chillers.csv'}") rows = [] for chiller in CHILLERS: r = morans_results[chiller] rows.append( { "chiller": chiller, "Plate_A_I": r["I_A"], "Plate_A_p": r["p_A"], "Plate_B_I": r["I_B"], "Plate_B_p": r["p_B"], } ) pd.DataFrame(rows).to_csv(OUTPUT_DIR / "morans_i_results.csv", index=False) print(f"Saved: {OUTPUT_DIR / 'morans_i_results.csv'}") summary_rows = [] for chiller in CHILLERS: for plate in PLATES: s = spatial_summary[chiller][plate] summary_rows.append( { "chiller": chiller, "plate_position": plate, "n_wells": s["n_wells"], "weighted_mean_C": s["weighted_mean"], "weighted_sigma_C": s["weighted_sigma"], "half_range_C": s["half_range"], "cmin_C": s["cmin"], "cmax_C": s["cmax"], } ) pd.DataFrame(summary_rows).to_csv(OUTPUT_DIR / "spatial_bias_summary.csv", index=False) print(f"Saved: {OUTPUT_DIR / 'spatial_bias_summary.csv'}") ci_rows = [] for chiller in CHILLERS: for plate in PLATES: bias = all_bias_C_A[chiller] if plate == "A" else all_bias_C_B[chiller] count = all_count_C_A[chiller] if plate == "A" else all_count_C_B[chiller] for idx in range(96): ci_rows.append( { "chiller": chiller, "plate_position": plate, "well_id": idx, "c_i_C": bias[idx], "n_events": int(count[idx]), } ) pd.DataFrame(ci_rows).to_csv(OUTPUT_DIR / "spatial_ci_fields.csv", index=False) print(f"Saved: {OUTPUT_DIR / 'spatial_ci_fields.csv'}") # --------------------------------------------------------------------------- # Entry point # --------------------------------------------------------------------------- if __name__ == "__main__": if not DATA_FILE.exists(): raise FileNotFoundError( f"{DATA_FILE} not found. Run build_freezing_temps_pooled.py first." ) df = pd.read_csv(DATA_FILE) df["experiment_id"] = pd.to_numeric(df["experiment_id"], errors="coerce") df["well_id"] = pd.to_numeric(df["well_id"], errors="coerce") df["freezing_temp_BT"] = pd.to_numeric(df["freezing_temp_BT"], errors="coerce") df = df.dropna(subset=["experiment_id", "well_id", "freezing_temp_BT"]) df["experiment_id"] = df["experiment_id"].astype(int) df["well_id"] = df["well_id"].astype(int) print(f"Loaded {len(df)} freezing events from {DATA_FILE}") print( df.groupby(["chiller", "plate_position"]) .agg( n_events=("freezing_temp_BT", "count"), n_experiments=("experiment_id", "nunique"), ) .to_string() ) print() results = process_pooled(df) if len(results["all_bias_C_A"]) == 2: plot_spatial_heatmaps_4panel( results["all_bias_C_A"], results["all_bias_C_B"], results["morans"], ) else: print("Need both chillers to produce the 4-panel figure.") save_results( results["all_bias_A"], results["all_bias_B"], results["all_bias_C_A"], results["all_bias_C_B"], results["all_count_C_A"], results["all_count_C_B"], results["morans"], results["spatial_summary"], )