fixed indices of ladders

src/atom.py

@@ -22,7 +22,7 @@ class Atom:
         self.coords = coordinates
 
 class Residue:
-    def __init__(self, atoms_list):
+    def __init__(self, atoms_list, indice):
         self.atoms = {}
         for atom in atoms_list:
             self.atoms[atom.atom_name] = atom
@@ -30,6 +30,7 @@ class Residue:
             self.res_name = atom.res_name
             self.chain_id = atom.chain_id
             self.insertion_code = atom.insertion_code
+            self.indice = indice
             
             
     def h_bond(self, res2):
@@ -120,11 +121,12 @@ class Residue:
             # finally add the strongest bridge at i and j pos
             if(strongest_bridge):
                 bridges[coord_bridge[0]] = (Bridge(strongest_bridge['btype'],
-                                                            strongest_bridge['ipos'],
-                                                            strongest_bridge['jpos']))
+                                                   strongest_bridge['ipos'],
+                                                   strongest_bridge['jpos']))
+                
                 bridges[coord_bridge[1]] = (Bridge(strongest_bridge['btype'],
-                                                            strongest_bridge['jpos'],
-                                                            strongest_bridge['ipos']))
+                                                   strongest_bridge['jpos'],
+                                                   strongest_bridge['ipos']))
         if(len(bridges)>0):
             return(bridges[coord_bridge[0]])
         else:
@@ -145,36 +147,6 @@ class Residue:
             return(self.get_turns(residues).turn_type, residues[i].resid)
         return(False)
 
-    def get_helix2(self, residues):
-        """
-        Return if there is an helix at a given residue,
-        as well as its type.
-        """
-        i = residues.index(self)
-        # if there are no turns or it is the first residue, skip
-        if i == 0:
-            return False
-        
-        if(i in turns.keys() and i-1 in turns.keys()):
-            print(turns[i].turn_type,"- HELIX at", residues[i].resid)
-            return(turns[i].turn_type, residues[i].resid)
-        return(False)
-    
-    def get_ladders(self, residues, ladders={}):
-        #ladders = {}
-        i = residues.index(self)
-        if i != 0:
-            if self.get_bridges(residues):
-                if (residues[i-1].get_bridges(residues)):
-                    local_bridge = self.get_bridges(residues)
-                    consec_bridge = residues[i-1].get_bridges(residues)
-                    if local_bridge.bridge_type == consec_bridge.bridge_type:
-                        print("ladder", consec_bridge.res_num, local_bridge.res_num)
-                        ladders[i] = {'start':i-1,
-                                      'end':i,
-                                      'bridges':[consec_bridge, local_bridge]}
-        return ladders
-
     def get_ladder(self, residues):
         #ladders = {}
         i = residues.index(self)
@@ -189,76 +161,112 @@ class Residue:
                                   'end':local_bridge.res_num,
                                   'bridges':[consec_bridge, local_bridge]}
                         return ladder
-        return False
+        return False   
 
-    def get_sheets(self, residues):
-        """
-        Bridges between ladders.
-        Check if 1 bridge between one ladder and one or more other ladders.
-        Iterate over all residues of one ladder and check if bridge with other residues
-        of the other ladders.
-        """
-        sheets = {}
-        sheet_start = residues.index(self)
-        j=sheet_start
-        k=2
-        if(self.get_ladder(residues)):
-            local_ladder = self.get_ladder(residues)
-            while j < len(residues)-2:
-                k = 2
-                while residues[j+k].get_ladder(residues):
-                    if(residues[j+k].get_bridges(residues)):
-                        j = j+k     
-                    k+=1
-                if(k>):
-                    last = j
-                j+=1
-                print(j)
-        if(k>2):                
-            sheet_end = last      
-                        
-            print("SHEET", sheet_start, sheet_end)       
+def get_bridges(residues):
+    bridges = {}
+    bridge = {}
+    strongest_bridge = {}
+    for i in range(1,len(residues)-4):
+        E_min = 0
+        for j in range(i+2,len(residues)-1):
+            # select triplet with the minimal energy 
 
