Merge pull request #1412 from haddocking/alascan-ligand

Alascan with non-standard ligand/molecules/residues

Author: VGPReys

CHANGELOG.md

@@ -1,5 +1,6 @@
 # Changelog
 
+- 2025-10-30: Added possibility to use alascan with ligands - Issue #1411
 - 2025-10-22: Allow the definition of chain combinations to be used for scoring - Issue #1414
 - 2025-09-11: Added `grid` mode
 - 2025-09-11: Corrected antibody-antigen notebook - Issue #1383

integration_tests/golden_data/protlig_complex_1.pdb

@@ -27,6 +27,30 @@ def alascan_module():
         yield alascan
 
 
+@pytest.fixture
+def alascan_module_protlig():
+    """Return a default alascan module."""
+    with tempfile.TemporaryDirectory() as tmpdir:
+        alascan = AlascanModule(
+            order=0, path=Path(tmpdir), initial_params=DEFAULT_ALASCAN_CONFIG
+        )
+        alascan.params["int_cutoff"] = 3.5
+        alascan.params["output_mutants"] = True
+        # Copy parameters and toplogy of the ligand
+        shutil.copy(
+            Path(GOLDEN_DATA, "ligand.top"),
+            Path(alascan.path, "ligand.top"),
+            )
+        shutil.copy(
+            Path(GOLDEN_DATA, "ligand.param"),
+            Path(alascan.path, "ligand.param"),
+            )
+        # Set the parameters to point the file
+        alascan.params["ligand_param_fname"] = Path(alascan.path, "ligand.param")
+        alascan.params["ligand_top_fname"] = Path(alascan.path, "ligand.top")
+        yield alascan
+
+
 class MockPreviousIO:
     def __init__(self, path):
         self.path = path
@@ -70,6 +94,26 @@ def output(self):
         return None
 
 
+class MockPreviousIO_protlig:
+    def __init__(self, path):
+        self.path = path
+
+    def retrieve_models(self, individualize: bool = False):
+        fname = "protlig_complex_1.pdb"
+        shutil.copy(
+            Path(GOLDEN_DATA, fname),
+            Path(self.path, fname),
+        )
+        model_list = [
+            PDBFile(file_name=fname, path=self.path),
+        ]
+
+        return model_list
+
+    def output(self):
+        return None
+
+
 def test_alascan_default(alascan_module, mocker):
     """Test the alascan module."""
     alascan_module.previous_io = MockPreviousIO(path=alascan_module.path)
@@ -148,3 +192,47 @@ def test_alascan_mutation_resiudes():
     config_allowed_resiudes = set(default_config["scan_residue"]["choices"])
     script_allowed_resiudes = set(list(RES_CODES.keys()))
     assert config_allowed_resiudes == script_allowed_resiudes
+
+
+def test_alascan_with_ligand_topar(alascan_module_protlig, mocker):
+    """Test the use of alascan in presence of a ligand."""
+    alascan_module_protlig.previous_io = MockPreviousIO_protlig(path=alascan_module_protlig.path)
+    alascan_module_protlig.run()
+
+    expected_csv = Path(alascan_module_protlig.path, "scan_protlig_complex_1.tsv")
+    expected_clt_csv = Path(alascan_module_protlig.path, "scan_clt_unclustered.tsv")
+
+    assert expected_csv.exists(), f"{expected_csv} does not exist"
+    assert expected_clt_csv.exists(), f"{expected_clt_csv} does not exist"
+
+    # List mutated files
+    mutated_filepaths = list(Path(alascan_module_protlig.path).glob("protlig_complex_1-*.pdb"))
+    assert len(mutated_filepaths) >= 1
+
+    # Loop over files
+    for mutated_fpath in mutated_filepaths:
+        # Make sure the ligand is in it
+        file_content = mutated_fpath.read_text()
+        assert file_content.count("G39") > 20
+
+
+def test_alascan_without_ligand_topar(alascan_module, mocker):
+    """Test the use of alascan in presence of a ligand without topo/param."""
+    alascan_module.previous_io = MockPreviousIO_protlig(path=alascan_module.path)
+    alascan_module.run()
+
+    expected_csv = Path(alascan_module.path, "scan_protlig_complex_1.tsv")
+    expected_clt_csv = Path(alascan_module.path, "scan_clt_unclustered.tsv")
+
+    assert expected_csv.exists(), f"{expected_csv} does not exist"
+    assert expected_clt_csv.exists(), f"{expected_clt_csv} does not exist"
+
+    # List mutated files
+    mutated_filepaths = list(Path(alascan_module.path).glob("protlig_complex_1-*.pdb"))
+    assert len(mutated_filepaths) >= 1
+
+    # Loop over files
+    for mutated_fpath in mutated_filepaths:
+        # Make sure the ligand is not in it
+        file_content = mutated_fpath.read_text()
+        assert file_content.count("G39") == 0

