Merge pull request #307 from Acellera/cyclic_peptides_bb_mods

Cyclic peptide support

pdb2pqr/biomolecule.py

@@ -501,19 +501,34 @@ def assign_termini(self, chain, neutraln=False, neutralc=False):
         if len(chain.residues) == 0:
             text = f'Error: chain "{chain.chain_id}" has 0 residues!'
             raise IndexError(text)
-        # Set the N-Terminus/ 5' Terminus
+
         res0 = chain.residues[0]
+        reslast = chain.residues[-1]
+        # Check if chain is cyclic. Amide distance ranges between 1.325 - 1.346
+        if "N" in res0.map and "C" in reslast.map:
+            dist = util.distance(res0.map["N"].coords, reslast.map["C"].coords)
+            if dist < 1.35:
+                # If the chain is cyclic, don't apply termini.
+                return
+
+        # Set the N-Terminus/ 5' Terminus
         if isinstance(res0, aa.Amino):
             res0.is_n_term = True
-            if isinstance(res0, aa.PRO) or neutraln:
+            # If N is bonded to more than one heavy atom switch to neutral-nterm
+            heavy_n_bonds = []
+            if "N" in res0.map:
+                heavy_n_bonds = [
+                    a for a in res0.map["N"].bonds if a.name[0] != "H"
+                ]
+
+            if neutraln or len(heavy_n_bonds) > 1:
                 self.apply_patch("NEUTRAL-NTERM", res0)
             else:
                 self.apply_patch("NTERM", res0)
         elif isinstance(res0, na.Nucleic):
             res0.is5term = True
             self.apply_patch("5TERM", res0)
         # Set the C-Terminus/ 3' Terminus
-        reslast = chain.residues[-1]
         if isinstance(reslast, aa.Amino):
             reslast.is_c_term = True
             if neutralc:
@@ -752,6 +767,9 @@ def apply_name_scheme(self, forcefield_):
                                 # Remove the C/N prefix to keep the protonation state of the residue
                                 # in the terminal residues
                                 rname = rname[1:]
+                            elif rname.startswith("NEUTRAL-"):
+                                # Remove the NEUTRAL-C and NEUTRAL-N prefixes to keep protonation state
+                                rname = rname[9:]
                             else:
                                 # This is the old code which will overwrite the protonation state
                                 # of the residue but have wrong hydrogens. Not sure if needed but

tests/core_test.py

@@ -118,3 +118,27 @@ def test_protonated_terminals(input_pdb, expected_pqr, tmp_path):
         tmp_path=tmp_path,
         compare_resnames=True,
     )
+
+
+@pytest.mark.parametrize(
+    "input_pdb, expected_pqr",
+    [
+        pytest.param(
+            "5vav_cyclic_peptide.pdb",
+            "5vav_cyclic_peptide_out.pqr",
+            id="Cyclic peptide",
+        )
+    ],
+)
+def test_cyclic_peptide(input_pdb, expected_pqr, tmp_path):
+    """Tests for cyclic peptide protonation."""
+    args = "--log-level=INFO --ff=AMBER --ffout AMBER"
+    output_pqr = Path(input_pdb).stem + ".pqr"
+    common.run_pdb2pqr(
+        args=args,
+        input_pdb=common.DATA_DIR / input_pdb,
+        output_pqr=output_pqr,
+        expected_pqr=common.DATA_DIR / expected_pqr,
+        tmp_path=tmp_path,
+        compare_resnames=True,
+    )