diff --git a/src/ops_model/post_process/combination/pca_optimization/__init__.py b/src/ops_model/post_process/combination/pca_optimization/__init__.py index b6eb363..8b3cb6e 100644 --- a/src/ops_model/post_process/combination/pca_optimization/__init__.py +++ b/src/ops_model/post_process/combination/pca_optimization/__init__.py @@ -177,6 +177,7 @@ _compute_and_plot_embeddings, _score_consistency, _score_distinctiveness, + _score_ebi_plus, ) from ops_model.post_process.combination.pca_optimization.phase1 import ( pca_sweep_pooled_signal, diff --git a/src/ops_model/post_process/combination/pca_optimization/aggregation.py b/src/ops_model/post_process/combination/pca_optimization/aggregation.py index 7e16242..bddaf56 100644 --- a/src/ops_model/post_process/combination/pca_optimization/aggregation.py +++ b/src/ops_model/post_process/combination/pca_optimization/aggregation.py @@ -180,9 +180,11 @@ def _score_single_reporter_metrics( "distinctiveness", "corum", "chad", + "ebi_plus", "distinctiveness_all", "corum_all", "chad_all", + "ebi_plus_all", ) } # Per-gene / per-complex DataFrames (raw mAP scores per entity for @@ -194,12 +196,14 @@ def _score_single_reporter_metrics( result["distinctiveness_df"] = None result["corum_df"] = None result["chad_df"] = None + result["ebi_plus_df"] = None try: from ops_utils.analysis.map_scores import ( phenotypic_activity_assesment, phenotypic_distinctivness, phenotypic_consistency_corum, phenotypic_consistency_manual_annotation, + phenotypic_ebi_plus, ) g_norm = normalize_guide_adata(g_raw.copy(), norm_method) @@ -225,6 +229,16 @@ def _score_single_reporter_metrics( result["distinctiveness_all"] = result["distinctiveness"] result["distinctiveness_df"] = dist_df + ebi_plus_df, ebi_plus_ratio = phenotypic_ebi_plus( + g_norm, + plot_results=False, + null_size=null_size, + distance=distance, + ) + result["ebi_plus"] = float(ebi_plus_ratio) + result["ebi_plus_all"] = result["ebi_plus"] + result["ebi_plus_df"] = ebi_plus_df + e_norm = aggregate_to_level( g_norm, "gene", preserve_batch_info=False, subsample_controls=False ) @@ -315,6 +329,7 @@ def _load_per_unit_blocks(per_unit_dir, norm_method, _logger, distance="cosine") "distinctiveness": reporter_metrics.get("distinctiveness_df"), "corum": reporter_metrics.get("corum_df"), "chad": reporter_metrics.get("chad_df"), + "ebi_plus": reporter_metrics.get("ebi_plus_df"), } report_rows.append( @@ -332,14 +347,17 @@ def _load_per_unit_blocks(per_unit_dir, norm_method, _logger, distance="cosine") "distinctiveness": reporter_metrics["distinctiveness"], "corum": reporter_metrics["corum"], "chad": reporter_metrics["chad"], + "ebi_plus": reporter_metrics["ebi_plus"], "distinctiveness_all": reporter_metrics["distinctiveness_all"], "corum_all": reporter_metrics["corum_all"], "chad_all": reporter_metrics["chad_all"], + "ebi_plus_all": reporter_metrics["ebi_plus_all"], } ) _logger.info( f" {sig}: {g.n_obs} guides x {g.n_vars} PCs @ {g.uns.get('pca_threshold', '?')} | " f"act={reporter_metrics['activity']:.1%} dist={reporter_metrics['distinctiveness']:.1%} " + f"ebi+={reporter_metrics['ebi_plus']:.1%} " f"corum={reporter_metrics['corum']:.1%} chad={reporter_metrics['chad']:.1%} | " f"all: dist={reporter_metrics['distinctiveness_all']:.1%} " f"corum={reporter_metrics['corum_all']:.1%} chad={reporter_metrics['chad_all']:.1%}" @@ -477,6 +495,10 @@ def _pivot(metric_key, key_col, value_col, out_csv): "distinctiveness", "perturbation", "mean_average_precision", overlay_dir / "gene_reporter_distinctiveness_raw.csv", ) + _pivot( + "ebi_plus", "ebi_group", "mean_average_precision", + overlay_dir / "group_reporter_ebi_plus_raw.csv", + ) _pivot( "corum", "complex_num", "mean_average_precision", overlay_dir / "complex_reporter_corum_consistency.csv", diff --git a/src/ops_model/post_process/combination/pca_optimization/embeddings.py b/src/ops_model/post_process/combination/pca_optimization/embeddings.py index c1f8234..547bec0 100644 --- a/src/ops_model/post_process/combination/pca_optimization/embeddings.py +++ b/src/ops_model/post_process/combination/pca_optimization/embeddings.py @@ -349,6 +349,68 @@ def _score_distinctiveness( return None, 0.0 +def _score_ebi_plus( + adata_guide, + activity_map, + r, + total_feats, + plots_dir, + metrics_dir, + plt, + _logger, + distance="cosine", + suffix="", +): + """Run EBI+ scoring (guide-level, complex-grouped), save CSV and plot. + + EBI+ is a guide-level metric like distinctiveness: every guide's positive + group is its gene's EBI complex (else the gene itself), so complex members + collapse into one group and non-complex genes reduce to distinctiveness. + Returns ``(ebi_plus_map, ebi_plus_ratio)`` or ``(None, 0.0)`` on failure. + Computed across all geneKOs (``suffix`` only affects output filenames). + """ + label = "all geneKOs" + if activity_map is None: + return None, 0.0 + try: + from ops_utils.analysis.map_scores import phenotypic_ebi_plus + + _logger.info(f"Running EBI+ ({label})...") + ebi_plus_map, ebi_plus_ratio = phenotypic_ebi_plus( + adata_guide, + plot_results=False, + null_size=100_000, + distance=distance, + ) + ebi_plus_map.to_csv( + metrics_dir / f"phenotypic_ebi_plus{suffix}.csv", index=False + ) + _logger.info(f" EBI+ ({label}): {ebi_plus_ratio:.1%}") + + try: + fig, ax = plt.subplots(figsize=(8, 7)) + plot_map_scatter( + ax, ebi_plus_map, f"EBI+ ({label})", ebi_plus_ratio, show_ntc=False + ) + fig.suptitle( + f"EBI+ ({label}) — {total_feats} features", + fontsize=13, fontweight="bold", + ) + fig.tight_layout() + fig.savefig( + plots_dir / f"map_ebi_plus{suffix}.png", dpi=150, bbox_inches="tight" + ) + plt.close(fig) + _logger.info(f" Saved plots/map_ebi_plus{suffix}.png") + except Exception as exc: + _logger.warning(f" EBI+ plot failed: {exc}") + + return ebi_plus_map, ebi_plus_ratio + except Exception as exc: + _logger.error(f" EBI+ ({label}) failed: {exc}") + return None, 0.0 + + def _score_consistency( adata_gene, activity_map, diff --git a/src/ops_model/post_process/combination/pca_optimization/phase2.py b/src/ops_model/post_process/combination/pca_optimization/phase2.py index 9ba8f70..7c4f4d6 100644 --- a/src/ops_model/post_process/combination/pca_optimization/phase2.py +++ b/src/ops_model/post_process/combination/pca_optimization/phase2.py @@ -54,6 +54,7 @@ _compute_and_plot_embeddings, _score_consistency, _score_distinctiveness, + _score_ebi_plus, ) from ops_model.post_process.combination.pca_optimization.sweep_core import ( _init_sweep_logger, @@ -212,6 +213,17 @@ def aggregate_channels( _logger, distance=distance, ) + ebi_plus_map, ebi_plus_ratio = _score_ebi_plus( + adata_guide, + activity_map, + r, + total_feats, + plots_dir, + metrics_dir, + plt, + _logger, + distance=distance, + ) corum_map, corum_ratio, chad_map, chad_ratio, ebi_map, ebi_ratio = _score_consistency( adata_gene, activity_map, @@ -302,6 +314,7 @@ def aggregate_channels( metric_maps={ "Activity": activity_map, "Distinctiveness": dist_map, + "EBI+": ebi_plus_map, "EBI": ebi_map, "CHAD": chad_map, "CORUM": corum_map, @@ -949,6 +962,17 @@ def apply_second_pass_pca( _logger, distance=distance, ) + ebi_plus_map, ebi_plus_ratio = _score_ebi_plus( + adata_guide, + activity_map, + r, + total_feats, + plots_dir, + metrics_dir, + plt, + _logger, + distance=distance, + ) corum_map, corum_ratio, chad_map, chad_ratio, ebi_map, ebi_ratio = _score_consistency( adata_gene, activity_map, @@ -1007,6 +1031,7 @@ def apply_second_pass_pca( metric_maps={ "Activity": activity_map, "Distinctiveness": dist_map, + "EBI+": ebi_plus_map, "EBI": ebi_map, "CHAD": chad_map, "CORUM": corum_map, diff --git a/src/ops_model/post_process/combination/titration/titration.py b/src/ops_model/post_process/combination/titration/titration.py index b9e5f17..27f338b 100644 --- a/src/ops_model/post_process/combination/titration/titration.py +++ b/src/ops_model/post_process/combination/titration/titration.py @@ -59,7 +59,7 @@ DOWNSAMPLE_RATIO = 0.75 # multiply cell count by this each step MIN_CELLS = 5_000 # stop titrating below this NULL_SIZE = 10_000 # smaller null for speed (per-reporter) -METRICS = ("activity", "distinctiveness", "corum", "chad", "ebi") +METRICS = ("activity", "distinctiveness", "corum", "chad", "ebi", "ebi_plus") SCALES = ("linear", "log2", "log10") # x-axis scale variants to save # Shared plot styling / labels (used by compare_titration_versions and below) @@ -69,6 +69,7 @@ "corum": "mediumpurple", "chad": "darkorange", "ebi": "crimson", + "ebi_plus": "deeppink", } TITRATION_RATIO_LABELS = { "activity": "% Active", @@ -76,6 +77,7 @@ "corum": "% CORUM consistent", "chad": "% CHAD consistent", "ebi": "% EBI consistent", + "ebi_plus": "% EBI+ significant", } TITRATION_MAP_LABELS = { "activity": "Activity mAP", @@ -83,6 +85,7 @@ "corum": "CORUM mAP", "chad": "CHAD mAP", "ebi": "EBI mAP", + "ebi_plus": "EBI+ mAP", } SCALE_LABEL_SHORT = {"linear": "linear", "log2": "log₂", "log10": "log₁₀"} @@ -330,6 +333,7 @@ def _score_all_metrics( phenotypic_consistency_corum, phenotypic_consistency_ebi, phenotypic_consistency_manual_annotation, + phenotypic_ebi_plus, ) result = { @@ -343,6 +347,8 @@ def _score_all_metrics( "chad_map_mean": math.nan, "ebi_ratio": math.nan, "ebi_map_mean": math.nan, + "ebi_plus_ratio": math.nan, + "ebi_plus_map_mean": math.nan, } try: @@ -388,6 +394,23 @@ def _score_all_metrics( except Exception as exc: _logger.warning(f" Distinctiveness scoring failed: {exc}") + try: + ebi_plus_map, ebi_plus_ratio = phenotypic_ebi_plus( + g_copairs, + plot_results=False, + null_size=NULL_SIZE, + distance=distance, + ) + # EBI+ groups by complex, not perturbation, so subset_targets (a set of + # perturbation names) does not cleanly filter it — report the full ratio. + result["ebi_plus_ratio"] = float(ebi_plus_ratio) + if ebi_plus_map is not None and "mean_average_precision" in ebi_plus_map.columns: + result["ebi_plus_map_mean"] = float( + ebi_plus_map["mean_average_precision"].mean() + ) + except Exception as exc: + _logger.warning(f" EBI+ scoring failed: {exc}") + try: e_norm = aggregate_to_level( g_copairs, "gene", preserve_batch_info=False, subsample_controls=False @@ -503,6 +526,7 @@ def _run_titration_points( "corum_ratio", "corum_map_mean", "chad_ratio", "chad_map_mean", "ebi_ratio", "ebi_map_mean", + "ebi_plus_ratio", "ebi_plus_map_mean", ] rows = [] @@ -598,6 +622,7 @@ def _run_titration_points( f"corum={scores['corum_ratio']:.1%}±{scores.get('corum_ratio_sem', 0):.1%} " f"chad={scores['chad_ratio']:.1%}±{scores.get('chad_ratio_sem', 0):.1%} " f"ebi={scores['ebi_ratio']:.1%}±{scores.get('ebi_ratio_sem', 0):.1%} " + f"ebi+={scores['ebi_plus_ratio']:.1%}±{scores.get('ebi_plus_ratio_sem', 0):.1%} " f"({time.time() - t_step:.0f}s)" ) return pd.DataFrame(rows)