Repository: marrink-lab/vermouth-martinize

Issue:
Martini 2.2 topologies for alpha-helical peptides are still incorrect
**Context**
This issue is related to Issue #560. Thank you for addressing the main problems mentioned in the issue. However, not all problems of the opened issue have been fixed.

**Short description**
While bead types for charged termini backbone beads of alpha-helical peptides in Martini 2.2 have been fixed (thank you for that!), the `vermouth-martinize` script treats the backbone beads _near_ the termini incorrectly.

**Detailed description**
The current dev version of `vermouth-martinize` (version 0.9.7_dev178) does not assign correct bead types to backbone beads of amino acids that are _near_ the peptide termini.

Using the command `martinize2 -f l20.pdb -o topol.top -x l20_martini.pdb -ss HHHHHHHHHHHHHHHHHHHH -p backbone -ff martini22`, `martinize2` generates the following `[ atoms ]` block:
```
[ atoms ]
 1 Qd   1 LEU BB   1   1
 2 AC1  1 LEU SC1  2 0.0
 3 N0   2 LEU BB   3 0.0
 4 AC1  2 LEU SC1  4 0.0
 5 N0   3 LEU BB   5 0.0
 6 AC1  3 LEU SC1  6 0.0
 7 N0   4 LEU BB   7 0.0
 8 AC1  4 LEU SC1  8 0.0
 9 N0   5 LEU BB   9 0.0
10 AC1  5 LEU SC1 10 0.0
11 N0   6 LEU BB  11 0.0
12 AC1  6 LEU SC1 12 0.0
13 N0   7 LEU BB  13 0.0
14 AC1  7 LEU SC1 14 0.0
15 N0   8 LEU BB  15 0.0
16 AC1  8 LEU SC1 16 0.0
17 N0   9 LEU BB  17 0.0
18 AC1  9 LEU SC1 18 0.0
19 N0  10 LEU BB  19 0.0
20 AC1 10 LEU SC1 20 0.0
21 N0  11 LEU BB  21 0.0
22 AC1 11 LEU SC1 22 0.0
23 N0  12 LEU BB  23 0.0
24 AC1 12 LEU SC1 24 0.0
25 N0  13 LEU BB  25 0.0
26 AC1 13 LEU SC1 26 0.0
27 N0  14 LEU BB  27 0.0
28 AC1 14 LEU SC1 28 0.0
29 N0  15 LEU BB  29 0.0
30 AC1 15 LEU SC1 30 0.0
31 N0  16 LEU BB  31 0.0
32 AC1 16 LEU SC1 32 0.0
33 N0  17 LEU BB  33 0.0
34 AC1 17 LEU SC1 34 0.0
35 N0  18 LEU BB  35 0.0
36 AC1 18 LEU SC1 36 0.0
37 N0  19 LEU BB  37 0.0
38 AC1 19 LEU SC1 38 0.0
39 Qa  20 LEU BB  39  -1
40 AC1 20 LEU SC1 40 0.0
```

Meanwhile, using `martinize.py` script version 2.6_3 (`python martinize.py -f l20.pdb -o topol.top -x l20_martini.pdb -ss HHHHHHHHHHHHHHHHHHHH -p backbone`) generates this atom block:
```
[ atoms ]
    1    Qd     1   LEU    BB     1  1.0000 ; 1
    2   AC1     1   LEU   SC1     2  0.0000 ; 1
    3    Nd     2   LEU    BB     3  0.0000 ; 1
    4   AC1     2   LEU   SC1     4  0.0000 ; 1
    5    Nd     3   LEU    BB     5  0.0000 ; 1
    6   AC1     3   LEU   SC1     6  0.0000 ; 1
    7    Nd     4   LEU    BB     7  0.0000 ; 1
    8   AC1     4   LEU   SC1     8  0.0000 ; 1
    9    N0     5   LEU    BB     9  0.0000 ; H
   10   AC1     5   LEU   SC1    10  0.0000 ; H
   11    N0     6   LEU    BB    11  0.0000 ; H
   12   AC1     6   LEU   SC1    12  0.0000 ; H
   13    N0     7   LEU    BB    13  0.0000 ; H
   14   AC1     7   LEU   SC1    14  0.0000 ; H
   15    N0     8   LEU    BB    15  0.0000 ; H
   16   AC1     8   LEU   SC1    16  0.0000 ; H
   17    N0     9   LEU    BB    17  0.0000 ; H
   18   AC1     9   LEU   SC1    18  0.0000 ; H
   19    N0    10   LEU    BB    19  0.0000 ; H
   20   AC1    10   LEU   SC1    20  0.0000 ; H
   21    N0    11   LEU    BB    21  0.0000 ; H
   22   AC1    11   LEU   SC1    22  0.0000 ; H
   23    N0    12   LEU    BB    23  0.0000 ; H
   24   AC1    12   LEU   SC1    24  0.0000 ; H
   25    N0    13   LEU    BB    25  0.0000 ; H
   26   AC1    13   LEU   SC1    26  0.0000 ; H
   27    N0    14   LEU    BB    27  0.0000 ; H
   28   AC1    14   LEU   SC1    28  0.0000 ; H
   29    N0    15   LEU    BB    29  0.0000 ; H
   30   AC1    15   LEU   SC1    30  0.0000 ; H
   31    N0    16   LEU    BB    31  0.0000 ; H
   32   AC1    16   LEU   SC1    32  0.0000 ; H
   33    Na    17   LEU    BB    33  0.0000 ; 2
   34   AC1    17   LEU   SC1    34  0.0000 ; 2
   35    Na    18   LEU    BB    35  0.0000 ; 2
   36   AC1    18   LEU   SC1    36  0.0000 ; 2
   37    Na    19   LEU    BB    37  0.0000 ; 2
   38   AC1    19   LEU   SC1    38  0.0000 ; 2
   39    Qa    20   LEU    BB    39 -1.0000 ; 2
   40   AC1    20   LEU   SC1    40  0.0000 ; 2
```