integration_tests/test_alascan.py

@@ -149,7 +149,7 @@ def main(
     from contextlib import suppress
     from pathlib import Path
 
-    from haddock import log
+    from haddock import log, EmptyPath
     from haddock.gear.haddockmodel import HaddockModel
     from haddock.gear.yaml2cfg import read_from_yaml_config
     from haddock.gear.zerofill import zero_fill
@@ -168,23 +168,29 @@ def main(
     ems_dict = default_emscoring.copy()
     n_warnings = 0
     for param, value in kwargs.items():
+        # Check if the parameter name is in the emscoring module ones
         if param not in default_emscoring:
             sys.exit(
                 f"* ERROR * Parameter {param!r} is not a "
                 f"valid `emscoring` parameter.{os.linesep}"
-                f"Valid emscoring parameters are: {', '.join(sorted(default_emscoring))}"
+                "Valid emscoring parameters are: "
+                f"{', '.join(sorted(default_emscoring))}"
                 )
+        # Compare the user-given value to the default one
         if value != default_emscoring[param]:
             print(
-                f"* ATTENTION * Value ({value}) of parameter {param} different from default ({default_emscoring[param]})"
-            )  # noqa:E501
+                f"* ATTENTION * Value ({value}) of parameter {param} "
+                f"different from default ({default_emscoring[param]})"
+            )
             # get the type of default value
             default_type = type(default_emscoring[param])
-            # convert the value to the same type
+            # cast the value to the same type
             if default_type == bool:
                 if value.lower() not in ["true", "false"]:
                     sys.exit(f"* ERROR * Boolean parameter {param} should be True or False")
                 value = value.lower() == "true"
+            elif param.endswith("_fname"):
+                value = EmptyPath() if str(value) == "" else Path(value).resolve()
             else:
                 value = default_type(value)
             ems_dict[param] = value
@@ -210,9 +216,14 @@ def main(
         # create a copy of the input pdb
         input_pdb_copy = Path(tmp.name)
         shutil.copy(input_pdb, input_pdb_copy)
-    
-        params = {
-            "topoaa": {"molecules": [input_pdb_copy]},
+
+        # Setting up a full workflow set of parameters
+        workflow_params = {
+            "topoaa": {
+                "molecules": [input_pdb_copy],
+                "ligand_param_fname": ems_dict["ligand_param_fname"],
+                "ligand_top_fname": ems_dict["ligand_top_fname"],
+                },
             "emscoring": ems_dict,
         }
 
@@ -221,13 +232,13 @@ def main(
         # run workflow
         with working_directory(run_dir):
             workflow = WorkflowManager(
-                workflow_params=params,
+                workflow_params=workflow_params,
                 start=0,
                 run_dir=run_dir,
             )
-
             workflow.run()
 
+    # Build expected pdb filepath
     minimized_mol = Path(run_dir, "1_emscoring", "emscoring_1.pdb")
     haddock_score_component_dic = HaddockModel(minimized_mol).energies
 
@@ -237,6 +248,7 @@ def main(
     air = haddock_score_component_dic["air"]
     bsa = haddock_score_component_dic["bsa"]
 
+    # Compute the haddock score
     # emscoring is equivalent to itw
     haddock_score_itw = (
         ems_dict["w_vdw"] * vdw

src/haddock/clis/cli_score.py

@@ -118,7 +118,7 @@ def _run(self):
                 mutation_res=self.params["scan_residue"],
                 model=model,
                 params=self.params,
-                library_mode = False
+                library_mode=False,
                 )
             for model in models
             ]
@@ -202,4 +202,4 @@ def _run(self):
             # Send models to the next step, no operation is done on them
             self.output_models = models
 
-        self.export_io_models()
+        self.export_io_models()

src/haddock/modules/analysis/alascan/__init__.py

@@ -67,4 +67,21 @@ output_mutants:
   short: Dump the mutated, energy-minimized PDB files.
   long: Dump the mutated, energy-minimized PDB files. As the number of mutants can be very large, this option is allowed only when a single model is provided in input.
   group: analysis
-  explevel: easy
+  explevel: easy
+ligand_param_fname:
+  default: ''
+  type: file
+  title: Custom ligand parameter file
+  short: Ligand parameter file in CNS format
+  long: Ligand parameter file in CNS format, for any ligand/residues/molecules not supported by default by HADDOCK.
+  group: 'force field'
+  explevel: easy
+ligand_top_fname:
+  default: ''
+  type: file
+  title: Custom ligand topology file
+  short: Ligand topology file in CNS format
+  long: Ligand topology file in CNS format containing the ligand topologies
+    (atoms, masses, charges, bond definitions...) for any ligand not supported by default by HADDOCK
+  group: 'force field'
+  explevel: easy