/* * Project: MoleCuilder * Description: creates and alters molecular systems * Copyright (C) 2010 University of Bonn. All rights reserved. * Please see the LICENSE file or "Copyright notice" in builder.cpp for details. */ /* * PdbParser.cpp * * Created on: Aug 17, 2010 * Author: heber */ // include config.h #ifdef HAVE_CONFIG_H #include #endif #include "Helpers/MemDebug.hpp" #include "Helpers/Assert.hpp" #include "Helpers/Log.hpp" #include "Helpers/toString.hpp" #include "Helpers/Verbose.hpp" #include "World.hpp" #include "atom.hpp" #include "bond.hpp" #include "element.hpp" #include "molecule.hpp" #include "periodentafel.hpp" #include "Descriptors/AtomIdDescriptor.hpp" #include "Parser/PdbParser.hpp" #include #include #include #include using namespace std; /** * Constructor. */ PdbParser::PdbParser() { knownTokens["ATOM"] = PdbKey::Atom; knownTokens["TER"] = PdbKey::Filler; knownTokens["END"] = PdbKey::EndOfFile; knownTokens["CONECT"] = PdbKey::Connect; knownTokens["REMARK"] = PdbKey::Remark; knownTokens[""] = PdbKey::EndOfFile; } /** * Destructor. */ PdbParser::~PdbParser() { additionalAtomData.clear(); atomIdMap.clear(); } /** Parses the initial word of the given \a line and returns the token type. * * @param line line to scan * @return token type */ enum PdbKey::KnownTokens PdbParser::getToken(string &line) { // look for first space const size_t space_location = line.find(' '); const size_t tab_location = line.find('\t'); size_t location = space_location < tab_location ? space_location : tab_location; string token; if (location != string::npos) { //DoLog(1) && (Log() << Verbose(1) << "Found space at position " << space_location << std::endl); token = line.substr(0,space_location); } else { token = line; } //DoLog(1) && (Log() << Verbose(1) << "Token is " << token << std::endl); if (knownTokens.count(token) == 0) return PdbKey::NoToken; else return knownTokens[token]; return PdbKey::NoToken; } /** * Loads atoms from a PDB-formatted file. * * \param PDB file */ void PdbParser::load(istream* file) { string line; size_t linecount = 0; enum PdbKey::KnownTokens token; molecule *newmol = World::getInstance().createMolecule(); newmol->ActiveFlag = true; bool NotEndOfFile = true; // TODO: Remove the insertion into molecule when saving does not depend on them anymore. Also, remove molecule.hpp include World::getInstance().getMolecules()->insert(newmol); while (NotEndOfFile) { std::getline(*file, line, '\n'); // extract first token token = getToken(line); //DoLog(1) && (Log() << Verbose(1) << " Recognized token of type : " << token << std::endl); switch (token) { case PdbKey::Atom: readAtomDataLine(line, newmol); break; case PdbKey::Remark: break; case PdbKey::Connect: readNeighbors(line); break; case PdbKey::Filler: break; case PdbKey::EndOfFile: NotEndOfFile = false; break; default: // TODO: put a throw here DoeLog(2) && (eLog() << Verbose(2) << "Unknown token: '" << line << "'" << std::endl); //ASSERT(0, "PdbParser::load() - Unknown token in line "+toString(linecount)+": "+line+"."); break; } NotEndOfFile = NotEndOfFile && (file->good()); linecount++; } } /** * Saves the World's current state into as a PDB file. * * \param file where to save the state */ void PdbParser::save(ostream* file) { DoLog(0) && (Log() << Verbose(0) << "Saving changes to pdb." << std::endl); { // add initial remark *file << "REMARK created by molecuilder on "; time_t now = time((time_t *)NULL); // Get the system time and put it into 'now' as 'calender time' // ctime ends in \n\0, we have to cut away the newline std::string time(ctime(&now)); size_t pos = time.find('\n'); if (pos != 0) *file << time.substr(0,pos); else *file << time; *file << endl; } // we distribute new atom numbers, hence clear map beforehand atomIdMap.clear(); { vector AtomList = World::getInstance().getAllAtoms(); std::vector elementNo(MAX_ELEMENTS,1); char name[MAXSTRINGSIZE]; // write ATOMs int AtomNo = 1; // serial number starts at 1 in pdb int MolNo = 1; // residue number starts at 1 in pdb for (vector::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) { const size_t Z = (*atomIt)->getType()->getAtomicNumber(); sprintf(name, "%2s%02d",(*atomIt)->getType()->getSymbol().c_str(), elementNo[Z]); elementNo[Z] = (elementNo[Z]+1) % 100; // confine to two digits const molecule *mol = (*atomIt)->getMolecule(); if (mol == NULL) { // for homeless atoms, MolNo = -1 is reserved MolNo = -1; } else { MolNo = mol->getId(); } saveLine(file, *atomIt, name, AtomNo, MolNo); setAtomId((*atomIt)->getId(), AtomNo); AtomNo++; } // write CONECTs for (vector::iterator atomIt = AtomList.begin(); atomIt != AtomList.end(); atomIt++) { writeNeighbors(file, 4, *atomIt); } } // END *file << "END" << endl; } /** * Writes one line of PDB-formatted data to the provided stream. * * \param stream where to write the line to * \param *currentAtom the atom of which information should be written * \param *name name of atom, i.e. H01 * \param AtomNo serial number of atom * \param ResidueNo number of residue */ void PdbParser::saveLine(ostream* file, const atom* currentAtom, const char *name, const int AtomNo, const int ResidueNo) { *file << "ATOM "; *file << setw(6) << AtomNo; /* atom serial number */ *file << setw(1) << " "; *file << setfill(' ') << left << setw(4) << name << right; /* atom name */ *file << setw(1) << " "; *file << setfill(' ') << setw(3) << ((currentAtom->getMolecule() != NULL) ? currentAtom->getMolecule()->getName().substr(0,3) : "-"); /* residue name */ *file << setw(1) << " "; *file << setfill(' ') << setw(1) << (char)('a'+(unsigned char)(AtomNo % 26)); /* letter for chain */ *file << setw(4) << ResidueNo; /* residue sequence number */ *file << setw(4) << " "; for (int i=0;iat(i); /* positional coordinate in Angstroem */ } *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getValence(); /* occupancy */ *file << setw(6) << setprecision(2) << showpoint << (double)currentAtom->getType()->getNoValenceOrbitals(); /* temperature factor */ *file << noshowpoint; *file << setw(6) << " "; *file << setw(4) << "0"; *file << setfill(' ') << setw(2) << currentAtom->getType()->getSymbol(); *file << setw(2) << "0"; *file << endl; } /** * Writes the neighbor information of one atom to the provided stream. * * \param *file where to write neighbor information to * \param MaxnumberOfNeighbors of neighbors * \param *currentAtom to the atom of which to take the neighbor information */ void PdbParser::writeNeighbors(ostream* file, int MaxnumberOfNeighbors, atom* currentAtom) { if (!currentAtom->ListOfBonds.empty()) { *file << "CONECT"; *file << setw(5) << getSerial(currentAtom->getId()); int MaxNo = 0; for(BondList::iterator currentBond = currentAtom->ListOfBonds.begin(); currentBond != currentAtom->ListOfBonds.end(); ++currentBond) { if (MaxNo < MaxnumberOfNeighbors) { *file << setw(5) << getSerial((*currentBond)->GetOtherAtom(currentAtom)->getId()); } MaxNo++; } *file << endl; } } /** Retrieves a value from PdbParser::atomIdMap. * \param atomid key * \return value */ size_t PdbParser::getSerial(const size_t atomid) const { ConvertTo toSize_t; ASSERT(additionalAtomData.find(atomid) != additionalAtomData.end(), "PdbParser::getSerial: atomid "+toString(atomid)+" not present in Map."); const PdbAtomInfoContainer &atomInfo = additionalAtomData.at(atomid); return toSize_t(atomInfo.get(PdbKey::serial)); } /** Retrieves a value from PdbParser::atomIdMap. * \param atomid key * \return value */ size_t PdbParser::getAtomId(const size_t atomid) const { ASSERT(atomIdMap.find(atomid) != atomIdMap.end(), "PdbParser::getAtomId: atomid not present in Map."); return (atomIdMap.find(atomid)->second); } /** Sets an entry in PdbParser::atomIdMap. * \param localatomid key * \param atomid value * \return true - key not present, false - value present */ void PdbParser::setAtomId(const size_t localatomid, const size_t atomid) { pair::iterator, bool > inserter; inserter = atomIdMap.insert( make_pair(localatomid, atomid) ); ASSERT(inserter.second, "PdbParser::setAtomId: atomId already present in Map."); } /** Parse an ATOM line from a PDB file. * * Reads one data line of a pdstatus file and interprets it according to the * specifications of the PDB 3.2 format: http://www.wwpdb.org/docs.html * * A new atom is created and filled with available information, non- * standard information is placed in additionalAtomData at the atom's id. * * \param line to parse as an atom * \param newmol molecule to add parsed atoms to */ void PdbParser::readAtomDataLine(std::string &line, molecule *newmol = NULL) { vector::iterator it; stringstream lineStream; atom* newAtom = World::getInstance().createAtom(); additionalAtomData[newAtom->getId()] = *(new PdbAtomInfoContainer); PdbAtomInfoContainer &atomInfo = additionalAtomData[newAtom->getId()]; string word; ConvertTo toSize_t; double tmp; lineStream << line; atomInfo.set(PdbKey::serial, line.substr(6,5)); std::pair< std::set::const_iterator, bool> Inserter = SerialSet.insert(toSize_t(atomInfo.get(PdbKey::serial))); ASSERT(Inserter.second, "PdbParser::readAtomDataLine() - ATOM contains entry with serial " +atomInfo.get(PdbKey::serial)+" already present!"); // assign hightest+1 instead, but then beware of CONECT entries! Another map needed! // if (!Inserter.second) { // const size_t id = (*SerialSet.rbegin())+1; // SerialSet.insert(id); // atomInfo.set(PdbKey::serial, toString(id)); // DoeLog(2) && (eLog() << Verbose(2) // << "Serial " << atomInfo.get(PdbKey::serial) << " already present, " // << "assigning " << toString(id) << " instead." << std::endl); // } // check whether serial exists, if so, assign next available // DoLog(2) && (Log() << Verbose(2) << "Split line:" // << line.substr(6,5) << "|" // << line.substr(12,4) << "|" // << line.substr(16,1) << "|" // << line.substr(17,3) << "|" // << line.substr(21,1) << "|" // << line.substr(22,4) << "|" // << line.substr(26,1) << "|" // << line.substr(30,8) << "|" // << line.substr(38,8) << "|" // << line.substr(46,8) << "|" // << line.substr(54,6) << "|" // << line.substr(60,6) << "|" // << line.substr(76,2) << "|" // << line.substr(78,2) << std::endl); setAtomId(toSize_t(atomInfo.get(PdbKey::serial)), newAtom->getId()); atomInfo.set(PdbKey::name, line.substr(12,4)); atomInfo.set(PdbKey::altloc, line.substr(16,1)); atomInfo.set(PdbKey::resName, line.substr(17,3)); atomInfo.set(PdbKey::chainID, line.substr(21,1)); atomInfo.set(PdbKey::resSeq, line.substr(22,4)); atomInfo.set(PdbKey::iCode, line.substr(26,1)); PdbAtomInfoContainer::ScanKey(tmp, line.substr(30,8)); newAtom->set(0, tmp); PdbAtomInfoContainer::ScanKey(tmp, line.substr(38,8)); newAtom->set(1, tmp); PdbAtomInfoContainer::ScanKey(tmp, line.substr(46,8)); newAtom->set(2, tmp); atomInfo.set(PdbKey::occupancy, line.substr(54,6)); atomInfo.set(PdbKey::tempFactor, line.substr(60,6)); atomInfo.set(PdbKey::charge, line.substr(78,2)); PdbAtomInfoContainer::ScanKey(word, line.substr(76,2)); newAtom->setType(World::getInstance().getPeriode()->FindElement(word)); if (newmol != NULL) newmol->AddAtom(newAtom); // printAtomInfo(newAtom); } /** Prints all PDB-specific information known about an atom. * */ void PdbParser::printAtomInfo(const atom * const newAtom) const { const PdbAtomInfoContainer &atomInfo = additionalAtomData.at(newAtom->getId()); // operator[] const does not exist DoLog(1) && (Log() << Verbose(1) << "We know about atom " << newAtom->getId() << ":" << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tserial is " << atomInfo.get(PdbKey::serial) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tname is " << atomInfo.get(PdbKey::name) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\taltloc is " << atomInfo.get(PdbKey::altloc) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tresName is " << atomInfo.get(PdbKey::resName) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tchainID is " << atomInfo.get(PdbKey::chainID) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tresSeq is " << atomInfo.get(PdbKey::resSeq) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tiCode is " << atomInfo.get(PdbKey::iCode) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tx is " << newAtom->getPosition() << std::endl); DoLog(1) && (Log() << Verbose(1) << "\toccupancy is " << atomInfo.get(PdbKey::occupancy) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\ttempFactor is " << atomInfo.get(PdbKey::tempFactor) << std::endl); DoLog(1) && (Log() << Verbose(1) << "\telement is '" << *(newAtom->getType()) << "'" << std::endl); DoLog(1) && (Log() << Verbose(1) << "\tcharge is " << atomInfo.get(PdbKey::charge) << std::endl); } /** * Reads neighbor information for one atom from the input. * * \param line to parse as an atom */ void PdbParser::readNeighbors(std::string &line) { const size_t length = line.length(); std::list ListOfNeighbors; ConvertTo toSize_t; // obtain neighbours // show split line for debugging string output; ASSERT(length >=16, "PdbParser::readNeighbors() - CONECT entry has not enough entries: "+line+"!"); // output = "Split line:|"; // output += line.substr(6,5) + "|"; const size_t id = toSize_t(line.substr(6,5)); for (size_t index = 11; index <= 26; index+=5) { if (index+5 <= length) { // output += line.substr(index,5) + "|"; const size_t otherid = toSize_t(line.substr(index,5)); ListOfNeighbors.push_back(otherid); } else { break; } } // DoLog(2) && (Log() << Verbose(2) << output << std::endl); // add neighbours atom *_atom = World::getInstance().getAtom(AtomById(getAtomId(id))); for (std::list::const_iterator iter = ListOfNeighbors.begin(); iter != ListOfNeighbors.end(); ++iter) { // DoLog(1) && (Log() << Verbose(1) << "Adding Bond (" << getAtomId(id) << "," << getAtomId(*iter) << ")" << std::endl); atom * const _Otheratom = World::getInstance().getAtom(AtomById(getAtomId(*iter))); _atom->addBond(_Otheratom); } } /** * Replaces atom IDs read from the file by the corresponding world IDs. All IDs * IDs of the input string will be replaced; expected separating characters are * "-" and ",". * * \param string in which atom IDs should be adapted * * \return input string with modified atom IDs */ //string PdbParser::adaptIdDependentDataString(string data) { // // there might be no IDs // if (data == "-") { // return "-"; // } // // char separator; // int id; // stringstream line, result; // line << data; // // line >> id; // result << atomIdMap[id]; // while (line.good()) { // line >> separator >> id; // result << separator << atomIdMap[id]; // } // // return result.str(); // return ""; //} bool PdbParser::operator==(const PdbParser& b) const { bool status = true; World::AtomComposite atoms = World::getInstance().getAllAtoms(); for (World::AtomComposite::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) { if ((additionalAtomData.find((*iter)->getId()) != additionalAtomData.end()) && (b.additionalAtomData.find((*iter)->getId()) != b.additionalAtomData.end())) { const PdbAtomInfoContainer &atomInfo = additionalAtomData.at((*iter)->getId()); const PdbAtomInfoContainer &OtheratomInfo = b.additionalAtomData.at((*iter)->getId()); status = status && (atomInfo.get(PdbKey::serial) == OtheratomInfo.get(PdbKey::serial)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in serials!" << std::endl); status = status && (atomInfo.get(PdbKey::name) == OtheratomInfo.get(PdbKey::name)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in names!" << std::endl); status = status && (atomInfo.get(PdbKey::altloc) == OtheratomInfo.get(PdbKey::altloc)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in altlocs!" << std::endl); status = status && (atomInfo.get(PdbKey::resName) == OtheratomInfo.get(PdbKey::resName)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resNames!" << std::endl); status = status && (atomInfo.get(PdbKey::chainID) == OtheratomInfo.get(PdbKey::chainID)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in chainIDs!" << std::endl); status = status && (atomInfo.get(PdbKey::resSeq) == OtheratomInfo.get(PdbKey::resSeq)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in resSeqs!" << std::endl); status = status && (atomInfo.get(PdbKey::iCode) == OtheratomInfo.get(PdbKey::iCode)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in iCodes!" << std::endl); status = status && (atomInfo.get(PdbKey::occupancy) == OtheratomInfo.get(PdbKey::occupancy)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in occupancies!" << std::endl); status = status && (atomInfo.get(PdbKey::tempFactor) == OtheratomInfo.get(PdbKey::tempFactor)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in tempFactors!" << std::endl); status = status && (atomInfo.get(PdbKey::charge) == OtheratomInfo.get(PdbKey::charge)); if (!status) DoeLog(1) && (eLog() << Verbose(1) << "Mismatch in charges!" << std::endl); } } return status; }