pdb_file_parser.hpp

Go to the documentation of this file.

// -*- mode: c++; indent-tabs-mode: nil; -*-
//
//The Biomolecule Toolkit (BTK) is a C++ library for use in the
//modeling, analysis, and design of biological macromolecules.
//Copyright (C) 2006, Tim Robertson <kid50@users.sourceforge.net>
//
//This program 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 2.1 of the License, or (at
//your option) any later version.
//
//This program 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 program; if not, write to the Free Software Foundation,
//Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

#ifndef BTK_PDB_FILE_PARSER_HPP
#define BTK_PDB_FILE_PARSER_HPP

#include <ios>
#include <sstream>
#include <fstream>
#include <cctype>
#include <string>
#include <algorithm>

#include <btk/core/math/btk_vector.hpp>
#include <btk/core/common/exceptions.hpp>
#include <btk/core/io/logging.hpp>
#include <btk/core/io/default_type_system.hpp>

namespace BTK {
namespace EXCEPTIONS {

//
// Possible exception classes for PDBSystem class.
//
NEW_BTK_EXCEPTION_TYPE(BTKFileOpenFailure,BTKRuntimeError);
NEW_BTK_EXCEPTION_TYPE(BTKFileReadFailure,BTKRuntimeError);
NEW_BTK_EXCEPTION_TYPE(BTKCorruptPDBFile,BTKRuntimeError);

} // namespace EXCEPTIONS

namespace IO {

using BTK::EXCEPTIONS::BTKInvalidArgument;
using BTK::EXCEPTIONS::BTKFileOpenFailure;
using BTK::EXCEPTIONS::BTKFileReadFailure;
using BTK::EXCEPTIONS::BTKCorruptPDBFile;

template <typename TypeSystem = DefaultTypeSystem>
class PDBFileParser
{
public:
  typedef TypeSystem type_system;
  
  typedef typename type_system::element_dictionary element_dictionary;
  typedef typename type_system::atom_dictionary atom_dictionary;
  typedef typename type_system::monomer_dictionary monomer_dictionary;

  typedef typename element_dictionary::id_type element_id_type;
  typedef typename atom_dictionary::id_type atom_id_type;
  typedef typename monomer_dictionary::id_type monomer_id_type; 

  PDBFileParser(type_system const & ts = type_system(),
                bool use_dynamic_types = true) : 
    _dictionary(ts), _dynamic_types(true) {}
 
  PDBFileParser(std::string const & filename,
                int which_model = 0,
                type_system const & ts = type_system(),
                bool use_dynamic_types = true) :
    _filename(filename), _cur_model(0), _dictionary(ts), 
    _dynamic_types(use_dynamic_types)
  {
    if (!open(filename,which_model)) {
      BTK_THROW_MSG(BTKFileReadFailure,
                    "Couldn't find structure in PDB file");
    }
  }

  PDBFileParser(PDBFileParser const & src) :
    _filename(src._filename), _header(src._header), _cur_model(src._cur_model), 
    _dictionary(src._dictionary), _dynamic_types(src._dynamic_types)
  {
    open_file();

    if (_cur_model > 0 && !find_model(_cur_model)) {
      BTK_THROW_MSG(BTKFileReadFailure,
                    "Couldn't find requested model in PDB file");
    }
  }

  virtual ~PDBFileParser() { close(); }

  bool open(std::string filename,
            int which_model = 0);
  
  void close() 
  {
    _pdb_fs.close();
    _pdb_fs.clear();
    _cur_model = 0;
    _cur_line = "";
  }

  std::string const & get_header() const { return _header; }

  std::string const & get_filename() const { return _filename; }
  
  template <typename PolymerSystemType>
  bool get_model(int which,
                 PolymerSystemType & model_out,
                 bool load_hetatoms = true);

  template <typename PolymerSystemType>
  bool get_current_model(PolymerSystemType & model_out, 
                         bool load_hetatoms = true)
  {
    return get_model(_cur_model,model_out,load_hetatoms);
  }

