Fingerprinting only showing Van der Waals interactions

[diaghil]

I am trying to see if my fingerprint displays any interactions other than VdWContact. My code is pretty much referencing the documentation directly, and I know the PDB file contains other kinds of interactions.

u = mda.Universe(args.pdb, args.dcd, refresh_offsets=True)
lig_def = args.lig_asl
prot_def = args.prot_asl
lig = u.select_atoms(lig_def)
prot = u.select_atoms(prot_def)
lmol = plf.Molecule.from_mda(lig)
pmol = plf.Molecule.from_mda(prot)
fp = plf.Fingerprint()
fp.run(u.trajectory, lig, prot)

After removing HOH residues from the PDB input file and making sure none of the atom IDs in the PDB were faulty, I'm still failing to see a full range of results, even without the trajectory file being passed as an argument.

GLY13.A:VdWContact,VAL14.A:VdWContact,GLY15.A:VdWContact,LYS16.A:VdWContact,ALA18.A:VdWContact,PHE28.A:VdWContact,ASP30.A:VdWContact,LYS117.A:VdWContact,LEU120.A:VdWContact,SER145.A:VdWContact,ALA146.A:VdWContact,LYS147.A:VdWContact
0.609,0.014,0.864,0.633,0.038,0.248,0.048,0.192,0.034,0.172,0.118,0.008

[cbouy]

Hi @diaghil, does your PDB file contain explicit hydrogen atoms for both the protein and the ligand?

If it does, can you also check if there are bonds defined in the PDB file, and if so is it all bonds or only a subset? If it's a subset, you can add the guess_bonds=True parameter in mda.Universe and rerun your script.

[cbouy]

@diaghil I haven't received anything (also not in my spam folder), could you double check?

[diaghil]

@cbouy Let me know if you received my new email.

[cbouy]

I did, will have a look in the coming days

[diaghil]

Please let me know.

[cbouy]

Okay the problem was indeed with only having partial bonds, but I managed to solve it by calling guess_bonds as in here:

import prolif as plf
import MDAnalysis as mda
from MDAnalysis.topology.tables import vdwradii

u = mda.Universe("file.pdb")
lig = u.select_atoms("resname UNK")

# I restrict the protein selection to 20 around the ligand but for the trajectory analysis you could
# only keep the `not group ligand` for the selection in case the ligand moves around a lot
prot = u.select_atoms("(byres around 20 group ligand) and not group ligand", ligand=lig)

# guess bonds for the salts
prot.guess_bonds(
    vdwradii={
        "Na": vdwradii["NA"],
        "Mg": vdwradii["MG"],
        "Cl": vdwradii["CL"],
    })

# run
fp = plf.Fingerprint()
fp.run(u.trajectory, lig, prot)
df = fp.to_dataframe()

and you get those interaction types:

df.groupby(level="interaction", axis=1).sum()

Frame	Anionic	Cationic	HBAcceptor	HBDonor	Hydrophobic	PiStacking	VdWContact
0	2	2	13	5	14	1	35

Tell me if this works for you,
Cédric