From a242afd6621e5fc879cb4aa6a6f163c4accdf389 Mon Sep 17 00:00:00 2001 From: Thomas Forest Date: Wed, 18 Sep 2019 03:34:06 +0200 Subject: [PATCH] fixed indices of ladders --- src/atom.py | 241 ++++++++++++++++++++++++++--------------------- src/dssp.py | 20 ++-- src/pdb.py | 13 ++- src/structure.py | 4 +- 4 files changed, 156 insertions(+), 122 deletions(-) diff --git a/src/atom.py b/src/atom.py index e869a2a..92a78da 100644 --- a/src/atom.py +++ b/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']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 diff --git a/src/dssp.py b/src/dssp.py index bc4c156..955df8d 100644 --- a/src/dssp.py +++ b/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)) + diff --git a/src/pdb.py b/src/pdb.py index 19615da..cb1b216 100644 --- a/src/pdb.py +++ b/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): diff --git a/src/structure.py b/src/structure.py index 0659596..ceb7b08 100644 --- a/src/structure.py +++ b/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):