Answers.Body: <p>OK, I tried to answer half of your question. I made a test answer showing how to get the delta-phi , delta-psi between two structure aligned with <strong>PyMOL</strong>, here my code , not sure it is the fastest/cleanest way  but only one I could figure out.
  If you align more structures to a reference you will end up with a raw_alignment that looks like:</p>
  <pre><code>[('3hr8', 1988), ('1xmv', 1764), ('1xp8', 1776)]
  [('3hr8', 1997), ('1xmv', 1773), ('1xp8', 1785)]
  [('3hr8', 2001), ('1xmv', 1777), ('1xp8', 1789)]
  [('3hr8', 2017), ('1xp8', 1808)]
  [('3hr8', 2028), ('1xmv', 1804), ('1xp8', 1817)]
  [('3hr8', 2049), ('1xp8', 1831)]
  [('3hr8', 2058), ('1xmv', 1829), ('1xp8', 1839)]
  [('3hr8', 2093), ('1xmv', 1860)]
  [('3hr8', 2134), ('1xmv', 1906)]
  [('3hr8', 2141), ('1xmv', 1914)]
  [('3hr8', 2205), ('1xmv', 1972)]
  [('3hr8', 2217), ('1xmv', 1978)]
  [('3hr8', 2229), ('1xmv', 1990)]
  [('3hr8', 2295), ('1xmv', 2021)]
  [('3hr8', 2345), ('1xp8', 2081)]
  [('3hr8', 2352), ('1xp8', 2089)]
  [('3hr8', 2361), ('1xmv', 2085), ('1xp8', 2094)]
  [('3hr8', 2372), ('1xmv', 2099), ('1xp8', 2106)]
  [('3hr8', 2380), ('1xp8', 2114)]
  [('3hr8', 2389), ('1xmv', 2116), ('1xp8', 2119)]
  [('3hr8', 2397), ('1xmv', 2124)]
  </code></pre>
  <p>so you'll have to use <code>cmd.iterate</code> over it selecting all the possible couples of reference-target structures and get their CA atoms coupled-indexes extract the corresponding resi (residue number/index) and use <code>phi_psi</code> command to get the two dihedral of interest, like it is done below:</p>
  <pre><code>import pymol
  
  
  from pymol import cmd
  
  from pymol import stored
  
  import time
  
  
  pymol.finish_launching()
  
  
  # get two structures
  # fetch from remote
  cmd.fetch('3hr8 ')
  cmd.fetch('3hr8 1xp8')
  
  
  # load from locale
  # cmd.load('3hr8.cif' , '3hr8')
  # cmd.load('1xp8.cif', '1xp8')
  
  time.sleep(3)
  
  cmd.remove('hetatm')
  
  
  ### Super and get raw alignment ###
  cmd.do('super 1xp8, 3hr8, object = aln')
  
  print('super ended !!!!')
  
  
  ### get raw alignment ###
  raw_aln = cmd.get_raw_alignment('aln')
  
  cmd.do('set seq_view , 1')
  
  # print residue pairs (atom index)
  for idx1, idx2 in raw_aln:
      print('%s`%d -&gt; %s`%d' % tuple(idx1 + idx2))
  
  
  ## To print residue numbers instead of atom indices:
  def myfunc(resi,model,index,name):
  
      if name == 'CA':
          
          idx2resi[model, index] = resi
          
      else :
          
          idx2resi[model, index] = None
          
      
  idx2resi = {}
  cmd.iterate('aln', 'myfunc(resi,model,index,name)', space={'idx2resi': idx2resi , 'myfunc': myfunc})
  
  
  # print(len(idx2resi), len(raw_aln))
  
  errors = []
  
  # print residue pairs (residue number)
  for idx1, idx2 in raw_aln:
      
      # print(idx1 , idx2)
  
      if idx2resi[idx1] and idx2resi[idx2] :
          
              a = idx1[0]+'...'+idx2resi[idx1]
              
              b = idx2[0]+'...'+idx2resi[idx2]
              
              a_die = cmd.phi_psi(f'{idx1[0]} and resi {idx2resi[idx1]}' )
              
              if a_die :
                  
                  # print(a_die , ' OK')
                  pass
                  
              else :
                  
                  a_die[idx1] = (0,000000 , 0,00000000)
                  
                  # print('\n\n error set to 0 !!!!!!!!' , a_die)
                  
                  errors.append((idx1 , idx2resi[idx1]))
              
              b_die = cmd.phi_psi(f'{idx2[0]} and resi {idx2resi[idx2]}' )
              
              if b_die :
                  
                  # print(b_die , ' OK')
                  pass
              
              else :
                  
                  b_die[idx2] = (0,000000 , 0,00000000)
                  
                  # print('\n\n error set to 0 !!!!!!!!' , b_die)
                  
                  errors.append((idx2, idx2resi[idx2]))
              
              # print(a, b)
              
              # print(a_die, b_die)
  
              stored.a_die_resi = str()
              
              stored.b_die_resi = str()
              
              
              cmd.iterate(f&quot;{idx1[0]} and resi {idx2resi[idx1]}&quot; , &quot;stored.a_die_resi = resn&quot; )
              
              cmd.iterate(f&quot;{idx2[0]} and resi {idx2resi[idx2]}&quot; , &quot;stored.b_die_resi = resn&quot; )
          
              # print('%s -&gt; %s' % (a , b ), ' phi_psi :' , a_die)
              
              s = (f&quot;{a} {stored.a_die_resi} ---&gt; phi_psi: {a_die[idx1]}&quot;
                   
              f&quot;   {b} {stored.b_die_resi} ---&gt; phi_psi: {b_die[idx2]}&quot;
              
              f&quot;   DELTA : {abs(a_die[idx1][0] -b_die[idx2][0])}  , {abs(a_die[idx1][1] -b_die[idx2][1])}&quot;)
              
              print(s)
              
  
  
  print('\n\n Errors : \n ', errors)
  </code></pre>
  <p>If for same reason (gap in structure) one of the residues doesnt return phi_psi its value gets assigned to 0 and is appended in the error list, probably would be better to remove it from the results, but took already a lot of time to answer this ... so ...).</p>
  <p>Last bit of output:</p>
  <pre><code>3hr8...313 VAL ---&gt; phi_psi: (-64.3411865234375, -44.5438232421875)   1xp8...319 ILE ---&gt; phi_psi: (-64.2035903930664, -27.23161506652832)   DELTA : 0.13759613037109375  , 17.31220817565918
  3hr8...314 GLN ---&gt; phi_psi: (-61.05774688720703, -39.42118453979492)   1xp8...320 ALA ---&gt; phi_psi: (-75.51972198486328, -28.401968002319336)   DELTA : 14.46197509765625  , 11.019216537475586
  3hr8...315 PHE ---&gt; phi_psi: (-56.17056655883789, -43.013572692871094)   1xp8...321 TYR ---&gt; phi_psi: (-71.33622741699219, -41.687713623046875)   DELTA : 15.165660858154297  , 1.3258590698242188
  3hr8...316 LEU ---&gt; phi_psi: (-62.212162017822266, -40.03855895996094)   1xp8...322 ILE ---&gt; phi_psi: (-69.78271484375, -35.62084197998047)   DELTA : 7.570552825927734  , 4.417716979980469
  3hr8...317 LYS ---&gt; phi_psi: (-62.3416633605957, -42.122344970703125)   1xp8...323 ALA ---&gt; phi_psi: (-59.84389114379883, -19.88322639465332)   DELTA : 2.497772216796875  , 22.239118576049805
  3hr8...318 ASP ---&gt; phi_psi: (-77.59044647216797, -13.573332786560059)   1xp8...324 GLU ---&gt; phi_psi: (-95.1947250366211, -7.630105972290039)   DELTA : 17.604278564453125  , 5.9432268142700195
  
  
   Errors : 
    [(('1xp8', 1386), '223')]
  </code></pre>
  <p>here a nice pics of <strong>PyMOL</strong> showing the correspondence between aligned/superposed atoms :</p>
  <p><a href="https://i.sstatic.net/Gg75e8QE.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Gg75e8QE.png" alt="enter image description here" /></a></p>
  <p>try playing with the alignment algorithm and different selections for target and reference and see if you get any different results:</p>
  <p><code>cmd.super( 'mobile' , 'target' , object='aln') </code></p>
  <p>vs</p>
  <p><code>cmd.align( 'mobile' , 'target' , object='aln') </code></p>
  <pre><code>cmd.select('target' , '/3hr8//A//CA')
  cmd.select('mobile' , '/1xp8//A//CA')
  cmd.super( 'mobile' , 'target'  , cycles=0, transform = 0, object='aln') 
  </code></pre>
  <p>vs</p>
  <pre><code>cmd.select('target' , '3hr8')
  cmd.select('mobile' , '1xp8')
  cmd.super( 'mobile' , 'target'  , cycles=0, transform = 0, object='aln') 
  </code></pre>
  <p>and cycles parameter :</p>
  <blockquote>
  <p>cycles = int: maximum number of outlier rejection cycles {default: 5}</p>
  </blockquote>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-05-24T17:40:37.303Z
