diff --git a/src/dssp.py b/src/dssp.py index ae4d35f..b3709be 100644 --- a/src/dssp.py +++ b/src/dssp.py @@ -18,3 +18,4 @@ else: turns = get_turns(pdb_file.residues) get_helix(pdb_file.residues, turns) get_bends(pdb_file.residues) + get_chirality(pdb_file.residues) diff --git a/src/structure.py b/src/structure.py index 8d48c02..51d2a0c 100644 --- a/src/structure.py +++ b/src/structure.py @@ -1,4 +1,5 @@ import math +import numpy as np class Structure: def __init__(self, res): @@ -30,8 +31,8 @@ def get_turns(residues): for j in range(3,6): if(i+j1: + ladders[bridges[i].res_num] = k-1 + 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", residues[i].resid, angle) + bends[residues[i].resid] = angle return(bends) -def vecteur_deux_points (a, b): +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 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 @@ -171,7 +166,7 @@ def dot_product(v1,v2): return(dot_product) def position_vector(c): - vector = vecteur_deux_points([0,0,0],c) + vector = vector_from_pos([0,0,0],c) return vector def vectors_substr(v1, v2): @@ -188,11 +183,59 @@ def vector_angles(v1,v2): 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 +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("angle", residues[i].resid, "+", angle) + + if(angle<=0 and angle > -180): + sign="-" + print("angle", residues[i].resid, "-", angle) + + chirality[residues[i].resid] = sign + + return chirality +