PyMol: Output Polar Contacts Between Chains to Text File

Question

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.

Answer 1

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.