-    def get_sheets2(self):
-        """
-        Bridges between ladders.
-        Check if 1 bridge between one ladder and one or more other ladders.
-        Iterate over all residues of one ladder and check if bridge with other residues
-        of the other ladders.
-        """
-        sheets = {}
-        for ladder in ladders:
-            for ladd2 in ladders:
-                for bridge in ladders[ladder]['bridges']:
-                    if bridge.res_partner in res_list(ladders[ladd2]):
-                        print("ladder",ladders[ladder]['start'], ladders[ladder]['end'],"bridge",bridge.res_num, bridge.res_partner,
-                        "ladder 2",ladders[ladd2]['start'], ladders[ladd2]['end'])
-                #print("stop ladder 2")
-            print("stop ladder 1")            
+            if(residues[i-1].h_bond(residues[j])<-0.5
+               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':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"para"}
+                   
+            if(residues[j-1].h_bond(residues[i])<-0.5
+               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':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"para"}
 
-            
+            if(residues[i].h_bond(residues[j])<-0.5
+               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':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"anti"}
+                    
+            if(residues[i-1].h_bond(residues[j+1])<-0.5
+               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':residues[i].resid,
+                          'jpos':residues[j].resid,
+                          'i':residues[i].indice,
+                          'j':residues[j].indice,
+                          'btype':"anti"}
+
+            if(bridge):
+                if(bridge['res1']+bridge['res2']<E_min):
+                    E_min = bridge['res1']+bridge['res2']
+                    strongest_bridge = bridge
+                    coord_bridge = [i,j]
+                    bridge = {}
+        # finally add the strongest bridge at i and j pos
+        if(strongest_bridge):
+            bridges[strongest_bridge['i']] = (Bridge(strongest_bridge['btype'],
+                                                     strongest_bridge['ipos'],
+                                                     strongest_bridge['jpos'],
+                                                     [strongest_bridge['i'],
+                                                      strongest_bridge['j']]))
+            bridges[strongest_bridge['j']] = (Bridge(strongest_bridge['btype'],
+                                               strongest_bridge['jpos'],
+                                               strongest_bridge['ipos'],
+                                               [strongest_bridge['i'],
+                                               strongest_bridge['j']]))
+    if(len(bridges)>0):
+        return(bridges)
+    else:
+        return(False)
+
+def get_ladders(bridges, residues):
+    ladders = {}
+    i = 1
+    while i < len(residues):
+        k = 1
+        if i in bridges.keys():
+            temp_bridges = [bridges[i]]
+            while ((i+k in bridges.keys()) and
+                   (bridges[i].bridge_type == bridges[i+k].bridge_type)):
+                temp_bridges.append(bridges[i+k])
+                k+=1
+        if k>1:
+            #print("ladder", bridges[i].res_num, bridges[i+k-1].res_num)
+            ladders[i] = {'start':bridges[i].res_num,
+                          'end':bridges[i+k-1].res_num,
+                          'bridges':temp_bridges,
+                          'i':i,
+                          'j':i+k-1}
+            i+=k-1
+        else:
+            i+=1
+    return ladders
+
+def connected_ladders(ladd_1, ladd_2):
+    links = []
+    for bridge in ladd_1['bridges']:
+        if bridge.res_partner in res_list(ladd_2):
+            return([ladd_1['i'], ladd_1['j'], bridge.i, bridge.j,
+                    ladd_2['i'], ladd_2['j']])
+
+    return False
 
-    def get_ladders2(self, bridges, residues):
-        ladders = {}
-        i = residues.index(self)
-        if i != 0:
-            k = 1
-            if i in bridges.keys():
-                temp_bridges = [bridges[i]]
-                while ((i+k in bridges.keys()) and
-                       (bridges[i].bridge_type == bridges[i+k].bridge_type)):
-                    temp_bridges.append(bridges[i+k])
-                    k+=1
-            if k>1:
-                print("ladder", bridges[i].res_num, bridges[i+k-1].res_num)
-                ladders[i] = {'start':bridges[i].res_num,
-                              'end':bridges[i+k-1].res_num,
-                              'bridges':temp_bridges}
-                i+=k-1
-            else:
-                i+=1
-        return ladders
-    
 def get_sheets(ladders):
     """
     Bridges between ladders.
@@ -266,15 +274,30 @@ def get_sheets(ladders):
     Iterate over all residues of one ladder and check if bridge with other residues
     of the other ladders.
     """
+    ladds = [ elem for elem in ladders.values() ]
     sheets = {}
-    for ladder in ladders:
-        for ladd2 in ladders:
-            for bridge in ladders[ladder]['bridges']:
-                if bridge.res_partner in res_list(ladders[ladd2]):
-                    print("ladder",ladders[ladder]['start'], ladders[ladder]['end'],"bridge",bridge.res_num, bridge.res_partner,
-                    "ladder 2",ladders[ladd2]['start'], ladders[ladd2]['end'])
-            #print("stop ladder 2")
-        print("stop ladder 1")            
+    i = 0
+    while i < len(ladds):
+        j = i+1
+        ladd = ladds[i]
+        ladd1 = ladds[i]
+        ladd2 = ladds[j]
+
+        while j < len(ladds):
+            print(connected_ladders(ladd1, ladd2))
+            if connected_ladders(ladd1, ladd2)!=False:
+                ladd1 = ladd2
+                #print(connected_ladders(ladd1, ladd2)[4])
+                ladd2 = ladders[connected_ladders(ladd1, ladd2)[4]]
+            j+=1
+        print(ladd['i'], ladd2['i'])
+        i+=1
+    # for ladder in ladders:
+    #     for ladd2 in ladders:
+    #        if connected_ladders(ladders[ladder], ladders[ladd2]):
+    #            pass
+               #print("ladder",ladders[ladder]['i'], ladders[ladder]['j'],"bridge",bridge.i, bridge.j,
+               #      "ladder 2",ladders[ladd2]['i'], ladders[ladd2]['j'])
 
 def res_list(ladder):
     # TODO : method in ladder class

src/dssp.py

@@ -24,11 +24,11 @@ else:
 
     #turns = get_turns(pdb_file.residues)
     #get_helix(pdb_file.residues, turns)
-    get_bends(pdb_file.residues)
+    #get_bends(pdb_file.residues)
     # get_chirality(pdb_file.residues)
-    # bridges = get_bridges(pdb_file.residues)
-    # ladders = get_ladders(bridges, pdb_file.residues)
-    # get_sheets(ladders)  
+    bridges = get_bridges(pdb_file.residues)
+    ladders = get_ladders(bridges, pdb_file.residues)
+    get_sheets(ladders)  
 
     # print("NBRIDGES",len(bridges))
 
@@ -47,13 +47,15 @@ else:
 
     #get_phi_psi(residues)
     #print(residues[2].atoms, residues[0].resid, residues[1].resid)
-    print(get_TCO(residues[2],residues[3]))
+    #print(get_TCO(residues[2],residues[3]))
 
-    #turns = {}
-    ladders = {}
     for i,res in enumerate(residues):
         # res.get_turns(residues, turns)
-        #res.get_helix(residues)
+        res.get_helix(residues)
         #res.get_bridges(residues)
-        res.get_ladders(residues, ladders)
+        #res.get_ladders(residues, ladders)
         #res.get_sheets(residues)
+
+for i,ladder in enumerate(ladders.values()):
+    print(chr(65+i))
+    

src/pdb.py

@@ -50,13 +50,20 @@ class PDBFile:
                 # get the current indice of atom
                 i = self.atoms.index(atom)
                 # if this is a brand new residue
+
+                if(len(self.atoms)>1 and
+                   (atom.res_seq_nb == self.atoms[i-1].res_seq_nb
+                    and atom.insertion_code!=self.atoms[i-1].insertion_code)):
+                    self.residues.append(Residue(temp_atoms, len(self.residues)+1))
+                    temp_atoms=[]
+                
                 if(len(self.atoms)>1
-                   and atom.res_seq_nb != self.atoms[i-1].res_seq_nb):
-                    self.residues.append(Residue(temp_atoms))
+                   and (atom.res_seq_nb != self.atoms[i-1].res_seq_nb)):
+                    self.residues.append(Residue(temp_atoms, len(self.residues)+1))
                     temp_atoms=[]
                 temp_atoms.append(atom)
         # last residue
-        self.residues.append(Residue(temp_atoms))
+        self.residues.append(Residue(temp_atoms, len(self.residues)))
         # hydrogens should represent in average 50% of total atoms...
         # We use 30% threshold...
         if(count_h/len(temp_atoms)<0.30):

src/structure.py

@@ -14,11 +14,13 @@ class Turn(Structure):
                         
 class Bridge(Structure):
 
-    def __init__(self, bridge_type, res_num, res_partner):
+    def __init__(self, bridge_type, res_num, res_partner, indices):
         self.bridge_type = bridge_type
         self.res_num = res_num
         self.res_partner = res_partner
         Structure.res = res_num
+        self.i = indices[0]
+        self.j = indices[1]
 
 class Helix(Structure):