Martinize2 fails on C-terminus canonicalization via KeyError due to dropped OXT atom on ARG 340

[Kostalab]

There seems to be an inconsistency with my martinize2 commands. the issue is with input pdbs.

I am encountering a KeyError when running martinize2. Our IT department helped me troubleshoot the problem, and it appears the script is failing during "canonicalization modifications" for a C-terminal residue.

The script drops the OXT atom (atom 1800) from my C-terminal ARG 340, which causes a KeyError later in the process.

Here is the summary from our IT support (email thread):

...Me:

For some reason, when I do my initial conversion from atomic detail level pdb structure to coarse grain structure, some of my conformers dont
work, and they break. as a perfect example;
conformer 1, BREAKS.
conformer 2 works perfectly. ... why?

Do you happen to know why this is happening? it isn't my script's fault (probably). I think this happened when I ran it outside of my script, which leads me to believe this is martinize2's fault. the coarse-grainifier process. Can it be fixed?

fyi; the flow is :
*martinize2 ---> BREAKS AFTER THIS FOR SOME REASON <----
*insane
*sed commands and edits
*em
*nvt
*npt
*md for 0.8 ns
*md_2 for 10 ns

...IT help:

Hi,

I have spent quite some time troubleshooting your issue, I am pretty certain it's not our system causing the issue.

Since I lack domain knowledge in this area, I can only make some guesses based on the logic of the source code. The error occurs during canonicalization modifications:

unknown-input - Could not identify the modifications for residues ['LEU338', 'ARG340'], involving atoms ['1800-OXT']

For some reason, the following line does not work:

ATOM 1800 OXT ARG B 340 150.701 -41.918 80.450 0.00 0.00 B O

But this one does:

ATOM 1800 OXT ARG B 340 152.014 -46.316 81.585 0.00 0.00 B O

With this error, atom 1800 is dropped from the input, and later you will encounter the key error because it was looking for this atom.

I am not sure what condition causes such modification to occur, but if you know what the correct modification should be, you might try passing it in
with the "-modify" option. If you are unsure as well, then you may need to contact the developer of martinize.

This seems to be a bug, as the script's behavior appears to be dependent on the atom's coordinates, which it shouldn't be.

**Could you please help me figure out why my input pdbs are inconsistently working.

Thank you.**