  template <typename PolymerSystemType>
  bool get_next_model(PolymerSystemType & model_out,
                      bool load_hetatoms = true)
  {
    return get_model(_cur_model + 1,model_out,load_hetatoms);
  }

  type_system dictionary() {return _dictionary;}
  type_system const dictionary() const { return _dictionary; }
  
protected:

  std::string eat_whitespace(std::string const & s) const;

  void open_file();
  void rewind_file();
  void save_header();
  bool find_model(int which);
  bool find_next_atom(bool find_hetatoms);

  template <typename AtomType>
  AtomType parse_atom(std::string const & atom_record,
                      bool is_hetatom = false);


  template <typename PolymerSystemType>
  void parse_current_model(PolymerSystemType & model_out,
                           bool load_hetatoms);

  template <typename AtomIterator, typename PolymerSystemType>
  void save_hetatoms(AtomIterator het_begin, AtomIterator het_end,
                     PolymerSystemType & model_out);
  
private:
  std::string _filename;
  std::ifstream _pdb_fs;
  std::string _header;
  int _cur_model;
  std::string _cur_line;

  type_system _dictionary;
  const bool _dynamic_types;
};

} // namespace IO
} // namespace BTK

// Implementations

template <typename TypeSystem>
template <typename PolymerSystemType>
bool
BTK::IO::PDBFileParser<TypeSystem>:: 
get_model(int which,
          PolymerSystemType & model_out,
          bool load_hetatoms)
{
  if (find_model(which)) {
    parse_current_model(model_out,load_hetatoms);
    return true;
  } else 
    return false;
}

template <typename TypeSystem>
bool 
BTK::IO::PDBFileParser<TypeSystem>::
open(std::string filename,
     int which_model)
{
  close();
  
  open_file();
  save_header();

  if (!_pdb_fs) {
    BTK_THROW_MSG(BTKFileReadFailure,"Couldn't read past PDB header");
  }
  
  // first, check what model we're in
  if (_cur_line.substr(0,5) == "MODEL") {
    _cur_model = std::atoi(_cur_line.substr(10).c_str());
  }
  
  // if the user requested a model, and it isn't the one we're in,
  // we need to find the requested model.
  if (which_model > 0 && _cur_model != which_model) {
    if (!find_model(which_model)) return false;
  } 
  
  return true;
}

template <typename TypeSystem>
std::string 
BTK::IO::PDBFileParser<TypeSystem>::
eat_whitespace(std::string const & s) const
{
  std::string tmp(s);
  
  tmp.erase(std::remove_if(tmp.begin(),tmp.end(),
                           static_cast<int (*)(int)>(std::isspace)),
            tmp.end());
  
  return tmp;
}

template <typename TypeSystem>
void 
BTK::IO::PDBFileParser<TypeSystem>::
open_file()
{
  if (!_filename.empty()) {
    _pdb_fs.open(_filename.c_str());
    
    if (_pdb_fs.fail()) {
        BTK_THROW(BTKFileOpenFailure);
    }
  } else {
    BTK_THROW_MSG(BTKInvalidArgument,
                  "Empty filename passed to PDBFileParser");
  }
  
  _cur_model = 0;
  _cur_line = "";
}

template <typename TypeSystem>
void 
BTK::IO::PDBFileParser<TypeSystem>::
rewind_file()
{
  _pdb_fs.clear();
  _pdb_fs.seekg(0,std::ios::beg);
  if (!_pdb_fs) {
    BTK_THROW_MSG(BTKFileReadFailure,
                  "Rewind seek failed in rewind_file()");
  }
  _cur_model = 0;
  _cur_line = "";
}

