From e80fdb57ab18e6bcded6b71ea55a21d63cdae30c Mon Sep 17 00:00:00 2001 From: Thomas Forest Date: Fri, 13 Sep 2019 19:09:15 +0200 Subject: [PATCH] hydrogen interaction formula + turns and bridges first implementation --- src/atom.py | 34 +++++++++++++++++++++ src/dssp.py | 18 +++++++++++ src/pdb.py | 79 ++++++++++++++++++++++++++++-------------------- src/structure.py | 49 ++++++++++++++++++++++++++++++ 4 files changed, 147 insertions(+), 33 deletions(-) create mode 100644 src/dssp.py create mode 100644 src/structure.py diff --git a/src/atom.py b/src/atom.py index 949141f..9337f02 100644 --- a/src/atom.py +++ b/src/atom.py @@ -1,4 +1,12 @@ +import math 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, coordinates): self.atom_id = atom_id @@ -15,5 +23,31 @@ class Residue: 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 + def h_bond(self, res2): + if("H" not in res2.atoms.keys()): + return(False) + # elementary charge in Coulomb + #e = 1.602176634 * 10*math.exp(-19) + # dimensionnal factor f + f = 332 + q1 = 0.42 + q2 = 0.20 + # r, distance between O-N atoms, in angströms + r_ON = self.atoms["O"].dist_atoms(res2.atoms["N"]) + # r, distance between C-H atoms, in angströms + r_CH = self.atoms["C"].dist_atoms(res2.atoms["H"]) + # r, distance between O-H atoms, in angströms + r_OH = self.atoms["O"].dist_atoms(res2.atoms["H"]) + # r, 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 + if(E<-0.5): + return(True) + else: + return(False) diff --git a/src/dssp.py b/src/dssp.py new file mode 100644 index 0000000..ab0c5e3 --- /dev/null +++ b/src/dssp.py @@ -0,0 +1,18 @@ +import pdb +import atom +import sys +from structure import * + +if(len(sys.argv)<2): + print("Not enough arguments! Run with --help to learn more about proper" + "call structure and parameters.") +else: + pdb_file = pdb.PDBFile(sys.argv[1]) + # print(pdb_file.residues[15].atoms["C"].coord_x) + #print(pdb_file.residues[2].atoms["N"].res_seq_nb, pdb_file.residues[2].atoms["N"].coord_x, pdb_file.residues[2].atoms["N"].coord_y, pdb_file.residues[2].atoms["N"].coord_z) + #print(pdb_file.residues[2].atoms["H"].coord_x) + #print(pdb_file.residues[2].h_bond(pdb_file.residues[40])) + #print(get_turns(pdb_file.residues)) + #print(pdb_file.residues[27].h_bond(pdb_file.residues[28])) + get_bridges(pdb_file.residues) + diff --git a/src/pdb.py b/src/pdb.py index 2ee6b6a..9e9e2c5 100644 --- a/src/pdb.py +++ b/src/pdb.py @@ -1,33 +1,33 @@ # custom imports from atom import * +# import libs import sys -import pymol -#import collections +import pymol + class PDBFile: - def getContent(self, filename): + def get_content(self, filename): with open(filename) as f: return(f.readlines()) - def getHeader(self): - #Metadata = collections.namedtuple("Metadata", ["header", "compound", "source", "author"]) - Metadata = {} + def get_header(self): + metadata = {} for line in self.rawLines: # no need to continue if meta are complete - if(len(Metadata) <4): + if(len(metadata) <4): if(line[0:10] == "HEADER "): - Metadata['header']=line + metadata['header']=line elif(line[0:10] == "COMPND 2"): - Metadata['compound']=line + metadata['compound']=line elif(line[0:10] == "SOURCE 2"): - Metadata['source']=line + metadata['source']=line elif(line[0:10] == "AUTHOR "): - Metadata['author']=line + metadata['author']=line else: # if meta are complete, stop parsing break - return(Metadata) + return(metadata) - def getAtoms(self, filename): + def get_atoms(self, filename): self.atoms = [] self.residues = [] temp_atoms = [] @@ -56,10 +56,12 @@ class PDBFile: temp_atoms.append(atom) # last residue self.residues.append(Residue(temp_atoms)) - # hydrogens should represent in average 50% of total atoms... We use 30% threshold... + # 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") + print("Need to add hydrogens ! If you want the modified PDB file, " + "please use the -o output.pdb argument") self.add_hydrogens(filename) def check_hydrogens(self, atoms): @@ -71,25 +73,36 @@ class PDBFile: pymol.cmd.load(filename) pymol.cmd.select("nitrogens",'name n') pymol.cmd.h_add("nitrogens") - pymol.stored.pos = [] - pymol.cmd.iterate_state(1, "hydrogens", 'stored.pos.append([name,resi,x,y,z])') + pymol.stored.pos = {} + pymol.cmd.iterate_state(1, "hydrogens", + 'stored.pos[resi] = []; ' \ + 'stored.pos[resi].append(name,resi,x,y,z) '\ + 'if resi not in stored.pos.keys() ' \ + 'else ' \ + 'stored.pos[resi].append([name,resi,x,y,z])') if(output_pdb!=None): pymol.cmd.save(output_file) + for resi in self.residues: + if(str(resi.resid) in pymol.stored.pos.keys()): + # iterate over residue H atoms + for hydrogen in pymol.stored.pos[str(resi.resid)]: + atom = Atom(atom_id = len(self.atoms)+1, + atom_name = "H", + res_name = resi.res_name, + chain_id = resi.chain_id, + res_seq_nb = resi.resid, + coordinates = [hydrogen[2], + hydrogen[3], + hydrogen[4]]) + # add hydrogen to atoms list + self.atoms.append(atom) + resi.atoms[atom.atom_name] = atom return(pymol.stored.pos) - - def __init__(self, filename, output_pdb=None): - self.rawLines = self.getContent(filename) - self.Metadata = self.getHeader() - self.getAtoms(filename) - # for elem in self.Metadata : - # print(self.Metadata[elem], end="") - -if __name__ == "__main__": - if(len(sys.argv)<2): - print("Not enough arguments! Run with --help to learn more about proper" - "call structure and parameters.") - else: - pdbFile = PDBFile(sys.argv[1]) - print(pdbFile.residues[15].atoms["C"].coord_x) - print(pdbFile.add_hydrogens(sys.argv[1])[3]) + def __init__(self, filename, output_pdb=None): + self.rawLines = self.get_content(filename) + self.metadata = self.get_header() + self.get_atoms(filename) + # for elem in self.metadata : + # print(self.metadata[elem], end="") + diff --git a/src/structure.py b/src/structure.py new file mode 100644 index 0000000..a6fecb5 --- /dev/null +++ b/src/structure.py @@ -0,0 +1,49 @@ + +class Turn: + + def __init__(self, turn_type): + self.turn_type = turn_type + +class Bridge: + + def __init__(self, bridge_type): + self.bridge_type = bridge_type + +class Helix: + + def __init__(self, residues): + self.residues = residues + +def get_turns(residues): + turns = [] + for i,res in enumerate(residues): + for j in range(3,6): + if(i+j=0 and i+1=0 and j+1=0 and i+1=0 and j+1