BioPython internal_coords module returns different dihedral angles for the (seemingly) same protein structure

Question

I am using the internal_coords module from BioPython in order to compare dihedral angles of two different conformations of the same protein. The conformations are processed by me in order to make sure that their residues are aligned and they always contain exactly the same atoms, so that in the end I expect to have the same number of dihedral angles in the two structures.

However, despite the match in the residues and atoms, I obtain a different number of dihedral angles for the two conformations.

You can find a reproducible code sample here.

from Bio.PDB import *
 
def extract_torsion_angles(structure):
    ic_chain = internal_coords.IC_Chain(structure)
    ic_chain.atom_to_internal_coordinates()
 
    d = ic_chain.dihedra
 
    """ 
    Update the keys of a dihedra dictionary to make searching easier.
 
    The keys of the dihedra dictionary (as obtained from 
    extract_torsion_angles()) are 4-tuples where each element is
    an AtomKey containing 5 fields (respos, icode, resname, atm,
    altloc, occ). We use the residue position, residue name and 
    atom name in order to create a new representation for each
    AtomKey. The final dihedra dictionary will have 4-tuple keys
    where each element of the tuple is one such string.
    """
 
    new_d = {}
    for key, value in sorted(d.items(), key=lambda x: x[0]):
        new_key = []
        for atom_key in key:
            # Keep: residue position + resiude name + atom name.
            atom_key_str = (atom_key.akl[0] + '_' + 
                atom_key.akl[2] + '_' + atom_key.akl[3])
 
            new_key.append(atom_key_str)
 
        new_d[tuple(new_key)] = value
 
    return new_d
 
 
f1 = 'protein1.pdb'
f2 = 'protein2.pdb'
 
parser = PDBParser(PERMISSIVE=1, QUIET=0)
s1 = parser.get_structure('', f1)
s2 = parser.get_structure('', f2)
 
c1 = list(s1.get_chains())
c2 = list(s2.get_chains())
 
print('Chains: ', len(c1), ' ', len(c2))  
 
r1 = list(c1[0].get_residues())
r2 = list(c2[0].get_residues())
 
print('Residues: ', len(r1), ' ', len(r2))
 
for i in range(len(r1)):
    resname1 = r1[i].get_resname()
    resname2 = r2[i].get_resname()
 
    if resname1 != resname2:
        print('Residues not matching')
        break
 
    atoms1 = list(r1[i].get_atoms())
    atoms2 = list(r2[i].get_atoms())
 
    atom_names1 = sorted([atom.get_id() for atom in atoms1])
    atom_names2 = sorted([atom.get_id() for atom in atoms2])
 
    if atom_names1 != atom_names2:
        print('Different atoms!')
 
d1 = extract_torsion_angles(s1)
d2 = extract_torsion_angles(s2)
 
print('Torsion angles: ', len(d1), ' ', len(d2))
 
diff_keys = set(d2.keys()) - set(d1.keys())
print(diff_keys)

These are the two PDB files protein1.pdb and protein2.pdb used.

As you can see, the number of chains, the number of residues, the residue names and the atom names in each residue are the same, yet there are 4 extra dihedral angles in the second protein.

My question is, am I missing something from my comparison? Is it unreasonable to expect the same dihedral angles for these two structures? (And if yes, why?). Is this the expected behaviour of the atom_to_internal_coordinates() function?

And extra: do you have any recommendations for other tools that can easily compute all dihedral angles (including the ones in side chains) for a protein?

By different you mean wildly different? Do the residues look like different when you superpose the protein in PyMOL (align) and look at the residues by eye —PyMOL does have the command get_dihedral, but the laziest way would be to click the four atoms (in edit mode) and it will show the dihedral. Do the values differ? — Matteo Ferla, Mar 10, 2022 at 8:50
I should also say, the two conformers should have some slight differences. Especially, surface sidechains and especially if the two protein are crystallised under different conditions... — Matteo Ferla, Mar 10, 2022 at 8:54
@MatteoFerla Thanks for your input! The problem in this case was not that the angles are different (which was indeed expected as the proteins were coming from different structures), but rather that BioPython just didn't output certain angles at all for one of the proteins. — CubeHead, Mar 11, 2022 at 15:03

CubeHead · Accepted Answer · 2022-03-11 15:07:00Z

6

The problem got solved after posting on the official BioPython repo.

Because of the processing (where some residues were removed in order to align the proteins), some consecutive residues were too far away from each other and thus the bond length was too large. The solution was to increase the maximum peptide bond allowed, as follows:

internal_coords.IC_Chain.MaxPeptideBond = 5.0
ic_chain = internal_coords.IC_Chain(structure, verbose=True)

answered Mar 11, 2022 at 15:07

CubeHead

4152 silver badges8 bronze badges

2

$\begingroup$ Thanks for posting back @cubehead. If you mark the post as accepted that would be cool $\endgroup$
– M__ ♦
Mar 11, 2022 at 15:11
$\begingroup$ Just out of curiosity: were these naughty bonds Asp-Pro or Glu-Pro in loops (self-cleavage and isoglutamate)? Or the more classic case of bad density fitting of loops? $\endgroup$
– Matteo Ferla
Mar 11, 2022 at 16:02

Add a comment |

pippo1980 · Accepted Answer · 2023-11-18 21:32:42Z

OK , visualizing both proteins in PyMOL , using set cartoon_gap_cutoff = 0 ,

I get :

Apparently the two conformations , have more than one difference in various loops.

I modify the code above , thanks to your hint :

internal_coords.IC_Chain.MaxPeptideBond = 5.0
ic_chain = internal_coords.IC_Chain(structure, verbose=True)

My code :

import Bio

print('########## Biopython  VERSION ##########################################')
print('\n                BIOPYTHON_VERSION : ', Bio.__version__)
print('\n########################################################################\n\n')


from Bio.PDB import *


def extract_torsion_angles(structure):
    
    print('\nSTRUCTURE ---> ' , structure)
    
    ic_chain = internal_coords.IC_Chain(structure, verbose=True)
    
    ic_chain.atom_to_internal_coordinates()
 
    d = ic_chain.dihedra
 
    """ 
    Update the keys of a dihedra dictionary to make searching easier.
 
    The keys of the dihedra dictionary (as obtained from 
    extract_torsion_angles()) are 4-tuples where each element is
    an AtomKey containing 5 fields (respos, icode, resname, atm,
    altloc, occ). We use the residue position, residue name and 
    atom name in order to create a new representation for each
    AtomKey. The final dihedra dictionary will have 4-tuple keys
    where each element of the tuple is one such string.
    """
 
    new_d = {}
    for key, value in sorted(d.items(), key=lambda x: x[0]):
        new_key = []
        for atom_key in key:
            # Keep: residue position + resiude name + atom name.
            atom_key_str = (atom_key.akl[0] + '_' + 
                atom_key.akl[2] + '_' + atom_key.akl[3])
 
            new_key.append(atom_key_str)
 
        # new_d[tuple(new_key)] = value
        
        new_d[tuple(new_key)] = tuple((value , structure.get_full_id()))
 
    return new_d
 
 
f1 = '18760-BIOINFORMATICS-protein1.pdb'
f2 = '18760-BIOINFORMATICS-protein2.pdb'

 
parser = PDBParser(PERMISSIVE=1, QUIET=0)
s1 = parser.get_structure('protein1', f1)
s2 = parser.get_structure('protein2', f2)
 
c1 = list(s1.get_chains())
c2 = list(s2.get_chains())
 
print('Chains: ', len(c1), ' ', len(c2))  
 
r1 = list(c1[0].get_residues())
r2 = list(c2[0].get_residues())
 
print('Residues: ', len(r1), ' ', len(r2))
 
for i in range(len(r1)):
    resname1 = r1[i].get_resname()
    resname2 = r2[i].get_resname()
 
    if resname1 != resname2:
        print('Residues not matching')
        break
 
    atoms1 = list(r1[i].get_atoms())
    atoms2 = list(r2[i].get_atoms())
 
    atom_names1 = sorted([atom.get_id() for atom in atoms1])
    atom_names2 = sorted([atom.get_id() for atom in atoms2])
 
    if atom_names1 != atom_names2:
        print('Different atoms!')
 


def get_no_tors_diff(structure1 , structure2) :
    
    pept_bond = 1.4
    
    while pept_bond < 10.0 :
        
        pept_bond = round(pept_bond , 2)
        
        print('\n\n Peptide Bond Lenght = ' , pept_bond)
        
        internal_coords.IC_Chain.MaxPeptideBond = pept_bond
        
        d1 = extract_torsion_angles(s1)
        d2 = extract_torsion_angles(s2)
         
        print('Torsion angles: ', len(d1), ' ', len(d2))
         
        diff_keys = sorted((set(d1.keys()) ^ set(d2.keys())) , key = lambda tup : tup[0].split('_')[0])

        print('diff_keys lenght : ', len(diff_keys))

        print(diff_keys , type(diff_keys))

        for i in diff_keys :
            
            print(i , i[0].split('_')[0] , type(i))



        diff_keys_values = {}

        for i in diff_keys :
            
            if i in d1.keys() :
                
                diff_keys_values[i] = d1[i]
                
            
            if i in d2.keys()  :
                
                diff_keys_values[i] = d2[i]
                
            else:
                
                pass


        print('diff_keys_values lenght : ', len(diff_keys_values))
            
        for i in diff_keys_values :
            
            print(i , diff_keys_values[i][1])
    
            
        print('\n\n')
        
        if len(diff_keys) == 0 :
            
            print('\n\n Peptide Bond Lenght without differences = ' , pept_bond)
            
            return pept_bond
        
        pept_bond += 0.1
        
        
    print('\n\n Cannot get rid of differences up to 10.0A in Peptide Bond ')
    
    return None
    


print('get_no_tors_diff(structure1 , structure2)' , get_no_tors_diff(s1 , s2))

First and last bit of output :

Chains:  1   1
Residues:  313   313


 Peptide Bond Lenght =  1.4

STRUCTURE --->  <Structure id=protein1>
chain break at GLY  6  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  37  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  38  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  75  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  84  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  92  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PRO  163  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  164  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at HIS  198  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  224  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  246  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  251  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  279  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  285  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (1.4 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at PHE  1  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at GLY  6  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  37  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  38  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  75  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  84  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  92  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PRO  163  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  164  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at HIS  198  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  224  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  246  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  251  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  278  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  285  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (1.4 angstroms) exceeded
Torsion angles:  2504   2501
diff_keys lenght :  9
[('0_N_C', '1_F_N', '1_F_CA', '1_F_C'), ('0_N_CA', '0_N_C', '1_F_N', '1_F_CA'), ('0_N_N', '0_N_CA', '0_N_C', '1_F_N'), ('277_D_C', '278_A_N', '278_A_CA', '278_A_C'), ('277_D_N', '277_D_CA', '277_D_C', '278_A_N'), ('277_D_CA', '277_D_C', '278_A_N', '278_A_CA'), ('278_A_N', '278_A_CA', '278_A_C', '279_K_N'), ('278_A_C', '279_K_N', '279_K_CA', '279_K_C'), ('278_A_CA', '278_A_C', '279_K_N', '279_K_CA')] <class 'list'>
('0_N_C', '1_F_N', '1_F_CA', '1_F_C') 0 <class 'tuple'>
('0_N_CA', '0_N_C', '1_F_N', '1_F_CA') 0 <class 'tuple'>
('0_N_N', '0_N_CA', '0_N_C', '1_F_N') 0 <class 'tuple'>
('277_D_C', '278_A_N', '278_A_CA', '278_A_C') 277 <class 'tuple'>
('277_D_N', '277_D_CA', '277_D_C', '278_A_N') 277 <class 'tuple'>
('277_D_CA', '277_D_C', '278_A_N', '278_A_CA') 277 <class 'tuple'>
('278_A_N', '278_A_CA', '278_A_C', '279_K_N') 278 <class 'tuple'>
('278_A_C', '279_K_N', '279_K_CA', '279_K_C') 278 <class 'tuple'>
('278_A_CA', '278_A_C', '279_K_N', '279_K_CA') 278 <class 'tuple'>
diff_keys_values lenght :  9
('0_N_C', '1_F_N', '1_F_CA', '1_F_C') ('protein1',)
('0_N_CA', '0_N_C', '1_F_N', '1_F_CA') ('protein1',)
('0_N_N', '0_N_CA', '0_N_C', '1_F_N') ('protein1',)
('277_D_C', '278_A_N', '278_A_CA', '278_A_C') ('protein1',)
('277_D_N', '277_D_CA', '277_D_C', '278_A_N') ('protein1',)
('277_D_CA', '277_D_C', '278_A_N', '278_A_CA') ('protein1',)
('278_A_N', '278_A_CA', '278_A_C', '279_K_N') ('protein2',)
('278_A_C', '279_K_N', '279_K_CA', '279_K_C') ('protein2',)
('278_A_CA', '278_A_C', '279_K_N', '279_K_CA') ('protein2',)

.........................
.................
..........
......
 Peptide Bond Lenght =  4.6

STRUCTURE --->  <Structure id=protein1>
chain break at PHE  160  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.6 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at LEU  84  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.6 angstroms) exceeded
Torsion angles:  2546   2543
diff_keys lenght :  3
[('83_K_CA', '83_K_C', '84_L_N', '84_L_CA'), ('83_K_N', '83_K_CA', '83_K_C', '84_L_N'), ('83_K_C', '84_L_N', '84_L_CA', '84_L_C')] <class 'list'>
('83_K_CA', '83_K_C', '84_L_N', '84_L_CA') 83 <class 'tuple'>
('83_K_N', '83_K_CA', '83_K_C', '84_L_N') 83 <class 'tuple'>
('83_K_C', '84_L_N', '84_L_CA', '84_L_C') 83 <class 'tuple'>
diff_keys_values lenght :  3
('83_K_CA', '83_K_C', '84_L_N', '84_L_CA') ('protein1',)
('83_K_N', '83_K_CA', '83_K_C', '84_L_N') ('protein1',)
('83_K_C', '84_L_N', '84_L_CA', '84_L_C') ('protein1',)





 Peptide Bond Lenght =  4.7

STRUCTURE --->  <Structure id=protein1>
chain break at PHE  160  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.7 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at PHE  160  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.7 angstroms) exceeded
Torsion angles:  2546   2546
diff_keys lenght :  0
[] <class 'list'>
diff_keys_values lenght :  0





 Peptide Bond Lenght without differences =  4.7
get_no_tors_diff(structure1 , structure2) 4.7

Using internal_coords.IC_Chain.MaxPeptideBond = 4.7 both models have same Chain Breaks at resi 160, 245 e 296 and are described by 2546 dihedral angles.

My question now is how, if feasible , to force PyMOL to paint petide bonds of such lenght both in cartoons and lines/sticks mode.

Second question if Biopython has a tool to detect chain breaks in Structures other than just printing it out with ic_chain = internal_coords.IC_Chain(structure, verbose=True) , asking this because setting internal_coords.IC_Chain.MaxPeptideBond = 50 gets rid of all the chainbreaks but don't have a way to monitor their disappearing. Pics below show chain breaks at 4.7.

Stack Exchange Network

BioPython internal_coords module returns different dihedral angles for the (seemingly) same protein structure

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
python
proteins
protein-structure
pdb
or ask your own question.

Hot Network Questions

BioPython internal_coords module returns different dihedral angles for the (seemingly) same protein structure

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
python
proteins
protein-structure
pdb
or ask your own question.