attempt of NH-->O bonds discovery + correction of get_helix

src/dssp.py

@@ -3,6 +3,10 @@ import atom
 import sys
 from structure import *
 
+def print_dssp():
+
+    empt_line=" " * 145
+
 if(len(sys.argv)<2):
     print("Not enough arguments! Run with --help to learn more about proper"
           "call structure and parameters.")
@@ -15,10 +19,29 @@ else:
     #print(get_turns(pdb_file.residues))
     #print(pdb_file.residues[27].h_bond(pdb_file.residues[28]))
     #print(get_bridges(pdb_file.residues))
+
+
     turns = get_turns(pdb_file.residues)
     get_helix(pdb_file.residues, turns)
-    get_bends(pdb_file.residues)
-    get_chirality(pdb_file.residues)
+    # get_bends(pdb_file.residues)
+    # get_chirality(pdb_file.residues)
+    # bridges = get_bridges(pdb_file.residues)
+    # ladders = get_ladders(bridges)
+    # get_sheets(ladders)
+
+    # print("NBRIDGES",len(bridges))
+
     bridges = get_bridges(pdb_file.residues)
-    ladders = get_ladders(bridges)
-    get_sheets(ladders)
+
+    get_bonds(pdb_file.residues)
+
+    print(pdb_file.residues[1].h_bond(pdb_file.residues[2]))
+
+    residues = pdb_file.residues
+    e_min = 0
+    for res in residues:
+        ene = residues[1].h_bond(res)
+        if ene <= e_min:
+            e_min = ene
+            best_res = res
+    print(residues[1].resid, best_res.resid,e_min)

src/structure.py

@@ -34,7 +34,7 @@ def get_turns(residues):
             if(i+j<len(residues)):
                 if(res.h_bond(residues[i+j])<-0.5):
                     print(j,"TURN", residues[i].resid, residues[i+j].resid)
-                    turns[residues[i].resid] = Turn(j,residues[i].resid)
+                    turns[i] = Turn(j,residues[i].resid)
     return(turns)
 
 def get_bridges(residues):
@@ -53,11 +53,6 @@ def get_bridges(residues):
                           '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:
-                #    E_min = residues[i-1].h_bond(residues[j])
-                #    +residues[j].h_bond(residues[i+1])
-                #    bridge_type = "para"
                    
             if(residues[j-1].h_bond(residues[i])<-0.5
                and residues[i].h_bond(residues[j+1])<-0.5):
@@ -66,11 +61,6 @@ def get_bridges(residues):
                           '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:
-                #    E_min = residues[j-1].h_bond(residues[i])
-                #    +residues[i].h_bond(residues[j+1])
-                #    bridge_type = "para"
 
             if(residues[i].h_bond(residues[j])<-0.5
                and residues[j].h_bond(residues[i])<-0.5):
@@ -79,11 +69,6 @@ def get_bridges(residues):
                           '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:
-                #    E_min = residues[i].h_bond(residues[j])
-                #    +residues[j].h_bond(residues[i])
-                #    bridge_type = "anti"
                    
             if(residues[i-1].h_bond(residues[j+1])<-0.5
                and residues[j-1].h_bond(residues[i+1])<-0.5):
@@ -92,11 +77,7 @@ def get_bridges(residues):
                           '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:
-                 #   E_min = residues[i-1].h_bond(residues[j+1])
-                 #   +residues[j-1].h_bond(residues[i+1])
-                 #   bridge_type = "anti"
+
             if(bridge):
                 if(bridge['res1']+bridge['res2']<E_min):
                     E_min = bridge['res1']+bridge['res2']
@@ -107,14 +88,51 @@ def get_bridges(residues):
             bridges[strongest_bridge['ipos']] = (Bridge(strongest_bridge['btype'],
                                                         strongest_bridge['ipos'],
                                                         strongest_bridge['jpos']))
-            bridges[strongest_bridge['jpos']] = (Bridge(strongest_bridge['btype'],
-                                                        strongest_bridge['jpos'],
-                                                        strongest_bridge['ipos']))
+            # bridges[strongest_bridge['jpos']] = (Bridge(strongest_bridge['btype'],
+            #                                             strongest_bridge['jpos'],
+            #                                             strongest_bridge['ipos']))
     if(len(bridges)>0):
         return(bridges)
     else:
         return(False)
 
+# def get_bonds(residues):
+#     for i,res1 in enumerate(residues):
+#         E_min = 0
+#         strongest_bridge = []
+#         for j in range(-5,6):
+#             if(i+j < len(residues)):
+#                 res2 = residues[i+j]
+#                 if(res1.h_bond(res2)<-0.5) and (res1.h_bond(res2)<E_min):
+#                     E_min = res1.h_bond(res2)
+#                     strongest_bridge = [res1, res2]
+#         if(len(strongest_bridge)>0):
+#             diff = strongest_bridge[0].resid - strongest_bridge[1].resid
+#             if( abs (diff) > 1 and abs(diff)<=5):
+#                 print(diff)
+
+def get_bonds(residues):
+    bonds = {}
+    k = 0
+    for i,res in enumerate(residues):
+        E_min = 0
+        for j in range(-5,-2):
+            if(i+j<len(residues)):
+                if(residues[i+j].h_bond(res)<-0.5):
+                    k+=1
+                    #if(res.h_bond(residues[i+j])<E_min):
+                    E_min = res.h_bond(residues[i+j])
+                    res_a = res
+                    res_b = residues[i+j]
+                    if j not in bonds:
+                        bonds[j] = []
+                    bonds[j].append([res_a.resid, res_b.resid])
+                    #turns[residues[i].resid] = Turn(j,residues[i].resid)
+    for key in bonds.keys():
+        print(key, len(bonds[key]))
+
+    print("LH",k)
+    
 def get_helix(residues, turns):
     i = 1
     helix = []