since my files are not uploading I will paste a preview of my terminal pdbs:
conformer 1:
...
ATOM 3565 C GLN A 337 146.905 -42.235 79.587 1.00 0.52 A C
ATOM 3566 O GLN A 337 146.559 -43.401 79.410 1.00 0.52 A O
ATOM 3567 CB GLN A 337 145.889 -40.837 77.840 1.00 0.52 A C
ATOM 3568 CG GLN A 337 144.611 -40.085 77.422 1.00 0.52 A C
ATOM 3569 CD GLN A 337 144.293 -40.287 75.946 1.00 0.52 A C
ATOM 3570 OE1 GLN A 337 144.951 -40.986 75.183 1.00 0.52 A O
ATOM 3571 NE2 GLN A 337 143.172 -39.673 75.504 1.00 0.52 A N
ATOM 3572 N LEU A 338 148.171 -41.911 79.925 1.00 0.70 A N
ATOM 3573 CA LEU A 338 149.168 -42.919 80.248 1.00 0.70 A C
ATOM 3574 C LEU A 338 148.943 -43.606 81.591 1.00 0.70 A C
ATOM 3575 O LEU A 338 148.676 -44.807 81.634 1.00 0.70 A O
ATOM 3576 CB LEU A 338 150.606 -42.314 80.219 1.00 0.70 A C
ATOM 3577 CG LEU A 338 151.110 -41.830 78.835 1.00 0.70 A C
ATOM 3578 CD1 LEU A 338 152.580 -41.381 78.937 1.00 0.70 A C
ATOM 3579 CD2 LEU A 338 150.964 -42.896 77.731 1.00 0.70 A C
ATOM 3580 N THR A 339 148.358 -42.414 81.499 0.00 0.00 A N
ATOM 3581 CA THR A 339 149.116 -41.187 81.728 0.00 0.00 A C
ATOM 3582 C THR A 339 149.570 -41.073 83.178 0.00 0.00 A C
ATOM 3583 O THR A 339 150.703 -40.671 83.452 0.00 0.00 A O
ATOM 3584 CB THR A 339 148.290 -39.938 81.362 0.00 0.00 A C
ATOM 3585 OG1 THR A 339 147.935 -39.988 79.975 0.00 0.00 A O
ATOM 3586 CG2 THR A 339 149.095 -38.673 81.620 0.00 0.00 A C
ATOM 3587 H THR A 339 147.492 -42.278 81.299 0.00 0.00 A H
ATOM 3588 N ARG A 340 148.679 -41.430 84.101 0.00 0.00 A N
ATOM 3589 CA ARG A 340 148.981 -41.387 85.529 0.00 0.00 A C
ATOM 3590 C ARG A 340 150.073 -42.377 85.917 0.00 0.00 A C
ATOM 3591 O ARG A 340 151.008 -42.018 86.697 0.00 0.00 A O
ATOM 3592 CB ARG A 340 147.721 -41.665 86.352 0.00 0.00 A C
ATOM 3593 CG ARG A 340 146.652 -40.592 86.233 0.00 0.00 A C
ATOM 3594 CD ARG A 340 145.425 -40.934 87.062 0.00 0.00 A C
ATOM 3595 NE ARG A 340 144.385 -39.913 86.948 0.00 0.00 A N
ATOM 3596 CZ ARG A 340 143.203 -39.971 87.555 0.00 0.00 A C
ATOM 3597 NH1 ARG A 340 142.903 -41.008 88.324 0.00 0.00 A N
ATOM 3598 NH2 ARG A 340 142.324 -38.993 87.389 0.00 0.00 A N
ATOM 3599 H ARG A 340 147.848 -41.717 83.910 0.00 0.00 A H
ATOM 3600 OXT ARG A 340 150.047 -43.562 85.459 0.00 0.00 A O
...

conformer 2:
...
ATOM 3565 C GLN A 337 146.905 -42.235 79.587 1.00 0.52 A C
ATOM 3566 O GLN A 337 146.559 -43.401 79.410 1.00 0.52 A O
ATOM 3567 CB GLN A 337 145.889 -40.837 77.840 1.00 0.52 A C
ATOM 3568 CG GLN A 337 144.611 -40.085 77.422 1.00 0.52 A C
ATOM 3569 CD GLN A 337 144.293 -40.287 75.946 1.00 0.52 A C
ATOM 3570 OE1 GLN A 337 144.951 -40.986 75.183 1.00 0.52 A O
ATOM 3571 NE2 GLN A 337 143.172 -39.673 75.504 1.00 0.52 A N
ATOM 3572 N LEU A 338 148.171 -41.911 79.925 1.00 0.70 A N
ATOM 3573 CA LEU A 338 149.168 -42.919 80.248 1.00 0.70 A C
ATOM 3574 C LEU A 338 148.943 -43.606 81.591 1.00 0.70 A C
ATOM 3575 O LEU A 338 148.676 -44.807 81.634 1.00 0.70 A O
ATOM 3576 CB LEU A 338 150.606 -42.314 80.219 1.00 0.70 A C
ATOM 3577 CG LEU A 338 151.110 -41.830 78.835 1.00 0.70 A C
ATOM 3578 CD1 LEU A 338 152.580 -41.381 78.937 1.00 0.70 A C
ATOM 3579 CD2 LEU A 338 150.964 -42.896 77.731 1.00 0.70 A C
ATOM 3580 N THR A 339 147.620 -43.716 81.681 0.00 0.00 A N
ATOM 3581 CA THR A 339 146.749 -42.871 80.870 0.00 0.00 A C
ATOM 3582 C THR A 339 146.010 -41.910 81.794 0.00 0.00 A C
ATOM 3583 O THR A 339 145.387 -42.350 82.764 0.00 0.00 A O
ATOM 3584 CB THR A 339 145.749 -43.712 80.050 0.00 0.00 A C
ATOM 3585 OG1 THR A 339 146.468 -44.580 79.166 0.00 0.00 A O
ATOM 3586 CG2 THR A 339 144.843 -42.809 79.227 0.00 0.00 A C
ATOM 3587 H THR A 339 147.179 -44.287 82.219 0.00 0.00 A H
ATOM 3588 N ARG A 340 146.093 -40.622 81.474 0.00 0.00 A N
ATOM 3589 CA ARG A 340 145.315 -39.638 82.224 0.00 0.00 A C
ATOM 3590 C ARG A 340 143.829 -39.974 82.186 0.00 0.00 A C
ATOM 3591 O ARG A 340 143.448 -41.133 81.837 0.00 0.00 A O
ATOM 3592 CB ARG A 340 145.553 -38.232 81.671 0.00 0.00 A C
ATOM 3593 CG ARG A 340 146.962 -37.707 81.890 0.00 0.00 A C
ATOM 3594 CD ARG A 340 147.137 -36.317 81.304 0.00 0.00 A C
ATOM 3595 NE ARG A 340 148.492 -35.806 81.501 0.00 0.00 A N
ATOM 3596 CZ ARG A 340 148.930 -34.637 81.042 0.00 0.00 A C
ATOM 3597 NH1 ARG A 340 148.117 -33.848 80.353 0.00 0.00 A N
ATOM 3598 NH2 ARG A 340 150.181 -34.260 81.275 0.00 0.00 A N
ATOM 3599 H ARG A 340 146.596 -40.266 80.818 0.00 0.00 A H
ATOM 3600 OXT ARG A 340 142.971 -39.092 82.504 0.00 0.00 A O
...

