From 59b08e32f75b01102b21bb818787cbcda6b69280 Mon Sep 17 00:00:00 2001 From: Comp-era Date: Wed, 29 Jul 2026 07:34:35 +0000 Subject: [PATCH 1/8] fix: mode-aware scenario discovery + metric-correct melquiplot legend; rename script to lowercase Rename AnalyseBenchmarkResults.py -> analysebenchmarkresults.py. Accept a run1/ directory (default mode) or run1_analysis.tgz (-a mode) in scenario discovery so non-archive benchmarks are no longer rejected with 'no scenario directories found'. Build the melquiplot legend from get_reverse_bool(): distance metrics (irmsd/ilrmsd) keep 'lo <= M < hi' in Angstrom; dockq is unitless and higher-is-better, so it now reads 'lo < M <= hi' without the unit. --- ...kResults.py => analysebenchmarkresults.py} | 700 +++++++++++++++--- 1 file changed, 578 insertions(+), 122 deletions(-) rename analysis/{AnalyseBenchmarkResults.py => analysebenchmarkresults.py} (69%) diff --git a/analysis/AnalyseBenchmarkResults.py b/analysis/analysebenchmarkresults.py similarity index 69% rename from analysis/AnalyseBenchmarkResults.py rename to analysis/analysebenchmarkresults.py index a3ab1bb..e60c5de 100644 --- a/analysis/AnalyseBenchmarkResults.py +++ b/analysis/analysebenchmarkresults.py @@ -10,6 +10,19 @@ Usage: >python3 AnalyseBenchmarkResults.py + +Expected input structure +------------------------- +/ + / + / + run1/ + _caprieval/ + capri_ss.tsv + +(scenario is the OUTER directory, PDBid the INNER one -- this matches the +directory layout currently produced by haddock-runner. The previous +`//run1/` layout is no longer supported.) """ import argparse @@ -40,7 +53,7 @@ ) -__version__ = "1.1.1" # September 2025 +__version__ = "2.0.0" # July 2026 __author__ = ", ".join(( "BonvinLab", "Computational Structural Biology group", @@ -52,32 +65,22 @@ )) __dev__ = ( "Victor G.P. Reys", + "Shantanu Khatri", ) #################### # GLOBAL VARIABLES # #################### -# Define custom caprieval steps names -# NOTE: "Feel free to modify this `CAPRIEVAL_STEPS` dict content to fit your experiment" # noqa : E501 -# This dict must have: -# as keys -> Index of the caprieval stage (used to parse data) -# as values -> Name to give to this stage (used as legends plots) -# NOTE: e.g.: for the following run: [topoaa, rigidbody, caprieval, seletop, caprieval.1, flexref, caprieval.2, emref, caprieval.3, clustfcc, seletopclusts, caprieval.4] # noqa : E501 -# CAPRIEVAL_STEPS = { - # '02': 'rigidbody', - # '04': 'seletop 200', - # '06': 'flexref', - # '08': 'emref', - # '11': 'top 4 models per fcc clust', - # } -CAPRIEVAL_STEPS = { - '02': 'rigidbody', - '04': 'seletop 200', - '06': 'flexref', - '08': 'emref', - '11': 'top 4 models per fcc clust', - } +# Caprieval stage names are now auto-detected at runtime from the module each +# caprieval evaluated (see detect_stage_modules() / DETECTED_STAGE_NAMES), so +# they are correct per benchmark type without manual editing. +# +# This dict is only an OPTIONAL manual override / fallback, used solely for a +# stage that could not be auto-detected. Leave it empty to let auto-detection +# handle everything (undetected stages then fall back to `_caprieval`). +# To force a custom label for a specific stage, add e.g. '02': 'rigidbody'. +CAPRIEVAL_STEPS: dict[str, str] = {} # Set threshold of top X structures to take into account TOP_X_THRESHOLDS = (1, 5, 10, 20, 50, 100, 200, 500, 1000, ) @@ -123,16 +126,9 @@ "Missing": (-2, -0.5), }, }, - # FIXME: Optimize irmsd and dockq values for glycan + # NOTE: protein-glycan is assessed with ilrmsd (the default metric for this + # type); irmsd is not used. dockq follows the protein-peptide convention. "glycan": { - "irmsd": { - "High": (0, 0.5), - "Medium": (0.5, 1), - "Acceptable": (1, 2), - "Near-acceptable": (2, 3), - "Low": (3, 99999), - "Missing": (-2, -0.5), - }, "dockq": { "High": (0.895, 1), "Medium": (0.71, 0.895), @@ -150,6 +146,45 @@ "Missing": (-2, -0.5), }, }, + #irmsd and dockq values for dna + "dna": { + "irmsd": { + "High": (0, 1), + "Medium": (1, 2), + "Acceptable": (2, 4), + "Near-acceptable": (4, 6), + "Low": (6, 99999), + "Missing": (-2, -0.5), + }, + "dockq": { + "High": (0.8, 1), + "Medium": (0.6, 0.8), + "Acceptable": (0.5, 0.6), + "Near-acceptable": (0.4, 0.5), + "Low": (0, 0.4), + "Missing": (-2, -0.5), + }, + }, + # NOTE: protein-ligand (small molecule) is assessed with ilrmsd, same as + # glycan + "small_molecule": { + "ilrmsd": { + "High": (0, 1), + "Medium": (1, 2), + "Acceptable": (2, 3), + "Near-acceptable": (3, 4), + "Low": (4, 99999), + "Missing": (-2, -0.5), + }, + "dockq": { + "High": (0.895, 1), + "Medium": (0.71, 0.895), + "Acceptable": (0.43, 0.71), + "Near-acceptable": (0.35, 0.43), + "Low": (0, 0.35), + "Missing": (-2, -0.5), + }, + }, } PERFORMANCES_CLASSES = ALL_PERFORMANCES_CLASSES["protein"] @@ -159,7 +194,7 @@ "Medium": "lightgreen", "Acceptable": "lightblue", "Near-acceptable": "gainsboro", - "Low": "white", + "Low": "lightcoral", # was "white" (invisible on white background) "Missing": "dimgrey", } @@ -290,7 +325,7 @@ def gen_graph( # Orient X labels ax.set_xticks(ax.get_xticks()) ax.set_xticklabels(ax.get_xticklabels(), rotation=30, ha='right') - ylabel = "Nb. entries" if not percentage else "% sucess rate" + ylabel = "Nb. entries" if not percentage else "% success rate" ax.set_ylabel(ylabel) @@ -400,8 +435,18 @@ def adjacent_values( return lower_adjacent_value, upper_adjacent_value +# Auto-detected {stage_index: module_name}, filled at runtime by +# detect_stage_modules(). Takes precedence over the hardcoded CAPRIEVAL_STEPS. +DETECTED_STAGE_NAMES: dict[str, str] = {} + + def stage_name(cname: str) -> str: - """Try to return the user defined stage name, or return default. + """Return a human-readable name for a caprieval stage. + + Preference order: + 1. auto-detected module name (the module this caprieval evaluated), + 2. the hardcoded CAPRIEVAL_STEPS mapping, + 3. a `_caprieval` fallback. Parameters ---------- @@ -413,11 +458,55 @@ def stage_name(cname: str) -> str: name : str Name of the stage """ + if cname in DETECTED_STAGE_NAMES: + return DETECTED_STAGE_NAMES[cname] try: - name = CAPRIEVAL_STEPS[cname] + return CAPRIEVAL_STEPS[cname] except KeyError: - name = f"{cname}_caprieval" - return name + return f"{cname}_caprieval" + + +def detect_stage_modules( + dtmapper: dict, + read_from_archive: bool = False, + ) -> dict[str, str]: + """Auto-name caprieval stages by the module they evaluated. + + Every caprieval `capri_ss.tsv` lists model paths like + `../../10_seletopclusts/cluster_1_model_1.pdb`; the parent directory + (`10_seletopclusts`) is the module being evaluated. This reads one file + per stage (from the first scenario/target) and strips the numeric prefix + to get the module name (`seletopclusts`). Duplicate module names across + stages are disambiguated with their stage index, e.g. `flexref (07)`. + + Falls back silently (empty entry) for any stage it cannot read, so the + hardcoded CAPRIEVAL_STEPS / `_caprieval` naming still applies. + """ + mapping: dict[str, str] = {} + if not dtmapper: + return mapping + scen = sorted(dtmapper)[0] + for stage, pdbmap in dtmapper[scen].items(): + if not pdbmap: + continue + first_pdb = next(iter(pdbmap)) + try: + data = load_caprieval_data(pdbmap[first_pdb], read_from_archive) + model_path = next(iter(data)) + module_dir = os.path.basename(os.path.dirname(str(model_path))) + module = re.sub(r'^\d+_', '', module_dir) + if module: + mapping[stage] = module + except Exception: + continue # leave undetected; stage_name() falls back + # disambiguate repeated module names by appending the stage index + counts: dict[str, int] = {} + for m in mapping.values(): + counts[m] = counts.get(m, 0) + 1 + for stage in list(mapping): + if counts[mapping[stage]] > 1: + mapping[stage] = f"{mapping[stage]} ({stage})" + return mapping def gen_full_comparison_violins( @@ -537,16 +626,29 @@ def gen_full_comparison_violins( [so.replace('scenario-', '') for so in scenars_order], loc='outside lower center', ncols=4, - title="Screnarios", + title="Scenarios", ) plt.gca().set_ylim(bottom=0) - # adjust border to let annotations fit inside graph - fig.subplots_adjust(left=0.08, top=0.95, bottom=0.12, right=0.98) + # Adjust border to let annotations fit inside graph. The left margin + # needs to hold both the y-axis label and the row-title annotations + # (drawn at a small fixed point-offset from it). That fixed offset + # becomes proportionally huge on narrow figures (e.g. a single caprieval + # stage -> nb_steps=1), so scale the left fraction to preserve roughly + # a constant absolute margin instead of leaving it fixed. + fig_width_in = fig.get_size_inches()[0] + left_margin = max(0.08, min(0.35, 1.3 / fig_width_in)) + fig.subplots_adjust(left=left_margin, top=0.95, bottom=0.12, right=0.98) # save figure - plt.savefig(f"{basepath}_violins.png", format='png', dpi=DPI) + plt.savefig( + f"{basepath}_violins.png", + format='png', + dpi=DPI, + bbox_inches='tight', + pad_inches=0.3, + ) return @@ -576,7 +678,7 @@ def gen_full_comparison_melquiplots( """ # Clear pervious instances of matplotlib clear_plt() - + # Set progression variables all_generated_melquis: list[str] = [] processed = 0 @@ -673,10 +775,23 @@ def make_scenar_melquiplots( # Add figure title fig.suptitle(title, fontsize=16) # Add Legend + # Metric-aware label: distance metrics (irmsd/ilrmsd) are in Angstrom and + # lower-is-better (`lo <= M < hi`); dockq is unitless and higher-is-better + # (`lo < M <= hi`), so both the unit and the inequality direction flip. + higher_is_better = get_reverse_bool(perf_dtype) + unit = "" if higher_is_better else r"$\AA$" + + def _perf_label(perfclass: str) -> str: + lo, hi = dtype_perf_classes[perfclass] + metric = perf_dtype.upper() + if higher_is_better: + return rf"{perfclass} | {lo} < {metric} <= {hi}{unit}" + return rf"{perfclass} | {lo} <= {metric} < {hi}{unit}" + legend_data = [ ( plt.Rectangle((0, 0), 1, 1, fc=COLORS_MAPPER[perfclass]), - rf'{perfclass} | {dtype_perf_classes[perfclass][0]} <= {perf_dtype.upper()} < {dtype_perf_classes[perfclass][1]}$\AA$', # noqa : E501 + _perf_label(perfclass), ) for perfclass in PERF_ORDER ] @@ -695,7 +810,13 @@ def make_scenar_melquiplots( # save figure figpath = f"{basepath}_{title}_melquiplot.png" - plt.savefig(figpath, format='png', dpi=DPI) + plt.savefig( + figpath, + format='png', + dpi=DPI, + bbox_inches='tight', + pad_inches=0.3, + ) return figpath @@ -939,71 +1060,188 @@ def gen_full_comparison_barplots( title="performance classes", ) - # Adjust border to let annotations fit inside graph - fig.subplots_adjust(left=0.15, top=0.9, bottom=0.15, right=0.98) + # Adjust border to let annotations fit inside graph. Same reasoning as + # in `gen_full_comparison_violins`: scale the left margin so a narrow + # figure (e.g. a single caprieval stage -> nb_cols=1) still has enough + # absolute room for the y-axis label and row-title annotations, instead + # of clipping them. + fig_width_in = fig.get_size_inches()[0] + left_margin = max(0.15, min(0.35, 1.3 / fig_width_in)) + fig.subplots_adjust(left=left_margin, top=0.9, bottom=0.15, right=0.98) # Save figure - plt.savefig(f"{basepath}_capribarpolots.png", format='png', dpi=DPI) + plt.savefig( + f"{basepath}_capribarplots.png", + format='png', + dpi=DPI, + bbox_inches='tight', + pad_inches=0.3, + ) return def get_pdb_entries(basepath: str) -> list: - """Retrieve list of PDBid. + """Retrieve list of immediate subdirectory names. + + Generic directory lister -- reused both to list scenario names at the + top level of the benchmark directory, and to list PDBid/target names + inside each scenario directory. Kept its original name for API + stability, even though it's no longer PDBid-specific. Parameters ---------- basepath : str - Path to a scenario directory containing pdb entries. + Path to a directory whose immediate subdirectories should be listed. Return ------ - pdbids : str - List of pdb entries benchmarked in this scenario. + entries : list + List of subdirectory names found directly under `basepath`. """ - pdbids = [ + entries = [ pdbid_path.stem for pdbid_path in Path(basepath).glob("*/") if pdbid_path.is_dir() ] - return pdbids + return entries + +def _find_scenario_data_dir( + candidate_path: str, + _max_depth: int = 3, + ) -> Optional[str]: + """Find the directory that actually holds `/run1/` entries. -def get_scenarios_names(basepath: str) -> list: - """Retrieve list of scenario names. + Checks whether `candidate_path` directly contains one or more + ``/run1/`` subdirectories. If not, and `candidate_path` has + exactly one subdirectory, recurses into that single subdirectory (up to + `_max_depth` levels) to transparently see through extra wrapper + directories -- e.g. a scenario directory that ended up nested one level + too deep, such as ``///run1/`` instead of + ``//run1/``. Only drills down when there is a single, + unambiguous subdirectory to descend into, so it won't misfire on + directories that hold unrelated clutter (e.g. output files). Parameters ---------- - basepath : str - Path to the benchmark directory to analyse containing X scenarios. + candidate_path : str + Directory to inspect (should end with '/'). + _max_depth : int + Maximum number of wrapper levels to look through. Return ------ - scenarios_names : list - List of tested scenario names. + resolved : Optional[str] + Path to the directory directly containing `/run1/` entries, + or None if none could be found within `_max_depth` levels. """ - scenario_paths = glob.glob(f'{basepath}scenario*/') - scenarios_names = [sp.split('/')[-2] for sp in scenario_paths] - return scenarios_names + for entry in Path(candidate_path).glob("*/"): + # A genuine target holds either a `run1/` run directory (default, + # non-archive mode) or a `run1_analysis.tgz` archive (`-a` mode). + # Accept both so scenario discovery works regardless of which mode + # the benchmark was produced/analysed in. + if entry.is_dir() and ( + (entry / "run1").is_dir() + or (entry / "run1_analysis.tgz").is_file() + ): + return candidate_path + if _max_depth <= 0: + return None + subdirs = [entry for entry in Path(candidate_path).glob("*/") if entry.is_dir()] + if len(subdirs) == 1: + return _find_scenario_data_dir( + f"{subdirs[0]}/", + _max_depth=_max_depth - 1, + ) + return None -def _make_run_path(basepath: str, pdbid: str, scenario: str) -> str: - """Return the haddock3 run directory path. +def _looks_like_scenario_dir(scenario_path: str) -> bool: + """Heuristic check for whether a directory is a real scenario dir. - When scenario is empty the run layout is flat — modules live directly - in ``{basepath}{pdbid}/`` with no extra ``run1`` subdirectory. - When a scenario name is present the legacy structure - ``{basepath}{pdbid}/{scenario}/run1/`` is used. + A genuine scenario directory (as produced by haddock-runner) contains + at least one ``/run1/`` subdirectory, possibly behind a single + extra wrapper level (see `_find_scenario_data_dir`). This is used to + filter out unrelated top-level clutter that can end up alongside the + scenario directories -- most commonly a previous analysis output + directory (e.g. `Analysis/`) left inside the benchmark directory, or + stray/hidden directories -- so they aren't mistaken for scenarios. + + Parameters + ---------- + scenario_path : str + Candidate scenario directory path (should end with '/'). + + Return + ------ + bool + True if a `/run1/` pattern is found inside `scenario_path`, + directly or behind a single unambiguous wrapper level. + """ + return _find_scenario_data_dir(scenario_path) is not None + + +def _resolve_scenario_dir(basepath: str, scenario: str) -> str: + """Resolve a scenario name to its actual data directory on disk. + + Normally this is simply ``{basepath}{scenario}/``. However, some + benchmark directories end up with an extra wrapper level in between + (e.g. produced by manual reorganisation), such as:: + + ////run1/ + + instead of the expected:: + + ///run1/ + + This transparently drills down through such wrapper levels via + `_find_scenario_data_dir`. + + Parameters + ---------- + basepath : str + Root directory containing scenario directories. + scenario : str + Scenario name. + + Return + ------ + resolved : str + Path to the directory that actually contains `/run1/` + subdirectories for this scenario (or the original, unresolved + candidate if no unambiguous wrapper level could be found -- in + which case downstream existence checks will raise a clear error). """ - if scenario: - return f"{basepath}{pdbid}/{scenario}/run1/" - return f"{basepath}{pdbid}/" + candidate = f"{basepath}{scenario}/" + resolved = _find_scenario_data_dir(candidate) + return resolved if resolved is not None else candidate -def _make_archive_path(basepath: str, pdbid: str, scenario: str) -> str: - """Return the haddock3 run analysis archive path.""" - if scenario: - return f"{basepath}{pdbid}/{scenario}/run1_analysis.tgz" - return f"{basepath}{pdbid}/run1_analysis.tgz" +def get_scenario_entries(basepath: str) -> list: + """Retrieve list of scenario directory names at the top of `basepath`. + + Like `get_pdb_entries`, but filters out any top-level directory that + doesn't actually look like a scenario directory (see + `_looks_like_scenario_dir`). This protects against stray directories + such as a leftover `Analysis/` output folder, hidden directories, etc. + being mistaken for scenarios. + + Parameters + ---------- + basepath : str + Path to the benchmark directory to analyse. + + Return + ------ + scenarios : list + List of directory names under `basepath` that look like real + scenario directories. + """ + scenarios = [ + name for name in get_pdb_entries(basepath) + if _looks_like_scenario_dir(f'{basepath}{name}/') + ] + return scenarios def scenario_name_2_threshold(scenar_name: str) -> float: @@ -1073,7 +1311,7 @@ def hd3_module_2_stage(indexed_modulename: str) -> str: ---------- indexed_modulename : str Directory name of an indexed haddock3 module name. - + Return ------ @@ -1084,6 +1322,34 @@ def hd3_module_2_stage(indexed_modulename: str) -> str: return stage +def _make_run_path(scenario_basepath: str, pdbid: str) -> str: + """Return the haddock3 run directory path. + + Parameters + ---------- + scenario_basepath : str + Path to the scenario's own directory (already resolved -- may be + ``{basepath}{scenario}/`` in multi-scenario mode, or just + ``{basepath}`` itself in single-scenario mode). + pdbid : str + Target/PDBid name. + """ + return f"{scenario_basepath}{pdbid}/run1/" + + +def _make_archive_path(scenario_basepath: str, pdbid: str) -> str: + """Return the haddock3 run analysis archive path. + + Parameters + ---------- + scenario_basepath : str + Path to the scenario's own directory (already resolved). + pdbid : str + Target/PDBid name. + """ + return f"{scenario_basepath}{pdbid}/run1_analysis.tgz" + + def map_data( basepath: str, subset_scenarios: Optional[list[str]] = None, @@ -1091,16 +1357,27 @@ def map_data( ) -> dict[str, dict[str, dict[str, dict[str, Union[str, list[str, str]]]]]]: """Map data in one analysis dict to accessit easily. + Expects the benchmark directory to be laid out as produced by + haddock-runner:: + + ///run1/_caprieval/capri_ss.tsv + + i.e. scenario is the outer directory and PDBid the inner one. This + holds whether `basepath` contains a single scenario directory or + several -- there is no special-casing, the code below just handles + however many scenario directories are found. + Parameters ---------- basepath : str - Path where the benchmarking scenarios can be found. + Path to the root directory containing one or more scenario + directories. subset_scenarios : Optional[list[str]] List of scenario names on which to perform the analysis. read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False - + Return ------ dtmap : dict[str, dict[str, dict[str, dict[str, str]]]] @@ -1118,27 +1395,44 @@ def map_data( """ # Initiate mapper dtmap: dict[str, dict[str, dict[str, dict[str, Union[str, list[str, str]]]]]] = {} - # Gather all directories - pdbids = get_pdb_entries(basepath) - # Initiate scenarios names holder - all_scenarios: list[str] = [] + # `basepath` is always the root containing one or more scenario + # directories -- this naturally handles both a single scenario and + # multiple scenarios without any special-casing, since the loops below + # just iterate over however many scenarios are found. if subset_scenarios: - all_scenarios += subset_scenarios - # Search for all available scenarios + all_scenarios = subset_scenarios else: - # Gather all scenarios (deduplicated across all targets) - scenarios_set: set[str] = set() - for pdbid in pdbids: - scenarios_set.update(get_scenarios_names(f'{basepath}{pdbid}/')) - all_scenarios = sorted(scenarios_set) if scenarios_set else [''] + # Scenarios are the top-level directories of basepath (filtered to + # exclude stray/unrelated directories, e.g. a leftover `Analysis/` + # output folder sitting inside the benchmark directory). + all_scenarios = sorted(get_scenario_entries(basepath)) + assert all_scenarios, ( + f"[ERROR] no scenario directories found directly under `{basepath}`. " + "Expected layout is `///run1/` -- " + "this applies whether there is one scenario or several." + ) + + def scenario_dir(scenario: str) -> str: + """Resolve a scenario name to its actual directory on disk.""" + return _resolve_scenario_dir(basepath, scenario) + + # pdbids are gathered per scenario, deduplicated across all scenarios. + pdbids_set: set[str] = set() + for scenario in all_scenarios: + pdbids_set.update(get_pdb_entries(scenario_dir(scenario))) + pdbids = sorted(pdbids_set) + assert pdbids, ( + f"[ERROR] no target/PDBid directories found under any scenario in " + f"`{basepath}`." + ) # Make sure all pdbs have all scenarios - for pdbid in pdbids: - for scenario in all_scenarios: + for scenario in all_scenarios: + for pdbid in pdbids: # Default behavior if not read_from_archive: - scenar_rundir = _make_run_path(basepath, pdbid, scenario) + scenar_rundir = _make_run_path(scenario_dir(scenario), pdbid) assert os.path.exists(scenar_rundir), \ ( f"[ERROR] could not find scenario `{scenario}` " @@ -1146,13 +1440,13 @@ def map_data( ) # Case where we need to search data in the archive else: - analysis_archive = _make_archive_path(basepath, pdbid, scenario) + analysis_archive = _make_archive_path(scenario_dir(scenario), pdbid) assert os.path.exists(analysis_archive), \ ( f"[ERROR] could not find scenario `{scenario}` " f"archive for entry `{pdbid}` at: {analysis_archive}" ) - + all_caprieval_stages = [] # Add scenario data to data maper for scenario in all_scenarios: @@ -1160,9 +1454,9 @@ def map_data( for pdbid in pdbids: if not read_from_archive: # Generate scenario basepath - scenario_bp = _make_run_path(basepath, pdbid, scenario) + scenario_bp = _make_run_path(scenario_dir(scenario), pdbid) else: - scenario_bp = _make_archive_path(basepath, pdbid, scenario) + scenario_bp = _make_archive_path(scenario_dir(scenario), pdbid) # Retrieve caprieval stages caprieval_stages = get_caprieval_stages( scenario_bp, @@ -1170,17 +1464,20 @@ def map_data( ) all_caprieval_stages += caprieval_stages all_caprieval_stages = sorted(list(set(all_caprieval_stages))) - + assert all_caprieval_stages, ( + f"[ERROR] no caprieval stages found anywhere under `{basepath}`." + ) + # Make sure all stages are computed for all pdb in all scenarios... for scenario in all_scenarios: for pdbid in pdbids: - archive_path = _make_archive_path(basepath, pdbid, scenario) + archive_path = _make_archive_path(scenario_dir(scenario), pdbid) if read_from_archive: tararchive = tarfile.open(archive_path, "r:gz") for stage in all_caprieval_stages: if not read_from_archive: # Build caprieval tsv filepath - caprieval_tsv_path = f"{_make_run_path(basepath, pdbid, scenario)}{stage}_caprieval/capri_ss.tsv" # noqa : E501 + caprieval_tsv_path = f"{_make_run_path(scenario_dir(scenario), pdbid)}{stage}_caprieval/capri_ss.tsv" # noqa : E501 assert os.path.exists(caprieval_tsv_path), \ ( f"[ERROR] could not access CAPRIEVAL results file at: " @@ -1203,7 +1500,7 @@ def map_data( f" - Target `{pdbid}`" ) if read_from_archive: - tararchive.close() + tararchive.close() # Gather all data for scenario in all_scenarios: @@ -1213,10 +1510,10 @@ def map_data( for pdbid in pdbids: if not read_from_archive: # Build caprieval tsv filepath - caprieval_tsv_path = f"{_make_run_path(basepath, pdbid, scenario)}{stage}_caprieval/capri_ss.tsv" # noqa : E501 + caprieval_tsv_path = f"{_make_run_path(scenario_dir(scenario), pdbid)}{stage}_caprieval/capri_ss.tsv" # noqa : E501 else: caprieval_tsv_path = [ - _make_archive_path(basepath, pdbid, scenario), + _make_archive_path(scenario_dir(scenario), pdbid), f"run1_analysis/{stage}_caprieval_analysis/capri_ss.tsv", ] # Hold datapath @@ -1488,8 +1785,8 @@ def best_perfs( # Get best performing value best_perfomance = best_function(all_perfs) return best_perfomance - - + + def dtype_to_best_function(dtype: str) -> Callable: """Get best function. @@ -1505,7 +1802,7 @@ def dtype_to_best_function(dtype: str) -> Callable: """ _function = max if get_reverse_bool(dtype) else min return _function - + def get_reverse_bool(dt_type: str) -> bool: """Check sorting order. @@ -1689,6 +1986,7 @@ def gen_outputdir(outputdir: str) -> None: def set_output_path( benchmark_directory: str, outputpath: str, + label: Optional[str] = None, ) -> tuple[str, str]: """Initiate output paths and directory. @@ -1698,18 +1996,22 @@ def set_output_path( Path to the benchmark directory to analyse outputpath : str Path to directory where to store the results + label : Optional[str] + Custom name to use for output filenames and plot titles instead of + the benchmark directory's basename (e.g. useful when the directory + is named something generic like `test-analysis`). Return ------ basename : str - Name of the directory to analyse + Name used to prefix output files / title the plots. base_outputpath : str Base path where to write data """ # Initiate output directory gen_outputdir(outputpath) - # Get basename of analysis - basename = os.path.basename(os.path.dirname(benchmark_directory)) + # Get basename of analysis, or use the user-provided label if given + basename = label if label else os.path.basename(os.path.dirname(benchmark_directory)) # Generate baseoutput path base_outputpath = f'{outputpath}{basename}' return basename, base_outputpath @@ -1734,7 +2036,7 @@ def _vprint(silence: bool, msg: str) -> None: ################# def main( benchmark_directory: str, - outputpath: str = 'analysis/', + outputpath: str = 'Analysis/', scenarios: Optional[list[str]] = None, metric: str = "irmsd", silent: bool = False, @@ -1743,6 +2045,8 @@ def main( no_violinplots: bool = False, no_melquiplots: bool = False, read_from_archive: bool = False, + label: Optional[str] = None, + per_scenario_plots: bool = False, ) -> None: """Run the analysis procedure. @@ -1771,6 +2075,14 @@ def main( read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False + label : Optional[str], optional + Custom name to use for output filenames and plot titles instead of + the benchmark directory's basename, by default None + per_scenario_plots : bool, optional + In addition to the single combined comparison bar/violin plot + (which already shows every scenario as its own row), also + generate a standalone bar/violin plot for each individual + scenario, by default False """ # Pre-setting the STDOUT print function vprint = partial(_vprint, silent) @@ -1779,6 +2091,7 @@ def main( basename, base_outputpath = set_output_path( benchmark_directory, outputpath, + label=label, ) vprint(f"Setting the output directory path: `{outputpath}`") @@ -1792,6 +2105,16 @@ def main( read_from_archive=read_from_archive, ) + # Auto-detect a human name for each caprieval stage from the module it + # evaluated (overrides the hardcoded CAPRIEVAL_STEPS for plot titles). + global DETECTED_STAGE_NAMES + DETECTED_STAGE_NAMES = detect_stage_modules(dtmapper, read_from_archive) + if DETECTED_STAGE_NAMES: + vprint( + "- Auto-named caprieval stages: " + + ", ".join(f"{k}->{v}" for k, v in sorted(DETECTED_STAGE_NAMES.items())) + ) + # Initiate all scenarios perforamces mapper vprint( f"- Loading data from `{benchmark_directory}` " @@ -1816,13 +2139,18 @@ def main( # Write data as json write_json(all_scenar_perfs, f"{base_outputpath}_performances.json") + # Plot title: when the user passes -l/--label, use it verbatim as the + # title. Otherwise fall back to a plain generic title. (Note: `basename` + # is still used for output filenames, so it stays as-is above.) + plot_title = label if label else "Benchmark output" + # Draw general graph if not no_capriplots: vprint("- Generating Bar plots") gen_full_comparison_barplots( all_scenar_perfs, basepath=base_outputpath, - title=basename, + title=plot_title, progress=not silent, no_percentage=no_percentage, ) @@ -1832,27 +2160,64 @@ def main( gen_full_comparison_violins( all_scenar_perfs, basepath=base_outputpath, - title=basename, + title=plot_title, metric=metric, progress=not silent, ) if not no_melquiplots: - vprint("- Generating Melqui plots") - gen_full_comparison_melquiplots( + vprint("- Generating Melqui plots") + gen_full_comparison_melquiplots( all_scenar_melquis, perf_dtype=metric, basepath=base_outputpath, progress=not silent, ) + # Draw per-scenario graphs, so each scenario also gets its own + # standalone bar/violin plot, not just the combined comparison figure. + # Skipped when there's only a single scenario overall, since in that + # case the combined plot and the per-scenario plot would be identical + # (redundant duplicate file). + if per_scenario_plots and len(all_scenar_perfs) > 1: + for scenar_name, scenar_perfs in all_scenar_perfs.items(): + scenar_outputpath = f"{base_outputpath}_{scenar_name}" + scenar_title = f"{plot_title} - {scenar_name}" + if not no_capriplots: + vprint(f"- Generating Bar plot for scenario `{scenar_name}`") + gen_full_comparison_barplots( + {scenar_name: scenar_perfs}, + basepath=scenar_outputpath, + title=scenar_title, + progress=not silent, + no_percentage=no_percentage, + ) + if not no_violinplots: + vprint(f"- Generating Violin plot for scenario `{scenar_name}`") + gen_full_comparison_violins( + {scenar_name: scenar_perfs}, + basepath=scenar_outputpath, + title=scenar_title, + metric=metric, + progress=not silent, + ) + ################################### # COMMAND LINE ARGUMENTS HANDLERS # ################################### def must_end_with_slash(dirpath: str) -> str: """Make sure a directory paths is terminating by '/'. - + + NOTE: only appends the slash if `dirpath` already exists as a + directory. Suitable for `benchmark_directory`, which must already + exist. NOT suitable for an output path that hasn't been created yet + -- use `ensure_trailing_slash` for that instead, otherwise the missing + slash causes output filenames to get mashed together with whatever + prefix follows (e.g. `AnalysisHADDOCK24_ab_initio_...` instead of + `Analysis/HADDOCK24_ab_initio_...`), and can leave output files + sitting loose inside the benchmark directory itself. + Parameters ---------- dirpath : str @@ -1863,6 +2228,23 @@ def must_end_with_slash(dirpath: str) -> str: return dirpath +def ensure_trailing_slash(dirpath: str) -> str: + """Unconditionally make sure a path ends with '/'. + + Unlike `must_end_with_slash`, this does not require the path to + already exist -- appropriate for output directories, which are + created later via `gen_outputdir`. + + Parameters + ---------- + dirpath : str + Path to a directory (existing or not). + """ + if not dirpath.endswith('/'): + dirpath += '/' + return dirpath + + def parse_cmd_line() -> argparse.Namespace: """Parse command line argument. @@ -1876,22 +2258,64 @@ def parse_cmd_line() -> argparse.Namespace: args = _get_cmd_line_args() # Convert directory paths args.benchmark_directory = must_end_with_slash(args.benchmark_directory) - args.output_path = must_end_with_slash(args.output_path) + # Default output path (when -o/--output_path isn't given) lives INSIDE + # the benchmark directory itself (e.g. `/Analysis/`), + # not relative to wherever the script happens to be run from. If the + # user explicitly passes -o, that value is used as-is (relative to cwd, + # or absolute). + if args.output_path is None: + args.output_path = f"{args.benchmark_directory}Analysis" + args.output_path = ensure_trailing_slash(args.output_path) + # Type-aware default metric: glycan and small molecule are assessed with + # ilrmsd by default; every other type defaults to irmsd. An explicit -m + # always wins. + if args.metric is None: + args.metric = "ilrmsd" if args.type in ("glycan", "small_molecule") else "irmsd" # Set global variables DPI = args.dpi PERFORMANCES_CLASSES = ALL_PERFORMANCES_CLASSES[args.type] + # `ilrmsd` boundaries are only defined for protein-glycan. Reject any other + # combination early with a clear message rather than crashing on a KeyError. + if args.metric not in PERFORMANCES_CLASSES: + sys.exit( + f'INPUT ERROR: metric "{args.metric}" is not defined for type ' + f'"{args.type}". Available metrics for "{args.type}": ' + f'{", ".join(PERFORMANCES_CLASSES.keys())}. ' + f'(Note: `ilrmsd` is only for `-t glycan` and `-t small_molecule`.)' + ) # Check that directory exist if not os.path.exists(args.benchmark_directory): sys.exit( f'INPUT ERROR: path to directory to analyse ' f'"{args.benchmark_directory}" not found!' ) + # Echo the resolved configuration so every run states exactly which + # metric it is classifying by (and which are available for this type). + if not args.quiet: + default_note = " (default for this type)" if args.metric == ( + "ilrmsd" if args.type in ("glycan", "small_molecule") else "irmsd" + ) else " (explicitly set via -m)" + print( + f"[config] type = {args.type} | metric = {args.metric}" + f"{default_note} | metrics available for {args.type}: " + f"{', '.join(PERFORMANCES_CLASSES.keys())}" + ) return args def _get_cmd_line_args() -> argparse.Namespace: """Define and parse the command line arguments.""" - parser = argparse.ArgumentParser() + parser = argparse.ArgumentParser( + description=( + "Analyse a HADDOCK3 + haddock-runner benchmark. Expects " + "///run1/ layout " + "(one or more scenarios). Generates a combined bar plot, " + "violin plot, per-scenario melquiplots (zipped), and a JSON " + "summary, written by default to Analysis/ inside the " + "benchmark directory." + ), + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) parser.add_argument( "benchmark_directory", help=( @@ -1903,18 +2327,37 @@ def _get_cmd_line_args() -> argparse.Namespace: parser.add_argument( "-o", "--output_path", - help="Directory where to write output files.", + help=( + "Directory where to write output files. Defaults to " + "`Analysis/` inside the benchmark directory itself." + ), required=False, - default="analysis", + default=None, + type=str, + ) + parser.add_argument( + '-l', + '--label', + help=( + "Custom name to use for output filenames and plot titles, " + "instead of the benchmark directory's basename (useful when " + "the directory is named something generic, e.g. `test-analysis`)." + ), + required=False, + default=None, type=str, ) parser.add_argument( "-m", "--metric", - help="Performance metric to track.", + help=( + "Performance metric to track. If not given, defaults to `ilrmsd` " + "for `-t glycan` and `-t small_molecule`, and `irmsd` otherwise. " + "`ilrmsd` is only defined for glycan and small molecule." + ), required=False, - default="irmsd", - choices=("irmsd", "dockq", ), + default=None, + choices=("irmsd", "dockq", "ilrmsd", ), type=str, ) parser.add_argument( @@ -1932,7 +2375,7 @@ def _get_cmd_line_args() -> argparse.Namespace: help="Type of analysis to be conducted.", required=False, default="protein", - choices=("protein", "peptide", "glycan", ), + choices=("protein", "peptide", "glycan", "dna", "small_molecule", ), type=str, ) parser.add_argument( @@ -1961,6 +2404,17 @@ def _get_cmd_line_args() -> argparse.Namespace: action="store_true", default=False, ) + parser.add_argument( + '--per-scenario-plots', + help=( + "In addition to the single combined comparison plot " + "(all scenarios as rows in one bar/violin plot), also " + "generate a separate standalone bar/violin plot for each " + "individual scenario. Off by default." + ), + action="store_true", + default=False, + ) parser.add_argument( '-n', '--no-percentage', @@ -2028,6 +2482,8 @@ def maincli() -> None: no_violinplots=args.no_violinplots, no_melquiplots=args.no_melquiplots, read_from_archive=args.from_archive, + label=args.label, + per_scenario_plots=args.per_scenario_plots, ) From 0dd205ec3f6ccfeee2c2c5a3f05d80da3eb168cf Mon Sep 17 00:00:00 2001 From: comp-era Date: Thu, 30 Jul 2026 14:36:07 +0300 Subject: [PATCH 2/8] Plot tweaks: DPI 200, optional violin/melqui, cleaner legend, manual stage names --- analysis/analysebenchmarkresults.py | 348 +++++----------------------- 1 file changed, 52 insertions(+), 296 deletions(-) diff --git a/analysis/analysebenchmarkresults.py b/analysis/analysebenchmarkresults.py index e60c5de..989b2a5 100644 --- a/analysis/analysebenchmarkresults.py +++ b/analysis/analysebenchmarkresults.py @@ -1,16 +1,12 @@ """Analysis script dedicated to the HADDOCK3 + haddock-runner benchmarks. - Generates multiple-plots analysing different scenarios performances. - Barplots: Standard best performing model among top X from all targets. - Melquiplots: Per-target complex performances among top 200. - Violinplots: Performance distribution among top X from all targets. - Please modify the Global variable: CAPRIEVAL_STEPS to suite your needs. It is used to generate nice title to the caprieval steps. - Usage: >python3 AnalyseBenchmarkResults.py - Expected input structure ------------------------- / @@ -19,12 +15,14 @@ run1/ _caprieval/ capri_ss.tsv - (scenario is the OUTER directory, PDBid the INNER one -- this matches the directory layout currently produced by haddock-runner. The previous `//run1/` layout is no longer supported.) -""" +NOTE: bar plots are generated by default. Violin and melqui plots are OFF by +default -- enable them explicitly with --violinplots and/or --melquiplots. +Default figure DPI is 200. +""" import argparse import glob import json @@ -33,11 +31,9 @@ import sys import tarfile import zipfile - from functools import partial from pathlib import Path from typing import Callable, Optional, Union - # Try to load external libraries try: import numpy as np @@ -51,8 +47,6 @@ "Please make sure they are accessible with current environment.\n" "(e.g.: >pip install numpy matplotlib)\n" ) - - __version__ = "2.0.0" # July 2026 __author__ = ", ".join(( "BonvinLab", @@ -67,8 +61,6 @@ "Victor G.P. Reys", "Shantanu Khatri", ) - - #################### # GLOBAL VARIABLES # #################### @@ -81,12 +73,10 @@ # handle everything (undetected stages then fall back to `_caprieval`). # To force a custom label for a specific stage, add e.g. '02': 'rigidbody'. CAPRIEVAL_STEPS: dict[str, str] = {} - # Set threshold of top X structures to take into account TOP_X_THRESHOLDS = (1, 5, 10, 20, 50, 100, 200, 500, 1000, ) # Set number of entries to display in melquiplot MELQUIPLOT_NB_ENTRIES = 200 - # CAPRI performances classes # NOTE: for each class, we define the lower and upper limit ALL_PERFORMANCES_CLASSES = { @@ -187,17 +177,15 @@ }, } PERFORMANCES_CLASSES = ALL_PERFORMANCES_CLASSES["protein"] - # Add color mapper COLORS_MAPPER = { "High": "darkgreen", "Medium": "lightgreen", "Acceptable": "lightblue", "Near-acceptable": "gainsboro", - "Low": "lightcoral", # was "white" (invisible on white background) - "Missing": "dimgrey", + "Low": "white", + "Missing": "white", } - # Performance order PERF_ORDER = ( "High", @@ -207,10 +195,14 @@ "Low", "Missing", ) +# Performance classes shown in the plot legends. `Low` and `Missing` are +# intentionally omitted (both are drawn white / effectively invisible), so +# they don't clutter the legend. Bars are still stacked using PERF_ORDER. +LEGEND_PERF_ORDER = tuple( + p for p in PERF_ORDER if p not in ("Low", "Missing") + ) # DPI of the generated figures -DPI = 400 - - +DPI = 200 #################### # DEFINE FUNCTIONS # #################### @@ -227,7 +219,6 @@ def gen_graph( percentage: bool = True, ) -> None: """Plot a barplot on the provided axis `ax`. - Parameters ---------- ax : `matplotlib.pyplot.Axes` @@ -289,7 +280,6 @@ def gen_graph( tops = [f'Top{v}' for v in top] # initialize first bottom values with 0s bottom = np.zeros(len(high)) - # Loop over performances for label in PERF_ORDER: # point performances data @@ -327,15 +317,11 @@ def gen_graph( ax.set_xticklabels(ax.get_xticklabels(), rotation=30, ha='right') ylabel = "Nb. entries" if not percentage else "% success rate" ax.set_ylabel(ylabel) - - def clear_plt() -> None: """Clear all previous instances/data generated by matplotlib.""" plt.gca() plt.cla() plt.clf() - - def gen_violin( ax: plt.Axes, perf_data: list, @@ -343,10 +329,8 @@ def gen_violin( metric: str = "", ) -> None: """Draw a violinplot on the axis. - inspired from: https://matplotlib.org/stable/gallery/statistics/customized_violin.html - Parameters ---------- ax : `matplotlib.pyplot.Axes` @@ -363,11 +347,9 @@ def gen_violin( showmedians=False, showextrema=False, ) - # Get color ramp nbcolors = 10 if len(perf_data) <= 10 else 20 colorramp = mpl.colormaps[f'tab{nbcolors}'] - # Modify colors for vi, pc in enumerate(parts['bodies']): pc.set_facecolor(colorramp((vi + 0.5) / nbcolors)) @@ -388,7 +370,6 @@ def gen_violin( ax.vlines(inds, quartile1, quartile3, color='k', linestyle='-', lw=5) ax.vlines(inds, whiskers_min, whiskers_max, color='k', linestyle='-', lw=1) ax.set_ylabel(metric) - # Set labels if labels: ax.set_xticks(list(range(1, len(labels) + 1))) @@ -398,18 +379,14 @@ def gen_violin( ha='right', rotation_mode='anchor', ) - - def adjacent_values( vals: list[float], q1: float, q3: float, ) -> tuple[float, float]: """Find adjacent values. - Inspired from: https://matplotlib.org/stable/gallery/statistics/customized_violin.html - Parameters ---------- vals : list @@ -418,7 +395,6 @@ def adjacent_values( Value of the first quartil q3 : float Value of the 3rd quartil - Return ------ lower_adjacent_value : float @@ -433,52 +409,42 @@ def adjacent_values( lower_adjacent_value = q1 - (q3 - q1) * 1.5 lower_adjacent_value = np.clip(lower_adjacent_value, vals[0], q1) return lower_adjacent_value, upper_adjacent_value - - # Auto-detected {stage_index: module_name}, filled at runtime by # detect_stage_modules(). Takes precedence over the hardcoded CAPRIEVAL_STEPS. DETECTED_STAGE_NAMES: dict[str, str] = {} - - def stage_name(cname: str) -> str: """Return a human-readable name for a caprieval stage. - Preference order: - 1. auto-detected module name (the module this caprieval evaluated), - 2. the hardcoded CAPRIEVAL_STEPS mapping, + 1. an explicit CAPRIEVAL_STEPS entry (manual override wins), + 2. the auto-detected module name (the module this caprieval evaluated), 3. a `_caprieval` fallback. - Parameters ---------- cname : str Index of the caprieval stage - Returns ------- name : str Name of the stage """ + # A manually-specified name always wins, so CAPRIEVAL_STEPS can be used + # to override the auto-detected module name for any stage. + if cname in CAPRIEVAL_STEPS: + return CAPRIEVAL_STEPS[cname] if cname in DETECTED_STAGE_NAMES: return DETECTED_STAGE_NAMES[cname] - try: - return CAPRIEVAL_STEPS[cname] - except KeyError: - return f"{cname}_caprieval" - - + return f"{cname}_caprieval" def detect_stage_modules( dtmapper: dict, read_from_archive: bool = False, ) -> dict[str, str]: """Auto-name caprieval stages by the module they evaluated. - Every caprieval `capri_ss.tsv` lists model paths like `../../10_seletopclusts/cluster_1_model_1.pdb`; the parent directory (`10_seletopclusts`) is the module being evaluated. This reads one file per stage (from the first scenario/target) and strips the numeric prefix to get the module name (`seletopclusts`). Duplicate module names across stages are disambiguated with their stage index, e.g. `flexref (07)`. - Falls back silently (empty entry) for any stage it cannot read, so the hardcoded CAPRIEVAL_STEPS / `_caprieval` naming still applies. """ @@ -507,8 +473,6 @@ def detect_stage_modules( if counts[mapping[stage]] > 1: mapping[stage] = f"{mapping[stage]} ({stage})" return mapping - - def gen_full_comparison_violins( scenars_perfs: dict, basepath: str = "./", @@ -517,7 +481,6 @@ def gen_full_comparison_violins( progress: bool = True, ) -> None: """Combine all scenarios caprieval within same plot. - Parameters ---------- scenars_perfs : dict @@ -542,7 +505,6 @@ def gen_full_comparison_violins( nb_thresh = len(tops_order) # Compute total number of plots total_plots = nb_thresh * nb_steps - # Initate figures / axis fig, axes = plt.subplots( figsize=((nb_steps * 4) + 1, (nb_thresh * 3) + 1), @@ -552,7 +514,6 @@ def gen_full_comparison_violins( sharex=True, squeeze=0, # squeeze=0 allows to always return a 2d array ) - processed = 0 # Loop over thresholds for ri, topx in enumerate(tops_order): @@ -579,7 +540,6 @@ def gen_full_comparison_violins( labels, metric=metric, ) - # Add columns titles pad = 5 for ax, cname in zip(axes[0], steps_order): @@ -606,10 +566,8 @@ def gen_full_comparison_violins( ha='right', va='center', ) - # Add figure title fig.suptitle(title, fontsize=16) - # Get color ramp nbcolors = 10 if nb_scenar <= 10 else 20 colorramp = mpl.colormaps[f'tab{nbcolors}'] @@ -628,9 +586,7 @@ def gen_full_comparison_violins( ncols=4, title="Scenarios", ) - plt.gca().set_ylim(bottom=0) - # Adjust border to let annotations fit inside graph. The left margin # needs to hold both the y-axis label and the row-title annotations # (drawn at a small fixed point-offset from it). That fixed offset @@ -640,7 +596,6 @@ def gen_full_comparison_violins( fig_width_in = fig.get_size_inches()[0] left_margin = max(0.08, min(0.35, 1.3 / fig_width_in)) fig.subplots_adjust(left=left_margin, top=0.95, bottom=0.12, right=0.98) - # save figure plt.savefig( f"{basepath}_violins.png", @@ -650,8 +605,6 @@ def gen_full_comparison_violins( pad_inches=0.3, ) return - - def gen_full_comparison_melquiplots( scenars_perfs: dict, perf_dtype: str = "irmsd", @@ -659,7 +612,6 @@ def gen_full_comparison_melquiplots( progress: bool = True, ) -> str: """Generate multiple melquiplots for each scenario. - Parameters ---------- scenars_perfs : dict @@ -670,7 +622,6 @@ def gen_full_comparison_melquiplots( Prefix to give to archive, by default "benchmark_melquis" basepath : str, optional Where to write the files, by default "./" - Returns ------- archive_path : str @@ -678,12 +629,10 @@ def gen_full_comparison_melquiplots( """ # Clear pervious instances of matplotlib clear_plt() - # Set progression variables all_generated_melquis: list[str] = [] processed = 0 total_plots = len(scenars_perfs) - # Loop over scenarios for scenar_name, scenar_perfs in scenars_perfs.items(): processed += 1 @@ -697,7 +646,6 @@ def gen_full_comparison_melquiplots( basepath=basepath, ) all_generated_melquis.append(scenar_melqui_path) - # Generate archive of melqui plots archive_path = Path(f"{basepath}_melquiplots.zip") path = archive_path.parent @@ -708,15 +656,11 @@ def gen_full_comparison_melquiplots( for figure in all_generated_melquis: zipf.write(Path(figure).name) os.chdir(initdir) - # Remove all original files for generated_melqui in all_generated_melquis: os.remove(generated_melqui) - # Return archive path return archive_path - - def make_scenar_melquiplots( scenar_perfs: dict, perf_dtype: str = "irmsd", @@ -724,7 +668,6 @@ def make_scenar_melquiplots( basepath: str = "./", ) -> str: """Generate multiples melquiplots for each Caprieval steps of a scenario. - Parameters ---------- scenar_perfs : dict @@ -735,7 +678,6 @@ def make_scenar_melquiplots( Title to the figure, by default "melquiplot" basepath : str, optional Where to write the files, by default "./" - Returns ------- figpath : str @@ -754,7 +696,6 @@ def make_scenar_melquiplots( ) if not isinstance(axes, np.ndarray): axes = [axes] - # Loop over stages for si, (stage, stage_perfs) in enumerate(scenar_perfs.items()): # Point axis @@ -780,20 +721,18 @@ def make_scenar_melquiplots( # (`lo < M <= hi`), so both the unit and the inequality direction flip. higher_is_better = get_reverse_bool(perf_dtype) unit = "" if higher_is_better else r"$\AA$" - def _perf_label(perfclass: str) -> str: lo, hi = dtype_perf_classes[perfclass] metric = perf_dtype.upper() if higher_is_better: return rf"{perfclass} | {lo} < {metric} <= {hi}{unit}" return rf"{perfclass} | {lo} <= {metric} < {hi}{unit}" - legend_data = [ ( plt.Rectangle((0, 0), 1, 1, fc=COLORS_MAPPER[perfclass]), _perf_label(perfclass), ) - for perfclass in PERF_ORDER + for perfclass in LEGEND_PERF_ORDER ] legend_proxies, legend_labels = zip(*legend_data) # Add bars legend @@ -801,13 +740,11 @@ def _perf_label(perfclass: str) -> str: legend_proxies, legend_labels, loc='outside lower center', - ncols=len(PERF_ORDER), + ncols=len(LEGEND_PERF_ORDER), title="performance classes", ) - # adjust border to let annotations fit inside graph fig.subplots_adjust(left=0.02, top=0.95, bottom=0.06, right=0.98) - # save figure figpath = f"{basepath}_{title}_melquiplot.png" plt.savefig( @@ -818,15 +755,12 @@ def _perf_label(perfclass: str) -> str: pad_inches=0.3, ) return figpath - - def melquiplot( ax: plt.Axes, pdb_perfs: dict[str, list[str]], width: float = 0.3, ) -> None: """Draw a melquiplot on an sub-figure with provided input data. - Parameters ---------- ax : plt.Axes @@ -858,10 +792,8 @@ def melquiplot( linewidth=0, ) y_coord += 1 - max_stack_y = max(max_stack_y, y_coord) stack_h_labels_pos.append(x_coord + (width / 2)) - # Aesthetics ax.set_xlim((0, len(pdb_labels) * width)) ax.set_ylim((0, max_stack_y)) @@ -869,14 +801,11 @@ def melquiplot( ax.set_xticks(stack_h_labels_pos) ax.set_xticklabels(pdb_labels, rotation=45, fontsize='small') ax.xaxis.set_ticks_position('none') - - def melquiplot_original( ax: plt.Axes, pdb_perfs: dict[str, list[str]], ) -> None: """Draw a melquiplot on an sub-figure with provided input data. - Parameters ---------- ax : plt.Axes @@ -905,10 +834,8 @@ def melquiplot_original( linewidth=0, ) stack_v += 1 - max_stack_v = max(max_stack_v, stack_v) stack_h_labels_pos.append(entry_index - 0.501) - # Aesthetics ax.set_xlim((0.05, len(pdb_labels) + 0.05)) ax.set_ylim((0, max_stack_v)) @@ -916,8 +843,6 @@ def melquiplot_original( ax.set_xticks(stack_h_labels_pos) ax.set_xticklabels(pdb_labels, rotation=45, fontsize='small') ax.xaxis.set_ticks_position('none') - - def gen_full_comparison_barplots( scenars_perfs: dict, basepath: str = "./", @@ -926,7 +851,6 @@ def gen_full_comparison_barplots( no_percentage: bool = False, ) -> None: """Combine all scenarios caprieval within same plot. - Parameters ---------- scenars_perfs : dict @@ -947,13 +871,11 @@ def gen_full_comparison_barplots( # Compute total number of plots total_plots = nb_rows * nb_cols processed = 0 - # Initate figures / axis fig, axes = plt.subplots( figsize=((nb_cols * 4) + 1, (nb_rows * 3) + 4), nrows=nb_rows, ncols=nb_cols, ) - # Loop over rows for ri, rname in enumerate(rows_order): # Point row axe(s) @@ -994,10 +916,8 @@ def gen_full_comparison_barplots( topx, percentage=not no_percentage, ) - # Set padding between two plots pad = 5 # NOTE: nicely hardcoded ! - # Search for first row if len(rows_order) == 1: first_row = axes @@ -1005,7 +925,6 @@ def gen_full_comparison_barplots( first_row = axes[0] if not isinstance(first_row, np.ndarray): first_row = [first_row] - # Add columns titles for ax, cname in zip(first_row, cols_order): # Add column name @@ -1019,7 +938,6 @@ def gen_full_comparison_barplots( ha='center', va='baseline', ) - # Find all first row columns if len(rows_order) == 1: axrows = [axes] @@ -1029,7 +947,6 @@ def gen_full_comparison_barplots( axrows_firstcols = axrows else: axrows_firstcols = [ax[0] for ax in axrows] - # Add rows titles for ax, rname in zip(axrows_firstcols, rows_order): ax.annotate( @@ -1042,10 +959,8 @@ def gen_full_comparison_barplots( ha='right', va='center', ) - # Add figure title fig.suptitle(title, fontsize=16) - # Add bars legend fig.legend( [mpatches.FancyBboxPatch( @@ -1053,13 +968,12 @@ def gen_full_comparison_barplots( boxstyle=mpatches.BoxStyle("Round", pad=0.02), color=COLORS_MAPPER[perfclass], ) - for perfclass in PERF_ORDER], - PERF_ORDER, + for perfclass in LEGEND_PERF_ORDER], + LEGEND_PERF_ORDER, loc='outside lower center', - ncols=len(PERF_ORDER), + ncols=len(LEGEND_PERF_ORDER), title="performance classes", ) - # Adjust border to let annotations fit inside graph. Same reasoning as # in `gen_full_comparison_violins`: scale the left margin so a narrow # figure (e.g. a single caprieval stage -> nb_cols=1) still has enough @@ -1068,7 +982,6 @@ def gen_full_comparison_barplots( fig_width_in = fig.get_size_inches()[0] left_margin = max(0.15, min(0.35, 1.3 / fig_width_in)) fig.subplots_adjust(left=left_margin, top=0.9, bottom=0.15, right=0.98) - # Save figure plt.savefig( f"{basepath}_capribarplots.png", @@ -1078,21 +991,16 @@ def gen_full_comparison_barplots( pad_inches=0.3, ) return - - def get_pdb_entries(basepath: str) -> list: """Retrieve list of immediate subdirectory names. - Generic directory lister -- reused both to list scenario names at the top level of the benchmark directory, and to list PDBid/target names inside each scenario directory. Kept its original name for API stability, even though it's no longer PDBid-specific. - Parameters ---------- basepath : str Path to a directory whose immediate subdirectories should be listed. - Return ------ entries : list @@ -1104,14 +1012,11 @@ def get_pdb_entries(basepath: str) -> list: if pdbid_path.is_dir() ] return entries - - def _find_scenario_data_dir( candidate_path: str, _max_depth: int = 3, ) -> Optional[str]: """Find the directory that actually holds `/run1/` entries. - Checks whether `candidate_path` directly contains one or more ``/run1/`` subdirectories. If not, and `candidate_path` has exactly one subdirectory, recurses into that single subdirectory (up to @@ -1121,14 +1026,12 @@ def _find_scenario_data_dir( ``//run1/``. Only drills down when there is a single, unambiguous subdirectory to descend into, so it won't misfire on directories that hold unrelated clutter (e.g. output files). - Parameters ---------- candidate_path : str Directory to inspect (should end with '/'). _max_depth : int Maximum number of wrapper levels to look through. - Return ------ resolved : Optional[str] @@ -1154,11 +1057,8 @@ def _find_scenario_data_dir( _max_depth=_max_depth - 1, ) return None - - def _looks_like_scenario_dir(scenario_path: str) -> bool: """Heuristic check for whether a directory is a real scenario dir. - A genuine scenario directory (as produced by haddock-runner) contains at least one ``/run1/`` subdirectory, possibly behind a single extra wrapper level (see `_find_scenario_data_dir`). This is used to @@ -1166,12 +1066,10 @@ def _looks_like_scenario_dir(scenario_path: str) -> bool: scenario directories -- most commonly a previous analysis output directory (e.g. `Analysis/`) left inside the benchmark directory, or stray/hidden directories -- so they aren't mistaken for scenarios. - Parameters ---------- scenario_path : str Candidate scenario directory path (should end with '/'). - Return ------ bool @@ -1179,31 +1077,22 @@ def _looks_like_scenario_dir(scenario_path: str) -> bool: directly or behind a single unambiguous wrapper level. """ return _find_scenario_data_dir(scenario_path) is not None - - def _resolve_scenario_dir(basepath: str, scenario: str) -> str: """Resolve a scenario name to its actual data directory on disk. - Normally this is simply ``{basepath}{scenario}/``. However, some benchmark directories end up with an extra wrapper level in between (e.g. produced by manual reorganisation), such as:: - ////run1/ - instead of the expected:: - ///run1/ - This transparently drills down through such wrapper levels via `_find_scenario_data_dir`. - Parameters ---------- basepath : str Root directory containing scenario directories. scenario : str Scenario name. - Return ------ resolved : str @@ -1215,22 +1104,17 @@ def _resolve_scenario_dir(basepath: str, scenario: str) -> str: candidate = f"{basepath}{scenario}/" resolved = _find_scenario_data_dir(candidate) return resolved if resolved is not None else candidate - - def get_scenario_entries(basepath: str) -> list: """Retrieve list of scenario directory names at the top of `basepath`. - Like `get_pdb_entries`, but filters out any top-level directory that doesn't actually look like a scenario directory (see `_looks_like_scenario_dir`). This protects against stray directories such as a leftover `Analysis/` output folder, hidden directories, etc. being mistaken for scenarios. - Parameters ---------- basepath : str Path to the benchmark directory to analyse. - Return ------ scenarios : list @@ -1242,18 +1126,13 @@ def get_scenario_entries(basepath: str) -> list: if _looks_like_scenario_dir(f'{basepath}{name}/') ] return scenarios - - def scenario_name_2_threshold(scenar_name: str) -> float: """Gather threshold value from scenario name. - NOTE: function used for CPORT-ARCTIC3D-BM5 benchmark - Parameters ---------- scenar_name : str Name of a scenario. - Return ------ threshold : float @@ -1266,14 +1145,11 @@ def scenario_name_2_threshold(scenar_name: str) -> float: # Cast it to float threshold = float(dot_thresh) return threshold - - def get_caprieval_stages( basepath: str, read_from_archive: bool = False, ) -> list[str]: """Retrieve all caprieval stages inside a haddock3 run. - Parameters ---------- basepath : str @@ -1281,7 +1157,6 @@ def get_caprieval_stages( read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False - Return ------ caprieval_stages : list @@ -1302,17 +1177,12 @@ def get_caprieval_stages( for caprip in caprieval_paths ] return caprieval_stages - - def hd3_module_2_stage(indexed_modulename: str) -> str: """Split a haddock3 module name and retrieve it id. - Parameters ---------- indexed_modulename : str Directory name of an indexed haddock3 module name. - - Return ------ stage : str @@ -1320,11 +1190,8 @@ def hd3_module_2_stage(indexed_modulename: str) -> str: """ stage = indexed_modulename.split('_')[0] return stage - - def _make_run_path(scenario_basepath: str, pdbid: str) -> str: """Return the haddock3 run directory path. - Parameters ---------- scenario_basepath : str @@ -1335,11 +1202,8 @@ def _make_run_path(scenario_basepath: str, pdbid: str) -> str: Target/PDBid name. """ return f"{scenario_basepath}{pdbid}/run1/" - - def _make_archive_path(scenario_basepath: str, pdbid: str) -> str: """Return the haddock3 run analysis archive path. - Parameters ---------- scenario_basepath : str @@ -1348,25 +1212,19 @@ def _make_archive_path(scenario_basepath: str, pdbid: str) -> str: Target/PDBid name. """ return f"{scenario_basepath}{pdbid}/run1_analysis.tgz" - - def map_data( basepath: str, subset_scenarios: Optional[list[str]] = None, read_from_archive: bool = False, ) -> dict[str, dict[str, dict[str, dict[str, Union[str, list[str, str]]]]]]: """Map data in one analysis dict to accessit easily. - Expects the benchmark directory to be laid out as produced by haddock-runner:: - ///run1/_caprieval/capri_ss.tsv - i.e. scenario is the outer directory and PDBid the inner one. This holds whether `basepath` contains a single scenario directory or several -- there is no special-casing, the code below just handles however many scenario directories are found. - Parameters ---------- basepath : str @@ -1377,7 +1235,6 @@ def map_data( read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False - Return ------ dtmap : dict[str, dict[str, dict[str, dict[str, str]]]] @@ -1395,7 +1252,6 @@ def map_data( """ # Initiate mapper dtmap: dict[str, dict[str, dict[str, dict[str, Union[str, list[str, str]]]]]] = {} - # `basepath` is always the root containing one or more scenario # directories -- this naturally handles both a single scenario and # multiple scenarios without any special-casing, since the loops below @@ -1412,11 +1268,9 @@ def map_data( "Expected layout is `///run1/` -- " "this applies whether there is one scenario or several." ) - def scenario_dir(scenario: str) -> str: """Resolve a scenario name to its actual directory on disk.""" return _resolve_scenario_dir(basepath, scenario) - # pdbids are gathered per scenario, deduplicated across all scenarios. pdbids_set: set[str] = set() for scenario in all_scenarios: @@ -1426,7 +1280,6 @@ def scenario_dir(scenario: str) -> str: f"[ERROR] no target/PDBid directories found under any scenario in " f"`{basepath}`." ) - # Make sure all pdbs have all scenarios for scenario in all_scenarios: for pdbid in pdbids: @@ -1446,7 +1299,6 @@ def scenario_dir(scenario: str) -> str: f"[ERROR] could not find scenario `{scenario}` " f"archive for entry `{pdbid}` at: {analysis_archive}" ) - all_caprieval_stages = [] # Add scenario data to data maper for scenario in all_scenarios: @@ -1467,7 +1319,6 @@ def scenario_dir(scenario: str) -> str: assert all_caprieval_stages, ( f"[ERROR] no caprieval stages found anywhere under `{basepath}`." ) - # Make sure all stages are computed for all pdb in all scenarios... for scenario in all_scenarios: for pdbid in pdbids: @@ -1501,7 +1352,6 @@ def scenario_dir(scenario: str) -> str: ) if read_from_archive: tararchive.close() - # Gather all data for scenario in all_scenarios: dtmap[scenario] = {} @@ -1518,10 +1368,7 @@ def scenario_dir(scenario: str) -> str: ] # Hold datapath dtmap[scenario][stage][pdbid] = caprieval_tsv_path - return dtmap - - def analyse_scenario( scenario_dt: dict, _entries_thresholds: list, @@ -1530,7 +1377,6 @@ def analyse_scenario( read_from_archive: bool = False, ) -> tuple[dict[str, dict], dict[str, dict]]: """Process the analysis of a scenario. - Parameters ---------- scenario_dt : dict @@ -1546,7 +1392,6 @@ def analyse_scenario( read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False - Return ------ scenario_stages_perfs : dict @@ -1559,7 +1404,6 @@ def analyse_scenario( entries_thresholds.append(MELQUIPLOT_NB_ENTRIES) else: entries_thresholds = _entries_thresholds - # Initiate all stages performances holder scenario_stages_perfs = {} scenario_stages_pdb_perfs = {} @@ -1585,7 +1429,6 @@ def analyse_scenario( topx: [perf_to_class(v, dtype=perf_dtype) for v in perfs] for topx, perfs in top_x_perfs.items() } - # Summerize all pdb best performances stage_best_perfs = { n: [pdb_perf[n] for pdb_perf in pdb_best_perfs.values()] @@ -1596,30 +1439,24 @@ def analyse_scenario( n: perfs_to_classes(perfs, dtype=perf_dtype) for n, perfs in stage_best_perfs.items() } - # Hold stage performances scenario_stages_perfs[stage] = { "values": stage_best_perfs, "classes": stage_class_perfs, } scenario_stages_pdb_perfs[stage] = pdb_perfs - return scenario_stages_perfs, scenario_stages_pdb_perfs - - def perfs_to_classes( perfs: list, dtype: str = "irmsd", ) -> dict: """Convert performances values into classes. - Parameters ---------- perfs : list List of performance values. dtype : str Type of the performance. - Return ------ perfs_classes : dict @@ -1637,18 +1474,14 @@ def perfs_to_classes( for cls in PERFORMANCES_CLASSES[dtype].keys() } return perfs_classes - - def perf_to_class(value: float, dtype: str = "irmsd") -> str: """Cast a performance value into a performance class. - Parameters ---------- value : float Performance value. dtype : str Name of the data type. - Return ------ cls : str @@ -1663,16 +1496,12 @@ def perf_to_class(value: float, dtype: str = "irmsd") -> str: return cls # In case it is not found, return "Missing" return PERF_ORDER[-1] - - def dtype_to_boundary_function(dtype: str) -> Callable: """Get function to make a comparison between two boundaries. - Parameters ---------- dtype : str Name of the data type. - Return ------ _function : function @@ -1684,8 +1513,6 @@ def _include_lower(lowb: float, highb: float, value: float) -> bool: return lowb <= value < highb _function = _include_higher if get_reverse_bool(dtype) else _include_lower return _function - - def analyse_caprieval_performances( capireval_fpath: Union[str, list[str, str]], entries_thresholds: list, @@ -1694,7 +1521,6 @@ def analyse_caprieval_performances( read_from_archive: bool = False, ) -> tuple[dict[int, float], dict[int, list[float]]]: """Analyse a CAPRIeval step performances. - Parameters ---------- capireval_fpath : Union[str, list[str, str]] @@ -1708,7 +1534,6 @@ def analyse_caprieval_performances( read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False - Return ------ best_perfs_h : dict @@ -1722,7 +1547,6 @@ def analyse_caprieval_performances( key=lambda k: caprieval_data[k][sort_dtype], reverse=get_reverse_bool(sort_dtype), ) - # Get performances taking into account increasing number of entries best_perfs_h = { nb_entries: best_perfs( @@ -1741,11 +1565,8 @@ def analyse_caprieval_performances( ] for nb_entries in entries_thresholds } - # Return performances return best_perfs_h, top_x_perfs - - def best_perfs( caprieval_data: dict, order: list, @@ -1754,7 +1575,6 @@ def best_perfs( tolerance: float = 6, ) -> float: """Obtain best performing value. - Parameters ---------- caprieval_data : dict @@ -1768,7 +1588,6 @@ def best_perfs( tolerance : float Percentage of tolerated difference between number of accessible entries `len(order)` and number of queried `n` entries. - Return ------ best_perfomance : float @@ -1777,7 +1596,6 @@ def best_perfs( # Check if length of order << n # if (len(order) + (len(order) * tolerance / 100)) < n: # return -1 - # Obtain list of all performances all_perfs = [caprieval_data[entry][dtype] for entry in order[:n]] # Get best function @@ -1785,16 +1603,12 @@ def best_perfs( # Get best performing value best_perfomance = best_function(all_perfs) return best_perfomance - - def dtype_to_best_function(dtype: str) -> Callable: """Get best function. - Parameters ---------- dttype : str Name of the data type. - Return ------ _function : function @@ -1802,16 +1616,12 @@ def dtype_to_best_function(dtype: str) -> Callable: """ _function = max if get_reverse_bool(dtype) else min return _function - - def get_reverse_bool(dt_type: str) -> bool: """Check sorting order. - Parameters ---------- dt_type : str Name of the data type. - Return ------ _reversed : bool @@ -1819,15 +1629,12 @@ def get_reverse_bool(dt_type: str) -> bool: """ _reversed = True if dt_type in ('fnat', 'dockq', ) else False return _reversed - - def load_caprieval_data( tsvpath: Union[str, list[str, str]], read_from_archive: bool, sep: str = "\t", ) -> dict[str, dict[str, float]]: """Load caprieval data as dict. - Parameters ---------- tsvpath : str @@ -1837,7 +1644,6 @@ def load_caprieval_data( analysis archive instead of the run directory, by default False sep : str String used to separate data in the file. - Return ------ data : dict[str, dict[str, float]] @@ -1859,7 +1665,6 @@ def load_caprieval_data( file_content = tar_tsv.read().decode("utf-8") # Close archive tarin.close() - # Loop over file lines for i, _ in enumerate(file_content.split("\n")): # Split the line @@ -1889,11 +1694,8 @@ def load_caprieval_data( # Hold this guy data[complex_name] = complex_dt return data - - def write_json(data: dict, path: str) -> None: """Write a json file containing data. - Parameters ---------- data : dict / list @@ -1903,16 +1705,12 @@ def write_json(data: dict, path: str) -> None: """ with open(path, 'w', encoding='utf-8') as f: json.dump(data, f, indent=4) - - def load_json(path: str) -> Union[dict, list]: """Load a json file. - Parameters ---------- path : str Path to the .json file to load. - Return ------ data : dict or list @@ -1921,8 +1719,6 @@ def load_json(path: str) -> Union[dict, list]: with open(path, "r") as fin: data = json.load(fin) return data - - def get_data_mapper( basepath: str, overwrite: bool = False, @@ -1931,7 +1727,6 @@ def get_data_mapper( read_from_archive: bool = False, ) -> dict: """Retrieve or generate data mapper. - Parameters ---------- basepath : str @@ -1946,7 +1741,6 @@ def get_data_mapper( read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False - Return ------ dtmapper : dict @@ -1958,7 +1752,6 @@ def get_data_mapper( if not overwrite and os.path.exists(mapper_fpath): dtmapper = load_json(mapper_fpath) return dtmapper - # Gather data mapper dtmapper = map_data( basepath, @@ -1969,11 +1762,8 @@ def get_data_mapper( write_json(dtmapper, mapper_fpath) # Return data return dtmapper - - def gen_outputdir(outputdir: str) -> None: """Make sure output directory is created. - Parameters ---------- outputdir : str @@ -1981,15 +1771,12 @@ def gen_outputdir(outputdir: str) -> None: """ if not os.path.exists(outputdir): os.mkdir(outputdir) - - def set_output_path( benchmark_directory: str, outputpath: str, label: Optional[str] = None, ) -> tuple[str, str]: """Initiate output paths and directory. - Parameters ---------- benchmark_directory : str @@ -2000,7 +1787,6 @@ def set_output_path( Custom name to use for output filenames and plot titles instead of the benchmark directory's basename (e.g. useful when the directory is named something generic like `test-analysis`). - Return ------ basename : str @@ -2015,11 +1801,8 @@ def set_output_path( # Generate baseoutput path base_outputpath = f'{outputpath}{basename}' return basename, base_outputpath - - def _vprint(silence: bool, msg: str) -> None: """Print on screen - Parameters ---------- msg : str @@ -2029,8 +1812,6 @@ def _vprint(silence: bool, msg: str) -> None: """ if not silence: print(msg) - - ################# # Main function # ################# @@ -2042,14 +1823,13 @@ def main( silent: bool = False, no_percentage: bool = False, no_capriplots: bool = False, - no_violinplots: bool = False, - no_melquiplots: bool = False, + violinplots: bool = False, + melquiplots: bool = False, read_from_archive: bool = False, label: Optional[str] = None, per_scenario_plots: bool = False, ) -> None: """Run the analysis procedure. - Parameters ---------- benchmark_directory : str @@ -2068,10 +1848,10 @@ def main( no_capriplots : bool, optional Do not generate the usual CAPRI model quality bar plots, by default False - no_violinplots : bool, optional - Do not generate the violin plots, by default False - no_melquiplots : bool, optional - Do not generate the melquiplots, by default False + violinplots : bool, optional + Generate the violin plots. OFF by default. + melquiplots : bool, optional + Generate the melquiplots. OFF by default. read_from_archive : bool, optional When set to True, performs the search of capri_ss.tsv files from the analysis archive instead of the run directory, by default False @@ -2086,7 +1866,6 @@ def main( """ # Pre-setting the STDOUT print function vprint = partial(_vprint, silent) - # Initiate output paths basename, base_outputpath = set_output_path( benchmark_directory, @@ -2094,7 +1873,6 @@ def main( label=label, ) vprint(f"Setting the output directory path: `{outputpath}`") - # Gather data mapper vprint(f"- Searching for data in `{benchmark_directory}`") dtmapper = get_data_mapper( @@ -2104,17 +1882,20 @@ def main( subset_scenarios=scenarios, read_from_archive=read_from_archive, ) - # Auto-detect a human name for each caprieval stage from the module it # evaluated (overrides the hardcoded CAPRIEVAL_STEPS for plot titles). global DETECTED_STAGE_NAMES DETECTED_STAGE_NAMES = detect_stage_modules(dtmapper, read_from_archive) if DETECTED_STAGE_NAMES: + # Show the EFFECTIVE label per stage (a CAPRIEVAL_STEPS entry + # overrides the auto-detected name), so the log matches the plots. vprint( - "- Auto-named caprieval stages: " - + ", ".join(f"{k}->{v}" for k, v in sorted(DETECTED_STAGE_NAMES.items())) + "- Caprieval stage names: " + + ", ".join( + f"{k}->{stage_name(k)}" + for k in sorted(DETECTED_STAGE_NAMES) + ) ) - # Initiate all scenarios perforamces mapper vprint( f"- Loading data from `{benchmark_directory}` " @@ -2135,15 +1916,12 @@ def main( # Hold perforamnces all_scenar_perfs[scenar_name] = scenar_best_perfs all_scenar_melquis[scenar_name] = scenar_pdb_perfs - # Write data as json write_json(all_scenar_perfs, f"{base_outputpath}_performances.json") - # Plot title: when the user passes -l/--label, use it verbatim as the # title. Otherwise fall back to a plain generic title. (Note: `basename` # is still used for output filenames, so it stays as-is above.) - plot_title = label if label else "Benchmark output" - + plot_title = label if label else f"Benchmark output - {metric}" # Draw general graph if not no_capriplots: vprint("- Generating Bar plots") @@ -2154,8 +1932,7 @@ def main( progress=not silent, no_percentage=no_percentage, ) - - if not no_violinplots: + if violinplots: vprint("- Generating Violin plots") gen_full_comparison_violins( all_scenar_perfs, @@ -2164,8 +1941,7 @@ def main( metric=metric, progress=not silent, ) - - if not no_melquiplots: + if melquiplots: vprint("- Generating Melqui plots") gen_full_comparison_melquiplots( all_scenar_melquis, @@ -2173,7 +1949,6 @@ def main( basepath=base_outputpath, progress=not silent, ) - # Draw per-scenario graphs, so each scenario also gets its own # standalone bar/violin plot, not just the combined comparison figure. # Skipped when there's only a single scenario overall, since in that @@ -2192,7 +1967,7 @@ def main( progress=not silent, no_percentage=no_percentage, ) - if not no_violinplots: + if violinplots: vprint(f"- Generating Violin plot for scenario `{scenar_name}`") gen_full_comparison_violins( {scenar_name: scenar_perfs}, @@ -2201,14 +1976,11 @@ def main( metric=metric, progress=not silent, ) - - ################################### # COMMAND LINE ARGUMENTS HANDLERS # ################################### def must_end_with_slash(dirpath: str) -> str: """Make sure a directory paths is terminating by '/'. - NOTE: only appends the slash if `dirpath` already exists as a directory. Suitable for `benchmark_directory`, which must already exist. NOT suitable for an output path that hasn't been created yet @@ -2217,7 +1989,6 @@ def must_end_with_slash(dirpath: str) -> str: prefix follows (e.g. `AnalysisHADDOCK24_ab_initio_...` instead of `Analysis/HADDOCK24_ab_initio_...`), and can leave output files sitting loose inside the benchmark directory itself. - Parameters ---------- dirpath : str @@ -2226,15 +1997,11 @@ def must_end_with_slash(dirpath: str) -> str: if dirpath[-1] != '/' and Path(dirpath).is_dir(): dirpath += '/' return dirpath - - def ensure_trailing_slash(dirpath: str) -> str: """Unconditionally make sure a path ends with '/'. - Unlike `must_end_with_slash`, this does not require the path to already exist -- appropriate for output directories, which are created later via `gen_outputdir`. - Parameters ---------- dirpath : str @@ -2243,11 +2010,8 @@ def ensure_trailing_slash(dirpath: str) -> str: if not dirpath.endswith('/'): dirpath += '/' return dirpath - - def parse_cmd_line() -> argparse.Namespace: """Parse command line argument. - Return ------ args : argparse.Namespace @@ -2301,16 +2065,15 @@ def parse_cmd_line() -> argparse.Namespace: f"{', '.join(PERFORMANCES_CLASSES.keys())}" ) return args - - def _get_cmd_line_args() -> argparse.Namespace: """Define and parse the command line arguments.""" parser = argparse.ArgumentParser( description=( "Analyse a HADDOCK3 + haddock-runner benchmark. Expects " "///run1/ layout " - "(one or more scenarios). Generates a combined bar plot, " - "violin plot, per-scenario melquiplots (zipped), and a JSON " + "(one or more scenarios). Generates a combined bar plot by " + "default; violin plots and per-scenario melquiplots (zipped) " + "are optional (--violinplots / --melquiplots), plus a JSON " "summary, written by default to Analysis/ inside the " "benchmark directory." ), @@ -2383,7 +2146,7 @@ def _get_cmd_line_args() -> argparse.Namespace: '--dpi', help="DPI of the generated figures.", required=False, - default=400, + default=200, type=int, ) parser.add_argument( @@ -2393,14 +2156,14 @@ def _get_cmd_line_args() -> argparse.Namespace: default=False, ) parser.add_argument( - '--no-violinplots', - help="Do not generate violin plots", + '--violinplots', + help="Also generate violin plots (off by default).", action="store_true", default=False, ) parser.add_argument( - '--no-melquiplots', - help="Do not generate melqui plots", + '--melquiplots', + help="Also generate melqui plots (off by default).", action="store_true", default=False, ) @@ -2443,11 +2206,8 @@ def _get_cmd_line_args() -> argparse.Namespace: ) args = parser.parse_args() return args - - def welcome_msg(silent: bool) -> None: """Print welcome message - Parameters ---------- silent : bool @@ -2462,8 +2222,6 @@ def welcome_msg(silent: bool) -> None: f"{'#' * 80}\n" ) _vprint(silent, msg) - - ############################ # COMMAND LINE ENTRY POINT # ############################ @@ -2479,13 +2237,11 @@ def maincli() -> None: silent=args.quiet, no_percentage=args.no_percentage, no_capriplots=args.no_capriplots, - no_violinplots=args.no_violinplots, - no_melquiplots=args.no_melquiplots, + violinplots=args.violinplots, + melquiplots=args.melquiplots, read_from_archive=args.from_archive, label=args.label, per_scenario_plots=args.per_scenario_plots, ) - - if __name__ == "__main__": maincli() From 9d2cd3944ac9512d5ea22dfdced56392b50658ef Mon Sep 17 00:00:00 2001 From: Shantanu Khatri <88526089+Comp-era@users.noreply.github.com> Date: Thu, 30 Jul 2026 14:01:16 +0200 Subject: [PATCH 3/8] Fix stale references to old AnalyseBenchmarkResults.py filename The script was renamed to lowercase (analysebenchmarkresults.py) but analyse.sh, README.md, analysis/README.md, and CONTRIBUTING.md still pointed at the old capitalized name, which would break on case-sensitive filesystems. --- CONTRIBUTING.md | 2 +- README.md | 4 ++-- analyse.sh | 4 ++-- analysis/README.md | 12 ++++++------ analysis/analysebenchmarkresults.py | 2 +- 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 2c1eac7..e70d596 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -26,7 +26,7 @@ If you want to add an entirely new molecular system (e.g. protein-RNA, antibody- ### Analysis Improvements -The analysis pipeline lives in `analysis/`. Improvements to `AnalyseBenchmarkResults.py` — such as support for new metrics, additional plot types, or better output formatting — are welcome. If you change the output format, update `analysis/README.md` accordingly. +The analysis pipeline lives in `analysis/`. Improvements to `analysebenchmarkresults.py` — such as support for new metrics, additional plot types, or better output formatting — are welcome. If you change the output format, update `analysis/README.md` accordingly. ### Bug Fixes and Path Issues diff --git a/README.md b/README.md index f0c3dee..1b1c987 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ Benchmarking/ ├── setup.sh # Environment setup entry point ├── run.sh # Wrapper: activates env, execs haddock-runner -├── analyse.sh # Wrapper: runs analysis/AnalyseBenchmarkResults.py +├── analyse.sh # Wrapper: runs analysis/analysebenchmarkresults.py ├── USAGE.md # Full usage guide ├── versions.env # Pinned dataset commit SHAs + haddock-runner/haddock3 versions ├── scripts/ # Individual setup steps, orchestrated by setup.sh @@ -118,7 +118,7 @@ sequenceDiagram HR-->>User: exit - All jobs completed - Note over User,FS: External - AnalyseBenchmarkResults.py reads from work_dir/ + Note over User,FS: External - analysebenchmarkresults.py reads from work_dir/ ``` ## Benchmark Systems diff --git a/analyse.sh b/analyse.sh index 73fd465..cdf952c 100755 --- a/analyse.sh +++ b/analyse.sh @@ -1,6 +1,6 @@ #!/usr/bin/env bash # -# Wrapper around analysis/AnalyseBenchmarkResults.py: runs it with the +# Wrapper around analysis/analysebenchmarkresults.py: runs it with the # venv's python3 so numpy/matplotlib are there without activating anything. # # Usage: @@ -21,4 +21,4 @@ if [ ! -x "$VENV_PATH/bin/python3" ]; then exit 1 fi -exec "$VENV_PATH/bin/python3" "$REPO_ROOT/analysis/AnalyseBenchmarkResults.py" "$@" +exec "$VENV_PATH/bin/python3" "$REPO_ROOT/analysis/analysebenchmarkresults.py" "$@" diff --git a/analysis/README.md b/analysis/README.md index f6d0852..5f01e57 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -6,11 +6,11 @@ This directory contains the post-processing and visualisation script for HADDOCK **Per-run analysis** — HADDOCK3 automatically generates an `analysis/` folder inside each target's `run1/` directory as part of the docking workflow. For every `caprieval` step it writes a `N_caprieval_analysis/` subfolder (e.g. `2_caprieval_analysis/`) containing an HTML report (`report.html`), per-metric plots (`*_clt.html`), and the `capri_ss.tsv` evaluation file. The HTML reports must be served over HTTP to render correctly — run `python -m http.server --directory .` from the `analysis/` folder and open the printed URL in a browser. This folder is produced by the `caprieval` module and requires no manual step. -**Benchmark-wide analysis (manual)** — once all targets across a scenario have finished, `AnalyseBenchmarkResults.py` aggregates the `capri_ss.tsv` files from every target and produces overall performance plots and a JSON summary across the full dataset. This is what you run manually after the benchmark completes. +**Benchmark-wide analysis (manual)** — once all targets across a scenario have finished, `analysebenchmarkresults.py` aggregates the `capri_ss.tsv` files from every target and produces overall performance plots and a JSON summary across the full dataset. This is what you run manually after the benchmark completes. ## Script -### `AnalyseBenchmarkResults.py` +### `analysebenchmarkresults.py` **Version**: 1.1.1 **Author**: BonvinLab, Computational Structural Biology group, Utrecht University @@ -48,7 +48,7 @@ If HADDOCK3 was run with `gen_archive = true`, the analysis can be read directly ## Basic Usage ```bash -python3 analysis/AnalyseBenchmarkResults.py +python3 analysis/analysebenchmarkresults.py ``` Output files are written to an `analysis/` subdirectory by default. @@ -133,13 +133,13 @@ The keys must match the two-digit index prefix of the caprieval directories in t Analyse all scenarios for a protein-protein benchmark run: ```bash -python3 analysis/AnalyseBenchmarkResults.py results/protein-protein/ -t protein -m irmsd +python3 analysis/analysebenchmarkresults.py results/protein-protein/ -t protein -m irmsd ``` Analyse only two specific scenarios and suppress plots: ```bash -python3 analysis/AnalyseBenchmarkResults.py results/protein-protein/ \ +python3 analysis/analysebenchmarkresults.py results/protein-protein/ \ -s scenario-HADDOCK3_clustfcc \ --no-melquiplots ``` @@ -147,5 +147,5 @@ python3 analysis/AnalyseBenchmarkResults.py results/protein-protein/ \ Read results from compressed archives (when `gen_archive = true` was set in HADDOCK3): ```bash -python3 analysis/AnalyseBenchmarkResults.py results/protein-peptide/ -t peptide -a +python3 analysis/analysebenchmarkresults.py results/protein-peptide/ -t peptide -a ``` diff --git a/analysis/analysebenchmarkresults.py b/analysis/analysebenchmarkresults.py index 989b2a5..c54c74d 100644 --- a/analysis/analysebenchmarkresults.py +++ b/analysis/analysebenchmarkresults.py @@ -6,7 +6,7 @@ Please modify the Global variable: CAPRIEVAL_STEPS to suite your needs. It is used to generate nice title to the caprieval steps. Usage: ->python3 AnalyseBenchmarkResults.py +>python3 analysebenchmarkresults.py Expected input structure ------------------------- / From c84313cd688ed1c448440bf0c0ba63298f34cd26 Mon Sep 17 00:00:00 2001 From: Shantanu Khatri <88526089+Comp-era@users.noreply.github.com> Date: Thu, 30 Jul 2026 14:09:55 +0200 Subject: [PATCH 4/8] Update analysis/README.md to document v2.0.0 script behaviour Reflects the scenario-outer/PDBid-inner layout, type-aware default metric, auto-detected caprieval stage names, new -l/--label and --per-scenario-plots flags, and violin/melquiplots being opt-in. --- analysis/README.md | 116 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 89 insertions(+), 27 deletions(-) diff --git a/analysis/README.md b/analysis/README.md index 5f01e57..ca05502 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -12,11 +12,12 @@ This directory contains the post-processing and visualisation script for HADDOCK ### `analysebenchmarkresults.py` -**Version**: 1.1.1 -**Author**: BonvinLab, Computational Structural Biology group, Utrecht University +**Version**: 2.0.0 The script reads `capri_ss.tsv` files produced by the `caprieval` module of HADDOCK3, ranks models by their HADDOCK score, and computes the quality of the best-ranking model at a series of top-X thresholds (Top1, Top5, Top10, Top20, Top50, Top100, Top200, Top500, Top1000). Results are classified into CAPRI performance categories and visualised across all scenarios and pipeline stages. +Caprieval stage names are **auto-detected** from the module each `caprieval` evaluated, so column titles read `rigidbody`, `flexref`, `emref`, `seletopclusts`, etc. without manual configuration (see *Customising Caprieval Stage Labels*). + ## Requirements ```bash @@ -27,12 +28,12 @@ Both packages are already installed if you ran `setup.sh`. ## Expected Input Structure -The script expects a benchmark results directory in the format produced by `haddock-runner`. Each target is a subdirectory named by its PDB ID, and within each target there is one subdirectory per scenario, each containing a `run1/` directory with the HADDOCK3 output. +The script expects a benchmark results directory in the format produced by `haddock-runner`. The **scenario** is the outer directory and the **target (PDB ID)** is the inner one — each scenario contains one subdirectory per target, each with a `run1/` directory holding the HADDOCK3 output. ``` / - / - / + / + / run1/ 02_caprieval/ capri_ss.tsv @@ -43,7 +44,7 @@ The script expects a benchmark results directory in the format produced by `hadd ... ``` -If HADDOCK3 was run with `gen_archive = true`, the analysis can be read directly from the compressed analysis archives (`run1_analysis.tgz`) using the `--from-archive` flag. +Scenario discovery accepts a target that holds **either** a live `run1/` directory (default mode) **or** a `run1_analysis.tgz` archive (`--from-archive` mode), so it works whether or not the runs were archived. If HADDOCK3 was run with `gen_archive = true`, the analysis can be read directly from the compressed archives using the `--from-archive` (`-a`) flag. ## Basic Usage @@ -51,28 +52,40 @@ If HADDOCK3 was run with `gen_archive = true`, the analysis can be read directly python3 analysis/analysebenchmarkresults.py ``` -Output files are written to an `analysis/` subdirectory by default. +By default, output files are written to an `Analysis/` subdirectory **inside the benchmark directory itself**. Use `-o` to write elsewhere. ## All Options ``` positional arguments: benchmark_directory Path to the directory where the benchmark was run + (contains the scenario subdirectories) optional arguments: - -o, --output_path Directory to write output files (default: analysis/) - -m, --metric Performance metric: irmsd or dockq (default: irmsd) + -o, --output_path Directory to write output files + (default: Analysis/ inside the benchmark directory) + -l, --label Custom name used for output filenames and plot titles + instead of the benchmark directory's basename + -m, --metric Performance metric: irmsd, dockq, or ilrmsd. + Default is type-aware: ilrmsd for glycan and + small_molecule, irmsd otherwise -s, --scenario Analyse only specific scenario(s) by name (default: all) - -t, --type System type for CAPRI thresholds: protein, peptide, glycan (default: protein) - -d, --dpi DPI for output figures (default: 400) + -t, --type System type for CAPRI thresholds: + protein, peptide, glycan, dna, small_molecule + (default: protein) + -d, --dpi DPI for output figures (default: 200) -n, --no-percentage Report raw counts instead of percentages in bar plots -a, --from-archive Read capri_ss.tsv from run1_analysis.tgz archives -q, --quiet Suppress all stdout output --no-capriplots Skip CAPRI bar plot generation - --no-violinplots Skip violin plot generation - --no-melquiplots Skip melquiplot generation + --violinplots Also generate violin plots (OFF by default) + --melquiplots Also generate melquiplots (OFF by default) + --per-scenario-plots In addition to the combined comparison figure, also + write a standalone bar/violin plot per scenario ``` +> **Changed in 2.0.0:** bar plots are generated by default; **violin and melqui plots are now OFF by default** and must be enabled with `--violinplots` / `--melquiplots` (the old `--no-violinplots` / `--no-melquiplots` flags were removed). + ## Output Files Running the script produces the following files in the output directory: @@ -80,16 +93,30 @@ Running the script produces the following files in the output directory: | File | Description | |---|---| | `*_performances.json` | JSON summary of best model quality at each top-X threshold per scenario and caprieval stage | -| `*_capribarpolots.png` | Stacked bar plots showing the fraction of targets with High / Medium / Acceptable / Near-acceptable / Low quality models at each threshold, for each scenario and caprieval stage | -| `*_violins.png` | Violin plots showing the distribution of performance metric values across all targets, broken down by scenario and caprieval stage | -| `*_melquiplots.zip` | Archive of per-scenario melquiplots — one column per target, one row per caprieval stage, colour-coded by CAPRI quality class | -| `*_benchmark_mapper.json` | Internal mapping of scenario/target/stage to file paths (cached to speed up re-runs) | +| `*_capribarplots.png` | Stacked bar plots showing the fraction of targets with High / Medium / Acceptable / Near-acceptable quality models at each threshold, for each scenario and caprieval stage | +| `*_violins.png` | (only with `--violinplots`) Violin plots of the metric distribution across all targets, by scenario and caprieval stage | +| `*_melquiplots.zip` | (only with `--melquiplots`) Archive of per-scenario melquiplots — one column per target, one row per caprieval stage, colour-coded by CAPRI quality class | +| `*_benchmark_mapper.json` | Internal mapping of scenario/target/stage to file paths | + +The plot title is `Benchmark output - ` by default (e.g. `Benchmark output - irmsd`), or your `-l/--label` value if given. + +## Metric by system type + +The default metric depends on the system type; an explicit `-m` always wins. Requesting a metric a type does not define (e.g. `-t glycan -m irmsd`) exits with a clear error. + +| Type | irmsd | ilrmsd | dockq | Default | +|---|:---:|:---:|:---:|---| +| protein | ✓ | | ✓ | irmsd | +| peptide | ✓ | | ✓ | irmsd | +| dna | ✓ | | ✓ | irmsd | +| glycan | | ✓ | ✓ | ilrmsd | +| small_molecule | | ✓ | ✓ | ilrmsd | ## Performance Classes -Models are classified according to CAPRI criteria. The thresholds differ by system type and metric: +Models are classified according to CAPRI criteria. The thresholds differ by system type and metric. -**Protein-Protein / Protein-DNA (iRMSD)** +**Protein / DNA (iRMSD)** | Class | iRMSD range | |---|---| @@ -98,9 +125,8 @@ Models are classified according to CAPRI criteria. The thresholds differ by syst | Acceptable | 2 – 4 Å | | Near-acceptable | 4 – 6 Å | | Low | > 6 Å | -| Missing | No model generated | -**Protein-Peptide / Protein-Glycan (iRMSD)** +**Peptide (iRMSD)** | Class | iRMSD range | |---|---| @@ -110,11 +136,19 @@ Models are classified according to CAPRI criteria. The thresholds differ by syst | Near-acceptable | 2 – 3 Å | | Low | > 3 Å | -DockQ thresholds are also available (`--metric dockq`) with corresponding class boundaries defined in the script. +**Glycan / Small molecule (ilRMSD)** + +| Class | ilRMSD range | +|---|---| +| High | < 1 Å | +| Medium | 1 – 2 Å | +| Acceptable | 2 – 3 Å | +| Near-acceptable | 3 – 4 Å | +| Low | > 4 Å | ## Customising Caprieval Stage Labels -By default, caprieval stages are labelled by their module index (e.g. `02`, `04`). You can assign human-readable names by editing the `CAPRIEVAL_STEPS` dictionary near the top of the script: +Caprieval stages are auto-named at runtime from the module each one evaluated, so no configuration is needed in the common case. To **override** a label, edit the `CAPRIEVAL_STEPS` dictionary near the top of the script — a manual entry takes precedence over the auto-detected name: ```python CAPRIEVAL_STEPS = { @@ -122,11 +156,33 @@ CAPRIEVAL_STEPS = { '04': 'seletop 200', '06': 'flexref', '08': 'emref', - '11': 'top 4 models per fcc clust', + '12': 'seletopclusts', } ``` -The keys must match the two-digit index prefix of the caprieval directories in the HADDOCK3 run output (e.g. `02_caprieval/`). Any stage not listed here will be labelled `_caprieval` automatically. +The keys must be strings matching the zero-padded index prefix of the caprieval directories (e.g. `'06'` for `06_caprieval/`), exactly as printed in the `- Caprieval stage names:` log line. Any stage not listed keeps its auto-detected name (or falls back to `_caprieval` if it could not be detected). + +The label can be **any text you like** — you are not limited to the module name. For example, to give fully custom, descriptive titles: + +```python +CAPRIEVAL_STEPS = { + '02': 'Rigid-body docking', + '04': 'Top-200 selection', + '06': 'Flexible refinement', + '08': 'Energy minimisation', + '12': 'Final clustered models', +} +``` + +With this, the plot column titles (and the `- Caprieval stage names:` log) read your custom text instead of the auto-detected `rigidbody` / `flexref` / `emref` / `seletopclusts`. You can also override just a subset — list only the stages you want to rename and leave the rest to auto-detection: + +```python +CAPRIEVAL_STEPS = { + '12': 'Final clustered models', # only this one is renamed +} +``` + +> **Reminder:** the keys are the stage **indices** as they appear in *your* run, which can differ between benchmarks/pipelines (e.g. `seletopclusts` may be `11` in one run and `12` in another). Check the `- Caprieval stage names:` line the script prints and use those indices. ## Examples @@ -136,12 +192,12 @@ Analyse all scenarios for a protein-protein benchmark run: python3 analysis/analysebenchmarkresults.py results/protein-protein/ -t protein -m irmsd ``` -Analyse only two specific scenarios and suppress plots: +Analyse only one scenario and also enable the melquiplot: ```bash python3 analysis/analysebenchmarkresults.py results/protein-protein/ \ -s scenario-HADDOCK3_clustfcc \ - --no-melquiplots + --melquiplots ``` Read results from compressed archives (when `gen_archive = true` was set in HADDOCK3): @@ -149,3 +205,9 @@ Read results from compressed archives (when `gen_archive = true` was set in HADD ```bash python3 analysis/analysebenchmarkresults.py results/protein-peptide/ -t peptide -a ``` + +Analyse a glycan benchmark (defaults to the ilrmsd metric): + +```bash +python3 analysis/analysebenchmarkresults.py results/protein-glycan/ -t glycan -a +``` From 68034a534ec1829ecb0373b0291637548fb0b848 Mon Sep 17 00:00:00 2001 From: Alexandre Bonvin Date: Fri, 31 Jul 2026 11:54:07 +0200 Subject: [PATCH 5/8] Added ligand top/param to topoaa stage --- docking_benchmarks/protein_protein/unbound_ab-initio_cm.yaml | 2 ++ .../protein_protein/unbound_ab-initio_ranair.yaml | 2 ++ docking_benchmarks/protein_protein/unbound_true-interface.yaml | 2 ++ .../protein_protein/unbound_true-interface5A-cltsel.yaml | 2 ++ .../protein_protein/unbound_true-interface5A.yaml | 2 ++ 5 files changed, 10 insertions(+) diff --git a/docking_benchmarks/protein_protein/unbound_ab-initio_cm.yaml b/docking_benchmarks/protein_protein/unbound_ab-initio_cm.yaml index 34daa11..35c8ef9 100644 --- a/docking_benchmarks/protein_protein/unbound_ab-initio_cm.yaml +++ b/docking_benchmarks/protein_protein/unbound_ab-initio_cm.yaml @@ -15,6 +15,8 @@ scenarios: - name: run-unbound_ab-initio_cm workflow: topoaa: + ligand_param_fname: _ligand.param + ligand_top_fname: _ligand.top rigidbody: sampling: 10000 diff --git a/docking_benchmarks/protein_protein/unbound_ab-initio_ranair.yaml b/docking_benchmarks/protein_protein/unbound_ab-initio_ranair.yaml index 6951e1a..0d26370 100644 --- a/docking_benchmarks/protein_protein/unbound_ab-initio_ranair.yaml +++ b/docking_benchmarks/protein_protein/unbound_ab-initio_ranair.yaml @@ -15,6 +15,8 @@ scenarios: - name: run-unbound_ab-initio_ranair workflow: topoaa: + ligand_param_fname: _ligand.param + ligand_top_fname: _ligand.top rigidbody: sampling: 10000 diff --git a/docking_benchmarks/protein_protein/unbound_true-interface.yaml b/docking_benchmarks/protein_protein/unbound_true-interface.yaml index 497698c..7d42c21 100644 --- a/docking_benchmarks/protein_protein/unbound_true-interface.yaml +++ b/docking_benchmarks/protein_protein/unbound_true-interface.yaml @@ -15,6 +15,8 @@ scenarios: - name: run-unbound_true-interface workflow: topoaa: + ligand_param_fname: _ligand.param + ligand_top_fname: _ligand.top rigidbody: ambig_fname: _ambig.tbl diff --git a/docking_benchmarks/protein_protein/unbound_true-interface5A-cltsel.yaml b/docking_benchmarks/protein_protein/unbound_true-interface5A-cltsel.yaml index df16750..71b19b0 100644 --- a/docking_benchmarks/protein_protein/unbound_true-interface5A-cltsel.yaml +++ b/docking_benchmarks/protein_protein/unbound_true-interface5A-cltsel.yaml @@ -15,6 +15,8 @@ scenarios: - name: run-unbound_true-interface5A-cltsel workflow: topoaa: + ligand_param_fname: _ligand.param + ligand_top_fname: _ligand.top rigidbody: ambig_fname: _ambig5.tbl diff --git a/docking_benchmarks/protein_protein/unbound_true-interface5A.yaml b/docking_benchmarks/protein_protein/unbound_true-interface5A.yaml index f8416bd..3a88497 100644 --- a/docking_benchmarks/protein_protein/unbound_true-interface5A.yaml +++ b/docking_benchmarks/protein_protein/unbound_true-interface5A.yaml @@ -15,6 +15,8 @@ scenarios: - name: run-unbound_true-interface5A workflow: topoaa: + ligand_param_fname: _ligand.param + ligand_top_fname: _ligand.top rigidbody: ambig_fname: _ambig5.tbl From 17606e33a66924748a6899a9f45e9c223ff4d574 Mon Sep 17 00:00:00 2001 From: VGPReys Date: Fri, 31 Jul 2026 12:00:53 +0200 Subject: [PATCH 6/8] custom labels fpath --- analysis/README.md | 40 +++++------------ analysis/analysebenchmarkresults.py | 67 +++++++++++++++++++++------- analysis/custom-labels-examples.json | 7 +++ 3 files changed, 70 insertions(+), 44 deletions(-) create mode 100644 analysis/custom-labels-examples.json diff --git a/analysis/README.md b/analysis/README.md index ca05502..38dd85c 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -148,40 +148,24 @@ Models are classified according to CAPRI criteria. The thresholds differ by syst ## Customising Caprieval Stage Labels -Caprieval stages are auto-named at runtime from the module each one evaluated, so no configuration is needed in the common case. To **override** a label, edit the `CAPRIEVAL_STEPS` dictionary near the top of the script — a manual entry takes precedence over the auto-detected name: - -```python -CAPRIEVAL_STEPS = { - '02': 'rigidbody', - '04': 'seletop 200', - '06': 'flexref', - '08': 'emref', - '12': 'seletopclusts', -} -``` +Caprieval stages are auto-named at runtime from the module each one evaluated, so no configuration is required by default. +Nevertheless, the user can provide a JSON file containing custom labeling scheme for all (or only some) for the steps. +This can be done by using the `custom-labels` parameter, followed by the path to that file. -The keys must be strings matching the zero-padded index prefix of the caprieval directories (e.g. `'06'` for `06_caprieval/`), exactly as printed in the `- Caprieval stage names:` log line. Any stage not listed keeps its auto-detected name (or falls back to `_caprieval` if it could not be detected). +The keys must be strings matching the zero-padded index prefix of the caprieval directories (e.g. `"06"` for `06_caprieval/`), exactly as printed in the `- Caprieval stage names:` log line. Any stage not listed keeps its auto-detected name (or falls back to `_caprieval` if it could not be detected). The label can be **any text you like** — you are not limited to the module name. For example, to give fully custom, descriptive titles: - -```python -CAPRIEVAL_STEPS = { - '02': 'Rigid-body docking', - '04': 'Top-200 selection', - '06': 'Flexible refinement', - '08': 'Energy minimisation', - '12': 'Final clustered models', -} -``` - -With this, the plot column titles (and the `- Caprieval stage names:` log) read your custom text instead of the auto-detected `rigidbody` / `flexref` / `emref` / `seletopclusts`. You can also override just a subset — list only the stages you want to rename and leave the rest to auto-detection: - -```python -CAPRIEVAL_STEPS = { - '12': 'Final clustered models', # only this one is renamed +```json +{ + "02": "Rigid-body docking", + "04": "Top-200 selection", + "06": "Flexible refinement", + "08": "Energy minimisation", + "12": "Final clustered models" } ``` +**Important note:** in JSON, only use double-quotes (`"`) for strings ! > **Reminder:** the keys are the stage **indices** as they appear in *your* run, which can differ between benchmarks/pipelines (e.g. `seletopclusts` may be `11` in one run and `12` in another). Check the `- Caprieval stage names:` line the script prints and use those indices. ## Examples diff --git a/analysis/analysebenchmarkresults.py b/analysis/analysebenchmarkresults.py index c54c74d..9d77ef0 100644 --- a/analysis/analysebenchmarkresults.py +++ b/analysis/analysebenchmarkresults.py @@ -64,15 +64,6 @@ #################### # GLOBAL VARIABLES # #################### -# Caprieval stage names are now auto-detected at runtime from the module each -# caprieval evaluated (see detect_stage_modules() / DETECTED_STAGE_NAMES), so -# they are correct per benchmark type without manual editing. -# -# This dict is only an OPTIONAL manual override / fallback, used solely for a -# stage that could not be auto-detected. Leave it empty to let auto-detection -# handle everything (undetected stages then fall back to `_caprieval`). -# To force a custom label for a specific stage, add e.g. '02': 'rigidbody'. -CAPRIEVAL_STEPS: dict[str, str] = {} # Set threshold of top X structures to take into account TOP_X_THRESHOLDS = (1, 5, 10, 20, 50, 100, 200, 500, 1000, ) # Set number of entries to display in melquiplot @@ -203,6 +194,12 @@ ) # DPI of the generated figures DPI = 200 +# Name of the stages +CAPRIEVAL_STEPS: dict[str, str] = {} +# Auto-detected {stage_index: module_name}, filled at runtime by +# detect_stage_modules(). Takes precedence over the hardcoded CAPRIEVAL_STEPS. +DETECTED_STAGE_NAMES: dict[str, str] = {} + #################### # DEFINE FUNCTIONS # #################### @@ -409,9 +406,25 @@ def adjacent_values( lower_adjacent_value = q1 - (q3 - q1) * 1.5 lower_adjacent_value = np.clip(lower_adjacent_value, vals[0], q1) return lower_adjacent_value, upper_adjacent_value -# Auto-detected {stage_index: module_name}, filled at runtime by -# detect_stage_modules(). Takes precedence over the hardcoded CAPRIEVAL_STEPS. -DETECTED_STAGE_NAMES: dict[str, str] = {} +def load_custom_stage_labels(custom_labels_fpath: Optional[str]) -> dict[str, str]: + """Load user provided custom labeling scheme. + + Parameters + ---------- + custom_labels_fpath : Optional[str] + Index of the caprieval stage + Returns + ------- + custom_labels : dict[str, str] + Name of the stages + """ + if custom_labels_fpath and os.path.exists(custom_labels_fpath): + try: + with open(custom_labels_fpath, "r") as fin: + return json.load(fin) + except: + return {} + return {} def stage_name(cname: str) -> str: """Return a human-readable name for a caprieval stage. Preference order: @@ -532,7 +545,7 @@ def gen_full_comparison_violins( # Set labels on last row only labels = None if ri + 1 == nb_thresh: - labels = [so.replace('scenario-', '') for so in scenars_order] + labels = [so for so in scenars_order] # Write bars gen_violin( ax, @@ -1812,6 +1825,7 @@ def _vprint(silence: bool, msg: str) -> None: """ if not silence: print(msg) + ################# # Main function # ################# @@ -1827,6 +1841,7 @@ def main( melquiplots: bool = False, read_from_archive: bool = False, label: Optional[str] = None, + custom_labels: Optional[str] = None, per_scenario_plots: bool = False, ) -> None: """Run the analysis procedure. @@ -1858,6 +1873,8 @@ def main( label : Optional[str], optional Custom name to use for output filenames and plot titles instead of the benchmark directory's basename, by default None + custom_labels: Optional[str], optional + Path to file containing the custom labeling scheme for scenarios per_scenario_plots : bool, optional In addition to the single combined comparison bar/violin plot (which already shows every scenario as its own row), also @@ -1882,6 +1899,9 @@ def main( subset_scenarios=scenarios, read_from_archive=read_from_archive, ) + # Potentially load the custom labeling scheme provided by the user + global CAPRIEVAL_STEPS + CAPRIEVAL_STEPS = load_custom_stage_labels(custom_labels) # Auto-detect a human name for each caprieval stage from the module it # evaluated (overrides the hardcoded CAPRIEVAL_STEPS for plot titles). global DETECTED_STAGE_NAMES @@ -1976,6 +1996,7 @@ def main( metric=metric, progress=not silent, ) + ################################### # COMMAND LINE ARGUMENTS HANDLERS # ################################### @@ -2100,7 +2121,7 @@ def _get_cmd_line_args() -> argparse.Namespace: ) parser.add_argument( '-l', - '--label', + '--plot-label', help=( "Custom name to use for output filenames and plot titles, " "instead of the benchmark directory's basename (useful when " @@ -2126,7 +2147,11 @@ def _get_cmd_line_args() -> argparse.Namespace: parser.add_argument( '-s', '--scenario', - help="Name(s) of a specific scenario(s) to analyze. Can be multiple of them, separated by space. By default, all scenarios will be analysed together.", + help=( + "Name(s) of a specific scenario(s) to analyze. " + "Can be multiple of them, separated by space. " + "By default, all scenarios will be analysed together." + ), required=False, default=None, nargs="+", @@ -2204,6 +2229,14 @@ def _get_cmd_line_args() -> argparse.Namespace: action="store_true", default=False, ) + parser.add_argument( + "--custom-labels", + "-c", + help="Path to json file containing custom labeling of steps. ", + required=False, + default=None, + type=str, + ) args = parser.parse_args() return args def welcome_msg(silent: bool) -> None: @@ -2222,6 +2255,7 @@ def welcome_msg(silent: bool) -> None: f"{'#' * 80}\n" ) _vprint(silent, msg) + ############################ # COMMAND LINE ENTRY POINT # ############################ @@ -2240,7 +2274,8 @@ def maincli() -> None: violinplots=args.violinplots, melquiplots=args.melquiplots, read_from_archive=args.from_archive, - label=args.label, + label=args.plot_label, + custom_labels=args.custom_labels, per_scenario_plots=args.per_scenario_plots, ) if __name__ == "__main__": diff --git a/analysis/custom-labels-examples.json b/analysis/custom-labels-examples.json new file mode 100644 index 0000000..4994fd3 --- /dev/null +++ b/analysis/custom-labels-examples.json @@ -0,0 +1,7 @@ +{ + "02": "Rigid-body docking", + "04": "Top-200 selection", + "06": "Flexible refinement", + "08": "Energy minimisation", + "12": "Final clustered models" +} From 39c10215092d726603e02013a6e66cc1ffa771f2 Mon Sep 17 00:00:00 2001 From: Shantanu Khatri <88526089+Comp-era@users.noreply.github.com> Date: Fri, 31 Jul 2026 12:17:31 +0200 Subject: [PATCH 7/8] added example to run custom-labels --- analysis/README.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/analysis/README.md b/analysis/README.md index 38dd85c..b84a079 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -195,3 +195,8 @@ Analyse a glycan benchmark (defaults to the ilrmsd metric): ```bash python3 analysis/analysebenchmarkresults.py results/protein-glycan/ -t glycan -a ``` +Run analysis with custom-lables + +```bash +python3 analysis/analysebenchmarkresults.py results/protein-protein/ --custom-labels custom-labels-examples.json +``` \ No newline at end of file From 1415c09500769f0b569b5906461aaedc60e39a7a Mon Sep 17 00:00:00 2001 From: Shantanu Khatri <88526089+Comp-era@users.noreply.github.com> Date: Fri, 31 Jul 2026 12:19:58 +0200 Subject: [PATCH 8/8] minor fix - added custom-labels to all options in README.ms --- analysis/README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/analysis/README.md b/analysis/README.md index b84a079..2b39e32 100644 --- a/analysis/README.md +++ b/analysis/README.md @@ -82,6 +82,9 @@ optional arguments: --melquiplots Also generate melquiplots (OFF by default) --per-scenario-plots In addition to the combined comparison figure, also write a standalone bar/violin plot per scenario + --custom-labels Path to JSON file containing custom labeling of + caprieval steps (see *Customising Caprieval Stage + Labels*) ``` > **Changed in 2.0.0:** bar plots are generated by default; **violin and melqui plots are now OFF by default** and must be enabled with `--violinplots` / `--melquiplots` (the old `--no-violinplots` / `--no-melquiplots` flags were removed).