2
$\begingroup$

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:

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 (reference) , 1XMV and 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'))

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Community Bot
    Commented Feb 2 at 15:52
  • $\begingroup$ 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 $\endgroup$
    – pippo1980
    Commented May 22 at 17:02

1 Answer 1

0
$\begingroup$

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

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}

$\endgroup$

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.