[Kostalab]

here is a better look at my command:

martinize2 -f conformer_1.pdb -o topol_1.top -x coarsegrained_1.pdb -ff martini3IDP -elastic -ef 700 -el 0.3 -eu 0.85 -eunit 131:209 -dssp -p backbone -pf 500
    INFO - general - Read 2 molecules from PDB file conformer_1.pdb
    INFO - step - Guessing the bonds.
    INFO - general - 1 molecules after guessing bonds
    INFO - step - Repairing the graph.
    INFO - general - Applying modification N-ter to residue B-GLY1
    INFO - general - Applying modification + to residue B-GLY1
    INFO - general - Applying modification N-ter to residue A-GLY1
    INFO - step - Dealing with modifications.
    INFO - general - Identified the modifications ['N-ter', 'N-ter'] on residues ['GLY1', 'GLY1', 'GLY1', 'GLY1']
 WARNING - unknown-input - Could not identify the modifications for residues ['LEU338', 'THR339', 'ARG340'], involving atoms ['1800-OXT']
    INFO - general - Identified the modifications ['C-ter'] on residues ['ARG340']
Traceback (most recent call last):
  File "/home/kvlasaka/miniconda3/envs/md-intro-tutorial/bin//martinize2", line 1205, in <module>
    entry()
  File "/home/kvlasaka/miniconda3/envs/md-intro-tutorial/bin//martinize2", line 916, in entry
    system = pdb_to_universal(
             ^^^^^^^^^^^^^^^^^
  File "/home/kvlasaka/miniconda3/envs/md-intro-tutorial/bin//martinize2", line 157, in pdb_to_universal
    vermouth.CanonicalizeModifications().run_system(canonicalized)
  File "/home/kvlasaka/miniconda3/envs/md-intro-tutorial/lib/python3.12/site-packages/vermouth/processors/processor.py", line 37, in run_system
    mols.append(self.run_molecule(molecule))
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/kvlasaka/miniconda3/envs/md-intro-tutorial/lib/python3.12/site-packages/vermouth/processors/canonicalize_modifications.py", line 380, in run_molecule
    fix_ptm(molecule)
  File "/home/kvlasaka/miniconda3/envs/md-intro-tutorial/lib/python3.12/site-packages/vermouth/processors/canonicalize_modifications.py", line 363, in fix_ptm
    node = molecule.nodes[n_idx]
           ~~~~~~~~~~~~~~^^^^^^^
  File "/home/kvlasaka/miniconda3/envs/md-intro-tutorial/lib/python3.12/site-packages/networkx/classes/reportviews.py", line 196, in __getitem__
    return self._nodes[n]
           ~~~~~~~~~~~^^^
KeyError: 1799

[Kostalab]

I should also mention that my command above, for some reason incorrectly merges my two chains into one, which causes me to have to write another separate script that splits them apart (as I have 20,000 conformers to process). This error seems to be caused by the OXT atoms ;

Read 2 molecules from PDB file conformer_1.pdb
...
1 molecules after guessing bonds

it seems to be a guessing the bonds step that messes up.

lastly, I have tried using the -sep flag, and it doesnt seem to work on my pdbs.

[jochuan]

I also have a similar problem where in the step of guessing bonds martinize detects fewer chains than the total

[pckroon]

Hello hello!
Thanks for the extensive report (and kudos to your IT department for digging into this as well!).
I have seen issues before where Martini2 generates KeyErrors rather than printing a nice and friendly error message, but I've not been able to chase down the root cause of it so far.

This seems to be a bug, as the script's behavior appears to be dependent on the atom's coordinates, which it shouldn't be.

This is absolutely not a bug, but a feature. Martinize2 guesses bonds based on interatomic distances (actually using the same algorithm as VMD), and uses the found bonds/connections/structure to actually identify which atom is which. This reduces the reliance on atom names.

In your case it means that the coordinate of your OXT is suspect, in that it is close enough to LEU338 and THR339 to create bonds with them. This usually means you have some pretty bad clashes, and any resulting simulation will be non-physical.

Martinize2 has a few debug flags you can use to write out the PDB file with CONECT records as M2 sees it, which may help your further investigation. In particular, take a look at -write-graph (which creates a pdb file immediately after guessing bonds) and -write-repair (which creates a pdb file after identifying non-modification atoms).

Dropping the offending atom may be a viable workaround, but make sure to double check you end up with a charged C-terminus.

I'll try to make some time to dig into the code for this again with your conformers as debug cases, but the goal there would be to fix the KeyError, not change the fundamental behaviour.

[csbrasnett]

The other suggestion along these lines could be to do a short atomistic simulation beforehand to energy minimise the structure before martinizing. That might help correct the interatomic distances and fix the bond guessing

[Kostalab]

Thank you for your help! I didn't think of doing a short atomistic simulation beforehand, but that may be a good solution. I am leaning towards simply deleting the atom, because it ran perfectly fine, with my entire pathway - no problem when I did that. The last issue I have is that it still stitched together the side chain A, with the side chain B, which i believe is still a problem even after deleting atom 1800. The fix to that is to simply renumber and re-chain them either by residue number, or other sequence specific modes.

I am using a collaborators software (IDPConfgen) from Julie-Foreman-Kay's lab - to create conformers of my protein, and it uses the whole Lennard Jones potential when calculating the geometric steric clashes, as opposed to hard spheres. This allows for some small overlap between the spheres, which may contribute to the errors I am seeing (or at least the inconsistent steric problems I am facing).

[pckroon]

Thank you for your help! I didn't think of doing a short atomistic simulation beforehand, but that may be a good solution. I am leaning towards simply deleting the atom, because it ran perfectly fine, with my entire pathway - no problem when I did that. The last issue I have is that it still stitched together the side chain A, with the side chain B, which i believe is still a problem even after deleting atom 1800. The fix to that is to simply renumber and re-chain them either by residue number, or other sequence specific modes.

Is this directly after the guess bonds step? That prints the number of molecules (and therefor itp files) left. Martinize2 merges chains if there's a bond between them, as with a disulfide bridge.

I am using a collaborators software (IDPConfgen) from Julie-Foreman-Kay's lab - to create conformers of my protein, and it uses the whole Lennard Jones potential when calculating the geometric steric clashes, as opposed to hard spheres. This allows for some small overlap between the spheres, which may contribute to the errors I am seeing (or at least the inconsistent steric problems I am facing).

Martinize2 uses the van der Waals radius of the atoms. A bonds gets made if distance < (vdw[atom1] + vdw[atom2])*1.2, and neither atom1 nor atom2 are hydrogens, and atom1 and atom2 are in separate residues. You can tune the factor 1.2 with the -bonds-fudge cli option.

You can find the vdw radii here: https://github.com/marrink-lab/vermouth-martinize/blob/master/vermouth/processors/make_bonds.py#L35-L62