Multiple Residues > 1 ligand?

[ONISniperOne]

Hi! I'm importing a topology (pdb... prmtop is also available and works) and trajectory (nc) into ProLif

Amber splits my ligand into 2 residues even though they're attached to each other. Is there a way to keep them attached as one and also running intramolecular interactions within the ligand as well as intermolecular between ligand and protein at the same time?

To show the interaction within the ligand, I'm currently doing the below and treating half of the ligand like the protein when it isn't one. My whole ligand is resid 209 + resid 210

import MDAnalysis as mda
import prolif as plf

# Load topology and trajectory
u = mda.Universe('/home/analog_OTP8_3_0/cocomplex_nowat_noion.pdb', '/home/analog_OTP8_3_0/10.md.nc')
# Create selections for the ligand and protein
ligand_selection = u.select_atoms("resid 210")
protein_selection = u.select_atoms("resid 1:209 and byres around 15.0 group ligand", ligand=ligand_selection)

It shows the intramolecular Hbond within the ligand (between resid 209 and 210), but they don't seem to be physically connected. Since 209 portion is now being treated like the protein (toget intramolecular bonds), it's also colored gray instead of blue.

Yellow highlight where the VDW bond is where the two halves of the ligand are supposed to be linked. Blue circle is correct. It's the H bond between the 2 halves

# Create selections for the ligand and protein
ligand_selection = u.select_atoms("not protein")
# Equivalent of ligand_selection = u.select_atoms("resid 209:210")
protein_selection = u.select_atoms("protein and byres around 15.0 group ligand", ligand=ligand_selection)
# Equivalent of protein_selection = u.select_atoms("resid 1:208 and byres around 15.0 group ligand", ligand=ligand_selection)

The code above gives the below... if I put the ligand together (resid 209:210). The intramolecular Hbond in the ligand is missing, but at least the ligand is whole.

And ofc, the 3D one is actually missing the bond whiel the 2D one doesn't (yellow highlight)

[cbouy]

Hi,

You can adjust the residue information on the MDAnalysis object directly to make it as 1 residue again, which should solve most of your problems:

ligand_selection = u.select_atoms("resid 209:210")
ligand_selection.residues.resnames = ["LIG"] * len(ligand_selection.residues)
ligand_selection.residues.resids = [209] * len(ligand_selection.residues)

Now your ligand will be labelled LIG209 and has only a single residue.

To calculate intramolecular interactions in one go, you can just add the ligand group to your protein selection and it will calculate the ligand-ligand interactions as if it were 2 separate entities:

protein_selection = u.select_atoms(
    "(protein and byres around 15.0 group ligand) or (group ligand)", ligand=ligand_selection
)

It should be displayed correctly on the 3D plot, but for the 2D plot it will display them as if they were 2 separate things i.e. one ligand in full details and the same ligand but just as a residue label and the interactions between those 2 objects (so not directly on the atomic-detail ligand unfortunately).

Hope this helps, feel free to reopen the issue if this doesn't answer your question!

[ONISniperOne]

Hi,

You can adjust the residue information on the MDAnalysis object directly to make it as 1 residue again, which should solve most of your problems:

ligand_selection = u.select_atoms("resid 209:210")
ligand_selection.residues.resnames = ["LIG"] * len(ligand_selection.residues)
ligand_selection.residues.resids = [209] * len(ligand_selection.residues)

Now your ligand will be labelled LIG209 and has only a single residue.

To calculate intramolecular interactions in one go, you can just add the ligand group to your protein selection and it will calculate the ligand-ligand interactions as if it were 2 separate entities:

protein_selection = u.select_atoms(
    "(protein and byres around 15.0 group ligand) or (group ligand)", ligand=ligand_selection
)

It should be displayed correctly on the 3D plot, but for the 2D plot it will display them as if they were 2 separate things i.e. one ligand in full details and the same ligand but just as a residue label and the interactions between those 2 objects (so not directly on the atomic-detail ligand unfortunately).

Hope this helps, feel free to reopen the issue if this doesn't answer your question!

Thanks! It kinda works? The 3D one is perfect. The hbonds are fully present and the ligand is connected correctly.

As you have noted, the 2D is... off. Rip, still at least the 3D is pretty and that any intramolecular interaction is there. We can always pinpoint when it happens with

%matplotlib ipympl
fp.plot_barcode()