0
$\begingroup$

I am new to PyMOL but have a very specific task that I need to do.

I have a PDB structure file of protein homo-oligomer, and I want to use PyMol to determine polar contacts between chain A and all other chains.

However, I not only want to visualize these polar contacts (which is easy to do via PyMOL GUI), but I also want to be able to store/save pairs of interacting residues in a dictionary or whatever may be the most convenient data structure to do this. The end goal to create a sort of "interaction map" in a csv or text file where each interacting pair of residues has its own line.

I would ideally like to do this via a Python script where I call the PyMol API or through a PyMol script - whichever is easier.

I have seen a similar question on Stack Overflow but found it unhelpful because I not only want to select and visualize these polar contacts but also want to save/record them in a list or dictionary (and then to a text or csv file), which the above question does not address.

$\endgroup$
2
  • 1
    $\begingroup$ Please do not shotgun post to multiple sites. (Post on biostars: biostars.org/p/9484676) $\endgroup$
    – Ram RS
    Aug 12 at 14:07
  • $\begingroup$ You should give a try to ChimeraX's hbonds or contacts. $\endgroup$ Aug 14 at 20:22
2
$\begingroup$

So you can install PyMOL as a standalone, but you can also install it as a bona fide Python 3 module via conda:

conda install -c schrodinger -c conda-forge pymol-bundle

It can be installed in other ways —apt-get or brew or even compiled, but the latter is excruciatingly painful.

In your python notebook you can do:

import pymol2
with pymol2.PyMOL() as pymol:
    pymol.cmd.load('structure_file.pdb', 'foo') # file
    # pymol.cmd.fetch('PDB_CODE') # PDB code
    # pymol.cmd.read_pdbstr(foo_str, 'foo') # string
    ...

where the ellipsis can be any API cmd.

This context manager approach is to avoid the problem of namespace pollution as the session (e.g. loading structures without removing the old ones etc.) is the same but you don't see what is happening. Also PyMOL() is parallel the parallel mode (so mutagenesis wizard does not work), while in singleton mode (pymol2.SingletonPyMOL), you get issues if you try to re-initialise etc.

Anyway. To "select" all interacting atoms at the interface you can do:

n = pymol.cmd.select('chain A around 3')

Expand includes the selection. Around does not.

This adds a selection to the session. So what we actually want is a pythonic instance to play with. So:

atoms = pymol.cmd.get_model('chain A around 3')

where atoms is a list of chempy.Atom objects. With attributes such as resi, resn, chain etc. Do dir on an atom for more...

One attribute is coord, which contains a length 3 list of floats.

NB. The above selection gives not a full residue, whereas 'byres chain A around 3' does. So if you want to make a list of all interacting residues:

interactors = [{'chain': atom.chain, 'index': atom.resi, 'number': atom.resn} for atom in pymol.cmd.get_model('name CA and byres chain A around 3')]

An iteration can give the atom to atom distances.

Do note that this is all proximity, so you will have to screen out hydrogens and carbons to get actual polar contacts.

$\endgroup$
1
  • 1
    $\begingroup$ After some searching, I discovered the command 'cmd.find_pairs(selection1, selection2, mode, angle, distance)'. If you pass it two selections, it will return pairs of atoms within the given criteria. This seems to be better than using "around" because it takes into account atom orientation - which is important when inferring hydrogen bonds. $\endgroup$ Aug 25 at 16:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.