Proteins with metal ions in the binding site

[agerexx]

Hi,
I have a problem when trying to calculate the fingerprint with proteins with metal ions in the binding site.
For example, I have a protein with manganese in the binding site. I loaded the protein with

u = mda.Universe("protein.mol2")
for i, bond_order in enumerate(u._topology.bonds.order):
    # you may need to replace double bonds ("2") as well
    if bond_order == "ar":
        u._topology.bonds.order[i] = 1
u._topology.bonds._cache.pop("bd", None)
protein_mol = plf.Molecule.from_mda(u)

and the poses with

poses_path = "poses.mol2"
pose_iterable = plf.mol2_supplier(poses_path)

I calculated the fingerprint with

fp = plf.Fingerprint()
fp.run_from_iterable(pose_iterable, protein_mol)

But I get the following error:

RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback: 
"""
Traceback (most recent call last):
  File "c:\Users\hello\.venv\Lib\site-packages\prolif\interactions\interactions.py", line 470, in detect
    vdw = self._vdw_cache[elements]
          ~~~~~~~~~~~~~~~^^^^^^^^^^
KeyError: frozenset({'O', 'Mn'})

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "c:\Users\hello\.venv\Lib\site-packages\multiprocess\pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
                    ^^^^^^^^^^^^^^^^^^^
  File "c:\Users\hello\.venv\Lib\site-packages\prolif\parallel.py", line 234, in executor
    return cls.fp.generate(mol, cls.pmol, residues=cls.residues, metadata=True)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "c:\Users\hello\.venv\Lib\site-packages\prolif\fingerprint.py", line 373, in generate
    interactions = get_interactions(lres, pres)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "c:\Users\hello\.venv\Lib\site-packages\prolif\fingerprint.py", line 311, in metadata
    if (metadata := interaction(res1, res2, metadata=True))
                    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
...
    [871](file:///C:/Users/hello/.venv/Lib/site-packages/multiprocess/pool.py:871) if success:
    [872](file:///C:/Users/hello/.venv/Lib/site-packages/multiprocess/pool.py:872)     return value
--> [873](file:///C:/Users/hello/.venv/Lib/site-packages/multiprocess/pool.py:873) raise value

KeyError: 'Mn'

I got the same error when trying to calculate the fingerprint for a protein with cobalt in the binding site.

I also have some proteins with an heme group with iron in the binding site and in these cases the calculation of the fingerprint just run for a very long time without giving any output (the tqdm bar doesn't appear) or any error. I stopped the calculation after 20 minutes without results.
In all the above cases, if I remove the metal ions there's no problem with the fingerprint, but it would be sensible for me to calculate the fingerprint with the metal ions, since the proteins had been used as they are for the docking and the metals are actually involved in interactions with the ligands. Since PROLIF is able to detect also MetalDonor and MetalAcceptor interactions, I don't know why I'm getting those errors.
Any help would be very appreciated! Thank you :)

[cbouy]

Hi,

First of, the default fingerprint does not calculate coordination with metal ions, you need to explicitly specify Metal interactions if you want them, otherwise you might only capture ionic interactions here. I'd suggest switching to

fp = plf.Fingerprint([
    'MetalAcceptor', 'MetalDonor', 'Hydrophobic',
    'HBAcceptor', 'HBDonor', 'Cationic', 'Anionic',
    'CationPi', 'PiCation', 'PiStacking', 'VdWContact',
])

As for the error message you see, it's caused by the VdWContact interaction that does not have access to the van der Waals radii for all elements, they are only implemented for a subset. I use the van der Waals radii as defined by MDAnalysis here, and any element that is not in there you'll have to provide to the VdWContact interaction:

fp = plf.Fingerprint(
    [
        'MetalAcceptor', 'MetalDonor', 'Hydrophobic',
        'HBAcceptor', 'HBDonor', 'Cationic', 'Anionic',
        'CationPi', 'PiCation', 'PiStacking', 'VdWContact',
    ],
    parameters={
        "VdWContact": {
            "vdwradii": {"Mn": 2.45}
        }
    }
)

Lastly, for your protein with the heme group not even showing the progress bar, it likely means that MDAnalysis is still trying to convert your protein (in that case I assume you've selected everything in your file) to an RDKit molecule, and that process is known to be slow. I would advise passing a selection of atoms around your ligand instead of the entire universe. Now because the ligands are in a separate file it complexifies the process a bit but essentially we can calculate the centroid of your first pose (or maybe you have a reference ligand) and then select all residues within a specified radius (here 15A) of that centroid:

ref = pose_iterable[0]
centroid = " ".join(f"{x:.2f}" for x in ref.centroid)
protein_selection = u.select_atoms(f"byres point {centroid} 15.0")
protein_mol = plf.Molecule.from_mda(protein_selection)

Hopefully with all this you should be able to calculate your fingerprint and it will be a bit faster!

[agerexx]

Thank you for you fast response. It totally solved my problem in both cases!