1
$\begingroup$

I use cealign in PyMol for structure alignment. Instead of visualization, I want to return the alignment object to my python script for further analysis. Is there a function to return the object without I/O to a clustalw file?

edit:

By alignment object, I meant the sequence alignment in cealign.

$\endgroup$
3
  • $\begingroup$ Thats a good question. The answer should be 'yes' outside Pymol based on my own experience in Biopython PDB, but in Pymol I don't know. Someone will definitely know around here. $\endgroup$
    – M__
    Commented Nov 24, 2020 at 14:45
  • $\begingroup$ By alignment object do you mean the 4x4 rototranslation matrix applied to the mobile element or do you mean a MSA object? with PyMOL, using Python you can get both painfully (i.e. using transform=0, object ='rototrans' attributes on cealign). For the former, 3D operations are an utter nightmare (the Mad Hatter's tea party in Alice in Wonderland is a mockery of quaternion rotations... and not unduly!) —is this to move a ED map with gemmi? Because there are less maddenning ways. For the latter, tcoffee.crg.cat/apps/tcoffee/do:expresso is probably a better choice. $\endgroup$ Commented Nov 24, 2020 at 16:29
  • $\begingroup$ @MatteoFerla good remark. Let me edit the question. The alignment object I meant is the sequence alignment. Since I am performing RMSD analysis, MSA in PyMol would be more convenient. $\endgroup$
    – Simon
    Commented Nov 25, 2020 at 15:14

1 Answer 1

0
$\begingroup$

In the application, unfortunately the only option is doing

  • cealign 👾👾,👾👾, object =👾👾
  • Menu > Export Alignment... on your cealign 👾👾 object
  • and converting the Clustal format using one of the millions of tools online

get_fastastr

As you know cealign adds gaps, but cmd.get_fastastr('👾👾') strips them. There is no option to keep the gaps —annoyingly especially with missing loop info! The help has this to say:

 get_fastastr(selection='all', state=-1, quiet=1, key='', _self=<module 'pymol.cmd' from '/Users/matteo/miniconda3/lib/python3.7/site-packages/pymol/cmd.py'>)

Correction and dissection

So I had a quick check to see if there is an alternative route to get the sequence string with gaps. These are my observations using the Python3 module, sorry if they aren't polished up.

In the comments I said something incorrect. cmd.cealign returns a dictionary with the 4x4 rototranslational matrix —no need for an object. Using Python3 and two TIM barrels as an example:

import pymol2

with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1TIM', 'tim')
    pymol.cmd.remove('not chain A')
    pymol.cmd.fetch('5AHE', 'hisA')
    x = pymol.cmd.cealign('tim', 'hisA')
    print(x)
# {'alignment_length': 184, 'RMSD': 3.8800425577361755, 'rotation_matrix': [...]}

While the object does not contain any useful info from PyMol directly*:

    ...
    pymol.cmd.cealign('tim', 'hisA', object ='rototrans')
    pymol.cmd.get_session('rototrans')['names'][0][5]
    # useful info is generally in ['names'][0][5]

(*) The code is actually wrapped C++ code, so the object rototrans is in the C++ layer and can only be access by session or possibly by some functions in cmd.exporting. In fact, it can be exported to an aln string (ClustalW).

    aln = pymol.cmd.exporting.get_alnstr('rototrans')

In the PyMOL application, doing the cealign with an object (cmd.cealign('tim', 'hisA', object ='rototrans'), allows the Menu>Export Alignment...>'rototrans' which is saved as a ClustalW file, as you mention. Namely:

CLUSTAL

tim          APRKFFVGGN----WK-------MNGKRKSLGELIHTLDGAKLSADTEVVCGAPSI------Y
hisA         MIIPALDLIDGTVVRLHQGDYARQRDYGN---DPLPRLQDYAAQGAGVLHLVDLTGAKDPAKR
                           ..       ......   .....*..................      .

tim          LDFARQK---LDAKIGVAAQNCYKVPKGAFTGEISPAMIKDI-GAAWVILGHSERRHVFGESD
hisA         QIPLIKTLVAGVNVPVQVGGGVRTEE--------DVAALLKAGVARVVIGS---------TAV
             .......   ................        ..*..... .*..**..         ...

tim          ELIGQKVAHALAEG---LGVIACIGEKLDEREAGIT-----------EKVVFQETKAIADNVK
hisA         KSPDV-VKGWFERFGAQALVLALDV---RIDEHGTKQVAVSGWQENSGVSLEQLVETYLPVGL
             ..... *.......   ..*.*...   ...*.*..           .....*..........

tim          DWSKVVLAYEPVWAIGTGKTATP------QQAQEVHEKLRGWLKTHVSDAVAVQSRIIYGGSV
hisA         ----KHVLCTDI----------SRDGTLAGSNVSLYEEVCARYP---------QIAFQSSGGI
                 ........          .      .......*.......         *......*..

tim          TGGNCKELASQHDVDGFLVGGASLKPEFVDIINAKH------------               
hisA         GDIDDIAALRGTGVRGVIVGRALLEGKFTVKEAIQ?WQNVKGHHHHHH               
             .............*.*..**.*.*...*..                                 

However, there is no method in cmd.exporting that gives out just the sequence with gaps.

Whereas the cmd.exporting has the following... none applicable.

 print(dir(pymol.cmd.exporting))
['DEFAULT_ERROR', 'DEFAULT_SUCCESS', 'QuietException', 'SDF', 'SDFRec', 'Shortcut', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_cmd', '_feedback', '_get_dump_str', '_get_mtl_obj', '_raising', '_resn_to_aa', '_session_convert_legacy', '_unit2px', 'absolute_import', 'assign_atom_types', 'cPickle', 'cache', 'cache_action_dict', 'cache_action_sc', 'cmd', 'copy', 'copy_image', 'export_coords', 'fb_mask', 'fb_module', 'get_alnstr', 'get_bytes', 'get_ccp4str', 'get_cifstr', 'get_fastastr', 'get_maestr', 'get_mol2str', 'get_pdbstr', 'get_pqrstr', 'get_psestr', 'get_sdfstr', 'get_session', 'get_str', 'get_xyzstr', 'io', 'is_dict', 'is_error', 'is_list', 'is_ok', 'is_tuple', 'loadable', 'lock', 'multifilenamegen', 'multifilesave', 'multisave', 'os', 'png', 'print_function', 'pymol', 're', 'save', 'savefunctions', 'selector', 'sys', 'thread', 'traceback', 'unlock']
$\endgroup$
2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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