From 4412ec42a9a4ed8475c2bcfc9ac71e682b3b386e Mon Sep 17 00:00:00 2001 From: Valerij Talagayev Date: Mon, 18 May 2026 00:42:27 +0200 Subject: [PATCH 1/3] added peptide pharmacophore and figures --- .../openmmdl_analysis/openmmdlanalysis.py | 50 +- .../visualization/figures.py | 462 +++++++++++++++++- .../visualization/pharmacophore.py | 171 ++++--- .../test_openmmdlanalysis.py | 31 +- .../visualization/test_figures.py | 76 ++- .../visualization/test_pharmacophore.py | 129 ++++- 6 files changed, 838 insertions(+), 81 deletions(-) diff --git a/openmmdl/openmmdl_analysis/openmmdlanalysis.py b/openmmdl/openmmdl_analysis/openmmdlanalysis.py index 8d313635..644102b4 100644 --- a/openmmdl/openmmdl_analysis/openmmdlanalysis.py +++ b/openmmdl/openmmdl_analysis/openmmdlanalysis.py @@ -18,7 +18,7 @@ def _load_analysis_dependencies(): global Preprocessing, RMSDAnalyzer, InteractionAnalyzer, BindingModeProcesser global MarkovChainAnalysis global FigureHighlighter, LigandImageGenerator - global FigureArranger, FigureMerger + global FigureArranger, FigureMerger, PeptideBindingModeFigureGenerator global BarcodeGenerator, BarcodePlotter global TrajectorySaver, PharmacophoreGenerator, StableWaters global update_dict, update_values @@ -46,6 +46,7 @@ def _load_analysis_dependencies(): from openmmdl.openmmdl_analysis.visualization.figures import ( FigureArranger, FigureMerger, + PeptideBindingModeFigureGenerator, ) from openmmdl.openmmdl_analysis.visualization.barcodes import ( BarcodeGenerator, @@ -73,6 +74,9 @@ def pushd(path: str): _FLOAT_RE = re.compile(r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?") +_NUMPY_FLOAT_WRAPPER_RE = re.compile( + r"(?:np\.)?float\d*\(\s*([-+]?\d*\.?\d+(?:[eE][-+]?\d+)?)\s*\)" +) _TRUE_BOOL_VALUES = {"1", "true", "t", "yes", "y", "on"} _FALSE_BOOL_VALUES = {"0", "false", "f", "no", "n", "off"} @@ -104,6 +108,7 @@ def parse_xyz(val): if val is None or val == 0 or val == "0": return None s = str(val) + s = _NUMPY_FLOAT_WRAPPER_RE.sub(r"\1", s) nums = _FLOAT_RE.findall(s) if len(nums) < 3: return None @@ -639,6 +644,34 @@ def run_analysis(args) -> int: fig_type, ) generator.generate_image() + else: + matplotlib.use("Agg") + merged_image_paths = [] + + peptide_figure_generator = PeptideBindingModeFigureGenerator( + complex_pdb_file=complex_pdb, + peptide_chain_id=peptide, + rdkit_max_residues=12, + ) + + for binding_mode, values in columns_with_value_1.items(): + occurrence_count = top_10_nodes_with_occurrences[binding_mode] + occurrence_percent = 100 * occurrence_count / total_frames + + output_png = peptide_figure_generator.generate_binding_mode_figure( + binding_mode=binding_mode, + interaction_values=values, + occurrence_percent=occurrence_percent, + output_dir=".", + ) + + merged_image_paths.append(output_png) + + figure_arranger = FigureArranger( + merged_image_paths, + "all_binding_modes_arranged.png", + ) + figure_arranger.arranged_figure_generation() logger.info(f"\033[1m{schema} bindig mode: Binding mode figure generated\033[0m") except Exception: logger.warning("Ligand could not be recognized, use the -l option") @@ -646,7 +679,13 @@ def run_analysis(args) -> int: df_all = pd.read_csv("df_all.csv") - pham_generator = PharmacophoreGenerator(df_all, ligand) + pharmacophore_partner_name = ligand if ligand is not None else f"peptide_chain_{peptide}" + pham_generator = PharmacophoreGenerator( + df_all, + pharmacophore_partner_name, + occurrence_threshold=threshold, + total_frames=total_frames, + ) # get the top 10 bindingmodes with the most occurrences binding_modes = grouped_frames_treshold["Binding_fingerprint_treshold"].str.split("\n") all_binding_modes = [mode.strip() for sublist in binding_modes for mode in sublist] @@ -661,9 +700,12 @@ def run_analysis(args) -> int: "Percentage Occurrence": [], } if generate_representative_frame: - DM = rmsd_analyzer.calculate_distance_matrix( - f"protein or nucleic or resname {ligand} or resname {special_ligand}", n_frames=md_len - 1 + representative_selection = ( + f"protein or nucleic or chainID {peptide} or resname {special_ligand}" + if peptide is not None + else f"protein or nucleic or resname {ligand} or resname {special_ligand}" ) + DM = rmsd_analyzer.calculate_distance_matrix(representative_selection, n_frames=md_len - 1) modes_to_process = top_10_binding_modes.index for mode in tqdm(modes_to_process): result_dict["Binding Mode"].append(mode) diff --git a/openmmdl/openmmdl_analysis/visualization/figures.py b/openmmdl/openmmdl_analysis/visualization/figures.py index 26d461d7..93460bff 100644 --- a/openmmdl/openmmdl_analysis/visualization/figures.py +++ b/openmmdl/openmmdl_analysis/visualization/figures.py @@ -1,8 +1,36 @@ import os +import re + +import cairosvg +import MDAnalysis as mda +import matplotlib.pyplot as plt +from rdkit import Chem +from rdkit.Chem import AllChem, rdCoordGen, rdDepictor +from rdkit.Chem.Draw import rdMolDraw2D import pylab from PIL import Image +AA3_TO_1 = { + "ALA": "A", "ARG": "R", "ASN": "N", "ASP": "D", "CYS": "C", + "GLN": "Q", "GLU": "E", "GLY": "G", "HIS": "H", "ILE": "I", + "LEU": "L", "LYS": "K", "MET": "M", "PHE": "F", "PRO": "P", + "SER": "S", "THR": "T", "TRP": "W", "TYR": "Y", "VAL": "V", + "HID": "H", "HIE": "H", "HIP": "H", +} + +INTERACTION_COLORS = { + "hbond": "#4C78A8", + "hydrophobic": "#F58518", + "waterbridge": "#72B7B2", + "saltbridge": "#E45756", + "pistacking": "#B279A2", + "pication": "#9D755D", + "halogen": "#54A24B", + "metal": "#BAB0AC", +} + + class FigureMerger: """ Handles the creation and merging of binding mode figures with corresponding legends. @@ -188,14 +216,438 @@ def arranged_figure_generation(self): # Save the big figure big_figure.save(self.output_path, "PNG") - # Ensure target directory exists + # Ensure target directories exist target_dir = "Binding_Modes_Markov_States" + individual_dir = os.path.join(target_dir, "individual_figures") os.makedirs(target_dir, exist_ok=True) + os.makedirs(individual_dir, exist_ok=True) - # Move the file + # Move the arranged overview figure into the main binding-mode directory new_path = os.path.join(target_dir, os.path.basename(self.output_path)) - os.rename(self.output_path, new_path) + os.replace(self.output_path, new_path) - # Remove the individual image files + # Move the individual binding-mode figures into a dedicated subfolder for path in self.merged_image_paths: - os.remove(path) + individual_path = os.path.join(individual_dir, os.path.basename(path)) + if os.path.abspath(path) != os.path.abspath(individual_path): + os.replace(path, individual_path) + + +class PeptideBindingModeFigureGenerator: + def __init__(self, complex_pdb_file, peptide_chain_id, rdkit_max_residues=12): + self.complex_pdb_file = complex_pdb_file + self.peptide_chain_id = peptide_chain_id + self.rdkit_max_residues = rdkit_max_residues + + def _load_peptide_residues(self): + u = mda.Universe(self.complex_pdb_file) + peptide_atoms = u.select_atoms(f"chainID {self.peptide_chain_id}") + + residues = [] + for residue in peptide_atoms.residues: + resname3 = str(residue.resname).upper() + residues.append( + { + "resnum": int(residue.resid), + "resname3": resname3, + "one_letter": AA3_TO_1.get(resname3, "X"), + } + ) + + return residues + + def _extract_interaction_annotations(self, interaction_values): + """ + Parse peptide-mode binding-mode columns like: + 44GLUA_176LYS_LYS_PI_saltbridge + 81ASPA_182LYS_LYS_PI_saltbridge + 53SERA_180TYR_hydrophobic + 12GLUA_177LYS_Donor_hbond + + Returns + ------- + residue_to_interactions : dict[int, set[str]] + residue number -> interaction types + residue_to_details : dict[int, list[dict]] + residue number -> detailed interaction descriptions + """ + residue_to_interactions = {} + residue_to_details = {} + + for value in interaction_values: + parts = str(value).split("_") + if len(parts) < 3: + continue + + protein_token = parts[0] + peptide_token = parts[1] + interaction_type = parts[-1].lower() + extra_detail = "_".join(parts[2:-1]) if len(parts) > 3 else "" + + match = re.match(r"(\d+)([A-Za-z0-9]+)", peptide_token) + if match is None: + continue + + resnum = int(match.group(1)) + peptide_label = peptide_token + + residue_to_interactions.setdefault(resnum, set()).add(interaction_type) + residue_to_details.setdefault(resnum, []) + + entry = { + "peptide_label": peptide_label, + "protein_label": protein_token, + "interaction_type": interaction_type, + "extra_detail": extra_detail, + } + + if entry not in residue_to_details[resnum]: + residue_to_details[resnum].append(entry) + + return residue_to_interactions, residue_to_details + + def _hex_to_rgb_tuple(self, hex_color): + hex_color = hex_color.lstrip("#") + return tuple(int(hex_color[i:i + 2], 16) / 255 for i in (0, 2, 4)) + + def _primary_interaction_type(self, interaction_types): + priority = [ + "saltbridge", + "hbond", + "pication", + "pistacking", + "hydrophobic", + "waterbridge", + "halogen", + "metal", + ] + + for interaction_type in priority: + if interaction_type in interaction_types: + return interaction_type + + return sorted(interaction_types)[0] + + def _interaction_color_rgb(self, interaction_type): + return self._hex_to_rgb_tuple( + INTERACTION_COLORS.get(interaction_type, "#444444") + ) + + def _write_peptide_only_pdb(self, output_pdb): + u = mda.Universe(self.complex_pdb_file) + peptide_atoms = u.select_atoms(f"chainID {self.peptide_chain_id}") + peptide_atoms.write(output_pdb) + return output_pdb + + def _draw_rdkit_peptide( + self, + peptide_pdb, + residue_to_interactions, + residue_label_map, + output_svg, + ): + mol = Chem.MolFromPDBFile( + peptide_pdb, + sanitize=False, + removeHs=False, + proximityBonding=True, + ) + if mol is None: + raise ValueError(f"Could not create RDKit molecule from {peptide_pdb}") + + Chem.SanitizeMol(mol, catchErrors=True) + mol = Chem.RemoveHs(mol) + + # Prefer CoordGen for cleaner 2D peptide layouts; fall back to rdDepictor. + try: + rdCoordGen.AddCoords(mol) + except Exception: + rdDepictor.SetPreferCoordGen(True) + rdDepictor.Compute2DCoords(mol, canonOrient=True, clearConfs=True) + + highlight_atoms = [] + highlight_atom_colors = {} + note_atom_by_residue = {} + + for atom in mol.GetAtoms(): + residue_info = atom.GetPDBResidueInfo() + if residue_info is None: + continue + + resnum = residue_info.GetResidueNumber() + if resnum in residue_to_interactions: + atom_idx = atom.GetIdx() + highlight_atoms.append(atom_idx) + primary_interaction = self._primary_interaction_type( + residue_to_interactions[resnum] + ) + highlight_atom_colors[atom_idx] = self._interaction_color_rgb( + primary_interaction + ) + + atom_name = residue_info.GetName().strip() + if resnum not in note_atom_by_residue or atom_name == "CA": + note_atom_by_residue[resnum] = atom_idx + + for resnum, label_id in residue_label_map.items(): + atom_idx = note_atom_by_residue.get(resnum) + if atom_idx is not None: + mol.GetAtomWithIdx(atom_idx).SetProp("atomNote", str(label_id)) + + drawer = rdMolDraw2D.MolDraw2DSVG(1400, 500) + options = drawer.drawOptions() + options.padding = 0.08 + options.fixedBondLength = 32 + options.additionalAtomLabelPadding = 0.15 + options.minFontSize = 10 + options.maxFontSize = 18 + options.annotationFontScale = 1.2 + + drawer.DrawMolecule( + mol, + highlightAtoms=highlight_atoms, + highlightAtomColors=highlight_atom_colors, + ) + drawer.FinishDrawing() + + svg = drawer.GetDrawingText().replace("svg:", "") + with open(output_svg, "w") as f: + f.write(svg) + + def generate_binding_mode_figure( + self, + binding_mode, + interaction_values, + occurrence_percent, + output_dir=".", + ): + residues = self._load_peptide_residues() + residue_to_interactions, residue_to_details = self._extract_interaction_annotations( + interaction_values + ) + residue_label_map = self._build_residue_label_map(residue_to_details) + title = f"{binding_mode} ({occurrence_percent:.1f}%)" + + if len(residues) <= self.rdkit_max_residues: + peptide_pdb = os.path.join(output_dir, f"{binding_mode}_peptide_only.pdb") + svg_path = os.path.join(output_dir, f"{binding_mode}.svg") + main_png_path = os.path.join(output_dir, f"{binding_mode}_main.png") + legend_png_path = os.path.join(output_dir, f"{binding_mode}_legend.png") + merged_png_path = os.path.join(output_dir, f"{binding_mode}.png") + + self._write_peptide_only_pdb(peptide_pdb) + self._draw_rdkit_peptide( + peptide_pdb, + residue_to_interactions, + residue_label_map, + svg_path, + ) + cairosvg.svg2png(url=svg_path, write_to=main_png_path) + self._draw_interaction_legend( + residue_label_map, + residue_to_details, + legend_png_path, + title=title, + ) + self._merge_main_and_legend( + main_png_path, + legend_png_path, + merged_png_path, + ) + + return merged_png_path + + main_png_path = os.path.join(output_dir, f"{binding_mode}_main.png") + legend_png_path = os.path.join(output_dir, f"{binding_mode}_legend.png") + merged_png_path = os.path.join(output_dir, f"{binding_mode}.png") + + self._draw_sequence_diagram( + residues, + residue_to_interactions, + residue_label_map, + main_png_path, + title=title, + ) + + self._draw_interaction_legend( + residue_label_map, + residue_to_details, + legend_png_path, + title=title, + ) + self._merge_main_and_legend( + main_png_path, + legend_png_path, + merged_png_path, + ) + + return merged_png_path + + def _draw_sequence_diagram( + self, + residues, + residue_to_interactions, + residue_label_map, + output_png, + title=None, + ): + n_residues = len(residues) + fig_width = max(12, n_residues * 0.6) + + fig, ax = plt.subplots(figsize=(fig_width, 3.8)) + ax.hlines(y=0, xmin=0, xmax=max(n_residues - 1, 0), linewidth=1.5) + + for i, residue in enumerate(residues): + ax.text( + i, + 0, + residue["one_letter"], + ha="center", + va="center", + fontsize=12, + ) + ax.text( + i, + -0.45, + str(residue["resnum"]), + ha="center", + va="center", + fontsize=8, + ) + + labels = residue_to_interactions.get(residue["resnum"], set()) + if not labels: + continue + + interaction_type = self._primary_interaction_type(labels) + color = INTERACTION_COLORS.get(interaction_type, "#444444") + + ax.annotate( + "", + xy=(i, 0.18), + xytext=(i, 0.95), + arrowprops=dict( + arrowstyle="->", + lw=1.2, + color=color, + ), + ) + ax.text( + i, + 1.15, + str(residue_label_map[residue["resnum"]]), + ha="center", + va="bottom", + fontsize=10, + fontweight="bold", + color=color, + ) + + ax.text(-0.8, 0, "N-term", ha="right", va="center", fontsize=10) + ax.text(n_residues - 1 + 0.8, 0, "C-term", ha="left", va="center", fontsize=10) + + if title: + ax.set_title(title) + + ax.set_xlim(-1.5, n_residues + 0.5) + ax.set_ylim(-0.9, 1.9) + ax.axis("off") + + plt.tight_layout() + plt.savefig(output_png, dpi=300, bbox_inches="tight") + plt.close(fig) + + def _build_residue_label_map(self, residue_to_details): + """ + Assign compact integer labels (1, 2, 3, ...) to interacting peptide residues. + """ + interacting_residues = sorted(residue_to_details.keys()) + return {resnum: i + 1 for i, resnum in enumerate(interacting_residues)} + + def _draw_interaction_legend(self, residue_label_map, residue_to_details, output_png, title=None): + """ + Draw a side legend like: + 1. 176LYS + - saltbridge ↔ 44GLUA (LYS_PI) + - hbond ↔ 12GLUA (Donor) + """ + line_count = 2 + for resnum in sorted(residue_to_details): + line_count += 1 + len(residue_to_details[resnum]) + + fig_height = max(4, line_count * 0.35) + fig, ax = plt.subplots(figsize=(7.5, fig_height)) + ax.axis("off") + + y = 0.98 + if title: + ax.text(0.00, y, title, fontsize=12, fontweight="bold", va="top", transform=ax.transAxes) + y -= 0.08 + + for resnum in sorted(residue_to_details): + label_id = residue_label_map[resnum] + entries = residue_to_details[resnum] + peptide_label = entries[0]["peptide_label"] + + ax.text( + 0.00, + y, + f"{label_id}. {peptide_label}", + fontsize=11, + fontweight="bold", + va="top", + transform=ax.transAxes, + ) + y -= 0.06 + + for entry in entries: + interaction_type = entry["interaction_type"] + protein_label = entry["protein_label"] + extra_detail = entry["extra_detail"] + color = INTERACTION_COLORS.get(interaction_type, "#444444") + + detail_text = f"{interaction_type} {protein_label}" + if extra_detail: + detail_text += f" ({extra_detail})" + + ax.text( + 0.04, + y, + "■", + color=color, + fontsize=11, + va="top", + transform=ax.transAxes, + ) + ax.text( + 0.08, + y, + detail_text, + fontsize=10, + va="top", + transform=ax.transAxes, + ) + y -= 0.05 + + y -= 0.03 + + plt.tight_layout() + plt.savefig(output_png, dpi=300, bbox_inches="tight") + plt.close(fig) + + def _merge_main_and_legend(self, main_png, legend_png, output_png): + main_img = Image.open(main_png).convert("RGB") + legend_img = Image.open(legend_png).convert("RGB") + + total_width = main_img.width + legend_img.width + total_height = max(main_img.height, legend_img.height) + + merged = Image.new("RGB", (total_width, total_height), "white") + merged.paste(main_img, (0, 0)) + merged.paste(legend_img, (main_img.width, 0)) + merged.save(output_png) + + os.remove(main_png) + os.remove(legend_png) + + return output_png \ No newline at end of file diff --git a/openmmdl/openmmdl_analysis/visualization/pharmacophore.py b/openmmdl/openmmdl_analysis/visualization/pharmacophore.py index 6eba9e23..0127570e 100644 --- a/openmmdl/openmmdl_analysis/visualization/pharmacophore.py +++ b/openmmdl/openmmdl_analysis/visualization/pharmacophore.py @@ -3,7 +3,22 @@ import xml.etree.ElementTree as ET -COORD_RE = re.compile(r"\(\s*([-\d.]+)\s*,\s*([-\d.]+)\s*,\s*([-\d.]+)\s*\)") +FLOAT_RE = re.compile(r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?") +NUMPY_FLOAT_WRAPPER_RE = re.compile( + r"(?:np\.)?float\d*\(\s*([-+]?\d*\.?\d+(?:[eE][-+]?\d+)?)\s*\)" +) + + +def parse_coord_triplet(value): + if value is None or value == 0 or value == "0" or value == "skip": + return None + + text = NUMPY_FLOAT_WRAPPER_RE.sub(r"\1", str(value)) + nums = FLOAT_RE.findall(text) + if len(nums) < 3: + return None + + return [float(nums[0]), float(nums[1]), float(nums[2])] class PharmacophoreGenerator: @@ -27,13 +42,47 @@ class PharmacophoreGenerator: Dictionary containing interaction types and associated 3D coordinates with visualization metadata. """ - def __init__(self, df_all, ligand_name): + def __init__(self, df_all, ligand_name, occurrence_threshold=None, total_frames=None): self.df_all = df_all - self.ligand_name = ligand_name - self.complex_name = f"{ligand_name}_complex" - self.coord_pattern = re.compile(r"\(([\d.-]+), ([\d.-]+), ([\d.-]+)\)") + self.ligand_name = str(ligand_name) if ligand_name is not None else "peptide" + self.complex_name = f"{self.ligand_name}_complex" + self.occurrence_threshold = self._normalize_occurrence_threshold(occurrence_threshold) + self.total_frames = total_frames self.clouds = self._generate_clouds() + @staticmethod + def _normalize_occurrence_threshold(occurrence_threshold): + if occurrence_threshold is None: + return None + occurrence_threshold = float(occurrence_threshold) + if occurrence_threshold > 1: + occurrence_threshold /= 100 + return occurrence_threshold + + def _interaction_passes_occurrence_threshold(self, interaction): + if self.occurrence_threshold is None: + return True + if interaction not in self.df_all.columns: + return False + + df_hits = self.df_all[self.df_all[interaction] == 1] + if df_hits.empty: + return False + + if self.total_frames is not None: + total_frames = int(self.total_frames) + elif "FRAME" in self.df_all.columns: + total_frames = int(self.df_all["FRAME"].nunique()) + else: + total_frames = len(self.df_all) + + if total_frames <= 0: + return False + + occurrence_count = df_hits["FRAME"].nunique() if "FRAME" in df_hits.columns else len(df_hits) + return occurrence_count >= self.occurrence_threshold * total_frames + + def to_dict(self): """ Export the interaction cloud as a dictionary. @@ -178,41 +227,46 @@ def generate_md_pharmacophore_cloudcenters(self, output_filename, id_num=0): ) elif interaction == "pistacking": pharm = self._generate_pharmacophore_vectors(self.df_all.filter(regex=interaction).columns) - feature_id_counter += 1 - lig_loc = position[0] - prot_loc = position[1] + for feature_name, position in pharm.items(): + feature_id_counter += 1 + lig_loc = position[0] + prot_loc = position[1] - vector = np.array(lig_loc) - np.array(prot_loc) - normal_vector = vector / np.linalg.norm(vector) - x, y, z = normal_vector + vector = np.array(lig_loc) - np.array(prot_loc) + norm = np.linalg.norm(vector) + if norm == 0: + continue - plane = ET.SubElement( - pharmacophore, - "plane", - name=feature_type, - featureId=interaction, - optional="false", - disabled="false", - weight="1.0", - coreCompound=self.ligand_name, - id=f"feature{str(feature_id_counter)}", - ) - ET.SubElement( - plane, - "position", - x3=str(lig_loc[0]), - y3=str(lig_loc[1]), - z3=str(lig_loc[2]), - tolerance="0.9", - ) - ET.SubElement( - plane, - "normal", - x3=str(x), - y3=str(y), - z3=str(z), - tolerance="0.43633232", - ) + normal_vector = vector / norm + x, y, z = normal_vector + + plane = ET.SubElement( + pharmacophore, + "plane", + name=feature_type, + featureId=feature_name, + optional="false", + disabled="false", + weight="1.0", + coreCompound=self.ligand_name, + id=f"feature{str(feature_id_counter)}", + ) + ET.SubElement( + plane, + "position", + x3=str(lig_loc[0]), + y3=str(lig_loc[1]), + z3=str(lig_loc[2]), + tolerance="0.9", + ) + ET.SubElement( + plane, + "normal", + x3=str(x), + y3=str(y), + z3=str(z), + tolerance="0.43633232", + ) tree = ET.ElementTree(root) tree.write(f"{output_filename}.pml", encoding="UTF-8", xml_declaration=True) @@ -478,9 +532,9 @@ def _generate_clouds(self): for index, row in self.df_all.iterrows(): if row["LIGCOO"] != 0: - coord_match = self.coord_pattern.match(row["LIGCOO"]) - if coord_match: - x, y, z = map(float, coord_match.groups()) + coord = parse_coord_triplet(row["LIGCOO"]) + if coord is not None: + x, y, z = coord x, y, z = round(x, 3), round(y, 3), round(z, 3) interaction = row["INTERACTION"] if interaction == "hbond": @@ -492,9 +546,9 @@ def _generate_clouds(self): for index, row in self.df_all.iterrows(): if row["TARGETCOO"] != 0: - coord_match = self.coord_pattern.match(row["TARGETCOO"]) - if coord_match: - x, y, z = map(float, coord_match.groups()) + coord = parse_coord_triplet(row["TARGETCOO"]) + if coord is not None: + x, y, z = coord x, y, z = round(x, 3), round(y, 3), round(z, 3) if row["INTERACTION"] == "metal": interaction_coords["metal"].append([x, y, z]) @@ -544,6 +598,8 @@ def _generate_pharmacophore_centers(self, interactions): # Only rows where this interaction occurs if interaction not in self.df_all.columns: continue + if not self._interaction_passes_occurrence_threshold(interaction): + continue df_hits = self.df_all[self.df_all[interaction] == 1] if df_hits.empty: continue @@ -552,10 +608,10 @@ def _generate_pharmacophore_centers(self, interactions): sum_x = sum_y = sum_z = 0.0 for _, row in df_hits.iterrows(): - m = COORD_RE.search(str(row.get("LIGCOO", ""))) - if not m: + coord = parse_coord_triplet(row.get("LIGCOO", "")) + if coord is None: continue - x, y, z = map(float, m.groups()) + x, y, z = coord sum_x += x sum_y += y sum_z += z @@ -579,6 +635,8 @@ def _generate_pharmacophore_vectors(self, interactions): for interaction in interactions: if interaction not in self.df_all.columns: continue + if not self._interaction_passes_occurrence_threshold(interaction): + continue df_hits = self.df_all[self.df_all[interaction] == 1] if df_hits.empty: continue @@ -588,15 +646,15 @@ def _generate_pharmacophore_vectors(self, interactions): sum_a = sum_b = sum_c = 0.0 for _, row in df_hits.iterrows(): - lig_m = COORD_RE.search(str(row.get("LIGCOO", ""))) - prot_m = COORD_RE.search(str(row.get("PROTCOO", ""))) + lig_coord = parse_coord_triplet(row.get("LIGCOO", "")) + prot_coord = parse_coord_triplet(row.get("PROTCOO", "")) # Only count frames where BOTH ends are available - if not lig_m or not prot_m: + if lig_coord is None or prot_coord is None: continue - x, y, z = map(float, lig_m.groups()) - a, b, c = map(float, prot_m.groups()) + x, y, z = lig_coord + a, b, c = prot_coord sum_x += x sum_y += y @@ -633,13 +691,14 @@ def _generate_pharmacophore_centers_all_points(self, interactions): coord_pattern = re.compile(r"\(([\d.-]+), ([\d.-]+), ([\d.-]+)\)") pharmacophore = {} for interaction in interactions: + if not self._interaction_passes_occurrence_threshold(interaction): + continue pharmacophore_points = [] for index, row in self.df_all.iterrows(): if row[interaction] == 1: - coord_match = coord_pattern.match(row["LIGCOO"]) - if coord_match: - x, y, z = map(float, coord_match.groups()) - pharmacophore_points.append([x, y, z]) + coord = parse_coord_triplet(row["LIGCOO"]) + if coord is not None: + pharmacophore_points.append(coord) if pharmacophore_points: pharmacophore[interaction] = pharmacophore_points diff --git a/openmmdl/tests/openmmdl_analysis/test_openmmdlanalysis.py b/openmmdl/tests/openmmdl_analysis/test_openmmdlanalysis.py index 1608e92d..5914f1fd 100644 --- a/openmmdl/tests/openmmdl_analysis/test_openmmdlanalysis.py +++ b/openmmdl/tests/openmmdl_analysis/test_openmmdlanalysis.py @@ -2,7 +2,7 @@ import pytest -from openmmdl.openmmdl_analysis.openmmdlanalysis import parse_bool_flag +from openmmdl.openmmdl_analysis.openmmdlanalysis import parse_bool_flag, parse_xyz @pytest.mark.parametrize( @@ -61,3 +61,32 @@ def test_parse_bool_flag_rejects_invalid_values(value): assert str(exc_info.value) == ( "Boolean value expected. Use one of: true/false, yes/no, y/n, 1/0, on/off." ) + + +@pytest.mark.parametrize( + "value, expected", + [ + ("(1.0, 2.0, 3.0)", [1.0, 2.0, 3.0]), + ("[1.0, 2.0, 3.0]", [1.0, 2.0, 3.0]), + ( + "(np.float64(30.156), np.float64(41.233), np.float64(92.595))", + [30.156, 41.233, 92.595], + ), + ], +) +def test_parse_xyz_parses_plain_and_numpy_float_coordinates(value, expected): + assert parse_xyz(value) == expected + + +def test_parse_xyz_does_not_parse_float64_width_as_coordinate(): + parsed = parse_xyz( + "(np.float64(30.156), np.float64(41.233), np.float64(92.595))" + ) + + assert parsed == [30.156, 41.233, 92.595] + assert 64.0 not in parsed + + +@pytest.mark.parametrize("value", [None, 0, "0"]) +def test_parse_xyz_returns_none_for_missing_coordinates(value): + assert parse_xyz(value) is None \ No newline at end of file diff --git a/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py b/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py index 5a2542ec..35a926d6 100644 --- a/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py +++ b/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py @@ -3,7 +3,11 @@ import shutil from PIL import Image from pathlib import Path -from openmmdl.openmmdl_analysis.visualization.figures import FigureMerger, FigureArranger +from openmmdl.openmmdl_analysis.visualization.figures import ( + FigureArranger, + PeptideBindingModeFigureGenerator, +) + test_data_directory = Path( "openmmdl/tests/data/openmmdl_analysis/rdkit_figure_generation" @@ -181,6 +185,76 @@ def test_arranged_figure_generation(): assert output_path is not None +def test_peptide_interaction_annotations_are_parsed(): + generator = PeptideBindingModeFigureGenerator( + complex_pdb_file="dummy.pdb", + peptide_chain_id="B", + ) + + values = { + "44GLUA_176LYS_LYS_PI_saltbridge", + "81ASPA_182LYS_LYS_PI_saltbridge", + "53SERA_180TYR_hydrophobic", + "12GLUA_177LYS_Donor_hbond", + } + + residue_to_interactions, residue_to_details = generator._extract_interaction_annotations(values) + + assert residue_to_interactions == { + 176: {"saltbridge"}, + 182: {"saltbridge"}, + 180: {"hydrophobic"}, + 177: {"hbond"}, + } + + assert { + "peptide_label": "176LYS", + "protein_label": "44GLUA", + "interaction_type": "saltbridge", + "extra_detail": "LYS_PI", + } in residue_to_details[176] + + +def test_peptide_residue_label_map_is_stable_and_sorted(): + generator = PeptideBindingModeFigureGenerator( + complex_pdb_file="dummy.pdb", + peptide_chain_id="B", + ) + + residue_to_details = { + 182: [{"peptide_label": "182LYS"}], + 176: [{"peptide_label": "176LYS"}], + 180: [{"peptide_label": "180TYR"}], + } + + assert generator._build_residue_label_map(residue_to_details) == { + 176: 1, + 180: 2, + 182: 3, + } + + +def test_figure_arranger_moves_individual_figures_to_subfolder(tmp_path, monkeypatch): + monkeypatch.chdir(tmp_path) + + Image.new("RGB", (20, 20), "white").save("Binding_Mode_1.png") + Image.new("RGB", (20, 20), "white").save("Binding_Mode_2.png") + + arranger = FigureArranger( + merged_image_paths=["Binding_Mode_1.png", "Binding_Mode_2.png"], + output_path="all_binding_modes_arranged.png", + ) + + arranger.arranged_figure_generation() + + assert os.path.exists("Binding_Modes_Markov_States/all_binding_modes_arranged.png") + assert os.path.exists("Binding_Modes_Markov_States/individual_figures/Binding_Mode_1.png") + assert os.path.exists("Binding_Modes_Markov_States/individual_figures/Binding_Mode_2.png") + + assert not os.path.exists("Binding_Mode_1.png") + assert not os.path.exists("Binding_Mode_2.png") + + # Run the tests if __name__ == "__main__": pytest.main() diff --git a/openmmdl/tests/openmmdl_analysis/visualization/test_pharmacophore.py b/openmmdl/tests/openmmdl_analysis/visualization/test_pharmacophore.py index 57ba5dea..c84396af 100644 --- a/openmmdl/tests/openmmdl_analysis/visualization/test_pharmacophore.py +++ b/openmmdl/tests/openmmdl_analysis/visualization/test_pharmacophore.py @@ -7,7 +7,10 @@ from unittest.mock import patch, mock_open from io import StringIO -from openmmdl.openmmdl_analysis.visualization.pharmacophore import PharmacophoreGenerator +from openmmdl.openmmdl_analysis.visualization.pharmacophore import ( + PharmacophoreGenerator, + parse_coord_triplet, +) @pytest.fixture @@ -40,7 +43,8 @@ def test_init(sample_df): assert generator.df_all.equals(sample_df) assert generator.ligand_name == "test_ligand" assert generator.complex_name == "test_ligand_complex" - assert isinstance(generator.coord_pattern, re.Pattern) + assert generator.occurrence_threshold is None + assert generator.total_frames is None assert isinstance(generator.clouds, dict) @@ -198,18 +202,34 @@ def test_format_clouds(): assert formatted_clouds["acceptor"]["radius"] == 0.1 -def test_coordinate_parsing(): - """Test the coordinate pattern regex.""" - generator = PharmacophoreGenerator(pd.DataFrame(), "test_ligand") - - coord_str = "(1.234, 2.345, 3.456)" - match = generator.coord_pattern.match(coord_str) - - assert match is not None - x, y, z = map(float, match.groups()) - assert x == 1.234 - assert y == 2.345 - assert z == 3.456 +@pytest.mark.parametrize( + "coord_str, expected", + [ + ("(1.234, 2.345, 3.456)", [1.234, 2.345, 3.456]), + ("[1.234, 2.345, 3.456]", [1.234, 2.345, 3.456]), + ( + "(np.float64(30.156), np.float64(41.233), np.float64(92.595))", + [30.156, 41.233, 92.595], + ), + ], +) +def test_coordinate_parsing(coord_str, expected): + """Test robust coordinate parsing, including numpy scalar repr strings.""" + assert parse_coord_triplet(coord_str) == expected + + +@pytest.mark.parametrize("coord_str", [None, 0, "0", "skip", "not a coordinate"]) +def test_coordinate_parsing_returns_none_for_missing_values(coord_str): + assert parse_coord_triplet(coord_str) is None + + +def test_coordinate_parsing_does_not_parse_float64_width_as_coordinate(): + parsed = parse_coord_triplet( + "(np.float64(30.156), np.float64(41.233), np.float64(92.595))" + ) + + assert parsed == [30.156, 41.233, 92.595] + assert 64.0 not in parsed def test_empty_dataframe(): @@ -220,3 +240,84 @@ def test_empty_dataframe(): clouds = generator.clouds for interaction in clouds: assert len(clouds[interaction]["coordinates"]) == 0 + + +def test_init_uses_safe_name_for_peptide_pharmacophore(): + generator = PharmacophoreGenerator(pd.DataFrame(), None) + + assert generator.ligand_name == "peptide" + assert generator.complex_name == "peptide_complex" + + +def test_saltbridge_centers_parse_numpy_float_coordinates(): + df = pd.DataFrame( + { + "FRAME": [0, 1], + "INTERACTION": ["saltbridge", "saltbridge"], + "LIGCOO": [ + "(np.float64(30.156), np.float64(41.233), np.float64(92.595))", + "(np.float64(32.156), np.float64(43.233), np.float64(94.595))", + ], + "PROTCOO": ["0", "0"], + "TARGETCOO": ["0", "0"], + "PROTISDON": ["False", "False"], + "PROTISPOS": ["False", "False"], + "PI_saltbridge": [1, 1], + } + ) + + generator = PharmacophoreGenerator(df, "peptide_chain_B") + + assert generator._generate_pharmacophore_centers(["PI_saltbridge"]) == { + "PI_saltbridge": [31.156, 42.233, 93.595] + } + + +def test_occurrence_threshold_filters_low_frequency_features(): + df = pd.DataFrame( + { + "FRAME": list(range(10)), + "INTERACTION": ["hydrophobic"] * 10, + "LIGCOO": ["(1.0, 2.0, 3.0)"] * 10, + "PROTCOO": ["0"] * 10, + "TARGETCOO": ["0"] * 10, + "PROTISDON": ["False"] * 10, + "PROTISPOS": ["False"] * 10, + "hydrophobic": [1, 1, 1, 0, 0, 0, 0, 0, 0, 0], + } + ) + + generator = PharmacophoreGenerator( + df, + "peptide_chain_B", + occurrence_threshold=40, + total_frames=10, + ) + + assert generator._generate_pharmacophore_centers(["hydrophobic"]) == {} + + +def test_occurrence_threshold_keeps_features_at_cutoff(): + df = pd.DataFrame( + { + "FRAME": list(range(10)), + "INTERACTION": ["hydrophobic"] * 10, + "LIGCOO": ["(1.0, 2.0, 3.0)"] * 10, + "PROTCOO": ["0"] * 10, + "TARGETCOO": ["0"] * 10, + "PROTISDON": ["False"] * 10, + "PROTISPOS": ["False"] * 10, + "hydrophobic": [1, 1, 1, 1, 0, 0, 0, 0, 0, 0], + } + ) + + generator = PharmacophoreGenerator( + df, + "peptide_chain_B", + occurrence_threshold=40, + total_frames=10, + ) + + assert generator._generate_pharmacophore_centers(["hydrophobic"]) == { + "hydrophobic": [1.0, 2.0, 3.0] + } \ No newline at end of file From cc2b1ebb924607c48cabf5a0267486e4bc69bdb0 Mon Sep 17 00:00:00 2001 From: Valerij Talagayev Date: Mon, 18 May 2026 01:06:46 +0200 Subject: [PATCH 2/3] fix import --- openmmdl/tests/openmmdl_analysis/visualization/test_figures.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py b/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py index 35a926d6..54125f93 100644 --- a/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py +++ b/openmmdl/tests/openmmdl_analysis/visualization/test_figures.py @@ -4,6 +4,7 @@ from PIL import Image from pathlib import Path from openmmdl.openmmdl_analysis.visualization.figures import ( + FigureMerger, FigureArranger, PeptideBindingModeFigureGenerator, ) From 391ea787931776947868f951a72872786ea9ec49 Mon Sep 17 00:00:00 2001 From: Valerij Talagayev <82884038+talagayev@users.noreply.github.com> Date: Mon, 18 May 2026 02:15:40 +0200 Subject: [PATCH 3/3] Update CHANGELOG.md --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f93c3230..8efec368 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ talagayev - Added `openmmdl check` for testing correct installation (2026-05-03, PR#207) ### Fixed +- Fixed peptide pharmacophore error (2026-05-18, PR#209) - Fixed remaining water type recognition errors (2026-04-29, PR#205) - Fixed `OPC` water recognition error (2026-04-28, PR#204)