From 25057924658fb61d1b4a90e823d9e5c00e75e0b7 Mon Sep 17 00:00:00 2001 From: Thomas Forest Date: Tue, 17 Sep 2019 20:15:31 +0200 Subject: [PATCH] get turns, bridges, ladders and helix oop --- src/atom.py | 221 ++++++++++++++++++++++++++++++++-- src/dssp.py | 15 ++- src/geom.py | 66 +++++++++++ src/geom.py~ | 66 +++++++++++ src/structure.py | 302 +---------------------------------------------- 5 files changed, 359 insertions(+), 311 deletions(-) create mode 100644 src/geom.py create mode 100644 src/geom.py~ diff --git a/src/atom.py b/src/atom.py index 95abcb1..b96dbb1 100644 --- a/src/atom.py +++ b/src/atom.py @@ -1,4 +1,5 @@ import math +from structure import * class Atom: def dist_atoms(self, atom2): @@ -34,21 +35,225 @@ class Residue: def h_bond(self, res2): if("H" not in res2.atoms.keys()): return(False) - # elementary charge in Coulomb - #e = 1.602176634 * 10*math.exp(-19) - # dimensionnal factor f + # dimensionnal factor, in kcal/mole f = 332 # partial charges q1 = 0.42 q2 = 0.20 - # r, distance between O-N atoms, in angströms + # distance between O-N atoms, in angströms r_ON = self.atoms["O"].dist_atoms(res2.atoms["N"]) - # r, distance between C-H atoms, in angströms + # distance between C-H atoms, in angströms r_CH = self.atoms["C"].dist_atoms(res2.atoms["H"]) - # r, distance between O-H atoms, in angströms + # distance between O-H atoms, in angströms r_OH = self.atoms["O"].dist_atoms(res2.atoms["H"]) - # r, distance between C-N atoms, in angströms + # distance between C-N atoms, in angströms r_CN = self.atoms["C"].dist_atoms(res2.atoms["N"]) - # Electrostatic interaction energy, in kcal/mole + # electrostatic interaction energy, in kcal/mole E = q1*q2*(1/r_ON + 1/r_CH - 1/r_OH - 1/r_CN)*f return(E) + + def get_turns(self, residues, turns): + """ + Get all the turns from a specific residue. + """ + i = residues.index(self) + for j in range(3,6): + if(i+j= 1 and i < len(residues)-4): + E_min = 0 + for j in range(i+2,len(residues)-1): + # select triplet with the minimal energy + + 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, + '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, + '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, + '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, + 'btype':"anti"} + + if(bridge): + if(bridge['res1']+bridge['res2']0): + return(bridges[coord_bridge[0]]) + else: + return(False) + + def get_helix(self, residues, turns): + """ + 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 not turns or i == 0: + return False + if (i <= len(residues)): + 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): + 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]} + + 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. + 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") + +def res_list(ladder): + # TODO : method in ladder class + l=[] + for i in range(ladder['start'], ladder['end']): + l.append(i) + return(l) + +def get_bends(residues): + bends = {} + for i in range(2,len(residues)-2): + angle = math.degrees(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 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("chirality", residues[i].resid, "+", angle) + + if(angle<=0 and angle > -180): + sign="-" + print("chirality", residues[i].resid, "-", angle) + + chirality[residues[i].resid] = sign + + return chirality + + +def get_phi_psi(residues): + for i in range(len(residues)): + if(i==0): + phi = 360.0 + else: + phi = calc_dihedral(residues[i-1].atoms["C"].coords, + residues[i].atoms["N"].coords, + residues[i].atoms["CA"].coords, + residues[i].atoms["C"].coords) + if(i==len(residues)-1): + psi = 360.0 + else: + psi = calc_dihedral(residues[i].atoms["N"].coords, + residues[i].atoms["CA"].coords, + residues[i].atoms["C"].coords, + residues[i+1].atoms["N"].coords) + + print("ANGLES", i, phi, psi) + +def get_TCO(res1, res2): + CO_res1 = vector_from_pos(res1.atoms["C"].coords, + res1.atoms["O"].coords) + CO_res2 = vector_from_pos(res2.atoms["C"].coords, + res2.atoms["O"].coords) + angle = vector_angles(CO_res1, CO_res2) + return(math.cos(angle)) diff --git a/src/dssp.py b/src/dssp.py index 39e7be6..6dd0c2a 100644 --- a/src/dssp.py +++ b/src/dssp.py @@ -2,6 +2,7 @@ import pdb import atom import sys from structure import * +from atom import * def print_dssp(): @@ -21,8 +22,8 @@ else: #print(get_bridges(pdb_file.residues)) - turns = get_turns(pdb_file.residues) - get_helix(pdb_file.residues, turns) + #turns = get_turns(pdb_file.residues) + #get_helix(pdb_file.residues, turns) get_bends(pdb_file.residues) # get_chirality(pdb_file.residues) # bridges = get_bridges(pdb_file.residues) @@ -47,3 +48,13 @@ else: #get_phi_psi(residues) #print(residues[2].atoms, residues[0].resid, residues[1].resid) print(get_TCO(residues[2],residues[3])) + + turns = {} + ladders = {} + for i,res in enumerate(residues): + # res.get_turns(residues, turns) + # res.get_helix(residues, turns) + res.get_bridges(residues) + res.get_ladders(residues, ladders) + + diff --git a/src/geom.py b/src/geom.py new file mode 100644 index 0000000..f69065a --- /dev/null +++ b/src/geom.py @@ -0,0 +1,66 @@ +import math +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 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 = vector_from_pos([0,0,0],c) + return vector + +def vectors_substr(v1, v2): + return ([v1[0]-v2[0], + v1[1]-v2[1], + v1[2]-v2[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) + return rad_angle + +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 diff --git a/src/geom.py~ b/src/geom.py~ new file mode 100644 index 0000000..4e96251 --- /dev/null +++ b/src/geom.py~ @@ -0,0 +1,66 @@ + +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 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 = vector_from_pos([0,0,0],c) + return vector + +def vectors_substr(v1, v2): + return ([v1[0]-v2[0], + v1[1]-v2[1], + v1[2]-v2[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) + return rad_angle + +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 diff --git a/src/structure.py b/src/structure.py index 4614579..0659596 100644 --- a/src/structure.py +++ b/src/structure.py @@ -1,5 +1,6 @@ import math import numpy as np +from geom import * class Structure: def __init__(self, res): @@ -25,305 +26,4 @@ class Helix(Structure): self.residues = residues self.res_num = res_num Structure.res = res_num - -def get_turns(residues): - # TODO : prevent redondency of overlapping turns - turns = {} - for i,res in enumerate(residues): - for j in range(3,6): - if(i+j0): - 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)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=1): - helix_size=last_res_pos - residues[i].resid - print(turns[i].turn_type,"- HELIX at", residues[i].resid) - helix.append(Helix(temp_res,residues[i].resid)) - i = i+k - else: - i+=1 - return(helix) - -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: - ladders[i] = 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. - 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") - -def res_list(ladder): - # TODO : method in ladder class - l=[] - for i in range(ladder['start'], ladder['end']): - l.append(i) - return(l) - -def get_bends(residues): - bends = {} - for i in range(2,len(residues)-2): - angle = math.degrees(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 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 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 = vector_from_pos([0,0,0],c) - return vector - -def vectors_substr(v1, v2): - return ([v1[0]-v2[0], - v1[1]-v2[1], - v1[2]-v2[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) - return rad_angle - -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("chirality", residues[i].resid, "+", angle) - - if(angle<=0 and angle > -180): - sign="-" - print("chirality", residues[i].resid, "-", angle) - - chirality[residues[i].resid] = sign - - return chirality - - -def get_phi_psi(residues): - for i in range(len(residues)): - if(i==0): - phi = 360.0 - else: - phi = calc_dihedral(residues[i-1].atoms["C"].coords, - residues[i].atoms["N"].coords, - residues[i].atoms["CA"].coords, - residues[i].atoms["C"].coords) - if(i==len(residues)-1): - psi = 360.0 - else: - psi = calc_dihedral(residues[i].atoms["N"].coords, - residues[i].atoms["CA"].coords, - residues[i].atoms["C"].coords, - residues[i+1].atoms["N"].coords) - - print("ANGLES", i, phi, psi) - -def get_TCO(res1, res2): - CO_res1 = vector_from_pos(res1.atoms["C"].coords, - res1.atoms["O"].coords) - CO_res2 = vector_from_pos(res2.atoms["C"].coords, - res2.atoms["O"].coords) - angle = vector_angles(CO_res1, CO_res2) - return(math.cos(angle))