Browse Source

full output except STRUCTURE column

Thomas Forest 5 years ago
parent
commit
651f01f626
3 changed files with 119 additions and 34 deletions
  1. 77 2
      src/atom.py
  2. 40 31
      src/dssp.py
  3. 2 1
      src/geom.py

+ 77 - 2
src/atom.py View File

28
             self.atoms[atom.atom_name] = atom
28
             self.atoms[atom.atom_name] = atom
29
             self.resid = atom.res_seq_nb
29
             self.resid = atom.res_seq_nb
30
             self.res_name = atom.res_name
30
             self.res_name = atom.res_name
31
+            self.res_letter = self.get_amino_letter(atom.res_name)
31
             self.chain_id = atom.chain_id
32
             self.chain_id = atom.chain_id
32
             self.insertion_code = atom.insertion_code
33
             self.insertion_code = atom.insertion_code
33
             self.indice = indice
34
             self.indice = indice
34
             
35
             
35
-            
36
+    def get_amino_letter(self, res_name):
37
+        code3 = ['ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLU', 'GLN', 'GLY',
38
+                 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'PHE', 'PRO', 'SER',
39
+                 'THR', 'TRP', 'TYR', 'VAL']
40
+
41
+
42
+        code1 = ['A','R','N','D','C','E','Q','G','H','I','L','K','M','F','P',
43
+                 'S','T','W','Y','V']
44
+
45
+        return code1[code3.index(res_name)]
46
+
36
     def h_bond(self, res2):
47
     def h_bond(self, res2):
37
         if("H" not in res2.atoms.keys()):
48
         if("H" not in res2.atoms.keys()):
38
             return(False)
49
             return(False)
69
             return Turn(k,residues[i].resid)
80
             return Turn(k,residues[i].resid)
70
 
81
 
71
         return False
82
         return False
72
-    
83
+
84
+    def get_bends(self, residues):
85
+        i = residues.index(self)
86
+        if i >=2 and i <len(residues)-2:
87
+            angle = math.degrees(vector_angles(vectors_substr(position_vector(residues[i].atoms["CA"].coords),
88
+                                                              position_vector(residues[i-2].atoms["CA"].coords)),
89
+                                               vectors_substr(position_vector(residues[i+2].atoms["CA"].coords),
90
+                                                              position_vector(residues[i].atoms["CA"].coords))))
91
+           
92
+            return angle
93
+        else:
94
+            return 360.0
95
+
73
     def get_bridges(self, residues):
96
     def get_bridges(self, residues):
74
         bridges = {}
97
         bridges = {}
75
         bridge = {}
98
         bridge = {}
163
                         return ladder
186
                         return ladder
164
         return False   
187
         return False   
165
 
188
 
189
+    def get_tco(self, residues):
190
+        i = residues.index(self)
191
+        if(i!=0):
192
+            res2 = residues[i-1]
193
+            CO_res1 = vector_from_pos(self.atoms["C"].coords,
194
+                                      self.atoms["O"].coords)
195
+            CO_res2 = vector_from_pos(res2.atoms["C"].coords,
196
+                                      res2.atoms["O"].coords)
197
+            angle = vector_angles(CO_res1, CO_res2)
198
+        else:
199
+            angle = math.pi/2
200
+        return(math.cos(angle))
201
+
202
+    def get_chirality(self, residues):
203
+        i = residues.index(self)
204
+        if (i >=1 and i < len(residues)-2):
205
+            chirality = {}
206
+            angle = calc_dihedral(residues[i-1].atoms["CA"].coords,
207
+                                  residues[i].atoms["CA"].coords,
208
+                                  residues[i+1].atoms["CA"].coords,
209
+                                  residues[i+2].atoms["CA"].coords)
210
+
211
+            if(angle>0 and angle<=180):
212
+                sign="+"
213
+
214
+            if(angle<=0 and angle > -180):
215
+                sign="-"
216
+
217
+        else:
218
+            angle = 360.0
219
+            sign = ''
220
+
221
+        return [angle, sign]
222
+    
223
+    def get_phi_psi(self, residues):
224
+        i = residues.index(self)  
225
+        if(i==0):
226
+            phi = 360.0
227
+        else:
228
+            phi = calc_dihedral(residues[i-1].atoms["C"].coords,
229
+                                residues[i].atoms["N"].coords,
230
+                                residues[i].atoms["CA"].coords,
231
+                                residues[i].atoms["C"].coords)
232
+        if(i==len(residues)-1):
233
+            psi = 360.0
234
+        else:
235
+            psi = calc_dihedral(residues[i].atoms["N"].coords,
236
+                                residues[i].atoms["CA"].coords,
237
+                                residues[i].atoms["C"].coords,
238
+                                residues[i+1].atoms["N"].coords)
239
+
240
+        return((phi, psi))
166
 def get_bridges(residues):
241
 def get_bridges(residues):
167
     bridges = {}
242
     bridges = {}
168
     bridge = {}
243
     bridge = {}

+ 40 - 31
src/dssp.py View File

13
           "call structure and parameters.")
13
           "call structure and parameters.")
14
 else:
