From 9130250d3ffad8658792fd526fe5f81aad11dfc1 Mon Sep 17 00:00:00 2001 From: ziadbkh Date: Thu, 24 Jul 2025 14:39:06 +1000 Subject: [PATCH 1/7] add reports --- modules/local/generate_report.nf | 23 +++++++++++++++++++ subworkflows/local/run_bindcraft.nf | 8 +++++++ .../utils_nfcore_bindflow_pipeline/main.nf | 2 +- tower.yml | 7 ++++-- 4 files changed, 37 insertions(+), 3 deletions(-) create mode 100644 modules/local/generate_report.nf diff --git a/modules/local/generate_report.nf b/modules/local/generate_report.nf new file mode 100644 index 0000000..83eef42 --- /dev/null +++ b/modules/local/generate_report.nf @@ -0,0 +1,23 @@ +process GENERATE_REPORT { + tag "${meta.id}" + label 'process_single' + container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? + 'https://depot.galaxyproject.org/singularity/multiqc:1.21--pyhdfd78af_0' : + 'biocontainers/multiqc:1.21--pyhdfd78af_0' }" + conda "bioconda::multiqc=1.21" + + input: + tuple val(meta), val(info) + + output: + tuple val(meta), path ("*report.html"), emit: report + + when: + task.ext.when == null || task.ext.when + script: + def args = task.ext.args ?: '' + + """ + echo "${info}" > ${meta.id}_report.html + """ +} diff --git a/subworkflows/local/run_bindcraft.nf b/subworkflows/local/run_bindcraft.nf index f74ee25..d0c0ff8 100644 --- a/subworkflows/local/run_bindcraft.nf +++ b/subworkflows/local/run_bindcraft.nf @@ -12,6 +12,8 @@ include { JSONMANAGER } from '../../modules/local/jsonmanager' include { samplesheetToList } from 'plugin/nf-schema' include { BINDCRAFT } from '../../modules/local/bindcraft' include { RANKER } from '../../modules/local/ranker' +include { GENERATE_REPORT } from '../../modules/local/generate_report' + /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBWORKFLOW TO INITIALISE PIPELINE @@ -82,6 +84,11 @@ workflow RUN_BINDCRAFT { BINDCRAFT.out.stats.map{[["id": it[0].id], it[1]]}.groupTuple(), BINDCRAFT.out.accepted_ranked.map{[["id": it[0].id], it[1]]}.groupTuple() ) + + GENERATE_REPORT( + RANKER.out.stats.map{[it[0], "Testing report information"]} + ) + emit: input_json = JSONMANAGER.out.json @@ -93,6 +100,7 @@ workflow RUN_BINDCRAFT { output_dir = BINDCRAFT.out.output_dir stats = RANKER.out.stats ranked = RANKER.out.accepted_ranked + reports = GENERATE_REPORT.out.report versions = ch_versions } diff --git a/subworkflows/local/utils_nfcore_bindflow_pipeline/main.nf b/subworkflows/local/utils_nfcore_bindflow_pipeline/main.nf index cfa1525..716a1b1 100644 --- a/subworkflows/local/utils_nfcore_bindflow_pipeline/main.nf +++ b/subworkflows/local/utils_nfcore_bindflow_pipeline/main.nf @@ -223,5 +223,5 @@ def validateInputParameters() { // Print a warning when using Bindcraft // def checkBindcraftContainer() { - log.warn "You need to provide a valid path for the singularity/Docker container! or bindcraft needs to be avaialble on the system!" + error ("You need to provide a valid path for the singularity/Docker container! or bindcraft needs to be avaialble on the system!") } \ No newline at end of file diff --git a/tower.yml b/tower.yml index 787aedf..459e32c 100644 --- a/tower.yml +++ b/tower.yml @@ -1,5 +1,8 @@ reports: - multiqc_report.html: - display: "MultiQC HTML report" + "*_report.html": + display: "HTML report" samplesheet.csv: display: "Auto-created samplesheet with collated metadata and FASTQ paths" + "*_final_design_stats.csv": + display: "Stats and Ranking" + \ No newline at end of file From 0992960f753af474530ff8a969522ec94331bebc Mon Sep 17 00:00:00 2001 From: ziadbkh Date: Tue, 18 Nov 2025 16:52:54 +1100 Subject: [PATCH 2/7] adding free bindcraft --- conf/gadi.config | 3 + conf/modules.config | 2 +- docker/Dockerfile | 35 +++++++++++ docker/Dockerfile-1 | 105 ++++++++++++++++++++++++++++++++ modules/local/bindcraft/main.nf | 2 +- nextflow.config | 18 +++--- nextflow_schema.json | 2 - 7 files changed, 154 insertions(+), 13 deletions(-) create mode 100644 docker/Dockerfile create mode 100644 docker/Dockerfile-1 diff --git a/conf/gadi.config b/conf/gadi.config index fb6eeb0..5f3e694 100644 --- a/conf/gadi.config +++ b/conf/gadi.config @@ -17,6 +17,9 @@ params { use_dgxa100 = false } +env.OPENMM_PLATFORM_ORDER="CPU" +env.OPENMM_DEFAULT_PLATFORM="CPU" + process { executor = "pbspro" storage = "gdata/if89+scratch/${params.project}+gdata/${params.project}" diff --git a/conf/modules.config b/conf/modules.config index bc37ad5..c3caf69 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -31,7 +31,7 @@ process { withName: 'JSONMANAGER' { memory = 2.GB time = 1.h - ext.args2 = "${params.quote_char}" + //ext.args2 = "${params.quote_char}" } withName: 'RANKER' { diff --git a/docker/Dockerfile b/docker/Dockerfile new file mode 100644 index 0000000..3a3402f --- /dev/null +++ b/docker/Dockerfile @@ -0,0 +1,35 @@ +# Use the latest official Ubuntu image from Docker Hub +FROM mambaorg/micromamba:2.3.3 +LABEL maintainer="Ziad Al Bkhetan " +LABEL version="1.0.3" +LABEL description="A Docker image for running Free Bindcraft " +LABEL org.opencontainers.image.authors="Ziad Al Bkhetan" +LABEL org.opencontainers.image.licenses="MIT" + +# Prevent interactive prompts during package installation +ENV DEBIAN_FRONTEND=noninteractive +ARG MAMBA_DOCKERFILE_ACTIVATE=1 + +WORKDIR / +USER root +RUN mkdir -p /work && chown -R mambauser:mambauser /work +USER mambauser +WORKDIR /work +RUN micromamba install -y -n base -c conda-forge git wget gcc && micromamba clean -a -y +RUN git clone --branch v1.0.3 --depth 1 https://github.com/cytokineking/FreeBindCraft +WORKDIR FreeBindCraft + +# making it work for micromamba by replacing conda by micromamba and installing to the base env + +RUN sed -i '/\$pkg_manager create --name BindCraft/s/^/#/' install_bindcraft.sh +RUN sed -i '/conda env list | grep -w/s/^/#/' install_bindcraft.sh +RUN sed -i '/source ${CONDA_BASE}\/bin\/activate/s/^/#/' install_bindcraft.sh +RUN sed -i '/[[:space:]]*\[ "\$CONDA_DEFAULT_ENV" = "BindCraft" \]/s/^/#/' install_bindcraft.sh +RUN sed -i 's/\/micromamba/g' install_bindcraft.sh +RUN sed -i 's/\/conda-forge/g' install_bindcraft.sh + +RUN bash install_bindcraft.sh --cuda '12.4' --pkg_manager 'micromamba' --no-pyrosetta +ENV PATH="$PATH:/work/FreeBindCraft" +ENV PYTHONPATH="${PYTHONPATH:-}/work/FreeBindCraft" +# Set default command +CMD ["python", "-u", "bindcraft.py"] \ No newline at end of file diff --git a/docker/Dockerfile-1 b/docker/Dockerfile-1 new file mode 100644 index 0000000..6707219 --- /dev/null +++ b/docker/Dockerfile-1 @@ -0,0 +1,105 @@ +# This Dockerfile was adapted from the below to create a bindcraft container by Uwe Winter +# This Dockerfile was adapted for running on Azure CycleCloud by Felipe Ayora +# (bizdata.co.nz) + +# Copyright 2021 DeepMind Technologies Limited +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# Look for latest version that matches NVIDIA drivers on the Azure GPU machines: +# https://hub.docker.com/r/nvidia/cuda/tags?page=1&name=cudnn8-runtime-ubuntu18 + +# Build image: +# 1. Run: +# git clone https://github.com/deepmind/alphafold/releases/tag/vx.x.x +# cd ./alphafold +# 2. Run: +# sudo docker build -t -f ./path/to/this/Dockerfile . + +ARG CUDA=12.4.1 +FROM nvidia/cuda:${CUDA}-cudnn-runtime-ubuntu22.04 + +# Remove entrypoint messaging to stdout +RUN rm /opt/nvidia/entrypoint.d/* + +# FROM directive resets ARGS, so we specify again + +ARG CUDA_MAJOR=12 +ARG CUDACMDTOOLS=12.4 + +# Use bash to support string substitution. +SHELL ["/bin/bash", "-c"] + +RUN apt-get update \ + && DEBIAN_FRONTEND=noninteractive apt-get install -y --no-install-recommends \ + build-essential \ + cmake \ + cuda-command-line-tools-${CUDACMDTOOLS/./-} \ + wget \ + git \ + && rm -rf /var/lib/apt/lists/* \ + && apt-get autoremove -y \ + && apt-get clean + +# Install Miniconda package manager. +RUN wget -q -P /tmp \ + https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh \ + && bash /tmp/Miniconda3-latest-Linux-x86_64.sh -b -p /opt/conda \ + && rm /tmp/Miniconda3-latest-Linux-x86_64.sh + +# Install conda packages. +ENV PATH="/opt/conda/bin:$PATH" + +RUN conda tos accept --override-channels --channel https://repo.anaconda.com/pkgs/main +RUN conda tos accept --override-channels --channel https://repo.anaconda.com/pkgs/r + +COPY bindcraft /app/ +RUN chmod +x /app/install_bindcraft.sh + +WORKDIR /app/ +RUN /app/install_bindcraft.sh --cuda 12.4.1 --pkg_manager 'conda' + +# Install pip packages. + +# Add SETUID bit to the ldconfig binary so that non-root users can run it. +RUN chmod u+s /sbin/ldconfig.real + +# properly configure ldconfig and run it +RUN echo "/opt/conda/lib" > /etc/ld.so.conf.d/conda.conf \ + && echo "/opt/conda/envs/BindCraft/lib" >> /etc/ld.so.conf.d/conda.conf \ + && ldconfig + +# ENTRYPOINT does not support easily running multiple commands, so instead we +# write a shell script to wrap them up. +# We need to run `ldconfig` first to ensure GPUs are visible, due to some quirk +# with Debian. See https://github.com/NVIDIA/nvidia-docker/issues/1399 for +# details. +RUN printf '#!/bin/bash\n\ +ldconfig\n\ +__conda_setup="$('/opt/conda/bin/conda' '\''shell.bash'\'' '\''hook'\'' 2> /dev/null)"\n\ +if [ $? -eq 0 ]; then\n\ + eval "$__conda_setup"\n\ +else\n\ + if [ -f "/opt/conda/etc/profile.d/conda.sh" ]; then\n\ + . "/opt/conda/etc/profile.d/conda.sh"\n\ + else\n\ + export PATH="/opt/conda/bin:$PATH"\n\ + fi\n\ +fi\n\ +unset __conda_setup\n\ +conda activate BindCraft\n\ +python -u /app/bindcraft.py "$@"\n' > /app/run_bindcraft.sh \ + && chmod +x /app/run_bindcraft.sh + +# Removed entry point to facilitate running via Galaxy bash scripts +ENTRYPOINT ["/app/run_bindcraft.sh"] \ No newline at end of file diff --git a/modules/local/bindcraft/main.nf b/modules/local/bindcraft/main.nf index e1e6cda..7abf03e 100644 --- a/modules/local/bindcraft/main.nf +++ b/modules/local/bindcraft/main.nf @@ -22,7 +22,7 @@ process BINDCRAFT { def args = task.ext.args ?: '' """ - /app/run_bindcraft.sh \\ + python /work/FreeBindCraft/bindcraft.py \\ --settings ${target_file} \\ --filters ${filters} \\ --advanced ${advanced_settings} \\ diff --git a/nextflow.config b/nextflow.config index 5d113e9..07c20db 100644 --- a/nextflow.config +++ b/nextflow.config @@ -17,7 +17,7 @@ params { settings_advanced = null batches = 1 quote_char = "\"" - bindcraft_container = null + bindcraft_container = "australianbiocommons/freebindcraft:1.0.3" // MultiQC options multiqc_config = null @@ -186,14 +186,14 @@ env { } // Set bash options -process.shell = """\ -bash - -set -e # Exit if a tool returns a non-zero status/exit code -set -u # Treat unset variables and parameters as an error -set -o pipefail # Returns the status of the last command to exit with a non-zero status or zero if all successfully execute -set -C # No clobber - prevent output redirection from overwriting files. -""" +process.shell = [ + "bash", + "-C", // No clobber - prevent output redirection from overwriting files. + "-e", // Exit if a tool returns a non-zero status/exit code + "-u", // Treat unset variables and parameters as an error + "-o", // Returns the status of the last command to exit.. + "pipefail" // ..with a non-zero status or zero if all successfully execute +] // Disable process selector warnings by default. Use debug profile to enable warnings. nextflow.enable.configProcessNamesValidation = false diff --git a/nextflow_schema.json b/nextflow_schema.json index 61c0652..eed34fb 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -52,8 +52,6 @@ "properties": { "bindcraft_container": { "type": "string", - "format": "file-path", - "exists": true, "description": "Path to bindcraft container to be used during execution.", "help_text": "", "fa_icon": "fas fa-file-csv" From 7ec06b4a56ec573856e29a8d8f9b38d1fee733f5 Mon Sep 17 00:00:00 2001 From: ziadbkh Date: Mon, 12 Jan 2026 11:46:27 +1100 Subject: [PATCH 3/7] update docs --- docs/usage.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/usage.md b/docs/usage.md index 4522c2c..58ce3bc 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -14,7 +14,7 @@ `outdir`: the output directory where the results will be stored. -`bindcraft_container`: the path to bindcraft container that is needed run the workflow. +`bindcraft_container`: the path to bindcraft container that is needed run the workflow. By default, teh workflow uses n open source version of Bindcraft ([FreeBindCraft](https://github.com/cytokineking/FreeBindCraft)) that is is hosted at Docker hub ([freebindcraft:1.0.3](https://hub.docker.com/r/australianbiocommons/freebindcraft)). ### Optional From 817bb7fe84636c327f01430d2c7a7bdc53b2bda0 Mon Sep 17 00:00:00 2001 From: ziadbkh Date: Wed, 1 Apr 2026 12:48:42 +1100 Subject: [PATCH 4/7] adding reporting --- assets/bindcraft_reporting.qmd | 1306 +++++++++++++++++++++++++++ modules/local/bindcraft/main.nf | 4 + modules/local/generate_report.nf | 48 +- modules/local/reporting/main.nf | 0 subworkflows/local/run_bindcraft.nf | 12 +- workflows/bindflow.nf | 1 + 6 files changed, 1359 insertions(+), 12 deletions(-) create mode 100644 assets/bindcraft_reporting.qmd create mode 100644 modules/local/reporting/main.nf diff --git a/assets/bindcraft_reporting.qmd b/assets/bindcraft_reporting.qmd new file mode 100644 index 0000000..d537f7c --- /dev/null +++ b/assets/bindcraft_reporting.qmd @@ -0,0 +1,1306 @@ +--- +title: "BindCraft Analysis Report" +author: "Australian-Protein-Design-Initiative/nf-binder-design" +date: last-modified +format: + html: + code-fold: true + toc: true + toc-depth: 3 + theme: cosmo + embed-resources: true + page-layout: full +execute: + echo: false +jupyter: python3 +--- + +# Overview + +This report presents a summary of the BindCraft run results, including failure metrics (`failure_csv.csv`), MPNN sequence design scoring (`mpnn_design_stats.csv`), and scores for all trajectories (`trajectory_stats.csv`) and final accepted design statistics (`final_design_stats.csv`). + +```{python} +# | label: setup +# | message: false +# | warning: false + +import os +import pandas as pd +import numpy as np +import matplotlib.pyplot as plt +from pathlib import Path +import warnings + +warnings.filterwarnings("ignore") + +# Set plotting style +plt.style.use("fast") + +boxplot_jitter = 0.05 + + +# Helper function to safely load CSV files +def safe_load_csv(filename, description): + try: + if Path(filename).exists(): + df = pd.read_csv(filename) + # print(f"✓ Loaded {description}: {len(df)} rows, {len(df.columns)} columns") + return df + else: + print(f"⚠ File not found: {filename}") + return None + except Exception as e: + print(f"✗ Error loading {filename}: {e}") + return None + + +print(f"Run path: {os.getcwd()}") + +# Load all available data files +failure_df = safe_load_csv("failure_csv.csv", "Failure metrics") +final_design_df = safe_load_csv("final_design_stats.csv", "Final design statistics") +mpnn_design_df = safe_load_csv("mpnn_design_stats.csv", "MPNN design statistics") +trajectory_df = safe_load_csv("trajectory_stats.csv", "Trajectory statistics") + +# print(f"\nData availability summary:") +# failure_status = "Available" if failure_df is not None else "Not available" +# final_status = "Available" if final_design_df is not None else "Not available" +# mpnn_status = "Available" if mpnn_design_df is not None else "Not available" +# trajectory_status = "Available" if trajectory_df is not None else "Not available" + +# print(f"- Failure metrics: {failure_status}") +# print(f"- Final design stats: {final_status}") +# print(f"- MPNN design stats: {mpnn_status}") +# print(f"- Trajectory stats: {trajectory_status}") +``` + +```{python} +# | label: acceptance-summary +# | message: false + +from IPython.display import HTML +from pathlib import Path +import glob + +# Calculate acceptance rate from CSV data +total_trajectories = len(trajectory_df) if trajectory_df is not None else 0 +accepted_count = len(final_design_df) if final_design_df is not None else 0 + + +# Count PDB files in different directories +def count_pdb_files(directory_pattern): + """Count .pdb/.pdb.gz files in directories matching the pattern""" + count = 0 + for path in Path(".").glob(directory_pattern): + if path.is_dir(): + # The glob pattern "*.{pdb,pdb.gz}" is not valid in Python's pathlib/glob. + # Instead, count both .pdb and .pdb.gz files separately: + count += len(list(path.glob("*.pdb"))) + count += len(list(path.glob("*.pdb.gz"))) + return count + + +# Count files in different categories +rejected_count = count_pdb_files("./batches/*/results/Rejected") +relaxed_count = count_pdb_files("./batches/*/results/Trajectory/Relaxed") +lowconf_count = count_pdb_files("./batches/*/results/Trajectory/LowConfidence") +clashing_count = count_pdb_files("./batches/*/results/Trajectory/Clashing") + +# Calculate total from directory counts +total_from_dirs = relaxed_count + lowconf_count + clashing_count + rejected_count + +# Use directory counts if available, otherwise fall back to CSV data +if total_from_dirs > 0: + total_trajectories = total_from_dirs + +if relaxed_count > 0: + acceptance_rate = (accepted_count / relaxed_count) * 100 + html_content = f""" +
+

+ Accept rate: {acceptance_rate:.1f}% ({accepted_count} / {relaxed_count}) +

+ + + + + + + + + + + + + + + + + + + + + +
Relaxed{relaxed_count}
Rejected{rejected_count}
LowConfidence{lowconf_count}
Clashing{clashing_count}
Total Trajectories{total_trajectories}
+
+ """ + display(HTML(html_content)) +else: + html_content = """ +
+

+ No trajectory data available +

+
+ """ + display(HTML(html_content)) +``` + +# Design Attrition Summary + +These graphs summarize the step in the BindCraft workflow or structural property filter that resulted in a design being rejected. + +```{python} +# | label: failure-metrics-trajectory +# | fig-width: 8 +# | fig-height: 6 + +if failure_df is not None and len(failure_df) > 0: + # Define column ranges for the three plots + trajectory_cols = ['Trajectory_logits_pLDDT', 'Trajectory_softmax_pLDDT', 'Trajectory_one-hot_pLDDT', + 'Trajectory_final_pLDDT', 'Trajectory_Contacts', 'Trajectory_Clashes', 'Trajectory_WrongHotspot', + 'MPNN_score', 'MPNN_seq_recovery', 'pLDDT', 'pTM', 'i_pTM', 'pAE', 'i_pAE', 'i_pLDDT', 'ss_pLDDT'] + + # Get available columns that exist in the dataframe + available_trajectory_cols = [col for col in trajectory_cols if col in failure_df.columns] + + if len(available_trajectory_cols) > 0: + # Get values for trajectory columns + if len(failure_df) > 1: + trajectory_values = failure_df[available_trajectory_cols].sum().values + else: + trajectory_values = failure_df.iloc[0][available_trajectory_cols].values + + # Create trajectory plot + fig, ax = plt.subplots(figsize=(8, 6)) + bars = ax.bar(range(len(available_trajectory_cols)), trajectory_values, alpha=0.7, color="steelblue") + + # Add value labels on bars for non-zero values + for i, (bar, val) in enumerate(zip(bars, trajectory_values)): + if val > 0: + ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, + str(int(val)), ha='center', va='bottom', fontsize=8, fontweight='bold') + + ax.set_xlabel("Trajectory Metrics", fontsize=12) + ax.set_ylabel("Number of Failures", fontsize=12) + ax.set_title("BindCraft Trajectory Failure Metrics", fontsize=14, fontweight="bold") + ax.set_xticks(range(len(available_trajectory_cols))) + ax.set_xticklabels(available_trajectory_cols, rotation=45, ha="right", fontsize=9) + ax.grid(axis="y", alpha=0.3) + ax.set_axisbelow(True) + plt.tight_layout() + plt.show() + + +``` + +```{python} +# | label: failure-metrics-structure +# | fig-width: 8 +# | fig-height: 6 + +if failure_df is not None and len(failure_df) > 0: + # Get all columns from Unrelaxed_Clashes to the last InterfaceAAs_ column + all_cols = list(failure_df.columns) + start_idx = all_cols.index('Unrelaxed_Clashes') if 'Unrelaxed_Clashes' in all_cols else 0 + + # Find the last InterfaceAAs_ column + interface_cols = [col for col in all_cols if col.startswith('InterfaceAAs_')] + if interface_cols: + last_interface_col = max(interface_cols, key=lambda x: all_cols.index(x)) + end_idx = all_cols.index(last_interface_col) + 1 + else: + end_idx = len(all_cols) + + structure_cols = all_cols[start_idx:end_idx] + + if len(structure_cols) > 0: + # Get values for structure columns + if len(failure_df) > 1: + structure_values = failure_df[structure_cols].sum().values + else: + structure_values = failure_df.iloc[0][structure_cols].values + + # Filter out InterfaceAAs_ columns unless they have non-zero values + filtered_structure_cols = [] + filtered_structure_values = [] + for col, val in zip(structure_cols, structure_values): + if not col.startswith('InterfaceAAs_') or val > 0: + filtered_structure_cols.append(col) + filtered_structure_values.append(val) + + if len(filtered_structure_cols) > 0: + # Create structure plot + fig, ax = plt.subplots(figsize=(8, 6)) + bars = ax.bar(range(len(filtered_structure_cols)), filtered_structure_values, alpha=0.7, color="coral") + + # Add value labels on bars for non-zero values + for i, (bar, val) in enumerate(zip(bars, filtered_structure_values)): + if val > 0: + ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, + str(int(val)), ha='center', va='bottom', fontsize=8, fontweight='bold') + + ax.set_xlabel("Structure Metrics", fontsize=12) + ax.set_ylabel("Number of Failures", fontsize=12) + ax.set_title("BindCraft Structure Failure Metrics", fontsize=14, fontweight="bold") + ax.set_xticks(range(len(filtered_structure_cols))) + ax.set_xticklabels(filtered_structure_cols, rotation=45, ha="right", fontsize=9) + ax.grid(axis="y", alpha=0.3) + ax.set_axisbelow(True) + plt.tight_layout() + plt.show() + + +``` + +```{python} +# | label: failure-metrics-final +# | fig-width: 8 +# | fig-height: 6 + +if failure_df is not None and len(failure_df) > 0: + # Get columns from Hotspot_RMSD to the end + all_cols = list(failure_df.columns) + start_idx = all_cols.index('Hotspot_RMSD') if 'Hotspot_RMSD' in all_cols else len(all_cols) - 4 + final_cols = all_cols[start_idx:] + + if len(final_cols) > 0: + # Get values for final columns + if len(failure_df) > 1: + final_values = failure_df[final_cols].sum().values + else: + final_values = failure_df.iloc[0][final_cols].values + + # Create final plot + fig, ax = plt.subplots(figsize=(8, 6)) + bars = ax.bar(range(len(final_cols)), final_values, alpha=0.7, color="lightgreen") + + # Add value labels on bars for non-zero values + for i, (bar, val) in enumerate(zip(bars, final_values)): + if val > 0: + ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, + str(int(val)), ha='center', va='bottom', fontsize=8, fontweight='bold') + + ax.set_xlabel("Final Metrics", fontsize=12) + ax.set_ylabel("Number of Failures", fontsize=12) + ax.set_title("BindCraft Final Failure Metrics", fontsize=14, fontweight="bold") + ax.set_xticks(range(len(final_cols))) + ax.set_xticklabels(final_cols, rotation=45, ha="right", fontsize=9) + ax.grid(axis="y", alpha=0.3) + ax.set_axisbelow(True) + plt.tight_layout() + plt.show() + + +``` + +# Trajectory Summary + +Summary of all the design trajectories, accepted and rejected. + +```{python} +#| label: trajectory-boxplots +#| fig-width: 12 +#| fig-height: 8 + + if trajectory_df is not None and len(trajectory_df) > 0: + # Key trajectory metrics for box plots (excluding dG and Target_RMSD) + trajectory_metrics = ['pLDDT', 'pTM', 'i_pTM', 'pAE', 'i_pAE', 'ShapeComplementarity'] + available_traj_metrics = [col for col in trajectory_metrics if col in trajectory_df.columns] + + if available_traj_metrics: + # Create single box plot with all metrics (excluding dG and Target_RMSD) + fig, ax = plt.subplots(figsize=(8, 6)) + + # Prepare data for box plot + plot_data = [] + labels = [] + + for metric in available_traj_metrics: + data = trajectory_df[metric].dropna() + if len(data) > 0: + plot_data.append(data.values) + labels.append(metric) + + if plot_data: + # Create box plot + bp = ax.boxplot(plot_data, patch_artist=True, labels=labels) + + # Color the boxes + for patch in bp['boxes']: + patch.set_facecolor('lightgreen') + for median in bp['medians']: + median.set_color('red') + + # Add individual points with jitter + for i, data in enumerate(plot_data): + x_points = np.ones(len(data)) * (i + 1) # Position points at box locations + # Add jitter to x-coordinates to avoid overlap + jitter = np.random.normal(0, boxplot_jitter, len(data)) + x_points_jittered = x_points + jitter + ax.plot(x_points_jittered, data, 'o', color='darkgreen', alpha=0.6, markersize=4) + + ax.set_xlabel('Trajectory Metrics', fontsize=12) + ax.set_ylabel('Value', fontsize=12) + ax.set_title('Trajectory Analysis Metrics Distribution', fontsize=14, fontweight='bold') + ax.tick_params(axis='x', rotation=45) + ax.set_xticklabels(labels, rotation=45, ha='right') + ax.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + + # Separate box plot for Target_RMSD + if 'Target_RMSD' in trajectory_df.columns: + rmsd_data = trajectory_df['Target_RMSD'].dropna() + if len(rmsd_data) > 0: + fig, ax = plt.subplots(figsize=(8, 6)) + + # Create box plot for Target_RMSD + bp = ax.boxplot(rmsd_data.values, patch_artist=True, labels=['Target_RMSD']) + + # Color the box + bp['boxes'][0].set_facecolor('lightblue') + bp['medians'][0].set_color('red') + + # Add individual points with jitter + x_points = np.ones(len(rmsd_data)) + jitter = np.random.normal(0, boxplot_jitter, len(rmsd_data)) + x_points_jittered = x_points + jitter + ax.plot(x_points_jittered, rmsd_data.values, 'o', color='darkblue', alpha=0.6, markersize=4) + + ax.set_xlabel('Metric', fontsize=12) + ax.set_ylabel('Target RMSD Value', fontsize=12) + ax.set_title('Target RMSD Distribution', fontsize=14, fontweight='bold') + ax.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + + # Separate box plot for dG + if 'dG' in trajectory_df.columns: + dg_data = trajectory_df['dG'].dropna() + if len(dg_data) > 0: + fig, ax = plt.subplots(figsize=(8, 6)) + + # Create box plot for dG + bp = ax.boxplot(dg_data.values, patch_artist=True, labels=['dG']) + + # Color the box + bp['boxes'][0].set_facecolor('lightcoral') + bp['medians'][0].set_color('red') + + # Add individual points with jitter + x_points = np.ones(len(dg_data)) + jitter = np.random.normal(0, boxplot_jitter, len(dg_data)) + x_points_jittered = x_points + jitter + ax.plot(x_points_jittered, dg_data.values, 'o', color='darkred', alpha=0.6, markersize=4) + + ax.set_xlabel('Metric', fontsize=12) + ax.set_ylabel('dG Value', fontsize=12) + ax.set_title('dG Distribution', fontsize=14, fontweight='bold') + ax.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() +``` + +```{python} +# | label: trajectory-ranking +# | message: false + +if trajectory_df is not None and len(trajectory_df) > 0: + # Define columns to show in the table + table_cols = [ + "Design", + "i_pTM", + "i_pAE", + "pLDDT", + "pTM", + "dG", + "ShapeComplementarity", + "Target_RMSD", + "Length", + "pAE", + "i_pLDDT", + "ss_pLDDT", + "Unrelaxed_Clashes", + "Relaxed_Clashes", + "Binder_Energy_Score", + "Surface_Hydrophobicity", + "PackStat", + "dSASA", + "dG/dSASA", + "Interface_SASA_%", + "Interface_Hydrophobicity", + "n_InterfaceResidues", + "n_InterfaceHbonds", + "InterfaceHbondsPercentage", + "n_InterfaceUnsatHbonds", + "InterfaceUnsatHbondsPercentage", + "Interface_Helix%", + "Interface_BetaSheet%", + "Interface_Loop%", + "Binder_Helix%", + "Binder_BetaSheet%", + "Binder_Loop%", + ] + available_table_cols = [col for col in table_cols if col in trajectory_df.columns] + + if available_table_cols and "Design" in trajectory_df.columns: + # Create ranking table sorted by i_pTM + ranking_df = trajectory_df[available_table_cols].copy() + ranking_df = ranking_df.sort_values("i_pTM", ascending=False) + + # Display top 20 designs (or less if fewer available) + print("\nTop 20 Trajectory Designs (ranked by i_pTM):") + display_df = ranking_df.head(20) + + # Format numeric columns - Length as integer, others as float + numeric_cols = [col for col in available_table_cols if col != "Design"] + format_dict = {} + for col in numeric_cols: + if col == "Length": + format_dict[col] = "{:.0f}" # Integer format + else: + format_dict[col] = "{:.3f}" # Float format + + styled_table = ( + display_df.style.format(format_dict) + .background_gradient(cmap="RdYlGn", subset=["i_pTM"]) + .hide(axis="index") + ) + display(styled_table) +``` + + +# MPNN Scores + +Summary of scores for designs after sequence generation with ProteinMPNN. + +```{python} +#| label: mpnn-design-boxplots +#| fig-width: 12 +#| fig-height: 8 + + if mpnn_design_df is not None and len(mpnn_design_df) > 0: + # Key metrics for box plots (excluding dG metrics) + key_metrics = ['MPNN_score', 'MPNN_seq_recovery', 'Average_pLDDT', 'Average_pTM', + 'Average_i_pTM', 'Average_i_pAE', 'Average_ShapeComplementarity'] + available_metrics = [col for col in key_metrics if col in mpnn_design_df.columns] + + if available_metrics: + # Create single box plot with all metrics (excluding dG) + fig, ax = plt.subplots(figsize=(8, 6)) + + # Prepare data for box plot + plot_data = [] + labels = [] + + for metric in available_metrics: + data = mpnn_design_df[metric].dropna() + if len(data) > 0: + plot_data.append(data.values) + labels.append(metric) + + if plot_data: + # Create box plot + bp = ax.boxplot(plot_data, patch_artist=True, labels=labels) + + # Color the boxes + for patch in bp['boxes']: + patch.set_facecolor('lightblue') + for median in bp['medians']: + median.set_color('red') + + # Add individual points + for i, data in enumerate(plot_data): + x_points = np.ones(len(data)) * (i + 1) # Position points at box locations + ax.plot(x_points, data, 'o', color='darkblue', alpha=0.6, markersize=4) + + ax.set_xlabel('MPNN Metrics', fontsize=12) + ax.set_ylabel('Value', fontsize=12) + ax.set_title('MPNN Design Metrics Distribution', fontsize=14, fontweight='bold') + ax.tick_params(axis='x', rotation=45) + ax.set_xticklabels(labels, rotation=45, ha='right') + ax.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() + + # Separate box plot for Average_dG + if 'Average_dG' in mpnn_design_df.columns: + dg_data = mpnn_design_df['Average_dG'].dropna() + if len(dg_data) > 0: + fig, ax = plt.subplots(figsize=(8, 6)) + + # Create box plot for dG + bp = ax.boxplot(dg_data.values, patch_artist=True, labels=['Average_dG']) + + # Color the box + bp['boxes'][0].set_facecolor('lightcoral') + bp['medians'][0].set_color('red') + + # Add individual points + ax.plot(np.ones(len(dg_data)), dg_data.values, 'o', color='darkred', alpha=0.6, markersize=4) + + ax.set_xlabel('Metric', fontsize=12) + ax.set_ylabel('dG Value', fontsize=12) + ax.set_title('Average dG Distribution', fontsize=14, fontweight='bold') + ax.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() +``` + +```{python} +#| label: mpnn-design-ranking +#| message: false + +if mpnn_design_df is not None and len(mpnn_design_df) > 0: + # Define columns to show in the table + table_cols = ['Design', 'Average_i_pTM', 'Average_i_pAE', 'Average_pLDDT', 'Average_pTM', + 'Average_dG', 'Average_ShapeComplementarity', 'Length'] + available_table_cols = [col for col in table_cols if col in mpnn_design_df.columns] + + if available_table_cols and 'Design' in mpnn_design_df.columns: + # Create ranking table sorted by Average_i_pTM + ranking_df = mpnn_design_df[available_table_cols].copy() + ranking_df = ranking_df.sort_values('Average_i_pTM', ascending=False) + + # Display top 20 designs + print("\nTop MPNN Designs (ranked by Average_i_pTM):") + display_df = ranking_df.head(20) + + # Format numeric columns - Length as integer, others as float + numeric_cols = [col for col in available_table_cols if col != 'Design'] + format_dict = {} + for col in numeric_cols: + if col == 'Length': + format_dict[col] = '{:.0f}' # Integer format + else: + format_dict[col] = '{:.3f}' # Float format + + styled_table = display_df.style.format(format_dict) \ + .background_gradient(cmap='RdYlGn', subset=['Average_i_pTM']) \ + .hide(axis='index') + display(styled_table) +``` + + +# Accepted Design Statistics + +Summary of the final accepted designs (if available). + +```{python} +# | label: accepted-designs-check +# | message: false + +if final_design_df is None or len(final_design_df) == 0: + from IPython.display import HTML + html_content = """ +
+

+ No accepted designs +

+
+ """ + display(HTML(html_content)) +``` + + + +```{python} +#| label: final-design-main-boxplots +#| fig-width: 8 +#| fig-height: 6 + +if final_design_df is not None and len(final_design_df) > 0: + # Key final design metrics for box plots (excluding dG) + final_metrics = ['Average_pLDDT', 'Average_pTM', 'Average_i_pTM', 'Average_i_pAE', 'Average_ShapeComplementarity', + 'Average_ss_pLDDT', 'Average_Unrelaxed_Clashes', 'Average_Relaxed_Clashes', 'Average_Binder_Energy_Score', + 'Average_Surface_Hydrophobicity', 'Average_PackStat', 'Average_Interface_SASA_%', 'Average_Interface_Hydrophobicity', + 'Average_n_InterfaceResidues', 'Average_n_InterfaceHbonds', 'Average_InterfaceHbondsPercentage', + 'Average_n_InterfaceUnsatHbonds', 'Average_InterfaceUnsatHbondsPercentage', 'Average_Interface_Helix%', + 'Average_Interface_BetaSheet%', 'Average_Interface_Loop%', 'Average_Binder_Helix%', 'Average_Binder_BetaSheet%', + 'Average_Binder_Loop%', 'Average_Hotspot_RMSD', 'Average_Target_RMSD', 'Average_Binder_pLDDT', 'Average_Binder_pTM', + 'Average_Binder_pAE', 'Average_Binder_RMSD'] + available_final_metrics = [col for col in final_metrics if col in final_design_df.columns] + + if available_final_metrics: + # Main metrics (excluding dG, Binder_Energy_Score, Interface_Hydrophobicity, n_* metrics, Clashes, and RMSD) + main_metrics = [metric for metric in available_final_metrics + if metric not in ['Average_Binder_Energy_Score', 'Average_dG', 'Average_Interface_Hydrophobicity'] + and not metric.startswith('Average_n_') + and 'Clashes' not in metric and 'RMSD' not in metric] + + # Create main box plot + if main_metrics: + fig, ax = plt.subplots(figsize=(8, 6)) + plot_data = [] + labels = [] + + for metric in main_metrics: + data = final_design_df[metric].dropna() + if len(data) > 0: + # Convert percentage columns to proportions (0-1 scale) + if any(percent_suffix in metric for percent_suffix in ['%', 'Percentage']): + data = data / 100.0 + # Update label to indicate it's now a proportion + label = metric.replace('Percentage', 'Prop').replace('%', '_prop') + else: + label = metric + + plot_data.append(data.values) + labels.append(label) + + if plot_data: + # Create box plot + bp = ax.boxplot(plot_data, patch_artist=True, labels=labels) + + # Color the boxes + for patch in bp['boxes']: + patch.set_facecolor('lightsteelblue') + for median in bp['medians']: + median.set_color('red') + + # Add individual points + for i, data in enumerate(plot_data): + x_points = np.ones(len(data)) * (i + 1) # Position points at box locations + ax.plot(x_points, data, 'o', color='darkblue', alpha=0.6, markersize=4) + + ax.set_xlabel('Final Design Metrics', fontsize=12) + ax.set_ylabel('Value', fontsize=12) + ax.set_title('Final Design Metrics Distribution', fontsize=14, fontweight='bold') + ax.tick_params(axis='x', rotation=45) + ax.set_xticklabels(labels, rotation=45, ha='right') + ax.grid(True, alpha=0.3) + plt.tight_layout() + plt.show() +``` + +```{python} +# | label: final-design-tiled-boxplots +# | fig-width: 3 +# | fig-height: 6 + +if final_design_df is not None and len(final_design_df) > 0: + # Tiled box plots for Interface Hydrophobicity and n_* metrics + tiled_metrics = [] + if "Average_Interface_Hydrophobicity" in final_design_df.columns: + tiled_metrics.append("Average_Interface_Hydrophobicity") + tiled_metrics.extend( + [col for col in available_final_metrics if col.startswith("Average_n_")] + ) + + if tiled_metrics: + # Calculate grid dimensions + n_plots = len(tiled_metrics) + cols = min(4, n_plots) # Max 4 columns + rows = (n_plots + cols - 1) // cols # Ceiling division + + fig, axes = plt.subplots( + rows, cols, figsize=(2 * cols, 2 * rows) + ) # Same height as standalone + if n_plots == 1: + axes = [axes] + elif rows == 1: + axes = axes.reshape(1, -1) + elif cols == 1: + axes = axes.reshape(-1, 1) + + for i, metric in enumerate(tiled_metrics): + row = i // cols + col = i % cols + ax = axes[row, col] + + data = final_design_df[metric].dropna() + if len(data) > 0: + # Create box plot + bp = ax.boxplot(data.values, patch_artist=True, labels=[metric]) + + # Color the box + bp["boxes"][0].set_facecolor("lightcoral") + bp["medians"][0].set_color("red") + + # Add individual points + x_points = np.ones( + len(data) + ) # Create array of ones matching data length + ax.plot( + x_points, data.values, "o", color="darkred", alpha=0.6, markersize=4 + ) + + ax.grid(True, alpha=0.3) + ax.tick_params(axis="x", rotation=45) + + # Hide empty subplots + for i in range(n_plots, rows * cols): + row = i // cols + col = i % cols + axes[row, col].set_visible(False) + + plt.tight_layout() + plt.show() +``` + +```{python} +# | label: final-design-standalone-boxplots +# | fig-width: 3 +# | fig-height: 6 + +if final_design_df is not None and len(final_design_df) > 0: + # Side-by-side standalone box plots + standalone_plots = [] + + # Collect all standalone plots + if "Average_Binder_Energy_Score" in final_design_df.columns: + standalone_plots.append(("Average_Binder_Energy_Score", "")) + if "Average_dG" in final_design_df.columns: + standalone_plots.append(("Average_dG", "")) + + # Add Clashes metrics + clashes_metrics = [col for col in available_final_metrics if "Clashes" in col] + if clashes_metrics: + standalone_plots.append(("clashes_group", "")) + + # Add RMSD metrics + rmsd_metrics = [col for col in available_final_metrics if "RMSD" in col] + if rmsd_metrics: + standalone_plots.append(("rmsd_group", "")) + + if standalone_plots: + n_standalone = len(standalone_plots) + fig, axes = plt.subplots( + 1, n_standalone, figsize=(2 * n_standalone, 6) + ) # Even narrower width + if n_standalone == 1: + axes = [axes] + + for i, (plot_type, title) in enumerate(standalone_plots): + ax = axes[i] + + if plot_type == "clashes_group": + # Create box plot for all Clashes metrics + plot_data = [] + labels = [] + for metric in clashes_metrics: + data = final_design_df[metric].dropna() + if len(data) > 0: + plot_data.append(data.values) + labels.append(metric) + + if plot_data: + bp = ax.boxplot(plot_data, patch_artist=True, labels=labels) + for patch in bp["boxes"]: + patch.set_facecolor("lightcoral") + for median in bp["medians"]: + median.set_color("red") + + # Add individual points + for j, data in enumerate(plot_data): + x_points = np.ones(len(data)) * (j + 1) + ax.plot( + x_points, + data, + "o", + color="darkred", + alpha=0.6, + markersize=4, + ) + + ax.set_xticklabels(labels, rotation=45, ha="right") + + elif plot_type == "rmsd_group": + # Create box plot for all RMSD metrics + plot_data = [] + labels = [] + for metric in rmsd_metrics: + data = final_design_df[metric].dropna() + if len(data) > 0: + plot_data.append(data.values) + labels.append(metric) + + if plot_data: + bp = ax.boxplot(plot_data, patch_artist=True, labels=labels) + for patch in bp["boxes"]: + patch.set_facecolor("lightcoral") + for median in bp["medians"]: + median.set_color("red") + + # Add individual points + for j, data in enumerate(plot_data): + x_points = np.ones(len(data)) * (j + 1) + ax.plot( + x_points, + data, + "o", + color="darkred", + alpha=0.6, + markersize=4, + ) + + ax.set_xticklabels(labels, rotation=45, ha="right") + + else: + # Single metric box plot + data = final_design_df[plot_type].dropna() + if len(data) > 0: + bp = ax.boxplot(data.values, patch_artist=True, labels=[plot_type]) + bp["boxes"][0].set_facecolor("lightcoral") + bp["medians"][0].set_color("red") + + x_points = np.ones(len(data)) + ax.plot( + x_points, + data.values, + "o", + color="darkred", + alpha=0.6, + markersize=4, + ) + + ax.grid(True, alpha=0.3) + ax.tick_params(axis="x", rotation=45) + + plt.tight_layout() + plt.show() +``` + +```{python} +# | label: final-design-ranking +# | message: false + +if final_design_df is not None and len(final_design_df) > 0: + # Define columns to show in the table + table_cols = [ + "Design", + "Average_i_pTM", + "Average_i_pAE", + "Average_pLDDT", + "Average_pTM", + "Average_dG", + "Average_ShapeComplementarity", + "Average_Target_RMSD", + "Length", + "Average_ss_pLDDT", + "Average_Unrelaxed_Clashes", + "Average_Relaxed_Clashes", + "Average_Binder_Energy_Score", + "Average_Surface_Hydrophobicity", + "Average_PackStat", + "Average_Interface_SASA_%", + "Average_Interface_Hydrophobicity", + "Average_n_InterfaceResidues", + "Average_n_InterfaceHbonds", + "Average_InterfaceHbondsPercentage", + "Average_n_InterfaceUnsatHbonds", + "Average_InterfaceUnsatHbondsPercentage", + "Average_Interface_Helix%", + "Average_Interface_BetaSheet%", + "Average_Interface_Loop%", + "Average_Binder_Helix%", + "Average_Binder_BetaSheet%", + "Average_Binder_Loop%", + "Average_Hotspot_RMSD", + "Average_Binder_pLDDT", + "Average_Binder_pTM", + "Average_Binder_pAE", + "Average_Binder_RMSD", + ] + available_table_cols = [col for col in table_cols if col in final_design_df.columns] + + if available_table_cols and "Design" in final_design_df.columns: + # Create ranking table sorted by Average_i_pTM + ranking_df = final_design_df[available_table_cols].copy() + ranking_df = ranking_df.sort_values("Average_i_pTM", ascending=False) + + # Display all accepted designs (or less if fewer available) + print("\nFinal Accepted Designs (ranked by Average_i_pTM):") + display_df = ranking_df + + # Format numeric columns - Length as integer, others as float + numeric_cols = [col for col in available_table_cols if col != "Design"] + format_dict = {} + for col in numeric_cols: + if col == "Length": + format_dict[col] = "{:.0f}" # Integer format + else: + format_dict[col] = "{:.3f}" # Float format + + styled_table = ( + display_df.style.format(format_dict) + .background_gradient(cmap="RdYlGn", subset=["Average_i_pTM"]) + .hide(axis="index") + ) + display(styled_table) +``` + +# Hotspot Analysis + +Analysis of design scores across different target hotspot selections. + +```{python} +# | label: final-design-hotspot-heatmap +# | fig-width: 10 +# | fig-height: 8 + +if ( + final_design_df is not None + and len(final_design_df) > 0 + and "Target_Hotspot" in final_design_df.columns + and "Average_i_pTM" in final_design_df.columns +): + # Parse Target_Hotspot column + def parse_hotspots(hotspot_str): + """Parse hotspot string, handling both single values and quoted lists""" + if pd.isna(hotspot_str): + return [] + + # Remove quotes if present + hotspot_str = str(hotspot_str).strip().strip("\"'") + + # Check if it's a comma-separated list + if "," in hotspot_str: + return [h.strip() for h in hotspot_str.split(",")] + else: + return [hotspot_str] + + # Extract unique hotspots + all_hotspots = [] + for hotspot_str in final_design_df["Target_Hotspot"]: + hotspots = parse_hotspots(hotspot_str) + all_hotspots.extend(hotspots) + + unique_hotspots = sorted(list(set(all_hotspots))) + + # Check if all designs have identical hotspots + if len(unique_hotspots) <= 1: + print("All designs have identical target hotspots. Skipping hotspot heatmap.") + else: + # Create pivot table for heatmap + heatmap_data = [] + design_names = [] + + for idx, row in final_design_df.iterrows(): + hotspots = parse_hotspots(row["Target_Hotspot"]) + i_pTM = row["Average_i_pTM"] + + # Create a row for each hotspot in this design + for hotspot in hotspots: + if hotspot in unique_hotspots: + heatmap_data.append([hotspot, row["Design"], i_pTM]) + design_names.append(row["Design"]) + + if heatmap_data: + # Create DataFrame for heatmap + heatmap_df = pd.DataFrame( + heatmap_data, columns=["Hotspot", "Design", "i_pTM"] + ) + pivot_df = heatmap_df.pivot( + index="Design", columns="Hotspot", values="i_pTM" + ) + + # Fill NaN values with 0 or appropriate value + pivot_df = pivot_df.fillna(0) + + # Perform hierarchical clustering + from scipy.cluster.hierarchy import dendrogram, linkage + from scipy.spatial.distance import pdist + import seaborn as sns + + # Prepare data for clustering + X = pivot_df.values + + # Sort rows by average i_pTM from highest to lowest + if len(pivot_df) > 1: + # Calculate average i_pTM for each design + design_avg_iptm = pivot_df.mean(axis=1) + # Sort by average i_pTM (highest to lowest) + row_order = design_avg_iptm.sort_values(ascending=False).index + pivot_df_clustered = pivot_df.loc[row_order] + else: + pivot_df_clustered = pivot_df + + # Sort columns by letter prefix, then by number + def sort_hotspot_key(hotspot): + """Sort hotspots by letter prefix, then by number""" + import re + + # Extract letter prefix and number + match = re.match(r"([A-Za-z]+)(\d+)", hotspot) + if match: + letter_prefix = match.group(1) + number = int(match.group(2)) + return (letter_prefix, number) + else: + # Fallback for non-standard format + return (hotspot, 0) + + if len(pivot_df.columns) > 1: + # Sort columns using the custom key function + col_order = sorted(pivot_df.columns, key=sort_hotspot_key) + pivot_df_clustered = pivot_df_clustered.loc[:, col_order] + + # Create heatmap with dendrogram + from scipy.cluster.hierarchy import dendrogram, linkage + from scipy.spatial.distance import pdist + import seaborn as sns + + # Calculate height based on number of designs + height_per_design = 0.3 # inches per design + min_height = 6 # minimum height in inches + max_height = 20 # maximum height in inches + calculated_height = max( + min_height, min(max_height, len(pivot_df_clustered) * height_per_design) + ) + + # Create figure with proper height (narrower width) + fig = plt.figure(figsize=(8, calculated_height)) + + # Single subplot for heatmap + ax_heatmap = fig.add_subplot(111) + + # Create custom colormap: white to red + import matplotlib.colors as mcolors + from matplotlib.colors import LinearSegmentedColormap + + # Create custom colormap: white -> red + colors = ["white", "red"] + n_bins = 100 + custom_cmap = LinearSegmentedColormap.from_list( + "custom_white_red", colors, N=n_bins + ) + + # Create heatmap + sns.heatmap( + pivot_df_clustered, + annot=False, # Remove cell annotations + cmap=custom_cmap, + center=0.5, + cbar_kws={"label": "Average i_pTM"}, + xticklabels=True, + yticklabels=True, + ax=ax_heatmap, + ) + + # Customize x-axis labels - show on both top and bottom + ax_heatmap.set_xticklabels( + ax_heatmap.get_xticklabels(), rotation=45, ha="center", fontsize=8 + ) + ax_heatmap.set_yticklabels(ax_heatmap.get_yticklabels(), fontsize=8) + + # Add x-axis labels on top as well + ax_heatmap.xaxis.set_tick_params(labeltop=True, labelbottom=True) + + # Adjust layout to prevent title overlap + plt.subplots_adjust(top=0.9) # Leave more space for title + + plt.suptitle( + "Accepted Design ipTM by Target Hotspot", + fontsize=14, + fontweight="bold", + y=1.0, + ) + plt.tight_layout() + plt.show() +``` + +```{python} +# | label: trajectory-hotspot-heatmap +# | fig-width: 10 +# | fig-height: 10 + +if ( + trajectory_df is not None + and len(trajectory_df) > 0 + and "Target_Hotspot" in trajectory_df.columns + and "i_pTM" in trajectory_df.columns +): + # Parse Target_Hotspot column + def parse_hotspots(hotspot_str): + """Parse hotspot string, handling both single values and quoted lists""" + if pd.isna(hotspot_str): + return [] + + # Remove quotes if present + hotspot_str = str(hotspot_str).strip().strip("\"'") + + # Check if it's a comma-separated list + if "," in hotspot_str: + return [h.strip() for h in hotspot_str.split(",")] + else: + return [hotspot_str] + + # Extract unique hotspots + all_hotspots = [] + for hotspot_str in trajectory_df["Target_Hotspot"]: + hotspots = parse_hotspots(hotspot_str) + all_hotspots.extend(hotspots) + + unique_hotspots = sorted(list(set(all_hotspots))) + + # Check if all designs have identical hotspots + if len(unique_hotspots) <= 1: + print( + "All trajectories have identical target hotspots. Skipping trajectory hotspot heatmap." + ) + else: + # Create pivot table for heatmap + heatmap_data = [] + design_names = [] + + for idx, row in trajectory_df.iterrows(): + hotspots = parse_hotspots(row["Target_Hotspot"]) + i_pTM = row["i_pTM"] + + # Create a row for each hotspot in this design + for hotspot in hotspots: + if hotspot in unique_hotspots: + heatmap_data.append([hotspot, row["Design"], i_pTM]) + design_names.append(row["Design"]) + + if heatmap_data: + # Create DataFrame for heatmap + heatmap_df = pd.DataFrame( + heatmap_data, columns=["Hotspot", "Design", "i_pTM"] + ) + pivot_df = heatmap_df.pivot( + index="Design", columns="Hotspot", values="i_pTM" + ) + + # Fill NaN values with 0 or appropriate value + pivot_df = pivot_df.fillna(0) + + # Perform hierarchical clustering + from scipy.cluster.hierarchy import dendrogram, linkage + from scipy.spatial.distance import pdist + import seaborn as sns + + # Prepare data for clustering + X = pivot_df.values + + # Sort rows by average i_pTM from highest to lowest + if len(pivot_df) > 1: + # Calculate average i_pTM for each design + design_avg_iptm = pivot_df.mean(axis=1) + # Sort by average i_pTM (highest to lowest) + row_order = design_avg_iptm.sort_values(ascending=False).index + pivot_df_clustered = pivot_df.loc[row_order] + else: + pivot_df_clustered = pivot_df + + # Sort columns by letter prefix, then by number + def sort_hotspot_key(hotspot): + """Sort hotspots by letter prefix, then by number""" + import re + + # Extract letter prefix and number + match = re.match(r"([A-Za-z]+)(\d+)", hotspot) + if match: + letter_prefix = match.group(1) + number = int(match.group(2)) + return (letter_prefix, number) + else: + # Fallback for non-standard format + return (hotspot, 0) + + if len(pivot_df.columns) > 1: + # Sort columns using the custom key function + col_order = sorted(pivot_df.columns, key=sort_hotspot_key) + pivot_df_clustered = pivot_df_clustered.loc[:, col_order] + + # Create heatmap with dendrogram + from scipy.cluster.hierarchy import dendrogram, linkage + from scipy.spatial.distance import pdist + import seaborn as sns + + # Calculate height based on number of designs + height_per_design = 0.3 # inches per design + min_height = 6 # minimum height in inches + max_height = 20 # maximum height in inches + calculated_height = max( + min_height, min(max_height, len(pivot_df_clustered) * height_per_design) + ) + + # Create figure with proper height (narrower width) + fig = plt.figure(figsize=(10, calculated_height)) + + # Single subplot for heatmap + ax_heatmap = fig.add_subplot(111) + + # Create custom colormap: white to red + import matplotlib.colors as mcolors + from matplotlib.colors import LinearSegmentedColormap + + # Create custom colormap: white -> red + colors = ["white", "red"] + n_bins = 100 + custom_cmap = LinearSegmentedColormap.from_list( + "custom_white_red", colors, N=n_bins + ) + + # Create heatmap + sns.heatmap( + pivot_df_clustered, + annot=False, # Remove cell annotations + cmap=custom_cmap, + center=0.5, + cbar_kws={"label": "i_pTM"}, + xticklabels=True, + yticklabels=True, + ax=ax_heatmap, + ) + + # Customize x-axis labels - show on both top and bottom + ax_heatmap.set_xticklabels( + ax_heatmap.get_xticklabels(), rotation=90, ha="center", fontsize=8 + ) + ax_heatmap.set_yticklabels(ax_heatmap.get_yticklabels(), fontsize=8) + + # Add x-axis labels on top as well + ax_heatmap.xaxis.set_tick_params(labeltop=True, labelbottom=True) + + # Adjust layout to prevent title overlap + plt.subplots_adjust(top=0.9) # Leave more space for title + + plt.suptitle( + "All Trajectory ipTM scores by Target Hotspot", + fontsize=14, + fontweight="bold", + y=1.0, + ) + plt.tight_layout() + plt.show() +``` diff --git a/modules/local/bindcraft/main.nf b/modules/local/bindcraft/main.nf index 7abf03e..62e687f 100644 --- a/modules/local/bindcraft/main.nf +++ b/modules/local/bindcraft/main.nf @@ -12,6 +12,10 @@ process BINDCRAFT { tuple val(meta), path("*_output/Accepted/Ranked"), emit: accepted_ranked tuple val(meta), path("*_output/Accepted/*pdb"), emit: accepted tuple val(meta), path("*_output"), emit: output_dir + tuple val(meta), path("*_output/failure_csv.csv") , emit: failure_csv + tuple val(meta), path("*_output/final_design_stats.csv"), emit: final_design_stats + tuple val(meta), path("*_output/mpnn_design_stats.csv"), emit: mpnn_design_stats + tuple val(meta), path("*_output/trajectory_stats.csv"), emit: trajectory_stats path "versions.yml", emit: versions when: diff --git a/modules/local/generate_report.nf b/modules/local/generate_report.nf index 83eef42..3395357 100644 --- a/modules/local/generate_report.nf +++ b/modules/local/generate_report.nf @@ -1,23 +1,51 @@ process GENERATE_REPORT { - tag "${meta.id}" + tag "all-run" label 'process_single' - container "${ workflow.containerEngine == 'singularity' && !task.ext.singularity_pull_docker_container ? - 'https://depot.galaxyproject.org/singularity/multiqc:1.21--pyhdfd78af_0' : - 'biocontainers/multiqc:1.21--pyhdfd78af_0' }" - conda "bioconda::multiqc=1.21" + container 'ghcr.io/australian-protein-design-initiative/containers/nf-binder-design-utils:0.1.5' input: - tuple val(meta), val(info) - + // Stage all batch result directories under ./batches/{n}/ + // Note the {n} in this case is a list index and may not match the batch_id + // For the purposes of generating aggregate stats, the batch_id is not important + path ('batches/*') + path('failure_csv.csv') + path('final_design_stats.csv') + path('mpnn_design_stats.csv') + path('trajectory_stats.csv') + path (qmd_template) + output: - tuple val(meta), path ("*report.html"), emit: report - + path('bindcraft_report.html'), emit: report + when: task.ext.when == null || task.ext.when + script: def args = task.ext.args ?: '' + """ + export XDG_CACHE_HOME="./.cache" + export XDG_DATA_HOME="./.local/share" + export JUPYTER_RUNTIME_DIR="./.jupyter" + export XDG_RUNTIME_DIR="/tmp" + + quarto render ${qmd_template} \\ + --execute-dir \${PWD} \\ + --output - > bindcraft_report.html + cat <<-END_VERSIONS > versions.yml + "${task.process}": + quarto: \$(quarto --version) + END_VERSIONS + """ + + stub: """ - echo "${info}" > ${meta.id}_report.html + + touch bindcraft_report.html + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + quarto: \$(quarto --version) + END_VERSIONS """ } diff --git a/modules/local/reporting/main.nf b/modules/local/reporting/main.nf new file mode 100644 index 0000000..e69de29 diff --git a/subworkflows/local/run_bindcraft.nf b/subworkflows/local/run_bindcraft.nf index d0c0ff8..bf08e93 100644 --- a/subworkflows/local/run_bindcraft.nf +++ b/subworkflows/local/run_bindcraft.nf @@ -86,10 +86,18 @@ workflow RUN_BINDCRAFT { ) GENERATE_REPORT( - RANKER.out.stats.map{[it[0], "Testing report information"]} + BINDCRAFT.out.output_dir.map{it[1]}.collect(), + BINDCRAFT.out.failure_csv.map{it[1]}, + RANKER.out.stats.map{it[1]}, + BINDCRAFT.out.mpnn_design_stats + .map{it[1].text} + .collectFile( name: "mpnn_design_stats.csv" ), + BINDCRAFT.out.trajectory_stats + .map{it[1].text} + .collectFile( name: "trajectory_stats.csv" ), + Channel.fromPath("${projectDir}/assets/bindcraft_reporting.qmd").first() ) - emit: input_json = JSONMANAGER.out.json .flatten() diff --git a/workflows/bindflow.nf b/workflows/bindflow.nf index 98739c7..bb2e30c 100644 --- a/workflows/bindflow.nf +++ b/workflows/bindflow.nf @@ -34,6 +34,7 @@ workflow BINDFLOW { ch_batches, quote_char ) + // // Collate and save software versions // From 2c003f57dae22f58864f9df9af52465510fbf111 Mon Sep 17 00:00:00 2001 From: Ziad Al-Bkhetan Date: Tue, 10 Mar 2026 13:50:17 +1100 Subject: [PATCH 5/7] Update Dockerfile --- docker/Dockerfile | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/docker/Dockerfile b/docker/Dockerfile index 3a3402f..ebd2819 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -30,6 +30,13 @@ RUN sed -i 's/\/conda-forge/g' install_bindcraft.sh RUN bash install_bindcraft.sh --cuda '12.4' --pkg_manager 'micromamba' --no-pyrosetta ENV PATH="$PATH:/work/FreeBindCraft" -ENV PYTHONPATH="${PYTHONPATH:-}/work/FreeBindCraft" +USER root +RUN apt-get update && apt-get install -y --no-install-recommends procps && rm -rf /var/lib/apt/lists/* +USER mambauser +RUN micromamba install -y -n base -c conda-forge ffmpeg && \ + micromamba clean -a -y +ENV PATH=/opt/conda/bin:$PATH +ENV PYTHONPATH="/work/FreeBindCraft" + # Set default command -CMD ["python", "-u", "bindcraft.py"] \ No newline at end of file +CMD ["python", "-u", "bindcraft.py"] From 9f13a6a76a9e09732e33f6128ded286b0becf94d Mon Sep 17 00:00:00 2001 From: ziadbkh Date: Wed, 1 Apr 2026 12:52:25 +1100 Subject: [PATCH 6/7] increase wall time for bindcraft --- conf/modules.config | 1 + 1 file changed, 1 insertion(+) diff --git a/conf/modules.config b/conf/modules.config index c3caf69..ee73cc2 100644 --- a/conf/modules.config +++ b/conf/modules.config @@ -42,6 +42,7 @@ process { withName: 'BINDCRAFT' { container = "${params.bindcraft_container}" memory = 24.GB + time = 24.h errorStrategy = params.error_strategy ?: "terminate" afterScript = { def output_path = new File("${params.outdir}") From 9842f91585d1c5ecb37fafd919cb7aa2eeedc85e Mon Sep 17 00:00:00 2001 From: ziadbkh Date: Thu, 9 Apr 2026 11:54:17 +1000 Subject: [PATCH 7/7] update reporting --- subworkflows/local/run_bindcraft.nf | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/subworkflows/local/run_bindcraft.nf b/subworkflows/local/run_bindcraft.nf index bf08e93..dffe74a 100644 --- a/subworkflows/local/run_bindcraft.nf +++ b/subworkflows/local/run_bindcraft.nf @@ -84,10 +84,25 @@ workflow RUN_BINDCRAFT { BINDCRAFT.out.stats.map{[["id": it[0].id], it[1]]}.groupTuple(), BINDCRAFT.out.accepted_ranked.map{[["id": it[0].id], it[1]]}.groupTuple() ) - + GENERATE_REPORT( BINDCRAFT.out.output_dir.map{it[1]}.collect(), - BINDCRAFT.out.failure_csv.map{it[1]}, + BINDCRAFT.out.failure_csv + .map{it[1]} + .splitCsv( header: true ) + .collect() + .map{ + it.inject([:]) { acc, map -> + map.each { k, v -> + def num = v?.toString()?.isNumber() ? v.toBigDecimal() : 0 + acc[k] = (acc[k] ?: 0) + num + } + def keys = acc.keySet().toList() + def values = keys.collect { acc[it] } + keys.join(',') + "\n" + values.join(',') + } + }.collectFile( name: "failure_csv.csv") + , RANKER.out.stats.map{it[1]}, BINDCRAFT.out.mpnn_design_stats .map{it[1].text}