dssp.py

Go to the documentation of this file.

#------------------------------------------------------------------------------
# This file is part of the OpenStructure project <www.openstructure.org>
#
# Copyright (C) 2008-2020 by the OpenStructure authors
#
# This library is free software; you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation; either version 3.0 of the License, or (at your option)
# any later version.
# This library is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this library; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
#------------------------------------------------------------------------------
"""
Utility functions to load secondary structure information from DSSP files
and assign them to entities.

Authors: Pascal Benkert, Marco Biasini
"""

import os
import tempfile,subprocess

from ost import io,mol
from ost import settings


def _SkipHeader(stream):
    line=stream.readline()
    while line:
      if line.strip().find('#')==0:
        return True
      line=stream.readline()
    return False

def _Cleanup(temp_dssp_path, pdb_path):
  if os.path.exists(temp_dssp_path):
    os.remove(temp_dssp_path)
  if os.path.exists(pdb_path):
    os.remove(pdb_path)

def _CalcRelativeSA(residue_type, absolute_sa):
  solvent_max_list=[118,317,238,243,183,262,286,154,258,228,
                    243,278,260,271,204,234,206,300,303,216] #TODO: source?
  residue_indices = "ARNDCQEGHILKMFPSTWYV"
  # cysteine bridges are marked with lower-case letters by DSSP. We don't 
  # really care which cysteines are forming covalent bonds, so let's set the 
  # one-letter-code to "C".
  if residue_type.islower():
    residue_type='C'
  if residue_indices.find(residue_type)==-1:
    raise RuntimeError('residue %s is a non-standard residue' %(residue_type))  
  else:
    rel=float(absolute_sa)/(solvent_max_list[residue_indices.find(residue_type)])
  return rel
  
#
#def AssignDSSP(ent, pdb_path="", extract_burial_status=False, tmp_dir=None, 
#               dssp_bin=None):
#  """
#  Assign secondary structure states to peptide residues in the structure. This
#  function uses the DSSP command line program.
#
#  If you already have a DSSP output file and would like to assign the secondary 
#  structure states to an entity, use :func:`LoadDSSP`.
#  
#  :param ent: The entity for which the secondary structure should be calculated
#  :type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
#  :param extract_burial_status: If true, also extract burial status and store
#                                as float-property
#                                ``relative_solvent_accessibility`` at residue
#                                level
#  :param tmp_dir: If set, overrides the default tmp directory of the
#                  operating system
#  :param dssp_bin: The path to the DSSP executable
#  :raises: :class:`~ost.settings.FileNotFound` if the dssp executable is not 
#      in the path.
#  :raises: :class:`RuntimeError` when dssp is executed with errors
#  """
#
#  if not ent.IsValid():
#    raise ValueError('model entity is not valid')
#  if ent.atom_count==0:
#    raise ValueError('model entity does not contain any atoms')
#
#  dssp_abs_path=settings.Locate(['dsspcmbi','dssp','mkdssp'], env_name='DSSP_EXECUTABLE', 
#                                explicit_file_name=dssp_bin)
#  if os.path.isdir(dssp_abs_path):
#    raise RuntimeError('"%s" is a directory. Specify path to DSSP binary' % dssp_abs_path)
#  if not os.access(dssp_abs_path, os.X_OK):
#    raise RuntimeError('"%s" is not executable' % dssp_abs_path)
#    
#  pdb_path=tempfile.mktemp(suffix=".pdb",prefix="temp_entity", dir=tmp_dir)
#  io.SavePDB(ent, pdb_path)
#  temp_dssp_path=tempfile.mktemp(suffix=".out",prefix="dssp", dir=tmp_dir)  
#
#  cmd = [dssp_abs_path, pdb_path, temp_dssp_path]
#  s = subprocess.run(cmd, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
#  if s.returncode == 0:
#    try:
#      LoadDSSP(temp_dssp_path, ent, extract_burial_status)
#      _Cleanup(temp_dssp_path, pdb_path)
#      return ent
#    except:
#      pass # continue with dssp 4 fallback
#  
#  # either dssp execution failed or output file cannot be loaded
#  # both can be the result of using dssp 4 (https://github.com/PDB-REDO/dssp) 
#  # => try fallback
#
#  # dssp 4 desperately wants a CRYST1 record (REMOVE IF NOT NEEDED ANYMORE!)
#  # Be aware, Even more stuff is needed if you don't set the CLIBD_MON env var
#  # The dssp binding has therefore been deprecated and we mimic to "old" interface
#  # using OpenStructure internal functionality
#  with open(pdb_path, 'r') as fh:
#    file_content = fh.read()
#  with open(pdb_path, 'w') as fh:
#    fh.write('CRYST1    1.000    1.000    1.000  90.00  90.00  90.00 P 1           1          ' +
#             os.linesep + file_content)
#
#  # explicitely request dssp output format
#  cmd = [dssp_abs_path, '--output-format', 'dssp', pdb_path, temp_dssp_path]
#  s_fallback = subprocess.run(cmd, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
#  
#  if s_fallback.returncode == 0:
#    try:
#      LoadDSSP(temp_dssp_path, ent, extract_burial_status)
#      _Cleanup(temp_dssp_path, pdb_path)
#      return ent
#    except:
#      pass # continue with cleanup and raise error
#      
#  _Cleanup(temp_dssp_path, pdb_path)
#  raise RuntimeError('Failed to assign DSSP')


