Fix peptide name error

CHANGELOG.md

@@ -28,6 +28,7 @@ talagayev
 - Fixed `OPC` water recognition error (2026-04-28, PR#204)
 
 ### Changed
+- Changed `-n` to optional to fix peptide issues in `OpenMMDL Analysis` (2026-05-13, PR#208)
 - Changed imports to lightweight imports in `OpenMMDL Analysis` (2026-05-03, PR#206)
 
 ## Version 1.3.0

docs/openmmdl_analysis.rst

@@ -22,13 +22,12 @@ Mandatory:
 
     -t = topology file of the simulation (in .pdb format)
     -d = trajectory file of the simulation (in .dcd format)
-    -n = Ligand name (3 letter code in PDB)
-
 
 Optional:
 
 .. code-block:: text
 
+    -n = Ligand name (3 letter code in PDB)
     -l = Ligand in SDF format
     -b = binding mode threshold. Is used to remove interactions under the defined procentual occurence from the binding mode generation. The default is 40% (accepted values: 0-100)
     -df = Dataframe (use if the interactions were already calculated, default name would be "interactions_gathered.csv")

openmmdl/openmmdl_analysis/analysis/interactions.py

@@ -167,8 +167,14 @@ def _retrieve_plip_interactions_peptide(self, pdb_file):
         dict
             A dictionary of the binding sites and the interactions.
         """
+        config.PEPTIDES = [self.peptide]
+
         protlig = PDBComplex()
         protlig.load_pdb(pdb_file)  # load the pdb file
+
+        if not protlig.ligands:
+            raise ValueError(f"PLIP did not detect peptide chain {self.peptide!r} as a ligand")
+
         protlig.characterize_complex(protlig.ligands[-1])  # find ligands and analyze interactions
         sites = {}
         # loop over binding sites
@@ -285,9 +291,14 @@ def _process_frame(self, frame):
         pd.DataFrame
             A dataframe conatining the interaction data for the processed frame.
         """
-        atoms_selected = self.pdb_md.select_atoms(
-            f"protein or nucleic or resname {self.lig_name} or (resname HOH and around 10 resname {self.lig_name}) or resname {self.special}"
-        )
+        if self.peptide is not None:
+            atoms_selected = self.pdb_md.select_atoms(
+                f"protein or nucleic or (resname HOH and around 10 chainID {self.peptide}) or resname {self.special}"
+            )
+        else:
+            atoms_selected = self.pdb_md.select_atoms(
+                f"protein or nucleic or resname {self.lig_name} or (resname HOH and around 10 resname {self.lig_name}) or resname {self.special}"
+            )
         for num in self.pdb_md.trajectory[(frame) : (frame + 1)]:
             atoms_selected.write(f"processing_frame_{frame}.pdb")
         if self.peptide is None:
@@ -926,7 +937,8 @@ def _process_trajectory_prolif(self) -> pd.DataFrame:
             )
         else:
             # peptide treated as ligand chain; pocket = rest of protein
-            ligand_ag = self.pdb_md.select_atoms(f"chainID {self.peptide}")
+            lig_sel = f"chainID {self.peptide}"
+            ligand_ag = self.pdb_md.select_atoms(lig_sel)
             base = "protein or nucleic" if getattr(config, "DNARECEPTOR", False) else "protein"
             protein_ag = self.pdb_md.select_atoms(f"({base}) and not chainID {self.peptide}")
 

openmmdl/openmmdl_analysis/analysis/rmsd.py

@@ -70,16 +70,16 @@ def rmsd_for_atomgroups(self, fig_type, selection1, selection2=None):
 
         return rmsd_df
 
-    def rmsd_dist_frames(self, fig_type, lig, nucleic=False):
+    def rmsd_dist_frames(self, fig_type, lig_selection, nucleic=False):
         """
         Calculate the RMSD between all frames in a matrix.
 
         Parameters
         ----------
         fig_type : str
             Type of the figure to save (e.g., 'png', 'jpg').
-        lig : str
-            Ligand name saved in the above PDB file. Selection string for the MDAnalysis AtomGroup to be investigated, also used during alignment.
+        lig_selection : str
+            Selection string for the MDAnalysis AtomGroup to be investigated.
         nucleic : bool, optional
             Bool indicating if the receptor to be analyzed contains nucleic acids. Defaults to False.
 
@@ -94,8 +94,8 @@ def rmsd_dist_frames(self, fig_type, lig, nucleic=False):
             pairwise_rmsd_prot = diffusionmap.DistanceMatrix(self.universe, select="nucleic").run().dist_matrix
         else:
             pairwise_rmsd_prot = diffusionmap.DistanceMatrix(self.universe, select="protein").run().dist_matrix
-        pairwise_rmsd_lig = diffusionmap.DistanceMatrix(self.universe, f"resname {lig}").run().dist_matrix
-
+        pairwise_rmsd_lig = diffusionmap.DistanceMatrix(self.universe, select=lig_selection).run().dist_matrix
+        
         max_dist = max(np.amax(pairwise_rmsd_lig), np.amax(pairwise_rmsd_prot))
 
         fig, ax = plt.subplots(1, 2)

openmmdl/openmmdl_analysis/openmmdlanalysis.py

@@ -160,7 +160,7 @@ def run_analysis(args) -> int:
     topology = os.path.abspath(args.topology)
     trajectory = os.path.abspath(args.trajectory)
 
-    if args.ligand_name is None:
+    if args.ligand_name is None and args.peptide is None:
         logger.error("Ligand name is missing. Add the name of your ligand from your topology file")
         sys.exit(1)
     # set variables for analysis and preprocess input files
@@ -213,7 +213,7 @@ def run_analysis(args) -> int:
     preprocessor.process_pdb_file(topology)
     logger.info("\033[1mFiles are preprocessed\033[0m")
 
-    if ligand_sdf is None:
+    if peptide is None and ligand_sdf is None:
         preprocessor.extract_and_save_ligand_as_sdf(topology, "./ligand_prepared.sdf", ligand)
         ligand_sdf = "./ligand_prepared.sdf"
 
@@ -224,9 +224,14 @@ def run_analysis(args) -> int:
 
     # TODO maybe put this part into a function possibly in visualization_functions.py TrajectorySaver
     # Writing out the complex of the protein and ligand with water around 10A of the ligand
-    complex = pdb_md.select_atoms(
-        f"protein or nucleic or resname {ligand} or (resname HOH and around 10 resname {ligand}) or resname {special_ligand}"
-    )
+    if peptide is not None:
+        complex = pdb_md.select_atoms(
+            f"protein or nucleic or (resname HOH and around 10 chainID {peptide}) or resname {special_ligand}"
+        )
+    else:
+        complex = pdb_md.select_atoms(
+            f"protein or nucleic or resname {ligand} or (resname HOH and around 10 resname {ligand}) or resname {special_ligand}"
+        )
     complex.write("complex.pdb")
     preprocessor.renumber_atoms_in_residues("complex.pdb", "complex.pdb", ligand)
     preprocessor.process_pdb("complex.pdb", "complex.pdb")
@@ -282,25 +287,29 @@ def run_analysis(args) -> int:
 
     os.makedirs("RMSD", exist_ok=True)
     rmsd_analyzer = RMSDAnalyzer(f"{topology}", f"{trajectory}")
-    if receptor_nucleic:
+    if peptide is not None:
         rmsd_analyzer.rmsd_for_atomgroups(
             fig_type,
-            selection1="nucleicbackbone",
-            selection2=["nucleic", f"resname {ligand}"],
+            selection1="nucleicbackbone" if receptor_nucleic else "backbone",
+            selection2=["nucleic" if receptor_nucleic else "protein", f"chainID {peptide}"],
         )
         if frame_rmsd:
             pairwise_rmsd_prot, pairwise_rmsd_lig = rmsd_analyzer.rmsd_dist_frames(
-                fig_type, lig=f"{ligand}", nucleic=True
+                fig_type,
+                lig_selection=f"chainID {peptide}",
+                nucleic=receptor_nucleic,
             )
             logger.info("\033[1mRMSD calculated\033[0m")
-    elif peptide is not None:
+    elif receptor_nucleic:
         rmsd_analyzer.rmsd_for_atomgroups(
             fig_type,
-            selection1="backbone",
-            selection2=["protein", f"chainID {peptide}"],
+            selection1="nucleicbackbone",
+            selection2=["nucleic", f"resname {ligand}"],
         )
         if frame_rmsd:
-            pairwise_rmsd_prot, pairwise_rmsd_lig = rmsd_analyzer.rmsd_dist_frames(fig_type, lig=f"chainID {peptide}")
+            pairwise_rmsd_prot, pairwise_rmsd_lig = rmsd_analyzer.rmsd_dist_frames(
+                fig_type, lig_selection=f"resname {ligand}", nucleic=True
+            )
             logger.info("\033[1mRMSD calculated\033[0m")
     else:
         rmsd_analyzer.rmsd_for_atomgroups(
@@ -309,7 +318,9 @@ def run_analysis(args) -> int:
             selection2=["protein", f"resname {ligand}"],
         )
         if frame_rmsd:
-            pairwise_rmsd_prot, pairwise_rmsd_lig = rmsd_analyzer.rmsd_dist_frames(fig_type, lig=f"{ligand}")
+            pairwise_rmsd_prot, pairwise_rmsd_lig = rmsd_analyzer.rmsd_dist_frames(
+                fig_type, lig_selection=f"resname {ligand}"
+            )
             logger.info("\033[1mRMSD calculated\033[0m")
 
     if receptor_nucleic:

openmmdl/tests/openmmdl_analysis/analysis/test_bindingmodes.py

@@ -105,3 +105,34 @@ def test_interactions_all_includes_rare_contacts_for_ligand_schema(sample_intera
     # But "all interactions" must still contain them
     assert "ligand_5_Acceptor_hbond" in bm_lig_high_thresh.unique_data_all
     assert "ligand_7_hydrophobic" in bm_lig_high_thresh.unique_data_all
+
+
+def test_peptide_ligand_schema_does_not_require_ligand_name():
+    """Peptide binding-mode processing should work with ligand=None."""
+    interaction_list = pd.DataFrame(
+        [
+            {
+                "FRAME": 1,
+                "INTERACTION": "hbond",
+                "Prot_partner": "12GLYA",
+                "PROTISDON": True,
+                "RESNR_LIG": 3,
+                "RESTYPE_LIG": "GLY",
+            }
+        ]
+    )
+
+    bm_pep = BindingModeProcesser(
+        pdb_md=None,
+        ligand=None,
+        peptide="B",
+        special=None,
+        ligand_rings=None,
+        interaction_list=interaction_list,
+        threshold=0,
+        total_frames=10,
+        schema="ligand",
+    )
+
+    assert "ligand_3GLY_Acceptor_hbond" in bm_pep.unique_data
+    assert bm_pep.interaction_list.loc[0, "ligand_3GLY_Acceptor_hbond"] == 1

openmmdl/tests/openmmdl_analysis/analysis/test_interactions.py

@@ -0,0 +1,44 @@
+import sys
+from unittest.mock import MagicMock, patch
+
+import pytest
+from plip.basic import config
+
+from openmmdl.openmmdl_analysis.analysis.interactions import InteractionAnalyzer
+
+
+def test_process_trajectory_prolif_peptide_empty_selection_raises_clean_error(monkeypatch):
+    """Wrong peptide chain IDs should raise a clean ValueError, not UnboundLocalError."""
+    monkeypatch.setitem(sys.modules, "prolif", MagicMock())
+
+    analyzer = InteractionAnalyzer.__new__(InteractionAnalyzer)
+    analyzer.special = None
+    analyzer.peptide = "B"
+    analyzer.pdb_md = MagicMock()
+
+    empty_ligand_ag = MagicMock()
+    empty_ligand_ag.__len__.return_value = 0
+    receptor_ag = MagicMock()
+    receptor_ag.__len__.return_value = 10
+    analyzer.pdb_md.select_atoms.side_effect = [empty_ligand_ag, receptor_ag]
+
+    with pytest.raises(ValueError, match="ProLIF: ligand selection returned 0 atoms: chainID B"):
+        analyzer._process_trajectory_prolif()
+
+
+@patch("openmmdl.openmmdl_analysis.analysis.interactions.PDBComplex")
+def test_retrieve_plip_interactions_peptide_sets_config_peptides(mock_pdb_complex, monkeypatch):
+    """Peptide PLIP workers should not rely on inherited parent-process config."""
+    monkeypatch.setattr(config, "PEPTIDES", [], raising=False)
+
+    analyzer = InteractionAnalyzer.__new__(InteractionAnalyzer)
+    analyzer.peptide = "B"
+
+    protlig = mock_pdb_complex.return_value
+    protlig.ligands = []
+
+    with pytest.raises(ValueError, match="PLIP did not detect peptide chain 'B' as a ligand"):
+        analyzer._retrieve_plip_interactions_peptide("processing_frame_1.pdb")
+
+    assert config.PEPTIDES == ["B"]
+    protlig.load_pdb.assert_called_once_with("processing_frame_1.pdb")

openmmdl/tests/openmmdl_analysis/analysis/test_rmsd.py

@@ -101,6 +101,18 @@ def test_rmsd_for_atomgroups_with_multiple_selections(mock_savefig, mock_to_csv,
         assert len(result) == 3
         assert result.columns.tolist() == ["protein", "resname LIG"]
 
+@patch("matplotlib.pyplot.savefig")
+@patch("openmmdl.openmmdl_analysis.analysis.rmsd.diffusionmap.DistanceMatrix")
+def test_rmsd_dist_frames_uses_full_ligand_selection(mock_distance_matrix, mock_savefig, analyzer, test_dir):
+    """rmsd_dist_frames must pass peptide chain selections through unchanged."""
+    mock_distance_matrix.return_value.run.return_value.dist_matrix = np.array(
+        [[0.0, 1.0], [1.0, 0.0]]
+    )
+
+    analyzer.rmsd_dist_frames("png", lig_selection="chainID B")
+
+    assert mock_distance_matrix.call_args_list[0].kwargs["select"] == "protein"
+    assert mock_distance_matrix.call_args_list[1].kwargs["select"] == "chainID B"
 
 def test_calc_rmsd_2frames(analyzer):
     """Test RMSD calculation between two frames."""