diff --git a/src/__pycache__/atom.cpython-37.pyc b/src/__pycache__/atom.cpython-37.pyc deleted file mode 100644 index cf2b4a7..0000000 Binary files a/src/__pycache__/atom.cpython-37.pyc and /dev/null differ diff --git a/src/__pycache__/geom.cpython-37.pyc b/src/__pycache__/geom.cpython-37.pyc deleted file mode 100644 index 9e45ca8..0000000 Binary files a/src/__pycache__/geom.cpython-37.pyc and /dev/null differ diff --git a/src/__pycache__/pdb.cpython-37.pyc b/src/__pycache__/pdb.cpython-37.pyc deleted file mode 100644 index cc64307..0000000 Binary files a/src/__pycache__/pdb.cpython-37.pyc and /dev/null differ diff --git a/src/__pycache__/structure.cpython-37.pyc b/src/__pycache__/structure.cpython-37.pyc deleted file mode 100644 index ee79b27..0000000 Binary files a/src/__pycache__/structure.cpython-37.pyc and /dev/null differ diff --git a/src/chem.py b/src/chem.py new file mode 100644 index 0000000..37bc8ee --- /dev/null +++ b/src/chem.py @@ -0,0 +1,507 @@ +import math +#from structure import * +from geom 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 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, res_partner, indices): + self.bridge_type = bridge_type + self.res_num = res_num + self.res_partner = res_partner + Structure.res = res_num + self.i = indices[0] + self.j = indices[1] + +class Helix(Structure): + + def __init__(self, residues, res_num, helix_type): + self.residues = residues + self.res_num = res_num + Structure.res = res_num + self.helix_type = helix_type + + +class Residue: + def __init__(self, atoms_list, indice): + 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.res_letter = self.get_amino_letter(atom.res_name) + self.chain_id = atom.chain_id + self.insertion_code = atom.insertion_code + self.indice = indice + + def get_amino_letter(self, res_name): + code3 = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY', + 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER', + 'THR', 'TRP', 'TYR', 'VAL'] + + + code1 = ['A','R','N','D','C','E','Q','G','H','I','L','K','M','F','P', + 'S','T','W','Y','V'] + + return code1[code3.index(res_name)] + + 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=2 and i 70): + return [angle, 'S'] + return [angle, ''] + else: + return [360.0, ''] + + 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']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].indice) + return(self.get_turns(residues).turn_type, residues[i].indice) + return(False) + + 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_tco(self, residues): + i = residues.index(self) + if(i!=0): + res2 = residues[i-1] + CO_res1 = vector_from_pos(self.atoms["C"].coords, + self.atoms["O"].coords) + CO_res2 = vector_from_pos(res2.atoms["C"].coords, + res2.atoms["O"].coords) + angle = vector_angles(CO_res1, CO_res2) + else: + angle = math.pi/2 + return(math.cos(angle)) + + def get_chirality(self, residues): + i = residues.index(self) + if (i >=1 and i < 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="+" + + if(angle<=0 and angle > -180): + sign="-" + + else: + angle = 360.0 + sign = '' + + return [angle, sign] + + def get_phi_psi(self, residues): + i = residues.index(self) + 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) + + return((phi, psi)) +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':residues[i].resid, + 'jpos':residues[j].resid, + 'i':residues[i].indice, + 'j':residues[j].indice, + '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, + 'i':residues[i].indice, + 'j':residues[j].indice, + '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, + 'i':residues[i].indice, + 'j':residues[j].indice, + '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, + 'i':residues[i].indice, + 'j':residues[j].indice, + 'btype':"anti"} + + if(bridge): + if(bridge['res1']+bridge['res2']0): + return(bridges) + else: + return(False) + +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: + #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':i, + 'j':i+k-1} + i+=k-1 + else: + i+=1 + return ladders + +def connected_ladders(ladd_1, ladd_2): + links = [] + for bridge in ladd_1['bridges']: + if bridge.res_partner in res_list(ladd_2): + return ladd_2 + + return False + +def connected_ladders2(ladd_1, ladd_2): + links = [] + for bridge in ladd_1['bridges']: + if bridge.res_partner in res_list(ladd_2): + return([ladd_1['i'], ladd_1['j'], bridge.i, bridge.j, + ladd_2['i'], ladd_2['j']]) + #return ladd_2 + + return False + + +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. + """ + ladds = [ elem for elem in ladders.values() ] + sheets = {} + + corresp = {} + for ladd1 in ladds: + for ladd2 in ladds: + if connected_ladders(ladd1, ladd2)!=False: + corresp_list = [ elem for elem in corresp.keys() ] + if ladd1['i'] not in corresp_list and ladd2['i'] not in corresp_list: + ind = len(sheets.keys()) + sheets[ind] = [] + sheets[ind].append(ladd1) + sheets[ind].append(ladd2) + corresp[ladd1['i']] = ind + corresp[ladd2['i']] = ind + elif ladd2 not in corresp_list and ladd1 in corresp_list: + sheets[corresp[ladd1['i']]].append(ladd2) + corresp[ladd2['i']] = corresp[ladd1['i']] + elif ladd1 not in corresp_list and ladd2 in corresp_list: + sheets[corresp[ladd2['i']]].append(ladd1) + corresp[ladd1['i']] = corresp[ladd2['i']] + + return sheets + +def get_sheets2(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: + if connected_ladders(ladders[ladder], ladders[ladd2]): + bridge_i = connected_ladders2(ladders[ladder], ladders[ladd2])[2] + bridge_j = connected_ladders2(ladders[ladder], ladders[ladd2])[3] + print("ladder",ladders[ladder]['i'], ladders[ladder]['j'],"bridge",bridge_i, bridge_j, + "ladder 2",ladders[ladd2]['i'], ladders[ladd2]['j']) + + +def res_list(ladder): + # TODO : method in ladder class + l=[] + for i in range(ladder['i'], ladder['j']): + l.append(i) + return(l) + +def build_turns_patterns(residues): + turns_3 = {} + turns_4 = {} + turns_5 = {} + for i,res in enumerate(residues): + turn = residues[i].get_turns(residues) + if(turn): + for k in range(turn.turn_type): + if turn.turn_type == 3: + turns_3[i+1+k] = turn.turn_type + if turn.turn_type == 4: + turns_4[i+1+k] = turn.turn_type + if turn.turn_type == 5: + turns_5[i+1+k] = turn.turn_type + return[turns_3, turns_4, turns_5] + +def build_helix_patterns(residues): + helix_3 = {} + helix_4 = {} + helix_5 = {} + for i,res in enumerate(residues): + helix = residues[i].get_helix(residues) + if(helix): + helix_type = residues[i].get_helix(residues)[0] + helix_pos = residues[i].get_helix(residues)[1] + #print("TYPE", helix_type) + for k in range(helix_type): + if helix_type == 3: + helix_3[i+1+k] = "G" + if helix_type == 4: + helix_4[i+1+k] = "H" + if helix_type == 5: + helix_5[i+1+k] = "I" + #print(helix_3) + return[helix_3, helix_4, helix_5] + +def print_helix_pattern(residues, res, helix): + i = residues.index(res)+1 + if i in helix.keys(): + return (helix[i]) + else: + return(' ') + +def print_turn_pattern(residues, res, turns): + i = residues.index(res)+1 + if i in turns.keys() and not i-1 in turns.keys(): + return(">") + elif i in turns.keys() and i-1 in turns.keys(): + return(turns[i]) + elif i not in turns.keys() and i-1 in turns.keys(): + return("<") + else: + return(' ') + +def raw_ladders_print(ladders, ): + for ladd1 in ladders: + ladd_1 = ladders[ladd1] + for ladd2 in ladders: + ladd_2 = ladders[ladd2] + for bridge in ladd_1['bridges']: + if bridge.j in res_list(ladd_2): + print("ladder", ladd_1['i'],"-",ladd_1['j'], "|", bridge.i, "...", bridge.j, "| ladder", ladd_2['i'], ladd_2['j']) + diff --git a/src/atom.py b/src/chem.py~ similarity index 95% rename from src/atom.py rename to src/chem.py~ index c3db1b2..5b9c05a 100644 --- a/src/atom.py +++ b/src/chem.py~ @@ -1,5 +1,7 @@ import math -from structure import * +#from structure import * +from geom import * + class Atom: def dist_atoms(self, atom2): @@ -21,6 +23,35 @@ class Atom: self.coord_z = coordinates[2] self.coords = coordinates +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, res_partner, indices): + self.bridge_type = bridge_type + self.res_num = res_num + self.res_partner = res_partner + Structure.res = res_num + self.i = indices[0] + self.j = indices[1] + +class Helix(Structure): + + def __init__(self, residues, res_num, helix_type): + self.residues = residues + self.res_num = res_num + Structure.res = res_num + self.helix_type = helix_type + + class Residue: def __init__(self, atoms_list, indice): self.atoms = {} diff --git a/src/dssp.py b/src/dssp.py index 3f1c9a2..84dbb3c 100644 --- a/src/dssp.py +++ b/src/dssp.py @@ -1,104 +1,86 @@ import pdb -import atom import sys from structure import * -from atom import * +from chem import * +import argparse -def print_dssp(): +from datetime import date - empt_line=" " * 145 +def print_dssp(pdb_file, verbose=False, raw_ladd=False): + # gets export date + export_date = date.today().strftime("%d/%m/%Y") + + ################### OUTPUT #################### + # print DSSP-style formatted header of PDB file + print("################### DSSP OUTPUT ####################\n####", + "Exported on",export_date, + "\n####################################################") + for elem in pdb_file.get_header() : + print(pdb_file.get_header()[elem], end="") + # Get preliminary data for print loop + residues = pdb_file.residues + bridges = get_bridges(residues) + ladders = get_ladders(bridges, residues) + sheets = get_sheets(ladders) + helix = build_helix_patterns(residues) + # iterating over residues + # print header of table + print(' # RESIDUE AA STRUCTURE BP1 TCO KAPPA ALPHA ' \ + 'PHI PSI X-CA Y-CA Z-CA') + for i,res in enumerate(residues): + #res.get_turns(residues, turns) + kappa = res.get_bends(residues)[0] + bend_symbol = res.get_bends(residues)[1] + t_co = res.get_tco(residues) + alpha = res.get_chirality(residues)[0] + chirality = res.get_chirality(residues)[1] + phi = res.get_phi_psi(residues)[0] + psi = res.get_phi_psi(residues)[1] + x_ca = res.atoms["CA"].coord_x + y_ca = res.atoms["CA"].coord_y + z_ca = res.atoms["CA"].coord_z + turns = build_turns_patterns(residues) + turn_3 = print_turn_pattern(residues, res, turns[0]) + turn_4 = print_turn_pattern(residues, res, turns[1]) + turn_5 = print_turn_pattern(residues, res, turns[2]) + helix_3 = print_helix_pattern(residues, res, helix[0]) + helix_4 = print_helix_pattern(residues, res, helix[1]) + helix_5 = print_helix_pattern(residues, res, helix[2]) + if i in bridges.keys(): + bp1 = bridges[i].j + else: + bp1 = 0 + + print('{0: >5}{1: >5}{2: >2}{3: >2}{4: >2}{5: >2}{6: >2}{7: >2}'\ + '{8: >2}{9: >2}{10: >2}{11: >1}{12: >5}{13: >7}{14: >7}' \ + '{15: >7}{16: >7}{17: >7}{18: >7}{19: >7}'\ + '{20: >7}'.format(i+1, res.resid, + res.chain_id, + res.res_letter, + helix_3, helix_4, + helix_5, turn_3, + turn_4, turn_5, bend_symbol, + chirality, bp1,round(t_co, 3), + round(kappa,1), round(alpha, 1), + round(phi, 1), round(psi,1), + round(x_ca,1), round(y_ca,1), + round(z_ca,1))) + if(raw_ladd): + raw_ladders_print(ladders) + +parser = argparse.ArgumentParser() if(len(sys.argv)<2): print("Not enough arguments! Run with --help to learn more about proper" "call structure and parameters.") else: + parser.add_argument("file", type=str, help="input file in PDB format") + parser.add_argument("-v","--verbose", help="add ladders verbose output", + action="store_true") + args = parser.parse_args() pdb_file = pdb.PDBFile(sys.argv[1]) - - #turns = get_turns(pdb_file.residues) - #get_helix(pdb_file.residues, turns) - #get_bends(pdb_file.residues) - # get_chirality(pdb_file.residues) - - - # print("NBRIDGES",len(bridges)) - - # bridges = get_bridges(pdb_file.residues) - - # get_bonds(pdb_file.residues) - - - # e_min = 0 - # for res in residues: - # ene = residues[1].h_bond(res) - # if ene <= e_min: - # e_min = ene - # best_res = res - #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])) - - -# for i,ladder in enumerate(ladders.values()): -# print(chr(65+i)) - -# for ind, sheet in sheets.items(): -# print(ind, sheet) -# for ladder in sheet: -# print(ladder['start'], ladder['end']) - - -################### OUTPUT #################### -# print DSSP-style formatted header of PDB file -for elem in pdb_file.get_header() : - print(pdb_file.get_header()[elem], end="") - -# Get preliminary data for print loop -residues = pdb_file.residues -bridges = get_bridges(residues) -ladders = get_ladders(bridges, residues) -sheets = get_sheets(ladders) -helix = build_helix_patterns(residues) -# iterating over residues - - -for i,res in enumerate(residues): - #res.get_turns(residues, turns) - kappa = res.get_bends(residues)[0] - bend_symbol = res.get_bends(residues)[1] - t_co = res.get_tco(residues) - alpha = res.get_chirality(residues)[0] - chirality = res.get_chirality(residues)[1] - phi = res.get_phi_psi(residues)[0] - psi = res.get_phi_psi(residues)[1] - x_ca = res.atoms["CA"].coord_x - y_ca = res.atoms["CA"].coord_y - z_ca = res.atoms["CA"].coord_z - turns = build_turns_patterns(residues) - turn_3 = print_turn_pattern(residues, res, turns[0]) - turn_4 = print_turn_pattern(residues, res, turns[1]) - turn_5 = print_turn_pattern(residues, res, turns[2]) - helix_3 = print_helix_pattern(residues, res, helix[0]) - helix_4 = print_helix_pattern(residues, res, helix[1]) - helix_5 = print_helix_pattern(residues, res, helix[2]) - if i in bridges.keys(): - bp1 = bridges[i].j + if(args.verbose): + print_dssp(pdb_file, raw_ladd=True) else: - bp1 = 0 - print(i+1, res.resid, res.chain_id, res.res_letter, helix_3, helix_4, helix_5, turn_3, turn_4, turn_5, bend_symbol, chirality, bp1,round(t_co, 3), round(kappa,1), - round(alpha, 1), round(phi, 1), round(psi,1), round(x_ca,1), round(y_ca,1), round(z_ca,1)) - - -for ladd1 in ladders: - ladd_1 = ladders[ladd1] - for ladd2 in ladders: - ladd_2 = ladders[ladd2] - for bridge in ladd_1['bridges']: - if bridge.j in res_list(ladd_2): - print("ladder", ladd_1['i'],"-",ladd_1['j'], "|", bridge.i, "...", bridge.j, "| ladder", ladd_2['i'], ladd_2['j']) - -# for ind, sheet in sheets.items(): -# print(ind, sheet) -# for ladder in sheet: -# print(ladder['start'], ladder['end']) + print_dssp(pdb_file) + diff --git a/src/geom.py~ b/src/geom.py~ deleted file mode 100644 index 4e96251..0000000 --- a/src/geom.py~ +++ /dev/null @@ -1,66 +0,0 @@ - -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/pdb.py b/src/pdb.py index cb1b216..ce616aa 100644 --- a/src/pdb.py +++ b/src/pdb.py @@ -1,5 +1,5 @@ # custom imports -from atom import * +from chem import * # import libs import sys import pymol @@ -67,9 +67,7 @@ class PDBFile: # hydrogens should represent in average 50% of total atoms... # We use 30% threshold... if(count_h/len(temp_atoms)<0.30): - #if(output_pdb==None): - print("Need to add hydrogens ! If you want the modified PDB file, " - "please use the -o output.pdb argument") + #Need to add hydrogens self.add_hydrogens(filename) def check_hydrogens(self, atoms):