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} |
+
+
+
+
+ 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/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..ee73cc2 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' {
@@ -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}")
diff --git a/docker/Dockerfile b/docker/Dockerfile
new file mode 100644
index 0000000..ebd2819
--- /dev/null
+++ b/docker/Dockerfile
@@ -0,0 +1,42 @@
+# 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"
+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"]
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/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
diff --git a/modules/local/bindcraft/main.nf b/modules/local/bindcraft/main.nf
index e1e6cda..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:
@@ -22,7 +26,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/modules/local/generate_report.nf b/modules/local/generate_report.nf
new file mode 100644
index 0000000..3395357
--- /dev/null
+++ b/modules/local/generate_report.nf
@@ -0,0 +1,51 @@
+process GENERATE_REPORT {
+ tag "all-run"
+ label 'process_single'
+ container 'ghcr.io/australian-protein-design-initiative/containers/nf-binder-design-utils:0.1.5'
+
+ input:
+ // 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:
+ 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:
+ """
+
+ 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/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"
diff --git a/subworkflows/local/run_bindcraft.nf b/subworkflows/local/run_bindcraft.nf
index f74ee25..dffe74a 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
@@ -83,6 +85,34 @@ workflow RUN_BINDCRAFT {
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]}
+ .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}
+ .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()
@@ -93,6 +123,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
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
//