Create barcode for protein-protein residue-residue interaction from single-frame PDB

[biomolsim]

Hello, I am unable to follow through with the documentation to create a barcode for a single-frame PDB file. Could anyone please share a basic template to get me started? Thanks a lot.

[cbouy]

Hi, it's hard to help without knowing which step you are getting stuck on or what error message you are seeing, please provide more details

[biomolsim]

I am currently facing the following error for both structures within a single-frame PDB file. The n_residues function shows the correct number of residues for both of them, and also just supplying the name shows the entire structure.

My ultimate aim is to create a 2-by-2 plot that shows the presence of interaction across all the residues between the two structures, if any.

Code:
fp.run_from_iterable(lig, rest)

Error:
AttributeError: 'Residue' object has no attribute 'residues'

[cbouy]

How did you load lig and rest? Please provide a complete code snippet and the full error trace

[biomolsim]

Thanks @cbouy. Please find them below.

Code:

import MDAnalysis as mda
import prolif as plf

structure_file = "plf_frame_0.pdb"

u_read = mda.Universe(structure_file)

lig_sel = u_read.select_atoms("segid 3N1")
rest_sel = u_read.select_atoms("not (segid 3N1)")

lig_mol = plf.Molecule.from_mda(lig_sel)
rest_mol = plf.Molecule.from_mda(rest_sel)

fp = plf.Fingerprint()

fp.run_from_iterable(lig_mol, rest_mol)

Log:

/homes/user/.local/lib/python3.10/site-packages/MDAnalysis/topology/tables.py:52: DeprecationWarning: Deprecated in version 2.8.0
MDAnalysis.topology.tables has been moved to MDAnalysis.guesser.tables. This import point will be removed in MDAnalysis version 3.0.0
  warnings.warn(wmsg, category=DeprecationWarning)
