Repository: marrink-lab/vermouth-martinize

Issue:
DSSP error (duplicated nr. of residues)
I'm getting the following error. The protein has 156 residues. So it appears that at some step, they are being incorrectly duplciated.

```
Traceback (most recent call last):
  File "/lustre1/work/courtade/software/miniconda3/envs/py37/bin/martinize2", line 669, in <module>
    entry()
  File "/lustre1/work/courtade/software/miniconda3/envs/py37/bin/martinize2", line 545, in entry
    AnnotateDSSP(executable=args.dssp, savedir='.').run_system(system)
  File "/lustre1/work/courtade/software/miniconda3/envs/py37/lib/python3.7/site-packages/vermouth/processors/processor.py", line 38, in run_system
    mols.append(self.run_molecule(molecule))
  File "/lustre1/work/courtade/software/miniconda3/envs/py37/lib/python3.7/site-packages/vermouth/dssp/dssp.py", line 444, in run_molecule
    annotate_dssp(molecule, self.executable, self.savedir)
  File "/lustre1/work/courtade/software/miniconda3/envs/py37/lib/python3.7/site-packages/vermouth/dssp/dssp.py", line 277, in annotate_dssp
    annotate_residues_from_sequence(molecule, attribute, secstructs)
  File "/lustre1/work/courtade/software/miniconda3/envs/py37/lib/python3.7/site-packages/vermouth/dssp/dssp.py", line 390, in annotate_residues_from_sequence
    raise ValueError(msg.format(len(sequence), len(residues)))
ValueError: The sequence length does not match the number of residues. The sequence has 156 elements for 312 residues.
```
I tried with a different PDB file, it throws the same error.

Edit: The error was triggered by reading a PDB without hydrogens (incorrectly defined termini?). I managed to solve it by running pdb2gmx on the all-atom coordinates before martinizing

Discussion:

Patch:
diff --git a/vermouth/processors/repair_graph.py b/vermouth/processors/repair_graph.py
index b95db34..1e4ea66 100644
--- a/vermouth/processors/repair_graph.py
+++ b/vermouth/processors/repair_graph.py
@@ -193,7 +193,15 @@ def make_reference(mol):
         residue = residues.nodes[residx]['graph']
         reference = _get_reference_residue(residues.nodes[residx], mol.force_field)
         if 'mutation' in residues.nodes[residx]:
+            # We need to do this, since we need to propagate all common
+            # attributes in the input residue to any reconstructed atoms. 
+            # make_residue_graph already picks up on all the common attributes
+            # and sets those as residue/node attributes, but any mutation needs
+            # to be propagated as well, otherwise the newly constructed atoms
+            # will end up with the old residue name. This way, we can just add
+            # *all* attributes in the residue graph to the new atoms.
             resname = residues.nodes[residx]['mutation'][0]
+            residues.nodes[residx]['resname'] = resname
         add_element_attr(reference)
         add_element_attr(residue)
         # We are going to sort the nodes of reference and residue by atomname.
@@ -273,9 +281,12 @@ def make_reference(mol):
         # "unsort" the matches
         match = {old_ref_names[ref]: old_res_names[res] for ref, res in match.items()}
 
-        reference_graph.add_node(residx, chain=chain, reference=reference,
-                                 found=residue, resname=resname, resid=resid,
-                                 match=match)
+        # The residue graph has attributes which are common to all atoms in the
+        # input residue, so propagate those to new atoms. These attributes are
+        # things like residue name (see above in case of mutations), resid, 
+        # insertion code, etc.
+        reference_graph.add_node(residx, reference=reference, found=residue, 
+                                 match=match, **residues.nodes[residx])
     reference_graph.add_edges_from(residues.edges())
     return reference_graph
 


Test patch:
diff --git a/vermouth/tests/test_repair_graph.py b/vermouth/tests/test_repair_graph.py
index f73b01a..e902d91 100644
--- a/vermouth/tests/test_repair_graph.py
+++ b/vermouth/tests/test_repair_graph.py
@@ -111,10 +111,11 @@ def build_system_mod(force_field):
         (16, {'resid': 1, 'resname': 'GLU', 'atomname': 'HN', 'chain': 'A', 'element': 'H'}),
         # Residue 2 has missing atoms (HA1 and HA2, the hydrogen of the C
         # alpha; HN the hydrogen on the N-ter nitrogen and O, the oxygen at the
-        # C-ter)
-        (17, {'resid': 2, 'resname': 'GLY', 'atomname': 'N', 'chain': 'A', 'element': 'N'}),
-        (18, {'resid': 2, 'resname': 'GLY', 'atomname': 'CA', 'chain': 'A', 'element': 'C'}),
-        (19, {'resid': 2, 'resname': 'GLY', 'atomname': 'C', 'chain': 'A', 'element': 'C'}),
+        # C-ter). 'common' is an arbitrary attribute name with identical values,
+        # which should be propagated to reconstructed atoms.
+        (17, {'resid': 2, 'resname': 'GLY', 'atomname': 'N', 'chain': 'A', 'element': 'N', 'common': 'a'}),
+        (18, {'resid': 2, 'resname': 'GLY', 'atomname': 'CA', 'chain': 'A', 'element': 'C', 'common': 'a'}),
+        (19, {'resid': 2, 'resname': 'GLY', 'atomname': 'C', 'chain': 'A', 'element': 'C', 'common': 'a'}),
         # Residue 3 has only N and C specified, everything else is missing.
         (20, {'resid': 3, 'resname': 'GLY', 'atomname': 'N', 'chain': 'A', 'element': 'N'}),
         (21, {'resid': 3, 'resname': 'GLY', 'atomname': 'O', 'chain': 'A', 'element': 'O'}),
@@ -243,6 +244,21 @@ def test_name_canonicalized(canonicalized_graph, key, expected_names):
     assert molecule.nodes[key]['atomname'] in expected_names
 
 
+@pytest.mark.parametrize('resid, expected_attrs', [
+    (2, {'resid': 2, 'resname': 'GLY', 'chain': 'A', 'common': 'a'})
+])
+def test_common_attributes(repaired_graph, resid, expected_attrs):
+    """Test that attributes that are common to all nodes in a residue end up in
+    reconstructed nodes."""
+    mol = repaired_graph.molecules[0]
+    for node_idx in mol:
+        node = mol.nodes[node_idx]
+        if node.get('resid') == resid:
+            for key, val in expected_attrs.items():
+                assert key in node
+                assert node[key] == val
+
+
 def test_renaming(renamed_graph):
     for node in renamed_graph.molecules[0].nodes.values():
         if node['resid'] == 1: