Adds Chirality evaluation

src/dssp.py

@@ -18,3 +18,4 @@ else:
     turns = get_turns(pdb_file.residues)
     get_helix(pdb_file.residues, turns)
     get_bends(pdb_file.residues)
+    get_chirality(pdb_file.residues)

src/structure.py

@@ -1,4 +1,5 @@
 import math
+import numpy as np
 
 class Structure:
     def __init__(self, res):
@@ -30,8 +31,8 @@ def get_turns(residues):
         for j in range(3,6):
             if(i+j<len(residues)):
                 if(res.h_bond(residues[i+j])<-0.5):
-                    print("TURN", i+1, i+j+1)
-                    turns[i] = Turn(j,i)
+                    print(j,"TURN", residues[i].resid, residues[i+j].resid)
+                    turns[residues[i].resid] = Turn(j,residues[i].resid)
     return(turns)
 
 def get_bridges(residues):
@@ -47,8 +48,8 @@ def get_bridges(residues):
                and residues[j].h_bond(residues[i+1])<-0.5):
                 bridge = {'res1':residues[i-1].h_bond(residues[j]),
                           'res2':residues[j].h_bond(residues[i+1]),
-                          'ipos':i,
-                          'jpos':j,
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
                           'btype':"para"}
                 # if(residues[i-1].h_bond(residues[j])+
                 #    residues[j].h_bond(residues[i+1]))<E_min:
@@ -60,8 +61,8 @@ def get_bridges(residues):
                and residues[i].h_bond(residues[j+1])<-0.5):
                 bridge = {'res1':residues[j-1].h_bond(residues[i]),
                           'res2':residues[i].h_bond(residues[j+1]),
-                          'ipos':i,
-                          'jpos':j,
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
                           'btype':"para"}
                 # if(residues[j-1].h_bond(residues[i])+
                 #    residues[i].h_bond(residues[j+1]))<E_min:
@@ -73,8 +74,8 @@ def get_bridges(residues):
                and residues[j].h_bond(residues[i])<-0.5):
                 bridge = {'res1':residues[i].h_bond(residues[j]),
                           'res2':residues[j].h_bond(residues[i]),
-                          'ipos':i,
-                          'jpos':j,
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
                           'btype':"anti"}
                 # if(residues[i].h_bond(residues[j])+
                 #    residues[j].h_bond(residues[i]))<E_min:
@@ -86,8 +87,8 @@ def get_bridges(residues):
                and residues[j-1].h_bond(residues[i+1])<-0.5):
                 bridge = {'res1':residues[i-1].h_bond(residues[j+1]),
                           'res2':residues[j-1].h_bond(residues[i+1]),
-                          'ipos':i,
-                          'jpos':j,
+                          'ipos':residues[i].resid,
+                          'jpos':residues[j].resid,
                           'btype':"anti"}
                  # if(residues[i-1].h_bond(residues[j+1])+
                  #   residues[j-1].h_bond(residues[i+1]))<E_min:
@@ -128,40 +129,34 @@ def get_helix(residues, turns):
             i+=1
     return(helix)
 
+def get_ladder(bridges):
+    ladders = {}
+    for i in range(len(bridges)):
+        k = 1
+        while bridges[i] == bridges[i+k]:
+            k+=1
+        if k>1:
+            ladders[bridges[i].res_num] = k-1
+
 def get_bends(residues):
     bends = {}
     for i in range(2,len(residues)-2):
-        # print("angle",i,calculer_angles(residues[i-2].atoms["CA"].coords,
-        #                                 residues[i].atoms["CA"].coords,
-        #                                 residues[i+2].atoms["CA"].coords))
-        
         angle = vector_angles(vectors_substr(position_vector(residues[i].atoms["CA"].coords),
                                              position_vector(residues[i-2].atoms["CA"].coords)),
                               vectors_substr(position_vector(residues[i+2].atoms["CA"].coords),
                                              position_vector(residues[i].atoms["CA"].coords)))
         if(angle>70):
             print("angle", residues[i].resid, angle)