Answers.LastActivityDate: 2024-05-24T19:38:17.613Z
Answers.LastEditDate: 2024-05-24T19:38:17.613Z
Answers.UserDisplayName: pippo1980
Body: <p>I aligned multiple protein structures from the same protein family to a reference structure in pymol and now I want to calculate the Difference in Phi psi values for each 'aligned' residue pair between Reference Structure and each Target Structure, I know about extracting the phi psi values using phi_psi command and I have the sequence alignment but I don't know how to extract Delta Phi Psi values for every aligned residue between Reference and Target.</p>
  <p>Here the code for fetching and generate a multiple alignment using PyMOL <a href="https://pymolwiki.org/index.php/Super" rel="nofollow noreferrer">Super</a>:</p>
  <blockquote>
  <p>super aligns two selections. It does a sequence-independent (unlike align) structure-based dynamic programming alignment followed by a series of refinement cycles intended to improve the fit by eliminating pairing with high relative variability (just like align). super is more robust than align for proteins with low sequence similarity.</p>
  </blockquote>
  <p>between PDBs <a href="https://www.rcsb.org/structure/3hr8" rel="nofollow noreferrer">3HR8</a> (reference) , <a href="https://www.rcsb.org/structure/1XMV" rel="nofollow noreferrer">1XMV</a> and <a href="https://www.rcsb.org/structure/1XP8" rel="nofollow noreferrer">1XP8</a> :</p>
  <pre><code>import pymol
  from pymol import cmd
  import time
  
  
  list_pdbs = ['1xp8' , '3hr8', '1xmv']
  
  pymol.finish_launching()
  
  time.sleep(5)
  
  for _i in list_pdbs :
      try :
          
          cmd.fetch(_i)
          
      except :
          
          print('cant load ', _i)
                
  cmd.remove('hetatm')
  
  sel = 'name CA'
  ref = list_pdbs[1]
  method = 'super'
  
  cmd.do(f&quot;extra_fit {sel} , {ref} , {method}, object = aln_super_all&quot;)
  
  cmd.do('set seq_view , 1')
  
  cmd.select('aln_super_all')
  
  print(cmd.get_type('aln_super_all'))
  
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-02-02T15:52:04.670Z
Comments.Text: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2024-05-22T17:02:59.887Z
Comments.Text: and at list provide the first half of the data i.e. the proteins lets say one reference and two other from the family and how you aligned them
Comments.UserDisplayName: pippo1980
Tags: <pymol>
Title: How to Calculate Delta Phi Psi between aligned pair of residue between Reference and Target Protein Structure in Pymol
ThreadText: # How to Calculate Delta Phi Psi between aligned pair of residue between Reference and Target Protein Structure in Pymol
  
  (Asked by: username_0 on 2024-02-02)
  
  I aligned multiple protein structures from the same protein family to a reference structure in pymol and now I want to calculate the Difference in Phi psi values for each 'aligned' residue pair between Reference Structure and each Target Structure, I know about extracting the phi psi values using phi\_psi command and I have the sequence alignment but I don't know how to extract Delta Phi Psi values for every aligned residue between Reference and Target.
  
  Here the code for fetching and generate a multiple alignment using PyMOL [Super](https://pymolwiki.org/index.php/Super):
  
  > super aligns two selections. It does a sequence-independent (unlike align) structure-based dynamic programming alignment followed by a series of refinement cycles intended to improve the fit by eliminating pairing with high relative variability (just like align). super is more robust than align for proteins with low sequence similarity.
  
  between PDBs [3HR8](https://www.rcsb.org/structure/3hr8) (reference) , [1XMV](https://www.rcsb.org/structure/1XMV) and [1XP8](https://www.rcsb.org/structure/1XP8) :
  
  ```
  import pymol
  from pymol import cmd
  import time
  
  list_pdbs = ['1xp8' , '3hr8', '1xmv']
  
  pymol.finish_launching()
  
  time.sleep(5)
  
  for _i in list_pdbs :
      try :
          
          cmd.fetch(_i)
          
      except :
          
          print('cant load ', _i)
                
  cmd.remove('hetatm')
  
  sel = 'name CA'
  ref = list_pdbs[1]
  method = 'super'
  
  cmd.do(f"extra_fit {sel} , {ref} , {method}, object = aln_super_all")
  
  cmd.do('set seq_view , 1')
  
  cmd.select('aln_super_all')
  
  print(cmd.get_type('aln_super_all'))
  
  ```
  
  **Tags:** pymol
  
  > **Comments:**
  > **username_1**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
  > **username_2**: and at list provide the first half of the data i.e. the proteins lets say one reference and two other from the family and how you aligned them  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2024-05-24. Score: 0
  
  OK, I tried to answer half of your question. I made a test answer showing how to get the delta-phi , delta-psi between two structure aligned with **PyMOL**, here my code , not sure it is the fastest/cleanest way but only one I could figure out.
  If you align more structures to a reference you will end up with a raw\_alignment that looks like:
  
  ```
  [('3hr8', 1988), ('1xmv', 1764), ('1xp8', 1776)]
  [('3hr8', 1997), ('1xmv', 1773), ('1xp8', 1785)]
  [('3hr8', 2001), ('1xmv', 1777), ('1xp8', 1789)]
  [('3hr8', 2017), ('1xp8', 1808)]
  [('3hr8', 2028), ('1xmv', 1804), ('1xp8', 1817)]
  [('3hr8', 2049), ('1xp8', 1831)]
  [('3hr8', 2058), ('1xmv', 1829), ('1xp8', 1839)]
  [('3hr8', 2093), ('1xmv', 1860)]
  [('3hr8', 2134), ('1xmv', 1906)]
  [('3hr8', 2141), ('1xmv', 1914)]
  [('3hr8', 2205), ('1xmv', 1972)]
  [('3hr8', 2217), ('1xmv', 1978)]
  [('3hr8', 2229), ('1xmv', 1990)]
  [('3hr8', 2295), ('1xmv', 2021)]
  [('3hr8', 2345), ('1xp8', 2081)]
  [('3hr8', 2352), ('1xp8', 2089)]
  [('3hr8', 2361), ('1xmv', 2085), ('1xp8', 2094)]
  [('3hr8', 2372), ('1xmv', 2099), ('1xp8', 2106)]
  [('3hr8', 2380), ('1xp8', 2114)]
  [('3hr8', 2389), ('1xmv', 2116), ('1xp8', 2119)]
  [('3hr8', 2397), ('1xmv', 2124)]
  
  ```
  
  so you'll have to use `cmd.iterate` over it selecting all the possible couples of reference-target structures and get their CA atoms coupled-indexes extract the corresponding resi (residue number/index) and use `phi_psi` command to get the two dihedral of interest, like it is done below:
  
  ```
  import pymol
  
  from pymol import cmd
  
  from pymol import stored
  
  import time
  
  pymol.finish_launching()
  
  # get two structures
  # fetch from remote
  cmd.fetch('3hr8 ')
  cmd.fetch('3hr8 1xp8')
  
  # load from locale
  # cmd.load('3hr8.cif' , '3hr8')
  # cmd.load('1xp8.cif', '1xp8')
  
  time.sleep(3)
  
  cmd.remove('hetatm')
  
  ### Super and get raw alignment ###
  cmd.do('super 1xp8, 3hr8, object = aln')
  
  print('super ended !!!!')
  
  ### get raw alignment ###
  raw_aln = cmd.get_raw_alignment('aln')
  
  cmd.do('set seq_view , 1')
  
  # print residue pairs (atom index)
  for idx1, idx2 in raw_aln:
      print('%s`%d -> %s`%d' % tuple(idx1 + idx2))
  
  ## To print residue numbers instead of atom indices:
  def myfunc(resi,model,index,name):
  
      if name == 'CA':
          
          idx2resi[model, index] = resi
          
      else :
          
          idx2resi[model, index] = None
          
      
  idx2resi = {}
  cmd.iterate('aln', 'myfunc(resi,model,index,name)', space={'idx2resi': idx2resi , 'myfunc': myfunc})
  
  # print(len(idx2resi), len(raw_aln))
  
  errors = []
  
  # print residue pairs (residue number)
  for idx1, idx2 in raw_aln:
      
      # print(idx1 , idx2)
  
      if idx2resi[idx1] and idx2resi[idx2] :
          
              a = idx1[0]+'...'+idx2resi[idx1]
              
              b = idx2[0]+'...'+idx2resi[idx2]
              
              a_die = cmd.phi_psi(f'{idx1[0]} and resi {idx2resi[idx1]}' )
              
              if a_die :
                  
                  # print(a_die , ' OK')
                  pass
                  
              else :
                  
                  a_die[idx1] = (0,000000 , 0,00000000)
                  
                  # print('\n\n error set to 0 !!!!!!!!' , a_die)
                  
                  errors.append((idx1 , idx2resi[idx1]))
              
              b_die = cmd.phi_psi(f'{idx2[0]} and resi {idx2resi[idx2]}' )
              
              if b_die :
                  
                  # print(b_die , ' OK')
                  pass
              
              else :
                  
                  b_die[idx2] = (0,000000 , 0,00000000)
                  
                  # print('\n\n error set to 0 !!!!!!!!' , b_die)
                  
                  errors.append((idx2, idx2resi[idx2]))
              
              # print(a, b)
              
              # print(a_die, b_die)
  
              stored.a_die_resi = str()
              
              stored.b_die_resi = str()
              
              
              cmd.iterate(f"{idx1[0]} and resi {idx2resi[idx1]}" , "stored.a_die_resi = resn" )
              
              cmd.iterate(f"{idx2[0]} and resi {idx2resi[idx2]}" , "stored.b_die_resi = resn" )
          
              # print('%s -> %s' % (a , b ), ' phi_psi :' , a_die)
              
              s = (f"{a} {stored.a_die_resi} ---> phi_psi: {a_die[idx1]}"
                   
              f"   {b} {stored.b_die_resi} ---> phi_psi: {b_die[idx2]}"
              
              f"   DELTA : {abs(a_die[idx1][0] -b_die[idx2][0])}  , {abs(a_die[idx1][1] -b_die[idx2][1])}")
              
              print(s)
              
  
  print('\n\n Errors : \n ', errors)
  
  ```
  
  If for same reason (gap in structure) one of the residues doesnt return phi\_psi its value gets assigned to 0 and is appended in the error list, probably would be better to remove it from the results, but took already a lot of time to answer this ... so ...).
  
  Last bit of output:
  
  ```
  3hr8...313 VAL ---> phi_psi: (-64.3411865234375, -44.5438232421875)   1xp8...319 ILE ---> phi_psi: (-64.2035903930664, -27.23161506652832)   DELTA : 0.13759613037109375  , 17.31220817565918
  3hr8...314 GLN ---> phi_psi: (-61.05774688720703, -39.42118453979492)   1xp8...320 ALA ---> phi_psi: (-75.51972198486328, -28.401968002319336)   DELTA : 14.46197509765625  , 11.019216537475586
  3hr8...315 PHE ---> phi_psi: (-56.17056655883789, -43.013572692871094)   1xp8...321 TYR ---> phi_psi: (-71.33622741699219, -41.687713623046875)   DELTA : 15.165660858154297  , 1.3258590698242188
  3hr8...316 LEU ---> phi_psi: (-62.212162017822266, -40.03855895996094)   1xp8...322 ILE ---> phi_psi: (-69.78271484375, -35.62084197998047)   DELTA : 7.570552825927734  , 4.417716979980469
  3hr8...317 LYS ---> phi_psi: (-62.3416633605957, -42.122344970703125)   1xp8...323 ALA ---> phi_psi: (-59.84389114379883, -19.88322639465332)   DELTA : 2.497772216796875  , 22.239118576049805
  3hr8...318 ASP ---> phi_psi: (-77.59044647216797, -13.573332786560059)   1xp8...324 GLU ---> phi_psi: (-95.1947250366211, -7.630105972290039)   DELTA : 17.604278564453125  , 5.9432268142700195
  
   Errors : 
    [(('1xp8', 1386), '223')]
  
  ```
  
  here a nice pics of **PyMOL** showing the correspondence between aligned/superposed atoms :
  
  [![enter image description here](https://i.sstatic.net/Gg75e8QE.png)](https://i.sstatic.net/Gg75e8QE.png)
  
  try playing with the alignment algorithm and different selections for target and reference and see if you get any different results:
  
  `cmd.super( 'mobile' , 'target' , object='aln')`
  
  vs
  
  `cmd.align( 'mobile' , 'target' , object='aln')`
  
  ```
  cmd.select('target' , '/3hr8//A//CA')
  cmd.select('mobile' , '/1xp8//A//CA')
  cmd.super( 'mobile' , 'target'  , cycles=0, transform = 0, object='aln') 
  
  ```
  
  vs
  
  ```
  cmd.select('target' , '3hr8')
  cmd.select('mobile' , '1xp8')
  cmd.super( 'mobile' , 'target'  , cycles=0, transform = 0, object='aln') 
  
  ```
  
  and cycles parameter :
  
  > cycles = int: maximum number of outlier rejection cycles {default: 5}
  
  
  ---
  
  