template <typename TypeSystem>
void 
BTK::IO::PDBFileParser<TypeSystem>::
save_header()
{
  // The strategy here is brain-dead, unfortunately.
  // Until someone volunteers to write a full-fledged parser for all
  // of the possible PDB header fields, we're taking an exclusionary
  // stance: everything that *isn't* interesting, and is before the start
  // of the first chain, is considered "header"
  
  rewind_file();
  
  if (!_pdb_fs) {
    BTK_THROW_MSG(BTKFileReadFailure,
                  "Bad stream in parse_header()");
  }
  
  _header.clear();
  
  while (getline(_pdb_fs,_cur_line)) {
    if (_cur_line.substr(0,4) == "ATOM" ||
        _cur_line.substr(0,6) == "HETATM" ||
        _cur_line.substr(0,5) == "MODEL") return;
    
    _header += _cur_line;
    _header += '\n';
  }
  
  if (!_pdb_fs) {
    BTK_THROW_MSG(BTKFileReadFailure,
                  "Reached EOF in parse_header()");
  }
}

template <typename TypeSystem>
bool
BTK::IO::PDBFileParser<TypeSystem>:: 
find_model(int which)
{
  // First, check the state of the stream, and
  // rewind if necessary to find the requested model.

  if (_pdb_fs.good()) {
    // file stream good/readable
    if (which == 0) {
      rewind_file();
      return find_next_atom(false);
    } else if (which <= _cur_model) {
      rewind_file();
    }
  } else if (_pdb_fs.eof()) {
    // at end of file
    rewind_file();
  } else {
    // bad stream!
    BTK_THROW_MSG(BTKFileReadFailure,
                  "Bad stream in find_model()");
  }

  // Once here, we're guaranteed to be upstream of the requested model,
  // and we need to seek forward to find it.
  while (_pdb_fs.good()) {
    // find the next MODEL entry
    do {
      if (_cur_line.substr(0,5) == "MODEL") break;
    } while (getline(_pdb_fs,_cur_line));

    if (_pdb_fs.eof()) {
      // We've reached the file end w/o finding the model.
      return false;
    } else if (_pdb_fs.fail()) {
      // bad stream -- read failure!
      BTK_THROW(BTKFileReadFailure);
    } else {
      // found a model -- get the model number
      int num = std::atoi(_cur_line.substr(10).c_str());

      if (which == 0 || num == which) {
        _cur_model = num;
        return true;
      }
    }
  }

  // if here, we definitely haven't found the model.
  return false;
}

template <typename TypeSystem>
bool 
BTK::IO::PDBFileParser<TypeSystem>::
find_next_atom(bool find_hetatoms)
{
  while (_cur_line.substr(0,4) != "ATOM" &&
         (!find_hetatoms || _cur_line.substr(0,6) != "HETATM"))
    {
      getline(_pdb_fs,_cur_line);
      if (!_pdb_fs) return false;
    }

  return true;
}

