Add bends discovery from position vectors

src/atom.py

@@ -16,7 +16,8 @@ class Atom:
         self.res_seq_nb = res_seq_nb
         self.coord_x = coordinates[0]
         self.coord_y = coordinates[1]
-        self.coord_z = coordinates[2]             
+        self.coord_z = coordinates[2]
+        self.coords = coordinates
 
 class Residue:
     def __init__(self, atoms_list):

src/dssp.py

@@ -14,5 +14,7 @@ else:
     #print(pdb_file.residues[2].h_bond(pdb_file.residues[40]))
     #print(get_turns(pdb_file.residues))
     #print(pdb_file.residues[27].h_bond(pdb_file.residues[28]))
-    print(get_bridges(pdb_file.residues))
-
+    #print(get_bridges(pdb_file.residues))
+    turns = get_turns(pdb_file.residues)
+    get_helix(pdb_file.residues, turns)
+    get_bends(pdb_file.residues)

src/structure.py

@@ -1,27 +1,37 @@
+import math
 
-class Turn:
+class Structure:
+    def __init__(self, res):
+        self.res = res
 
-    def __init__(self, turn_type, res_num):
+class Turn(Structure):
+
+    def __init__(self, turn_type, res):
         self.turn_type = turn_type
+        Structure.res = res
                         
-class Bridge:
+class Bridge(Structure):
 
     def __init__(self, bridge_type, res_num):
         self.bridge_type = bridge_type
         self.res_num = res_num
+        Structure.res = res_num
 
-class Helix:
+class Helix(Structure):
 
-    def __init__(self, residues):
+    def __init__(self, residues, res_num):
         self.residues = residues
+        self.res_num = res_num
+        Structure.res = res_num
         
 def get_turns(residues):
-    turns = []
+    turns = {}
     for i,res in enumerate(residues):
         for j in range(3,6):
             if(i+j<len(residues)):
-                if(res.h_bond(residues[i+j]<-0.5)):
-                    turns.append(Turn(j,i))
+                if(res.h_bond(residues[i+j])<-0.5):
+                    print("TURN", i+1, i+j+1)
+                    turns[i] = Turn(j,i)
     return(turns)
 
 def get_bridges(residues):
@@ -101,6 +111,88 @@ def get_bridges(residues):
         return(False)
 
 def get_helix(residues, turns):
-    for i in range(residues):
-        for j in range(i+1,i+5):
-            if
+    i = 1
+    helix = []
+    while i <= len(residues):
+        if(i in turns.keys() and i-1 in turns.keys()):
+            k = 0
+            temp_res = []
+            while(i+k in turns):
+                k+=1
+                temp_res.append(residues[i])
+            if(k>2):
+                print(k,"- HELIX at", i)
+                helix.append(Helix(temp_res,i))
+            i = i+k
+        else:
+            i+=1
+    return(helix)
+
+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", i+1, angle)
+    return(bends)
+
+def vecteur_deux_points (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
+
+def dot_product(v1,v2):
+    dot_product = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
+    return(dot_product)
+
+def position_vector(c):
+    vector = vecteur_deux_points([0,0,0],c)
+    return vector
+
+def vectors_substr(v1, v2):
+    return ([v2[0]-v1[0],
+             v2[1]-v1[1],
+             v2[2]-v1[2]])
+
+def vector_angles(v1,v2):
+    dot_prod = dot_product(v1,v2)
+    norm_v1 = vector_norm(v1)
+    norm_v2 = vector_norm(v2)
+    term = dot_prod/(abs(norm_v1)*abs(norm_v2))
+    rad_angle = math.acos(term)
+    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