While `martinize2` now generates correct bead types for the backbone beads of the N-terminus and the C-terminus (Qd and Qa, respectively), it does not properly modify the backbone beads of the amino acids near the termini (4 amino acids from the N/C-terminus, **including the first/last residue, if it is uncharged**). The backbone beads of leucines near the N-terminus should be of type Nd, while backbone beads of leucines near the C-terminus should be of type Na. Currently, `martinize2` assigns all leucine backbone beads in helices the type N0.

Refer to article https://doi.org/10.1021/ct700324x, specifically Table 2 for further details.

`martinize.py` achieves the correct behavior by defining three additional alpha-helical "structures" denoted as 1, 2, and 3 and replacing the helical elements at the start and end of the helix with them. Excerpt from the code of the `martinize.py` script:

```
# Helix start and end regions are special and require assignment of
# specific types. The following pattern substitutions are applied
# (in the given order). A dot matches any other type.
# Patterns can be added to the dictionaries. This only makes sense
# if for each key in patterns there is a matching key in pattypes.
patterns = {
    "H": pat(".H. .HH. .HHH. .HHHH. .HHHHH. .HHHHHH. .HHHHHHH. .HHHH HHHH.")                #@#
}
pattypes = {
    "H": pat(".3. .33. .333. .3333. .13332. .113322. .1113222. .1111 2222.")                #@#
}
```

Note that structure `3` is used for overlapping segments (that are close both to the N-terminus and the C-terminus, i.e. it is used in short helices). `Nda` bead type is used for the backbone beads in structure 3.

Also note that Gly, Pro, and Ala should be assigned yet different bead types in these terminal segments (as shown in the Table 2 in the article linked above).

Discussion:

Patch:
diff --git a/vermouth/dssp/dssp.py b/vermouth/dssp/dssp.py
index 79f89ed..b9c0998 100644
--- a/vermouth/dssp/dssp.py
+++ b/vermouth/dssp/dssp.py
@@ -417,8 +417,15 @@ def convert_dssp_to_martini(sequence):
     cg_sequence = ''.join(ss_cg[secstruct] for secstruct in sequence)
     wildcard_sequence = ''.join('H' if secstruct == 'H' else '.'
                                 for secstruct in cg_sequence)
+    # Flank the sequence with dots. Otherwise in a sequence consisting of only
+    # H will not have a start or end. See also issue 566.
+    # This should not cause further issues, since '..' doesn't map to anything
+    wildcard_sequence = '.' + wildcard_sequence + '.'
     for pattern, replacement in patterns.items():
-        wildcard_sequence = wildcard_sequence.replace(pattern, replacement)
+        while pattern in wildcard_sequence:  # EXPENSIVE! :'(
+            wildcard_sequence = wildcard_sequence.replace(pattern, replacement)
+    # And remove the flanking dots again
+    wildcard_sequence = wildcard_sequence[1:-1]
     result = ''.join(
         wildcard if wildcard != '.' else cg
         for wildcard, cg in zip(wildcard_sequence, cg_sequence)


Test patch:
diff --git a/vermouth/tests/test_dssp.py b/vermouth/tests/test_dssp.py
index 89c52fe..28474ff 100644
--- a/vermouth/tests/test_dssp.py
+++ b/vermouth/tests/test_dssp.py
@@ -680,3 +680,18 @@ def test_cterm_atomnames():
     vermouth.processors.CanonicalizeModifications().run_system(system)
     dssp_out = dssp.run_dssp(system, executable=DSSP_EXECUTABLE)
     assert dssp_out == list("CC")
+
+
+@pytest.mark.parametrize('sequence, expected', [
+    ('H', '3'),
+    ('HH', '33'),
+    ('CHH', 'C33'),
+    ('HHHHHH', '113322'),
+    ('EHHHHHHC', 'E113322C'),
+    ('HHHHHHHHH', '1111H2222'),
+    ('CHHHHHHHHHC', 'C1111H2222C'),
+    ('CHHHHEHHHHC', 'C3333E3333C'),
+])
+def test_convert_dssp_to_martini(sequence, expected):
+    found = dssp.convert_dssp_to_martini(sequence)
+    assert expected == found