template <typename TypeSystem>
template <typename AtomType>
AtomType 
BTK::IO::PDBFileParser<TypeSystem>::
parse_atom(std::string const & atom_record,
           bool is_hetatom)
{
  using std::atoi;
  using std::atof;
  using std::string;

  int atom_num = atoi(atom_record.substr(6,5).c_str());
  int res_num = atoi(atom_record.substr(22,4).c_str());
  string atom_name = eat_whitespace(atom_record.substr(12,4));
  string res_name = eat_whitespace(atom_record.substr(17,3));
  char alt_loc = atom_record[16];
  char chain_id = atom_record[21];
  char insert_code = atom_record[26];
  BTK::MATH::BTKVector position(atof(atom_record.substr(30,8).c_str()),
                                atof(atom_record.substr(38,8).c_str()),
                                atof(atom_record.substr(46,8).c_str()));
  
  monomer_id_type res_type = 
    _dictionary.get_monomer_dictionary().get_type(res_name);

  if (res_type == BTK::UTILITY::TypeIDTraits<monomer_id_type>::unknown()) {
    std::ostringstream msg;
    msg << "Warning: unknown residue name " << res_name << " in PDB file.";
    BTK::IO::LOGGING::log_msg(msg.str(),
                              BTK::IO::LOGGING::IO & BTK::IO::LOGGING::WARNING);
  }

  atom_id_type atom_type = 
    _dictionary.get_atom_dictionary().get_type(atom_name);

  if (atom_type == BTK::UTILITY::TypeIDTraits<atom_id_type>::unknown()) {
    std::ostringstream msg;
    msg << "Warning: unknown atom name " << atom_name << " in PDB file.";
    BTK::IO::LOGGING::log_msg(msg.str(),
                              BTK::IO::LOGGING::IO & BTK::IO::LOGGING::WARNING);
  }

  // Fields after columb 46 are often missing or corrupted
  // (even though they are required), so we have to treat
  // them as special cases.
  double occupancy(0), b_factor(0);

  if (atom_record.size() >= 60) {
    occupancy = atof(atom_record.substr(54,6).c_str());
  }

  if (atom_record.size() >= 66) {
    b_factor = atof(atom_record.substr(60,6).c_str());
  }

  string segment_id;

  if (atom_record.size() >= 76) {
    segment_id = eat_whitespace(atom_record.substr(72,4));
  }

  string element_name;
  element_id_type element_type;

  if (atom_record.size() >= 78 &&
      (std::isalpha(atom_record[76]) || std::isalpha(atom_record[77]))) {
    // well, i'll be damned...a well-formed element field!
    element_name = eat_whitespace(atom_record.substr(76,2));
    element_type = _dictionary.get_element_dictionary().get_type(element_name);
  } else {
    // No element field found in ATOM record.
    // Attempt to deduce element type from atom name.

    // these are obvious element types, but are hard to
    // catch with the more general code below...
    if (atom_name == "ZN" || atom_name == "MG" || atom_name == "SE") {
      element_name = atom_name;
      element_type = _dictionary.get_element_dictionary().get_type(atom_name);
    } else if (is_hetatom && atom_name == "CA") {
      // calcium is a bit difficult, because the same atom name
      // is used for protein alpha carbons.  So we'll only
      // acknowledge it if it's a HETATM (the most likely use).
      element_name = atom_name;
      element_type = _dictionary.get_element_dictionary().get_type(atom_name);
    } else {
      // look for simple organic elements -- N, C, O, S, P, H
      unsigned i = atom_name.find_first_of("NCHOSP");
      element_name = atom_name.substr(i,1);
        
      if (i != std::string::npos) {
        element_type = _dictionary.get_element_dictionary().get_type(element_name);
      }
    }

    if (element_type == BTK::UTILITY::TypeIDTraits<element_id_type>::unknown()) {
      std::ostringstream msg;
      msg << "Warning: uknown element name " << element_name << " found in PDB file.";
      BTK::IO::LOGGING::log_msg(msg.str(),
                                BTK::IO::LOGGING::IO & 
                                BTK::IO::LOGGING::WARNING);
    }
  }

  string charge;

  if (atom_record.size() >= 80) {
    charge = eat_whitespace(atom_record.substr(78,2));
  }

  return AtomType(dictionary(),
                  is_hetatom,
                  atom_num,
                  atom_type,
                  alt_loc,
                  res_type,
                  chain_id,
                  res_num,
                  insert_code,
                  position,
                  occupancy,
                  b_factor,
                  segment_id,
                  element_type,
                  charge);
}


