123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346 |
- import math
- from structure import *
- class Atom:
-
- def dist_atoms(self, atom2):
- return(math.sqrt((self.coord_x-atom2.coord_x)**2 +
- (self.coord_y-atom2.coord_y)**2 +
- (self.coord_z-atom2.coord_z)**2))
-
-
- def __init__(self, atom_id, atom_name, res_name, chain_id,
- res_seq_nb, insertion_code, coordinates):
- self.atom_id = atom_id
- self.atom_name = atom_name
- self.res_name = res_name
- self.chain_id = chain_id
- self.res_seq_nb = res_seq_nb
- self.insertion_code = insertion_code
- self.coord_x = coordinates[0]
- self.coord_y = coordinates[1]
- self.coord_z = coordinates[2]
- self.coords = coordinates
-
- class Residue:
- def __init__(self, atoms_list):
- self.atoms = {}
- for atom in atoms_list:
- self.atoms[atom.atom_name] = atom
- self.resid = atom.res_seq_nb
- self.res_name = atom.res_name
- self.chain_id = atom.chain_id
- self.insertion_code = atom.insertion_code
-
-
- def h_bond(self, res2):
- if("H" not in res2.atoms.keys()):
- return(False)
- # dimensionnal factor, in kcal/mole
- f = 332
- # partial charges
- q1 = 0.42
- q2 = 0.20
- # distance between O-N atoms, in angströms
- r_ON = self.atoms["O"].dist_atoms(res2.atoms["N"])
- # distance between C-H atoms, in angströms
- r_CH = self.atoms["C"].dist_atoms(res2.atoms["H"])
- # distance between O-H atoms, in angströms
- r_OH = self.atoms["O"].dist_atoms(res2.atoms["H"])
- # distance between C-N atoms, in angströms
- r_CN = self.atoms["C"].dist_atoms(res2.atoms["N"])
- # 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):
- """
- Get all the turns from a specific residue.
- """
- turns = {}
- i = residues.index(self)
- k = 0
- for j in range(3,6):
- if(i+j<len(residues)):
- if(self.h_bond(residues[i+j])<-0.5):
- k = j
- if k != 0:
- #print(k,"TURN", residues[i].resid, residues[i+k].resid)
- return Turn(k,residues[i].resid)
-
- return False
-
- def get_bridges(self, residues):
- bridges = {}
- bridge = {}
- strongest_bridge = {}
- i = residues.index(self)
- if(i >= 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']<E_min):
- E_min = bridge['res1']+bridge['res2']
- strongest_bridge = bridge
- bridge = {}
- coord_bridge = [i,j]
- # 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']))
- bridges[coord_bridge[1]] = (Bridge(strongest_bridge['btype'],
- strongest_bridge['jpos'],
- strongest_bridge['ipos']))
- if(len(bridges)>0):
- return(bridges[coord_bridge[0]])
- else:
- return(False)
-
- def get_helix(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(self.get_turns(residues) and residues[i-1].get_turns(residues)):
- print(self.get_turns(residues).turn_type,"- HELIX at", residues[i].resid)
- 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)
- 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)
- ladder = {'start':consec_bridge.res_num,
- 'end':local_bridge.res_num,
- 'bridges':[consec_bridge, local_bridge]}
- return ladder
- 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_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")
-
-
-
- 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))
|