cadscore.py

#------------------------------------------------------------------------------
# This file is part of the OpenStructure project <www.openstructure.org>
#
# Copyright (C) 2008-2009 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
#------------------------------------------------------------------------------
"""
Wrapper for the CAD score.

References:

Olechnovic K, Kulberkyte E, Venclovas C., CAD-score: A new contact area 
difference-based function for evaluation of protein structural models
Proteins. 2012 Aug 30. [Epub ahead of print]

Authors: Valerio Mariani, Alessandro Barbato
"""

import subprocess, os, tempfile, platform
from ost import settings, io

def _SetupFiles(model,reference):
  # create temporary directory
  tmp_dir_name=tempfile.mkdtemp()
  dia = 'PDB'
  for chain in model.chains:
    if chain.name==" ":
      raise RuntimeError("One of the chains in the model has no name. Cannot calculate CAD score")
    if len(chain.name) > 1:
      dia = 'CHARMM'
      break;
    for res in chain.residues:
      if len(res.name) > 3:
        dia = 'CHARMM'
        break;
  io.SavePDB(model, os.path.join(tmp_dir_name, 'model.pdb'), dialect=dia)
  dia = 'PDB'
  for chain in reference.chains:
    if chain.name==" ":
      raise RuntimeError("One of the chains in the reference has no name. Cannot calculate CAD score")
    if len(chain.name) > 1:
      dia = 'CHARMM'
      break;
    for res in chain.residues:
      if len(res.name) > 3:
        dia = 'CHARMM'
        break;
  io.SavePDB(reference, os.path.join(tmp_dir_name, 'reference.pdb'),dialect=dia)
  return tmp_dir_name

def _CleanupFiles(dir_name):
  import shutil
  shutil.rmtree(dir_name)


class CADResult:
  """
  Holds the result of running CAD
  
  .. attribute:: globalAA
    
    The global CAD's atom-atom (AA) score
  
  .. attribute:: localAA
  
    Dictionary containing local CAD's atom-atom (AA) scores. 
    
    :type: dictionary (key: chain.resnum (e.g.: A.24), value: CAD local AA score (see CAD Documentation online)
  """
  def __init__(self, globalAA, localAA):    
    self.globalAA=globalAA
    self.localAA=localAA    

def _ParseCADGlobal(lines):
  interesting_lines=lines[1]
  aa=float(interesting_lines.split()[10])
  return aa

def _ParseCADLocal(lines):
  local_aa_dict={}
  for lin in lines[11:]:
    items=lin.split()
    chain=items[0]
    resnum=int(items[1])
    key=chain+'.'+str(resnum)
    aa=float(items[2])
    local_aa_dict[key]=aa
  return local_aa_dict

def _RunCAD(max_iter, tmp_dir):
  model_filename=os.path.join(tmp_dir, 'model.pdb')
  reference_filename=os.path.join(tmp_dir, 'reference.pdb')
  if platform.system() == "Windows":
    raise RuntimeError('CAD score not available on Windows')
  else:
    cad_calc_path=settings.Locate('CADscore_calc.bash')  
    cad_read_g_path=settings.Locate('CADscore_read_global_scores.bash')  
    cad_read_l_path=settings.Locate('CADscore_read_local_scores.bash')  
    command1="\"%s\" -o %i -m \"%s\" -t \"%s\" -D \"%s\"" %(cad_calc_path, max_iter, model_filename, reference_filename, os.path.join(tmp_dir,"cadtemp"))
    command2="\"%s\" -D \"%s\"" %(cad_read_g_path, os.path.join(tmp_dir,"cadtemp"))
    command3="\"%s\" -m \"%s\" -t \"%s\" -D \"%s\" -c AA" %(cad_read_l_path, model_filename, reference_filename,os.path.join(tmp_dir,"cadtemp"))
  ps1=subprocess.Popen(command1, shell=True, stdout=subprocess.PIPE)
  ps1.wait()
  ps2=subprocess.Popen(command2, shell=True, stdout=subprocess.PIPE)
  ps2.wait()
  lines=ps2.stdout.readlines()
  try:
    globalAA=_ParseCADGlobal(lines)
  except:
    raise RuntimeError("CAD calculation failed")
  ps3=subprocess.Popen(command3, shell=True, stdout=subprocess.PIPE)
  ps3.wait()
  lines=ps3.stdout.readlines()
  try:
    localAA=_ParseCADLocal(lines)
  except:
    raise RuntimeError("CAD calculation failed")

  return CADResult(globalAA,localAA)


def CADScore(model, reference, max_iter=300):
  """
  Calculates global and local atom-atom (AA) CAD Scores
  

  :param model: The model structure. 
  :type model: :class:`~ost.mol.EntityView` or :class:`~ost.mol.EntityHandle`
  :param reference: The reference structure
  :type model2: :class:`~ost.mol.EntityView` or :class:`~ost.mol.EntityHandle`
  :param max_iter: Optional. The maximum number of iteration for the tessellation algorithm before giving up. By default 300
  :returns: The result of the CAD score calculation
  :rtype: :class:`CADResult`
  
  :raises: :class:`~ost.settings.FileNotFound` if any of the CAD score exacutables could not be located.
  :raises: :class:`RuntimeError` if the calculation failed
  """
  tmp_dir_name=_SetupFiles(model, reference)
  result=_RunCAD(max_iter, tmp_dir_name)
  _CleanupFiles(tmp_dir_name)
  return result