template <typename TypeSystem>
template <typename PolymerSystemType>
void
BTK::IO::PDBFileParser<TypeSystem>::
parse_current_model(PolymerSystemType & model_out,
                    bool load_hetatoms)
{
  typedef typename PolymerSystemType::atom_type AtomType;
  typedef typename PolymerSystemType::monomer_type MonomerType;
  typedef typename PolymerSystemType::chain_type PolymerType;
  typedef typename PolymerSystemType::chain_iterator ChainIterator;

  std::vector<AtomType> tmp_hetatoms, tmp_atoms;
  bool must_find_endmdl = false;
  bool found_endmdl = false;

  // set a cached logger object, so that duplicate parse errors
  // are filtered.
  {
    BTK::IO::LOGGING::Logger::logger_ptr
      tmp(new BTK::IO::LOGGING::
          CachedLogger(BTK::IO::LOGGING::
                       Logger::instance()));

    BTK::IO::LOGGING::Logger::push_logger(tmp);
  }

  // before parsing, clear the system
  model_out.clear();

  // see if the first line is actually a MODEL record.
  // if so, we know that we need to parse until we find an ENDMDL.
  if (_cur_line.substr(0,5) == "MODEL") {
    _cur_model = std::atoi(_cur_line.substr(10).c_str());
    must_find_endmdl = true;
  }

  // find the first ATOM in the model
  if (!find_next_atom(load_hetatoms)) {
    BTK_THROW_MSG(BTKFileReadFailure,
                  "Couldn't find first atom in PDB model");
  }

  // get first chain id, res number and res name
  bool found_ter = false;
  char prev_chain_id = _cur_line[21];
  std::string prev_res_name = eat_whitespace(_cur_line.substr(17,3));
  unsigned prev_res_num = std::atoi(_cur_line.substr(22,4).c_str());

  // add the first chain
  ChainIterator cur_chain = 
    model_out.insert(model_out.system_begin(),PolymerType());
  cur_chain->set_chain_id(_cur_line[21]);

  // process each relevant line in the model.
  // (right now, "relevant" lines are ATOMs, HETATMs, TERs, and ENDMDLs).
  do {
    if (_cur_line.substr(0,4) == "ATOM") {
      bool add_new_residue = false;
      bool add_new_chain = false;

      // get res name and nunber from atom record
      std::string res_name = eat_whitespace(_cur_line.substr(17,3));
      unsigned res_num = std::atoi(_cur_line.substr(22,4).c_str());
      
      if (_cur_line[21] == prev_chain_id) { // new atom in the current chain
        
        if (res_num != prev_res_num) // residue number change
          add_new_residue = true;
        else if (res_name != prev_res_name) {
          // res name change w/o res num change -- what the heck is going on?
          std::ostringstream msg;
          msg << "PDB res name change w/o res num change in " 
              << _cur_line.substr(0,11) << " (cur name=" 
              << res_name << ", prev name=" << prev_res_name << ")";
          BTK_THROW_MSG(BTKCorruptPDBFile,msg.str().c_str());
        }

      } else {  // chain ID changed: new atom in a new chain
        
        if (!found_ter) { // chain ID changed w/o a TER record
          std::ostringstream msg;
          msg << "Chain id changed from "
              << prev_chain_id << " to " << _cur_line[21]
              << " w/o a TER record.\n";
          using namespace BTK::IO::LOGGING;
          Logger::instance().log(msg.str(),WARNING & IO);
        }
        
        add_new_residue = true;
        add_new_chain = true;
      }        

      // handle addition of new residue to current chain
      if (add_new_residue && tmp_atoms.size()) {
        cur_chain->insert(cur_chain->polymer_end(),
                          MonomerType(tmp_atoms.begin(),
                                      tmp_atoms.end(),
                                      tmp_atoms.begin()->res_type(),
                                      tmp_atoms.begin()->res_number()));
          
        tmp_atoms.clear();
        prev_res_name = res_name;
        prev_res_num = res_num;
      }

      // handle addition of new chain
      if (add_new_chain) {
        cur_chain = model_out.insert(model_out.system_end(),PolymerType());
        cur_chain->set_chain_id(_cur_line[21]);
        prev_chain_id = _cur_line[21];
        found_ter = false;
      }

      //  Parse current atom
      tmp_atoms.push_back(parse_atom<AtomType>(_cur_line));

    } else if (_cur_line.substr(0,6) == "HETATM") {
      // parse heterogen atom & save for later
      if (!load_hetatoms) continue;

      tmp_hetatoms.push_back(parse_atom<AtomType>(_cur_line,true));

    } else if (_cur_line.substr(0,3) == "TER") {
      // normal chain termination.
      found_ter = true;

    } else if (_cur_line.substr(0,6) == "ENDMDL") {
      // normal model termination
      found_endmdl = true;
      break;
    } else if (_cur_line.substr(0,3) == "END") {
      // nothing for now
    } else if (_cur_line.substr(0,6) == "CONECT") {
      // nothing for now
    } else if (_cur_line.size() > 0) {
      std::ostringstream msg; 
      msg << "Unknown record type in pdb file : "
          << _cur_line.substr(0,6) << '\n';

      using namespace BTK::IO::LOGGING;
      Logger::instance().log(msg.str(),WARNING & IO);
    }
  } while (getline(_pdb_fs,_cur_line));

  // Add final monomer, if any.
  if (tmp_atoms.size()) {
    cur_chain->insert(cur_chain->polymer_end(),
                      MonomerType(tmp_atoms.begin(),
                                  tmp_atoms.end(),
                                  tmp_atoms.begin()->res_type(),
                                  tmp_atoms.begin()->res_number()));
    tmp_atoms.clear();
  }

  if (must_find_endmdl && !found_endmdl) {
    using namespace BTK::IO::LOGGING;
    Logger::instance().log("Reached end of PDB before finding ENDMDL record",
                           WARNING & IO);
  }

  if (load_hetatoms) save_hetatoms(tmp_hetatoms.begin(),
                                   tmp_hetatoms.end(),
                                   model_out);

  // remove the cached logger
  BTK::IO::LOGGING::Logger::pop_logger();
}

