OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
dssp.py
Go to the documentation of this file.
1 #------------------------------------------------------------------------------
2 # This file is part of the OpenStructure project <www.openstructure.org>
3 #
4 # Copyright (C) 2008-2009 by the OpenStructure authors
5 #
6 # This library is free software; you can redistribute it and/or modify it under
7 # the terms of the GNU Lesser General Public License as published by the Free
8 # Software Foundation; either version 3.0 of the License, or (at your option)
9 # any later version.
10 # This library is distributed in the hope that it will be useful, but WITHOUT
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 # FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
13 # details.
14 #
15 # You should have received a copy of the GNU Lesser General Public License
16 # along with this library; if not, write to the Free Software Foundation, Inc.,
17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 #------------------------------------------------------------------------------
19 """
20 Utility functions to load secondary structure information from DSSP files
21 and assign them to entities.
22 
23 Authors: Pascal Benkert, Marco Biasini
24 """
25 
26 import os
27 import tempfile,subprocess
28 
29 from ost import io,mol
30 from ost import settings
31 
32 
33 def _SkipHeader(stream):
34  line=stream.readline()
35  while line:
36  if line.strip().find('#')==0:
37  return True
38  line=stream.readline()
39  return False
40 
41 def _Cleanup(pdb_path, temp_path, entity_saved):
42  if entity_saved and os.path.exists(pdb_path):
43  os.remove(pdb_path)
44  if os.path.exists(temp_path):
45  os.remove(temp_path)
46 
47 def _ExecuteDSSP(path, dssp_bin, temp_dir=None):
48  # use of mktemp is a safty problem (use mkstemp and provide file handle to
49  # subsequent process
50  temp_dssp_path=tempfile.mktemp(suffix=".out",prefix="dssp", dir=temp_dir)
51  dssp_abs_path=settings.Locate(['dsspcmbi','dssp','mkdssp'], env_name='DSSP_EXECUTABLE',
52  explicit_file_name=dssp_bin)
53  if os.path.isdir(dssp_abs_path):
54  raise RuntimeError('"%s" is a directory. Specify path to DSSP binary' % dssp_abs_path)
55  if not os.access(dssp_abs_path, os.X_OK):
56  raise RuntimeError('"%s" is not executable' % dssp_abs_path)
57 
58  ps=subprocess.Popen([dssp_abs_path, path, temp_dssp_path],
59  stderr=subprocess.PIPE)
60  err_lines=ps.stderr.readlines()
61  return temp_dssp_path
62 
63 
64 def _CalcRelativeSA(residue_type, absolute_sa):
65  solvent_max_list=[118,317,238,243,183,262,286,154,258,228,
66  243,278,260,271,204,234,206,300,303,216] #TODO: source?
67  residue_indices = "ARNDCQEGHILKMFPSTWYV"
68  # cysteine bridges are marked with lower-case letters by DSSP. We don't
69  # really care which cysteines are forming covalent bonds, so let's set the
70  # one-letter-code to "C".
71  if residue_type.islower():
72  residue_type='C'
73  if residue_indices.find(residue_type)==-1:
74  raise RuntimeError('residue %s is a non-standard residue' %(residue_type))
75  else:
76  rel=float(absolute_sa)/(solvent_max_list[residue_indices.find(residue_type)])
77  return rel
78 
79 
80 def AssignDSSP(ent, pdb_path="", extract_burial_status=False, tmp_dir=None,
81  dssp_bin=None):
82  """
83  Assign secondary structure states to peptide residues in the structure. This
84  function uses the DSSP command line program.
85 
86  If you already have a DSSP output file and would like to assign the secondary
87  structure states to an entity, use :func:`LoadDSSP`.
88 
89  :param ent: The entity for which the secondary structure should be calculated
90  :type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
91  :param extract_burial_status: If true, also extract burial status and store
92  as float-property
93  ``relative_solvent_accessibility`` at residue
94  level
95  :param tmp_dir: If set, overrides the default tmp directory of the
96  operating system
97  :param dssp_bin: The path to the DSSP executable
98  :raises: :class:`~ost.settings.FileNotFound` if the dssp executable is not
99  in the path.
100  :raises: :class:`RuntimeError` when dssp is executed with errors
101  """
102  entity_saved = False
103  # use of mktemp is a safty problem (use mkstemp and provide file handle to
104  # subsequent process
105  pdb_path=tempfile.mktemp(suffix=".pdb",prefix="temp_entity",
106  dir=tmp_dir)
107  io.SavePDB(ent, pdb_path)
108  entity_saved = True
109 
110  #TODO: exception handling (currently errors occuring here
111  # are handled in the parser LoadDSSP)
112  temp_dssp_path=_ExecuteDSSP(pdb_path, dssp_bin)
113  if not os.path.exists(temp_dssp_path):
114  raise RuntimeError('DSSP output file does not exist.')
115  # assign DSSP to entity
116  try:
117  LoadDSSP(temp_dssp_path, ent, extract_burial_status,
118  entity_saved)
119  except Exception, e:
120  # clean up
121  print "Exception in DSSP:", e
122  _Cleanup(pdb_path, temp_dssp_path, entity_saved)
123  raise RuntimeError(e)
124 
125  # clean up
126  #print pdb_path, temp_dssp_path
127  _Cleanup(pdb_path, temp_dssp_path, entity_saved)
128 
129  return ent
130 
131 
132 
133 def LoadDSSP(file_name, model, extract_burial_status=False,
134  entity_saved=False, calculate_relative_sa=True):
135  """
136  Loads DSSP output and assigns secondary structure states to the peptidic
137  residues.
138 
139  If you would like to run dssp *and* assign the secondary structure,
140  use :func:`AssignDSSP` instead.
141 
142  :param file_name: The filename of the DSSP output file
143  :param model: The entity to which the secondary structure states should be
144  assigned
145  :param extract_burial_status: If true also calculates burial status of
146  residues and assigns it to the burial_status string property.
147  :param calculate_relative_sa: If true also relative solvent accessibility and
148  and assigns it to the relative_solvent_accessibility float property of
149  the residue.
150  :param entity_save: Whether the entity was saved.
151  """
152  if not model.IsValid():
153  raise ValueError('model entity is not valid')
154  if model.atom_count==0:
155  raise ValueError('model entity does not contain any atoms')
156  stream=open(file_name)
157  if not _SkipHeader(stream):
158  stream.close()
159  raise RuntimeError('Ill-formatted DSSP file')
160 
161  for line in stream:
162  num=line[6:10].strip()
163  ins_code=line[10].strip()
164  chain_name=line[11]
165  solvent_accessibility=float(line[34:39].strip())
166  #solvent_accessibility=line[34:39].strip()
167  amino_acid=line[13]
168  #print line
169 
170 
171  if isinstance(model,mol.ChainView):
172  chain=model
173  else:
174  chain=model.FindChain(chain_name)
175 
176  if not chain.IsValid():
177  continue
178  if num=='':
179  continue
180  residue=None
181  try:
182  if ins_code == "":
183  residue=chain.FindResidue(mol.ResNum(int(num)))
184  else:
185  residue=chain.FindResidue(mol.ResNum(int(num),ins_code))
186 
187  # set property "burial status:
188  if extract_burial_status:
189  #set default (dummy) burial status for incomplete residues:
190  residue.SetStringProp("burial_status", 'X')
191 
192  #handle seleno-methionine appearing as amino acid 'X' in DSSP:
193  if residue.name=="MSE" and amino_acid=='X':
194  amino_acid='M'
195 
196  residue.SetFloatProp("solvent_accessibility",
197  solvent_accessibility)
198  if calculate_relative_sa:
199  relative_sa=_CalcRelativeSA(amino_acid,solvent_accessibility)
200  residue.SetFloatProp("relative_solvent_accessibility",
201  relative_sa)
202  if relative_sa < 0.25:
203  residue.SetStringProp("burial_status", 'b')
204  else:
205  residue.SetStringProp("burial_status", 'e')
206  except Exception, e:
207  print "ERROR:",e
208  continue
209 
210  rtype=line[16:17]
211  rt=mol.SecStructure.COIL
212  if rtype=='H':
213  rt=mol.SecStructure.ALPHA_HELIX
214  elif rtype=='E':
215  rt=mol.SecStructure.EXTENDED
216  elif rtype=='B':
217  rt=mol.SecStructure.BETA_BRIDGE
218  elif rtype=='S':
219  rt=mol.SecStructure.BEND
220  elif rtype=='T':
221  rt=mol.SecStructure.TURN
222  elif rtype=='I':
223  rt=mol.SecStructure.PI_HELIX
224  elif rtype=='G':
225  rt=mol.SecStructure.THREE_TEN_HELIX
226  # for corrupted DSSP files. Catch in calling routine:
227  if not residue.IsValid():
228  #Todo: if residues with occupancy 0 have been removed before
229  #using a selection statement, they are missed here
230  #IMPORTANT: asign DSSP before any selections
231  stream.close()
232  raise RuntimeError('Ill-formatted DSSP file: invalid residue')
233  residue.SetSecStructure(mol.SecStructure(rt))
234  stream.close()