hydrogen interaction formula + turns and bridges first implementation

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)

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)
+

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="")
+        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="")
         
-
-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])

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<len(residues)):
+                if(res.h_bond(residues[i+j])):
+                    turns.append(Turn(j))
+    return(turns)
+
+def get_bridges(residues):
+    # TODO : non-overlaping check to add
+    bridges = []
+    for i,res in enumerate(residues):
+        for j, res in enumerate(residues):
+            if(i-1>=0 and i+1<len(residues)
+               and j-1>=0 and j+1<len(residues)):
+                if((residues[i-1].h_bond(residues[j])
+                    and residues[j].h_bond(residues[i+1]))
+                   or(residues[j-1].h_bond(residues[i])
+                      and residues[i].h_bond(residues[j+1]))):
+                    bridges.append(Bridge("para"))
+                    
+            if(i-1>=0 and i+1<len(residues)
+               and j-1>=0 and j+1<len(residues)):
+                if((residues[i].h_bond(residues[j])
+                    and residues[j].h_bond(residues[i]))
+                   or(residues[i-1].h_bond(residues[j+1])
+                      and residues[j-1].h_bond(residues[i+1]))):
+                    bridges.append(Bridge("anti"))
+    return(bridges)
+
+def get_helix(turns):
+    pass