struct hetatom_compare {
  template <typename AtomType>
  bool operator()(AtomType const & a1,
                  AtomType const & a2) const
  {
    if (a1.chain_id() < a2.chain_id()) return true;
    else if (a1.chain_id() == a2.chain_id()) {
      if (a1.res_number() < a2.res_number()) return true;
      else return false;
    } else return false;
  }
};

template <typename TypeSystem>
template <typename AtomIterator, typename PolymerSystemType>
void
BTK::IO::PDBFileParser<TypeSystem>:: 
save_hetatoms(AtomIterator het_begin, AtomIterator het_end,
              PolymerSystemType & model_out)
{
  typedef typename PolymerSystemType::chain_type PolymerType;
  typedef typename PolymerSystemType::monomer_type MonomerType;
  typedef typename PolymerSystemType::chain_iterator ChainIterator;

  if (het_begin == het_end) return;

  // sort hetatoms into chain/residue groups.
  std::stable_sort(het_begin,
                   het_end,
                   hetatom_compare());

  AtomIterator i, j, end = het_end;
  int prev_res_num;
  char prev_chain_id;

  // add first hetatom chain
  ChainIterator cur_chain = model_out.insert(model_out.system_end(),
                                             PolymerType());
  cur_chain->set_chain_id(het_begin->chain_id());

  prev_chain_id = het_begin->chain_id();
  prev_res_num = het_begin->res_number();

  // add hetatoms, adding monomers and chains as necessary.
  for (i = het_begin, j = het_begin; j != end; ++j) {

    // create a new monomer, if necessary
    if (j->res_number() != prev_res_num) {
      cur_chain->insert(cur_chain->polymer_end(),
                        MonomerType(i,j,i->res_type(),
                                    i->res_number()));

      prev_res_num = j->res_number();
      i = j;
    }

    // create a new chain, if necessary
    if (j->chain_id() != prev_chain_id) {
      cur_chain = model_out.insert(model_out.system_end(),PolymerType());
      cur_chain->set_chain_id(j->chain_id());

      prev_chain_id = j->chain_id();
      i = j; // just in case -- probably not necessary
    }
  }

  // add the last hetatom monomer, if any.
  if (i != j) cur_chain->insert(cur_chain->polymer_end(),
                                MonomerType(i,j,i->res_type(),
                                            i->res_number()));
}

#endif // BTK_PDB_FILE_PARSER_HPP

---