14
 else:
15
     pdb_file = pdb.PDBFile(sys.argv[1])
15
     pdb_file = pdb.PDBFile(sys.argv[1])
16
-   # print(pdb_file.residues[15].atoms["C"].coord_x)
17
-    #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)
18
-    #print(pdb_file.residues[2].atoms["H"].coord_x)
19
-    #print(pdb_file.residues[2].h_bond(pdb_file.residues[40]))
20
-    #print(get_turns(pdb_file.residues))
21
-    #print(pdb_file.residues[27].h_bond(pdb_file.residues[28]))
22
-    #print(get_bridges(pdb_file.residues))
23
-
24
-
16
+  
25
     #turns = get_turns(pdb_file.residues)
17
     #turns = get_turns(pdb_file.residues)
26
     #get_helix(pdb_file.residues, turns)
18
     #get_helix(pdb_file.residues, turns)
27
     #get_bends(pdb_file.residues)
19
     #get_bends(pdb_file.residues)
28
     # get_chirality(pdb_file.residues)
20
     # get_chirality(pdb_file.residues)
29
-    bridges = get_bridges(pdb_file.residues)
30
-    ladders = get_ladders(bridges, pdb_file.residues)
31
-    sheets = get_sheets(ladders)  
21
+
32
 
22
 
33
     # print("NBRIDGES",len(bridges))
23
     # print("NBRIDGES",len(bridges))
34
 
24
 
36
 
26
 
37
     # get_bonds(pdb_file.residues)
27
     # get_bonds(pdb_file.residues)
38
 
28
 
39
-    residues = pdb_file.residues
40
-    e_min = 0
41
-    for res in residues:
42
-        ene = residues[1].h_bond(res)
43
-        if ene <= e_min:
44
-            e_min = ene
45
-            best_res = res
29
+
30
+    # e_min = 0
31
+    # for res in residues:
32
+    #     ene = residues[1].h_bond(res)
33
+    #     if ene <= e_min:
34
+    #         e_min = ene
35
+    #         best_res = res
46
     #print(residues[1].resid, best_res.resid,e_min)
36
     #print(residues[1].resid, best_res.resid,e_min)
47
 
37
 
48
     #get_phi_psi(residues)
38
     #get_phi_psi(residues)
49
     #print(residues[2].atoms, residues[0].resid, residues[1].resid)
39
     #print(residues[2].atoms, residues[0].resid, residues[1].resid)
50
     #print(get_TCO(residues[2],residues[3]))
40
     #print(get_TCO(residues[2],residues[3]))
51
 
41
 
52
-    for i,res in enumerate(residues):
53
-        # res.get_turns(residues, turns)
54
-        res.get_helix(residues)
55
-        #res.get_bridges(residues)
56
-        #res.get_ladders(residues, ladders)
57
-        #res.get_sheets(residues)
58
 
42
 
59
-for i,ladder in enumerate(ladders.values()):
60
-    print(chr(65+i))
43
+# for i,ladder in enumerate(ladders.values()):
44
+#     print(chr(65+i))
61
     
45
     
62
-for ind, sheet in sheets.items():
63
-    print(ind, sheet)
64
-    for ladder in sheet:
65
-        print(ladder['start'], ladder['end'])
46
+# for ind, sheet in sheets.items():
47
+#     print(ind, sheet)
48
+#     for ladder in sheet:
49
+#         print(ladder['start'], ladder['end'])
50
+
51
+
52
+################### OUTPUT ####################
53
+# print DSSP-style formatted header of PDB file
54
+for elem in pdb_file.get_header() : 
55
+    print(pdb_file.get_header()[elem], end="")
66
 
56
 
57
+# Get preliminary data for print loop
58
+residues = pdb_file.residues
59
+bridges = get_bridges(residues)
60
+ladders = get_ladders(bridges, residues)
61
+sheets = get_sheets(ladders)  
67
 
62
 
63
+# iterating over residues
64
+for i,res in enumerate(residues):
65
+    #res.get_turns(residues, turns)
66
+    #res.get_helix(residues)
67
+    kappa = res.get_bends(residues)
68
+    t_co = res.get_tco(residues)
69
+    alpha = res.get_chirality(residues)[0]
70
+    phi = res.get_phi_psi(residues)[0]
71
+    psi = res.get_phi_psi(residues)[1]
72
+    x_ca = res.atoms["CA"].coord_x
73
+    y_ca = res.atoms["CA"].coord_y
74
+    z_ca = res.atoms["CA"].coord_z
75
+    print(i+1, res.resid, res.chain_id, res.res_letter, round(t_co, 3), round(kappa,1),
76
+          round(alpha, 1), round(phi, 1), round(psi,1), round(x_ca,1), round(y_ca,1), round(z_ca,1))

+ 2 - 1
src/geom.py View File

1
-import math 
1
+import math
2
+import numpy as np
2
 def vector_from_pos(a, b):
3
 def vector_from_pos(a, b):
3
     xAB = b[0]-a[0]
4
     xAB = b[0]-a[0]
4
     yAB = b[1]-a[1]
5
     yAB = b[1]-a[1]