projet_court/src/structure.py

import math

class Structure:
    def __init__(self, res):
        self.res = res

class Turn(Structure):

    def __init__(self, turn_type, res):
        self.turn_type = turn_type
        Structure.res = res
                        
class Bridge(Structure):

    def __init__(self, bridge_type, res_num):
        self.bridge_type = bridge_type
        self.res_num = res_num
        Structure.res = res_num

class Helix(Structure):

    def __init__(self, residues, res_num):
        self.residues = residues
        self.res_num = res_num
        Structure.res = res_num
        
def get_turns(residues):
    turns = {}
    for i,res in enumerate(residues):
        for j in range(3,6):
            if(i+j<len(residues)):
                if(res.h_bond(residues[i+j])<-0.5):
                    print("TURN", i+1, i+j+1)
                    turns[i] = Turn(j,i)
    return(turns)

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 

            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':i,
                          'jpos':j,
                          'btype':"para"}
                # if(residues[i-1].h_bond(residues[j])+
                #    residues[j].h_bond(residues[i+1]))<E_min:
                #    E_min = residues[i-1].h_bond(residues[j])
                #    +residues[j].h_bond(residues[i+1])
                #    bridge_type = "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':i,
                          'jpos':j,
                          'btype':"para"}
                # if(residues[j-1].h_bond(residues[i])+
                #    residues[i].h_bond(residues[j+1]))<E_min:
                #    E_min = residues[j-1].h_bond(residues[i])
                #    +residues[i].h_bond(residues[j+1])
                #    bridge_type = "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':i,
                          'jpos':j,
                          'btype':"anti"}
                # if(residues[i].h_bond(residues[j])+
                #    residues[j].h_bond(residues[i]))<E_min:
                #    E_min = residues[i].h_bond(residues[j])
                #    +residues[j].h_bond(residues[i])
                #    bridge_type = "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':i,
                          'jpos':j,
                          'btype':"anti"}
                 # if(residues[i-1].h_bond(residues[j+1])+
                 #   residues[j-1].h_bond(residues[i+1]))<E_min:
                 #   E_min = residues[i-1].h_bond(residues[j+1])
                 #   +residues[j-1].h_bond(residues[i+1])
                 #   bridge_type = "anti"
            if(bridge):
                if(bridge['res1']+bridge['res2']<E_min):
                    E_min = bridge['res1']+bridge['res2']
                    strongest_bridge = bridge
                    bridge = {}
        # finally add the strongest bridge at i and j pos
        if(strongest_bridge):
            bridges.append(Bridge(strongest_bridge['btype'],
                                  strongest_bridge['ipos']))
            bridges.append(Bridge(strongest_bridge['btype'],
                                  strongest_bridge['jpos']))
    if(len(bridges)>0):
        return(bridges)
    else:
        return(False)

def get_helix(residues, turns):
    i = 1
    helix = []
    while i <= len(residues):
        if(i in turns.keys() and i-1 in turns.keys()):
            k = 0
            temp_res = []
            while(i+k in turns):
                k+=1
                temp_res.append(residues[i])
            if(k>2):
                print(k,"- HELIX at", residues[i].resid)
                helix.append(Helix(temp_res,residues[i].resid))
            i = i+k
        else:
            i+=1
    return(helix)

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)
    return(bends)

def vecteur_deux_points (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

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 = vecteur_deux_points([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)
    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