From 8532e4ecd70aa335cba47b497467db500c8c3c94 Mon Sep 17 00:00:00 2001 From: Thomas Forest Date: Tue, 17 Sep 2019 17:56:52 +0200 Subject: [PATCH] Working phi_psi --- src/dssp.py | 20 +++++++++++--------- src/structure.py | 46 +++++++++++++++++++++++++++++++++++++--------- 2 files changed, 48 insertions(+), 18 deletions(-) diff --git a/src/dssp.py b/src/dssp.py index 1d166b3..39e7be6 100644 --- a/src/dssp.py +++ b/src/dssp.py @@ -23,19 +23,17 @@ 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)) - bridges = get_bridges(pdb_file.residues) + # bridges = get_bridges(pdb_file.residues) - get_bonds(pdb_file.residues) - - print(pdb_file.residues[1].h_bond(pdb_file.residues[2])) + # get_bonds(pdb_file.residues) residues = pdb_file.residues e_min = 0 @@ -44,4 +42,8 @@ else: if ene <= e_min: e_min = ene best_res = res - print(residues[1].resid, best_res.resid,e_min) + #print(residues[1].resid, best_res.resid,e_min) + + #get_phi_psi(residues) + #print(residues[2].atoms, residues[0].resid, residues[1].resid) + print(get_TCO(residues[2],residues[3])) diff --git a/src/structure.py b/src/structure.py index 42c008b..4614579 100644 --- a/src/structure.py +++ b/src/structure.py @@ -144,8 +144,10 @@ def get_helix(residues, turns): while(i+k in turns): k+=1 temp_res.append(residues[i]) + last_res_pos = residues[i].resid+k if(k>=1): - print(k,"- HELIX at", residues[i].resid) + 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: @@ -167,8 +169,8 @@ def get_ladders(bridges, residues): 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} + 'end':bridges[i+k-1].res_num, + 'bridges':temp_bridges} i+=k-1 else: i+=1 @@ -201,10 +203,10 @@ def res_list(ladder): def get_bends(residues): bends = {} for i in range(2,len(residues)-2): - 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))) + 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 @@ -240,8 +242,7 @@ def vector_angles(v1,v2): norm_v2 = vector_norm(v2) term = dot_prod/(abs(norm_v1)*abs(norm_v2)) rad_angle = math.acos(term) - deg_angle = rad_angle*(180/math.pi) - return deg_angle + return rad_angle def calc_dihedral(u1, u2, u3, u4): """ Calculate dihedral angle method. From bioPython.PDB @@ -299,3 +300,30 @@ def get_chirality(residues): 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))