/homes/user/.local/lib/python3.10/site-packages/MDAnalysis/coordinates/PDB.py:453: UserWarning: 1 A^3 CRYST1 record, this is usually a placeholder. Unit cell dimensions will be set to None.
  warnings.warn("1 A^3 CRYST1 record,"
/homes/user/.local/lib/python3.10/site-packages/MDAnalysis/converters/RDKit.py:497: UserWarning: No `bonds` attribute in this AtomGroup. Guessing bonds based on atoms coordinates
  warnings.warn(
 0/? [00:00<?, ?it/s]
---------------------------------------------------------------------------
RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback: 
"""
Traceback (most recent call last):
  File "/homes/user/.local/lib/python3.10/site-packages/multiprocess/pool.py", line 125, in worker
    result = (True, func(*args, **kwds))
  File "/homes/user/.local/lib/python3.10/site-packages/prolif/parallel.py", line 234, in executor
    return cls.fp.generate(mol, cls.pmol, residues=cls.residues, metadata=True)
  File "/homes/user/.local/lib/python3.10/site-packages/prolif/fingerprint.py", line 365, in generate
    for lresid, lres in lig.residues.items():
AttributeError: 'Residue' object has no attribute 'residues'
"""

The above exception was the direct cause of the following exception:

AttributeError                            Traceback (most recent call last)
Cell In [2], line 25
     19 rest_mol = plf.Molecule.from_mda(rest_sel)
     23 fp = plf.Fingerprint()
---> 25 fp.run_from_iterable(lig_mol, rest_mol)

File ~/.local/lib/python3.10/site-packages/prolif/fingerprint.py:600, in Fingerprint.run_from_iterable(self, lig_iterable, prot_mol, residues, progress, n_jobs)
    598     residues = list(prot_mol.residues)
    599 if n_jobs != 1:
--> 600     return self._run_iter_parallel(
    601         lig_iterable=lig_iterable,
    602         prot_mol=prot_mol,
    603         residues=residues,
    604         progress=progress,
    605         n_jobs=n_jobs,
    606     )
    608 iterator = tqdm(lig_iterable) if progress else lig_iterable
    609 ifp = {}

File ~/.local/lib/python3.10/site-packages/prolif/fingerprint.py:633, in Fingerprint._run_iter_parallel(self, lig_iterable, prot_mol, residues, progress, n_jobs)
    624 ifp = {}
    626 with MolIterablePool(
    627     n_jobs,
    628     fingerprint=self,
   (...)
    631     tqdm_kwargs={"total": total, "disable": not progress},
    632 ) as pool:
--> 633     for i, ifp_data in enumerate(pool.process(lig_iterable)):
    634         ifp[i] = ifp_data
    636 self.ifp = ifp

File ~/.local/lib/python3.10/site-packages/tqdm/notebook.py:250, in tqdm_notebook.__iter__(self)
    248 try:
    249     it = super().__iter__()
--> 250     for obj in it:
    251         # return super(tqdm...) will not catch exception
    252         yield obj
    253 # NB: except ... [ as ...] breaks IPython async KeyboardInterrupt

File ~/.local/lib/python3.10/site-packages/tqdm/std.py:1181, in tqdm.__iter__(self)
   1178 time = self._time
   1180 try:
-> 1181     for obj in iterable:
   1182         yield obj
   1183         # Update and possibly print the progressbar.
   1184         # Note: does not call self.update(1) for speed optimisation.

File ~/.local/lib/python3.10/site-packages/multiprocess/pool.py:873, in IMapIterator.next(self, timeout)
    871 if success:
    872     return value
--> 873 raise value

AttributeError: 'Residue' object has no attribute 'residues'

[cbouy]

run_from_iterable expects the first argument to be an iterable of ligand poses, not an individual molecule, so just wrap your lig_mol in a list:

fp.run_from_iterable([lig_mol], rest_mol)

[biomolsim]

Thanks a lot, this worked. I would still like your help with the following questions regarding data input and output visualization.

(1) Is it possible to use a DCD file (NAMD generated) instead of a PDB file?

(2) Would the above code work even for a multi-frame PDB or DCD file?

(3) Could you please suggest how I could generate a residue-residue interaction map between these two selections, like the one in the attached image?

[cbouy]

1. DCD is supported as a coordinate file format by MDAnalysis, so I believe you'll need the topology file associated with it: https://userguide.mdanalysis.org/stable/formats/index.html.
2. At that stage it's best to use fp.run if you have a "trajectory" rather than loading each frame in a separate analysis, see this tutorial: https://prolif.readthedocs.io/en/stable/notebooks/md-protein-protein.html#protein-protein-md
3. ProLIF has something similar but not exactly as seen in the above tutorial section, i.e it plots each pair of residue per trajectory/pose frame, rather than residue against residue. But you can achieve something similar with a bit of data wrangling:

df = fp.to_dataframe()

# prepare df for residue map
resmap_df = (
    df
    # aggregate all frames together
    .sum(axis=0)
    # aggregate all interaction types together
    .T.groupby(["ligand", "protein"]).sum().astype(bool).T
    # pivot the "protein" level to its own column
    .unstack("protein", fill_value=False)
    # sort both rows and columns by residue id
    .sort_index(axis=0, key=lambda index: [plf.ResidueId.from_string(x) for x in index])
    .sort_index(axis=1, key=lambda index: [plf.ResidueId.from_string(x) for x in index])
)

# plot
fig, ax = plt.subplots(figsize=(8, 5), dpi=100)
im = ax.imshow(
    resmap_df.values,
    aspect="auto",
    interpolation="none",
    cmap="binary",
)
ax.yaxis.set_ticks(range(resmap_df.shape[0]), resmap_df.index)
ax.xaxis.set_ticks(range(resmap_df.shape[1]), resmap_df.columns, rotation=90)

[biomolsim]

Hi @cbouy, thank you very much for the detailed information and for taking the time to resolve my queries.

1. I was unable to use the DCD trajectory as I could not add the elements column to the file. At the end, I had to use a PDB file instead, but I replaced the keyword END between the frames with MODEL for hassle-free parsing with MDAnalysis. Moreover, as with DCD, I could not use a PSF file to provide the topology information because of the missing elements column. I had to use a PDB frame from the PDB trajectory file as the reference.

2. For me, this code did not work:

fp.run(u_read.trajectory, lig_sel, rest_sel)

I kept getting something similar to 0 atoms in the frame error at <N-1>/<N> stage of the code execution, where N refers to the total number of frames in the trajectory.

However, the following code worked:

fp.run(u_read.trajectory[:<N-1>], lig_sel, rest_sel)

3. This really helped me get the desired output. I faced an unknown command issue for plt and added the following line after #plot to resolve it.

import matplotlib.pyplot as plt

[cbouy]

FYI you can add elements or any other missing topology attribute required by the converter with

from MDAnalysis.topology.guessers import guess_types

guessed_elements = guess_types(universe.atoms.names)
universe.add_TopologyAttr('elements', guessed_elements)