+            bends[residues[i].resid] = angle
     return(bends)
 
-def vecteur_deux_points (a, b):
+def vector_from_pos(a, b):
     xAB = b[0]-a[0]
     yAB = b[1]-a[1]
     zAB = b[2]-a[2]
     coord_AB = [xAB,yAB,zAB]
     return coord_AB
 
-def produit_scalaire (a, b, c):
-    coord_BA = vecteur_deux_points(b, a)
-    coord_BC = vecteur_deux_points(b, c)
-    dot_product = coord_BA[0]*coord_BC[0] + coord_BA[1]*coord_BC[1] + coord_BA[2]*coord_BC[2]
-    return dot_product
-
-
-def norme_vecteur (a,b):
-    vecteur = vecteur_deux_points(a,b)
-    norme = math.sqrt(vecteur[0]**2 + vecteur[1]**2 + vecteur[2]**2)
-    return norme
-
 def vector_norm(v):
     norm = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
     return norm
@@ -171,7 +166,7 @@ def dot_product(v1,v2):
     return(dot_product)
 
 def position_vector(c):
-    vector = vecteur_deux_points([0,0,0],c)
+    vector = vector_from_pos([0,0,0],c)
     return vector
 
 def vectors_substr(v1, v2):
@@ -188,11 +183,59 @@ def vector_angles(v1,v2):
     deg_angle = rad_angle*(180/math.pi)
     return deg_angle
 
-def calculer_angles (a,b,c):
-    dot_product_ABC = produit_scalaire(a,b,c)
-    norme_BA = norme_vecteur (b,a)
-    norme_BC = norme_vecteur (b,c)
-    terme = dot_product_ABC/(abs(norme_BA)*abs(norme_BC))
-    angle_radian = math.acos(terme)
-    angle_degre = angle_radian*(180/math.pi)
-    return angle_degre
+def calc_dihedral(u1, u2, u3, u4):
+    """ Calculate dihedral angle method. From bioPython.PDB
+    (adapted to np.array)
+    Calculate the dihedral angle between 4 vectors
+    representing 4 connected points. The angle is in
+    [-pi, pi].
+
+    Adapted function of dihedral_angle_from_coordinates.py 
+    by Eric Alcaide.
+    Source : https://gist.github.com/EricAlcaide
+    URL : https://gist.github.com/EricAlcaide/30d6bfdb8358d3a57d010c9a501fda56
+    """
+    #convert coords to numpy arrays
+    u1 = np.array(u1)
+    u2 = np.array(u2)
+    u3 = np.array(u3)
+    u4 = np.array(u4)
+    
+    a1 = u2 - u1
+    a2 = u3 - u2
+    a3 = u4 - u3
+
+    v1 = np.cross(a1, a2)
+    v1 = v1 / (v1 * v1).sum(-1)**0.5
+    v2 = np.cross(a2, a3)
+    v2 = v2 / (v2 * v2).sum(-1)**0.5
+    porm = np.sign((v1 * a3).sum(-1))
+    rad = np.arccos((v1*v2).sum(-1) / ((v1**2).sum(-1) * (v2**2).sum(-1))**0.5)
+    if not porm == 0:
+        rad = rad * porm
+
+    deg_angle = rad*(180/math.pi)
+    return deg_angle
+
+
+def get_chirality(residues):
+
+    for i in range(1,len(residues)-2):
+        chirality = {}
+        angle = calc_dihedral(residues[i-1].atoms["CA"].coords,
+                              residues[i].atoms["CA"].coords,
+                              residues[i+1].atoms["CA"].coords,
+                              residues[i+2].atoms["CA"].coords)
+        
+        if(angle>0 and angle<=180):
+            sign="+"
+            print("angle", residues[i].resid, "+", angle)
+            
+        if(angle<=0 and angle > -180):
+            sign="-"
+            print("angle", residues[i].resid, "-", angle)
+            
+        chirality[residues[i].resid] = sign
+        
+    return chirality
+