def AssignDSSP(ent, pdb_path="", extract_burial_status=False, tmp_dir=None, 
               dssp_bin=None):
  """
  Assign secondary structure states to peptide residues in the structure. This
  function replaces the "old" AssignDSSP which relies on the DSSP command line
  program and uses OpenStructure internal functionality only. The sole purpose
  is to retain the "old" interface and you're adviced to directly use
  :func:`ost.mol.alg.AssignSecStruct` and :func:`ost.mol.alg.Accessibility`.  

  If you already have a DSSP output file and would like to assign the secondary 
  structure states to an entity, use :func:`LoadDSSP`.
  
  :param ent: The entity for which the secondary structure should be calculated
  :type ent: :class:`~ost.mol.EntityHandle` or :class:`~ost.mol.EntityView`
  :param extract_burial_status: If true, also extract burial status and store
                                as float-property
                                ``relative_solvent_accessibility`` at residue
                                level
  :param tmp_dir: If set, overrides the default tmp directory of the
                  operating system - deprecated, has no effect
  :param dssp_bin: The path to the DSSP executable - deprecated, has no effect
  """

  if not ent.IsValid():
    raise ValueError('model entity is not valid')
  if ent.atom_count==0:
    raise ValueError('model entity does not contain any atoms')

  mol.alg.AssignSecStruct(ent)
  if extract_burial_status:
    mol.alg.Accessibility(ent, algorithm=mol.alg.DSSP)
    # map float properties from Accessibility algorithm to the original
    # properties that have been set in the binding
    for r in ent.residues:
      if r.HasProp("asaAbs"):
        asa = round(r.GetFloatProp("asaAbs")) # original DSSP output is rounded
        asa_rel = _CalcRelativeSA(r.one_letter_code, asa) # there would be the
                                                          # asaRel property but
                                                          # it relates to the
                                                          # non-rounded asa
        r.SetFloatProp("solvent_accessibility", asa)
        r.SetFloatProp("relative_solvent_accessibility", asa_rel)
        if asa_rel < 0.25:
          r.SetStringProp("burial_status", 'b')
        else:
          r.SetStringProp("burial_status", 'e')


def LoadDSSP(file_name, model, extract_burial_status=False,
             entity_saved=False, calculate_relative_sa=True):
    """
    Loads DSSP output and assigns secondary structure states to the peptidic
    residues.

    If you would like to run dssp *and* assign the secondary structure,
    use :func:`AssignDSSP` instead.

    :param file_name: The filename of the DSSP output file
    :param model: The entity to which the secondary structure states should be
                  assigned
    :param extract_burial_status: If true also calculates burial status of
        residues and assigns it to the burial_status string property.
    :param calculate_relative_sa: If true also relative solvent accessibility and
        and assigns it to the relative_solvent_accessibility float property of
        the residue.
    :param entity_save: Whether the entity was saved.
    """
    if not model.IsValid():
      raise ValueError('model entity is not valid')
    if model.atom_count==0:
      raise ValueError('model entity does not contain any atoms')
    stream=open(file_name)
    if not _SkipHeader(stream):
      stream.close()
      raise RuntimeError('Ill-formatted DSSP file')

    for line in stream:
      num=line[6:10].strip()
      ins_code=line[10].strip()
      chain_name=line[11]
      solvent_accessibility=float(line[34:39].strip())
      #solvent_accessibility=line[34:39].strip()
      amino_acid=line[13]
      #print line


      if isinstance(model,mol.ChainView):
        chain=model
      else:
        chain=model.FindChain(chain_name)
        
      if not chain.IsValid():
        continue
      if num=='':
        continue
      residue=None
      try:
        if ins_code == "":
          residue=chain.FindResidue(mol.ResNum(int(num)))
        else:
          residue=chain.FindResidue(mol.ResNum(int(num),ins_code))

        # set property "burial status:
        if extract_burial_status:
          #set default (dummy) burial status for incomplete residues:
          residue.SetStringProp("burial_status", 'X')

          #handle seleno-methionine appearing as amino acid 'X' in DSSP:
          if residue.name=="MSE" and amino_acid=='X':
            amino_acid='M'

          residue.SetFloatProp("solvent_accessibility",
                                          solvent_accessibility)
          if calculate_relative_sa:
            relative_sa=_CalcRelativeSA(amino_acid,solvent_accessibility)
            residue.SetFloatProp("relative_solvent_accessibility",
                                            relative_sa)
            if relative_sa < 0.25:
              residue.SetStringProp("burial_status", 'b')
            else:
              residue.SetStringProp("burial_status", 'e')
      except Exception as e:
        print("ERROR:",e)
        continue

      rtype=line[16:17]
      rt=mol.SecStructure.COIL
      if rtype=='H':
        rt=mol.SecStructure.ALPHA_HELIX
      elif rtype=='E':
        rt=mol.SecStructure.EXTENDED
      elif rtype=='B':
        rt=mol.SecStructure.BETA_BRIDGE
      elif rtype=='S':
        rt=mol.SecStructure.BEND
      elif rtype=='T':
        rt=mol.SecStructure.TURN
      elif rtype=='I':
        rt=mol.SecStructure.PI_HELIX
      elif rtype=='G':
        rt=mol.SecStructure.THREE_TEN_HELIX
      # for corrupted DSSP files. Catch in calling routine:
      if not residue.IsValid():
        #Todo: if residues with occupancy 0 have been removed before
        #using a selection statement, they are missed here
        #IMPORTANT: asign DSSP before any selections
        stream.close()
        raise RuntimeError('Ill-formatted DSSP file: invalid residue')  
      residue.SetSecStructure(mol.SecStructure(rt))
    stream.close()