Re: [PyMOL] How to best create a sliding window (a selection of specific residues) and compare specific residue names+numbers in pymol/python?

[Thomas H. <tho...@sc...>]

Hi Anne,

This is not at all an obvious question, it's quite advanced! Impressive code snippet for a newbie :-)

For the sliding window question, I suggest to have a look at the implementation of "local_rms" from the psico package:

https://github.com/speleo3/pymol-psico/blob/master/psico/fitting.py

Example usage:

import psico.fitting
fetch 1akeA, async=0
fetch 4akeA, async=0
local_rms 1akeA, 4akeA, window=5

It calculates the RMSD of the entire sliding window, not only one residue. And it only considers the C-alpha atoms (remove the "and guide" selector which is passed to Matchmaker to include all atoms).

My question would be: What does the RMSD of a single residue really tell you? The number will reflect something in between the fit of the sliding window and the sidechain conformation of the center residue. Are you interested in both? Or only one of them?

Cheers,
  Thomas


> On Feb 22, 2019, at 6:28 PM, Anne Nierobisch <ann...@ke...> wrote:
> 
> Hi,
> 
> I need to calculate the RMSD for the same residue, e.g. 131Thr, from 2 pdb files for the same target. As I need a local alignment, I use a sliding window of 5 residues (the residue of interest is in the middle.)
> I have adapted the script from https://pymolwiki.org/index.php/RmsdByResidue (also see below my code below) by adding a sliding window, but I need advice on the following:
> 
> - I would need to check whether I am actually comparing the intended residues, e.g. 131Thr, 131Thr, or whether the same residues (e.g. His, Asp, His), with specific residue names and numbers   are in the selection I have chosen from both pdb files for my sliding window. 
> The reason:
> I have coded up the script below, but often it skips and doesn't compare two residues,  because the atom count is different between the residue from protein A and protein B. 
> (It is exactly double the number, I have already checked and excluded problems like occupancy and different conformations as potential causes.)
> 
> - I constructed a sliding window by selecting two residues before the residue of interest and two residues which follow said residue. (I know that residues can be missing, I am working on this.)
> Is there a better way for constructing a sliding window? I have not found such a method in pymol.
> 
> For anyone interested, I have attached a code snipplet below.
> I am sorry, if these seem like obvious questions, but I tried various approaches and I feel that I need a push in the right direction. I have the feeling that I am missing something fundamental. 
> 
> Many thanks for any suggestion!
> Anne (Newbie in Pymol)
> 
> ######################################
> 
> Code snipplet (Python/Pymol interface):
>   referenceProteinChain = cmd.fetch(pbdstructure1)
>   targetProteinChain = cmd.fetch(pdbstructure2)
>   sel = referenceProteinChain
>   list_atoms = [[133, 133]] # example list, I want to calculate the rmsd between residue 133 and residue 133 of two pdb structures
> 
> for j in list_atoms:
>    # I formulate my sliding window of 5 residues, my residue of interest is in the middle
>     ref_begin = int(j[0])-2
>     ref_end = int(j[0])+2
>     target_begin = int(j[1])-2
>     target_end = int(j[1])+2
> 
>     # I create a selection for the reference protein and the target protein 
>     cmd.create("ref_gzt", referenceProteinChain+" and polymer and not alt B and not Alt C and not alt D  and resi %s-%s" 
>     cmd.alter("ref_gzt", "chain='A'")
>     cmd.alter("ref_gzt", "segi=''")
>     cmd.create("target_gzt", targetProteinChain+" and polymer and  not alt B and not Alt C and not alt D  and resi %s-%s" %(target_begin,target_end) )
>     cmd.alter("target_gzt", "chain='A'")
>     cmd.alter("target_gzt", "segi=''")
> 
>     # I align my selected 5 residues for a local alignment window
>     cmd.align("target_gzt","ref_gzt",object="align", cycles =5) 
> 
> 
>    # select alpha carbon of in reference structure to loop over
>     calpha=cmd.get_model(sel+" and name CA and not alt B and not Alt C and not alt D  and resi %s-%s" %(ref_begin,ref_end) )
> 
>    ## here I loop over all residues in the sliding window and calculte the rmsd for my residues of interest.
>     for g in calpha.atom : I loop over the slinding window
>         if g.resi == str(j[0]): ## we select only our residue of interest within the sliding window
>             if cmd.count_atoms("ref_gzt and polymer and resi "+g.resi)==cmd.count_atoms("target_gzt and polymer and resi "+g.resi):
> 
>                 ## calculte the pairwise RMSD between the residues I specified in list_atoms 
>                 rmsdRes=cmd.rms_cur("ref_gzt and polymer and resi "+g.resi,"target_gzt and polymer and resi "+g.resi)

--
Thomas Holder
PyMOL Principal Developer
Schrödinger, Inc.