Commit 6c935996 by maarten

renaming intermediate backup

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@170 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent ca881b82
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
#pragma once #pragma once
#include "libcif/config.h" #include "cif++/Config.h"
#include <boost/filesystem/operations.hpp> #include <boost/filesystem/operations.hpp>
#include <boost/math/quaternion.hpp> #include <boost/math/quaternion.hpp>
...@@ -10,7 +10,7 @@ ...@@ -10,7 +10,7 @@
namespace libcif namespace libcif
{ {
enum atom_type : uint8 enum AtomType : uint8
{ {
Nn = 0, // Unknown Nn = 0, // Unknown
...@@ -143,9 +143,9 @@ enum atom_type : uint8 ...@@ -143,9 +143,9 @@ enum atom_type : uint8
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// atom_type_info // AtomTypeInfo
enum radius_type { enum RadiusType {
eRadiusCalculated, eRadiusCalculated,
eRadiusEmpirical, eRadiusEmpirical,
eRadiusCovalentEmpirical, eRadiusCovalentEmpirical,
...@@ -159,9 +159,9 @@ enum radius_type { ...@@ -159,9 +159,9 @@ enum radius_type {
eRadiusTypeCount eRadiusTypeCount
}; };
struct atom_type_info struct AtomTypeInfo
{ {
atom_type type; AtomType type;
std::string name; std::string name;
std::string symbol; std::string symbol;
float weight; float weight;
...@@ -169,36 +169,36 @@ struct atom_type_info ...@@ -169,36 +169,36 @@ struct atom_type_info
float radii[eRadiusTypeCount]; float radii[eRadiusTypeCount];
}; };
extern const atom_type_info kKnownAtoms[]; extern const AtomTypeInfo kKnownAtoms[];
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// atom_type_traits // AtomTypeTraits
class atom_type_traits class AtomTypeTraits
{ {
public: public:
atom_type_traits(atom_type a); AtomTypeTraits(AtomType a);
atom_type_traits(const std::string& symbol); AtomTypeTraits(const std::string& symbol);
atom_type type() const { return m_info->type; } AtomType type() const { return mInfo->type; }
std::string name() const { return m_info->name; } std::string name() const { return mInfo->name; }
std::string symbol() const { return m_info->symbol; } std::string symbol() const { return mInfo->symbol; }
float weight() const { return m_info->weight; } float weight() const { return mInfo->weight; }
bool is_metal() const { return m_info->metal; } bool isMetal() const { return mInfo->metal; }
static bool is_element(const std::string& symbol); static bool isElement(const std::string& symbol);
static bool is_metal(const std::string& symbol); static bool isMetal(const std::string& symbol);
float radius(radius_type type = eRadiusSingleBond) const float radius(RadiusType type = eRadiusSingleBond) const
{ {
if (type >= eRadiusTypeCount) if (type >= eRadiusTypeCount)
throw std::invalid_argument("invalid radius requested"); throw std::invalid_argument("invalid radius requested");
return m_info->radii[type] / 100.f; return mInfo->radii[type] / 100.f;
} }
private: private:
const struct atom_type_info* m_info; const struct AtomTypeInfo* mInfo;
}; };
} }
#pragma once
#include "cif++/Cif++.h"
void WritePDBFile(std::ostream& pdbFile, cif::File& cifFile);
// CIF parser // CIF Parser
#include "libcif/cif++.h" #include "cif++/Cif++.h"
#include <stack> #include <stack>
...@@ -9,10 +9,10 @@ namespace cif ...@@ -9,10 +9,10 @@ namespace cif
// -------------------------------------------------------------------- // --------------------------------------------------------------------
class cif_parser_error : public std::runtime_error class CifParserError : public std::runtime_error
{ {
public: public:
cif_parser_error(uint32 line_nr, const std::string& message); CifParserError(uint32 lineNr, const std::string& message);
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -28,38 +28,38 @@ enum CharTraitsMask: uint8 { ...@@ -28,38 +28,38 @@ enum CharTraitsMask: uint8 {
kAnyPrintMask = 1 << 3 kAnyPrintMask = 1 << 3
}; };
inline bool is_white(int ch) inline bool isWhite(int ch)
{ {
return std::isspace(ch) or ch == '#'; return std::isspace(ch) or ch == '#';
} }
inline bool is_ordinary(int ch) inline bool isOrdinary(int ch)
{ {
return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kOrdinaryMask) != 0; return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kOrdinaryMask) != 0;
} }
inline bool is_non_blank(int ch) inline bool isNonBlank(int ch)
{ {
return ch > 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kNonBlankMask) != 0; return ch > 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kNonBlankMask) != 0;
} }
inline bool is_text_lead(int ch) inline bool isTextLead(int ch)
{ {
return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kTextLeadMask) != 0; return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kTextLeadMask) != 0;
} }
inline bool is_any_print(int ch) inline bool isAnyPrint(int ch)
{ {
return ch == '\t' or return ch == '\t' or
(ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0); (ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
} }
inline bool is_unquoted_string(const char* s) inline bool isUnquotedString(const char* s)
{ {
bool result = is_ordinary(*s++); bool result = isOrdinary(*s++);
while (result and *s != 0) while (result and *s != 0)
{ {
result = is_non_blank(*s); result = isNonBlank(*s);
++s; ++s;
} }
return result; return result;
...@@ -67,16 +67,16 @@ inline bool is_unquoted_string(const char* s) ...@@ -67,16 +67,16 @@ inline bool is_unquoted_string(const char* s)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
std::tuple<std::string,std::string> split_tag_name(const std::string& tag); std::tuple<std::string,std::string> splitTagName(const std::string& tag);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// sac parser, analogous to SAX parser (simple api for xml) // sac Parser, analogous to SAX Parser (simple api for xml)
class sac_parser class SacParser
{ {
public: public:
sac_parser(std::istream& is); SacParser(std::istream& is);
virtual ~sac_parser() {} virtual ~SacParser() {}
enum CIFToken enum CIFToken
{ {
...@@ -108,30 +108,30 @@ class sac_parser ...@@ -108,30 +108,30 @@ class sac_parser
static const char* kValueName[]; static const char* kValueName[];
int get_next_char(); int getNextChar();
void retract(); void retract();
void restart(); void restart();
CIFToken get_next_token(); CIFToken getNextToken();
void match(CIFToken token); void match(CIFToken token);
void parse_file(); void parseFile();
void parse_global(); void parseGlobal();
void parse_data_block(); void parseDataBlock();
virtual void parse_save_frame(); virtual void parseSaveFrame();
void parse_dictionary(); void parseDictionary();
void error(const std::string& msg); void error(const std::string& msg);
// production methods, these are pure virtual here // production methods, these are pure virtual here
virtual void produce_datablock(const std::string& name) = 0; virtual void produceDatablock(const std::string& name) = 0;
virtual void produce_category(const std::string& name) = 0; virtual void produceCategory(const std::string& name) = 0;
virtual void produce_row() = 0; virtual void produceRow() = 0;
virtual void produce_item(const std::string& category, const std::string& item, const string& value) = 0; virtual void produceItem(const std::string& category, const std::string& item, const string& value) = 0;
protected: protected:
...@@ -153,60 +153,60 @@ class sac_parser ...@@ -153,60 +153,60 @@ class sac_parser
eStateValue = 300 eStateValue = 300
}; };
std::istream& m_data; std::istream& mData;
// parser state // Parser state
bool m_validate; bool mValidate;
uint32 m_line_nr; uint32 mLineNr;
bool m_bol; bool mBol;
int m_state, m_start; int mState, mStart;
CIFToken m_lookahead; CIFToken mLookahead;
std::string m_token_value; std::string mTokenValue;
CIFValueType m_token_type; CIFValueType mTokenType;
std::stack<int> m_buffer; std::stack<int> mBuffer;
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
class parser : public sac_parser class Parser : public SacParser
{ {
public: public:
parser(std::istream& is, file& f); Parser(std::istream& is, File& f);
virtual void produce_datablock(const std::string& name); virtual void produceDatablock(const std::string& name);
virtual void produce_category(const std::string& name); virtual void produceCategory(const std::string& name);
virtual void produce_row(); virtual void produceRow();
virtual void produce_item(const std::string& category, const std::string& item, const std::string& value); virtual void produceItem(const std::string& category, const std::string& item, const std::string& value);
protected: protected:
file& m_file; File& mFile;
datablock* m_data_block; Datablock* mDataBlock;
datablock::iterator m_cat; Datablock::iterator mCat;
row m_row; Row mRow;
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
class dict_parser : public parser class DictParser : public Parser
{ {
public: public:
dict_parser(validator& validator, std::istream& is); DictParser(Validator& validator, std::istream& is);
~dict_parser(); ~DictParser();
void load_dictionary(); void loadDictionary();
private: private:
virtual void parse_save_frame(); virtual void parseSaveFrame();
bool collect_item_types(); bool collectItemTypes();
void link_items(); void linkItems();
validator& m_validator; Validator& mValidator;
file m_file; File mFile;
struct dict_parser_data_impl* m_impl; struct DictParserDataImpl* mImpl;
bool m_collected_item_types = false; bool mCollectedItemTypes = false;
}; };
} }
...@@ -5,7 +5,7 @@ ...@@ -5,7 +5,7 @@
#include <vector> #include <vector>
#include <set> #include <set>
#include "libcif/config.h" #include "cif++/Config.h"
namespace cif namespace cif
{ {
...@@ -19,8 +19,8 @@ int icompare(const std::string& a, const std::string& b); ...@@ -19,8 +19,8 @@ int icompare(const std::string& a, const std::string& b);
bool iequals(const char* a, const char* b); bool iequals(const char* a, const char* b);
int icompare(const char* a, const char* b); int icompare(const char* a, const char* b);
void to_lower(std::string& s); void toLower(std::string& s);
std::string to_lower_copy(const std::string& s); std::string toLowerCopy(const std::string& s);
// To make life easier, we also define iless and iset using iequals // To make life easier, we also define iless and iset using iequals
...@@ -46,11 +46,11 @@ inline char tolower(char ch) ...@@ -46,11 +46,11 @@ inline char tolower(char ch)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
std::tuple<std::string,std::string> split_tag_name(const std::string& tag); std::tuple<std::string,std::string> splitTagName(const std::string& tag);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// custom wordwrapping routine // custom wordwrapping routine
std::vector<std::string> word_wrap(const std::string& text, unsigned int width); std::vector<std::string> wordWrap(const std::string& text, unsigned int width);
} }
// cif parsing library
#include "cif++/Cif++.h"
#include <boost/filesystem/path.hpp>
// the std regex of gcc is crashing....
#include <boost/regex.hpp>
#include <set>
namespace cif
{
struct ValidateCategory;
// --------------------------------------------------------------------
class ValidationError : public std::exception
{
public:
ValidationError(const std::string& msg) : mMsg(msg) {}
const char* what() const noexcept { return mMsg.c_str(); }
std::string mMsg;
};
// --------------------------------------------------------------------
enum DDL_PrimitiveType
{
ptChar, ptUChar, ptNumb
};
DDL_PrimitiveType mapToPrimitiveType(const std::string& s);
struct ValidateType
{
std::string mName;
DDL_PrimitiveType mPrimitiveType;
boost::regex mRx;
bool operator<(const ValidateType& rhs) const
{
return icompare(mName, rhs.mName) < 0;
}
// compare values based on type
// int compare(const std::string& a, const std::string& b) const
// {
// return compare(a.c_str(), b.c_str());
// }
int compare(const char* a, const char* b) const;
};
struct ValidateItem
{
std::string mTag;
bool mMandatory;
const ValidateType* mType;
cif::iset mEnums;
ValidateItem* mParent = nullptr;
std::set<ValidateItem*>
mChildren;
ValidateCategory* mCategory = nullptr;
std::set<ValidateItem*>
mForeignKeys;
void setParent(ValidateItem* parent);
bool operator<(const ValidateItem& rhs) const
{
return icompare(mTag, rhs.mTag) < 0;
}
bool operator==(const ValidateItem& rhs) const
{
return iequals(mTag, rhs.mTag);
}
void operator()(std::string value) const;
};
struct ValidateCategory
{
std::string mName;
std::vector<string> mKeys;
cif::iset mGroups;
cif::iset mMandatoryFields;
std::set<ValidateItem> mItemValidators;
bool operator<(const ValidateCategory& rhs) const
{
return icompare(mName, rhs.mName) < 0;
}
void addItemValidator(ValidateItem&& v);
const ValidateItem* getValidatorForItem(std::string tag) const;
const std::set<ValidateItem>& itemValidators() const
{
return mItemValidators;
}
};
// --------------------------------------------------------------------
class Validator
{
public:
friend class DictParser;
Validator();
~Validator();
Validator(const Validator& rhs) = delete;
Validator& operator=(const Validator& rhs) = delete;
Validator(Validator&& rhs);
Validator& operator=(Validator&& rhs);
void addTypeValidator(ValidateType&& v);
const ValidateType* getValidatorForType(std::string typeCode) const;
void addCategoryValidator(ValidateCategory&& v);
const ValidateCategory* getValidatorForCategory(std::string category) const;
void reportError(const std::string& msg);
std::string dictName() const { return mName; }
void dictName(const std::string& name) { mName = name; }
std::string dictVersion() const { return mVersion; }
void dictVersion(const std::string& version) { mVersion = version; }
private:
// name is fully qualified here:
ValidateItem* getValidatorForItem(std::string name) const;
std::string mName;
std::string mVersion;
bool mStrict = false;
// std::set<uint32> mSubCategories;
std::set<ValidateType> mTypeValidators;
std::set<ValidateCategory> mCategoryValidators;
};
}
// Lib for working with structures as contained in mmCIF and PDB files
#pragma once
#include <set>
#include <tuple>
#include <vector>
#include <map>
#include "libcif++/AtomType.h"
namespace libcif
{
// --------------------------------------------------------------------
// The chemical composition of the structure in an mmCIF file is
// defined in the class composition. A compositon consists of
// entities. Each Entity can be either a polymer, a non-polymer
// a macrolide or a water molecule.
// Entities themselves are made up of compounds. And compounds
// contain CompoundAtom records for each atom.
class Composition;
class Entity;
class Compound;
struct CompoundAtom;
// --------------------------------------------------------------------
// struct containing information about an atom in a chemical compound
// This information comes from the CCP4 monomer library.
struct CompoundAtom
{
std::string id;
AtomType typeSymbol;
std::string typeEnergy;
float partialCharge;
};
// --------------------------------------------------------------------
// a class that contains information about a chemical compound.
// This information is derived from the ccp4 monomer library by default.
// To create compounds, you'd best use the factory method.
class Compound
{
public:
Compound(const std::string& id, const std::string& name,
const std::string& group, std::vector<CompoundAtom>&& atoms,
std::map<std::tuple<std::string,std::string>,float>&& bonds)
: mId(id), mName(name), mGroup(group)
, mAtoms(std::move(atoms)), mBonds(std::move(bonds))
{
}
~Compound();
// factory method, create a Compound based on the three letter code
// (for amino acids) or the one-letter code (for bases) or the
// code as it is known in the CCP4 monomer library.
static const Compound* create(const std::string& id);
// this second factory method can create a Compound even if it is not
// recorded in the library. It will take the values from the CCP4 lib
// unless the value passed to this function is not empty.
static const Compound* create(const std::string& id, const std::string& name,
const std::string& type, const std::string& formula);
// add an additional path to the monomer library.
static void addMonomerLibraryPath(const std::string& dir);
// accessors
std::string id() const { return mId; }
std::string name() const { return mName; }
std::string type() const;
// std::string group() const { return mGroup; }
std::vector<CompoundAtom> atoms() const { return mAtoms; }
CompoundAtom getAtomById(const std::string& atomId) const;
bool atomsBonded(const std::string& atomId_1, const std::string& atomId_2) const;
float atomBondValue(const std::string& atomId_1, const std::string& atomId_2) const;
std::string formula() const;
float formulaWeight() const;
int charge() const;
bool isWater() const;
private:
// Entity& mEntity;
std::string mId;
std::string mName;
std::string mGroup;
std::vector<CompoundAtom> mAtoms;
std::map<std::tuple<std::string,std::string>,float> mBonds;
};
// --------------------------------------------------------------------
// an Entity. This is a base class for PolymerEntity and NonPolyEntity
// The latter can be either a regular non-polymer (residue), a macrolide or
// water.
class Entity
{
public:
Entity(const std::string& id, const std::string& type, const std::string& description);
virtual ~Entity();
std::string id() const;
std::string type() const;
std::string description() const;
virtual float formulaWeight() const = 0;
private:
std::string mId;
std::string mType;
std::string mDescription;
};
// --------------------------------------------------------------------
// A polymer Entity
class PolymerEntity : public Entity
{
public:
PolymerEntity(const std::string& id, const std::string& description);
~PolymerEntity();
std::string seqOneLetterCode(bool cannonical) const;
std::string pdbxStrandId() const;
virtual float formulaWeight() const;
class monomer
{
public:
friend class PolymerEntity;
size_t num() const; // sequence number
bool hetero() const; // whether this position contains alternate Compounds
const Compound& comp(size_t altNr) const; // the chemical Compound of this monomer
private:
monomer* mNext;
monomer* mAlt;
size_t mNum;
Compound* mComp;
};
class iterator : public std::iterator<std::forward_iterator_tag, const monomer>
{
public:
typedef std::iterator<std::forward_iterator_tag, const monomer> baseType;
typedef baseType::reference reference;
typedef baseType::pointer pointer;
iterator(monomer* monomer = nullptr)
: mCursor(monomer) {}
iterator(const iterator& rhs)
: mCursor(rhs.mCursor)
{
}
iterator& operator=(const iterator& rhs)
{
mCursor = rhs.mCursor;
return *this;
}
reference operator*() { return *mCursor; }
pointer operator->() { return mCursor; }
iterator& operator++() { mCursor = mCursor->mNext; return *this; }
iterator operator++(int)
{
iterator tmp(*this);
operator++();
return tmp;
}
bool operator==(const iterator& rhs) const { return mCursor == rhs.mCursor; }
bool operator!=(const iterator& rhs) const { return mCursor != rhs.mCursor; }
private:
monomer* mCursor;
};
iterator begin() const { return iterator(mSeq); }
iterator end() const { return iterator(); }
const monomer& operator[](size_t index) const;
private:
Entity& mEntity;
monomer* mSeq;
};
// --------------------------------------------------------------------
// nonPoly Entity
class NonPolyEntity : public Entity
{
public:
NonPolyEntity(const std::string& id, const std::string& type, const std::string& description);
~NonPolyEntity();
Compound& comp() const;
virtual float formulaWeight() const;
private:
Compound* mCompound;
};
}
#pragma once
#include "libcif++/cif++.h"
// --------------------------------------------------------------------
struct PDBRecord
{
PDBRecord* mNext;
uint32 mLineNr;
char mName[11];
size_t mVlen;
char mValue[0];
PDBRecord(uint32 lineNr, const std::string& name, const std::string& value);
~PDBRecord();
void* operator new(size_t);
void* operator new(size_t size, size_t vLen);
void operator delete(void* p);
bool is(const char* name) const;
char vC(size_t column);
std::string vS(size_t columnFirst, size_t columnLast = std::numeric_limits<size_t>::max());
int vI(int columnFirst, int columnLast);
std::string vF(size_t columnFirst, size_t columnLast);
};
// --------------------------------------------------------------------
void ReadPDBFile(std::istream& pdbFile, cif::file& cifFile);
...@@ -11,15 +11,15 @@ class Remark3Parser ...@@ -11,15 +11,15 @@ class Remark3Parser
public: public:
virtual ~Remark3Parser() {} virtual ~Remark3Parser() {}
static bool Parse(const std::string& expMethod, PDBRecord* r, cif::datablock& db); static bool parse(const std::string& expMethod, PDBRecord* r, cif::datablock& db);
virtual std::string Program(); virtual std::string program();
virtual std::string Version(); virtual std::string version();
protected: protected:
Remark3Parser(const std::string& name, const std::string& expMethod, PDBRecord* r, cif::datablock& db, Remark3Parser(const std::string& name, const std::string& expMethod, PDBRecord* r, cif::datablock& db,
const TemplateLine templatelines[], uint32 templateLineCount, std::regex program_version); const TemplateLine templatelines[], uint32 templateLineCount, std::regex programVersion);
virtual float Parse(); virtual float Parse();
std::string NextLine(); std::string NextLine();
...@@ -31,17 +31,17 @@ class Remark3Parser ...@@ -31,17 +31,17 @@ class Remark3Parser
virtual void Fixup() {} virtual void Fixup() {}
std::string m_name; std::string mName;
std::string m_expMethod; std::string mExpMethod;
PDBRecord* m_rec; PDBRecord* mRec;
cif::datablock m_db; cif::datablock mDb;
std::string m_line; std::string mLine;
std::smatch m_m; std::smatch mM;
uint32 m_state; uint32 mState;
const TemplateLine* m_template; const TemplateLine* mTemplate;
uint32 m_templateCount; uint32 mTemplateCount;
std::regex m_program_version; std::regex mProgramVersion;
}; };
#pragma once
#include <map>
#include <string>
#include <boost/filesystem/path.hpp>
extern const std::map<std::string,char> kAAMap, kBaseMap;
class PeptideDB
{
public:
static PeptideDB& Instance();
void PushDictionary(boost::filesystem::path dict);
void PopDictionary();
bool IsKnownPeptide(const std::string& res_name) const;
bool IsKnownBase(const std::string& res_name) const;
std::string GetNameForResidue(const std::string& res_name) const;
std::string GetFormulaForResidue(const std::string& res_name) const;
std::string Unalias(const std::string& res_name) const;
private:
PeptideDB();
~PeptideDB();
PeptideDB(const PeptideDB&) = delete;
PeptideDB& operator=(const PeptideDB&) = delete;
struct PeptideDBImpl* mImpl;
static PeptideDB* sInstance;
};
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
#pragma once #pragma once
#include <libcif/config.h> #include "cif++/Config.h"
#include <boost/filesystem/operations.hpp> #include <boost/filesystem/operations.hpp>
#include <boost/math/quaternion.hpp> #include <boost/math/quaternion.hpp>
...@@ -19,76 +19,76 @@ const long double ...@@ -19,76 +19,76 @@ const long double
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// point, a location with x, y and z coordinates as float. // Point, a location with x, y and z coordinates as float.
// This one is derived from a tuple<float,float,float> so // This one is derived from a tuple<float,float,float> so
// you can do things like: // you can do things like:
// //
// float x, y, z; // float x, y, z;
// tie(x, y, z) = atom.loc(); // tie(x, y, z) = atom.loc();
struct point : public std::tuple<float,float,float> struct Point : public std::tuple<float,float,float>
{ {
typedef std::tuple<float,float,float> base_type; typedef std::tuple<float,float,float> base_type;
point() : base_type(0.f, 0.f, 0.f) {} Point() : base_type(0.f, 0.f, 0.f) {}
point(float x, float y, float z) : base_type(x, y, z) {} Point(float x, float y, float z) : base_type(x, y, z) {}
point(const clipper::Coord_orth& pt): base_type(pt[0], pt[1], pt[2]) {} Point(const clipper::Coord_orth& pt): base_type(pt[0], pt[1], pt[2]) {}
point& operator=(const clipper::Coord_orth& rhs) Point& operator=(const clipper::Coord_orth& rhs)
{ {
x(rhs[0]); setX(rhs[0]);
y(rhs[1]); setY(rhs[1]);
z(rhs[2]); setZ(rhs[2]);
return *this; return *this;
} }
float& x() { return std::get<0>(*this); } float& getX() { return std::get<0>(*this); }
float x() const { return std::get<0>(*this); } float getX() const { return std::get<0>(*this); }
void x(float x) { std::get<0>(*this) = x; } void setX(float x) { std::get<0>(*this) = x; }
float& y() { return std::get<1>(*this); } float& getY() { return std::get<1>(*this); }
float y() const { return std::get<1>(*this); } float getY() const { return std::get<1>(*this); }
void y(float y) { std::get<1>(*this) = y; } void setY(float y) { std::get<1>(*this) = y; }
float& z() { return std::get<2>(*this); } float& getZ() { return std::get<2>(*this); }
float z() const { return std::get<2>(*this); } float getZ() const { return std::get<2>(*this); }
void z(float z) { std::get<2>(*this) = z; } void setZ(float z) { std::get<2>(*this) = z; }
point& operator+=(const point& rhs) Point& operator+=(const Point& rhs)
{ {
x() += rhs.x(); getX() += rhs.getX();
y() += rhs.y(); getY() += rhs.getY();
z() += rhs.z(); getZ() += rhs.getZ();
return *this; return *this;
} }
point& operator-=(const point& rhs) Point& operator-=(const Point& rhs)
{ {
x() -= rhs.x(); getX() -= rhs.getX();
y() -= rhs.y(); getY() -= rhs.getY();
z() -= rhs.z(); getZ() -= rhs.getZ();
return *this; return *this;
} }
point& operator*=(float rhs) Point& operator*=(float rhs)
{ {
x() *= rhs; getX() *= rhs;
y() *= rhs; getY() *= rhs;
z() *= rhs; getZ() *= rhs;
return *this; return *this;
} }
point& operator/=(float rhs) Point& operator/=(float rhs)
{ {
x() *= rhs; getX() *= rhs;
y() *= rhs; getY() *= rhs;
z() *= rhs; getZ() *= rhs;
return *this; return *this;
} }
float normalize() float normalize()
{ {
auto length = x() * x() + y() * y() + z() * z(); auto length = getX() * getX() + getY() * getY() + getZ() * getZ();
if (length > 0) if (length > 0)
{ {
length = std::sqrt(length); length = std::sqrt(length);
...@@ -99,133 +99,132 @@ struct point : public std::tuple<float,float,float> ...@@ -99,133 +99,132 @@ struct point : public std::tuple<float,float,float>
void rotate(const boost::math::quaternion<float>& q) void rotate(const boost::math::quaternion<float>& q)
{ {
boost::math::quaternion<float> p(0, x(), y(), z()); boost::math::quaternion<float> p(0, getX(), getY(), getZ());
p = q * p * boost::math::conj(q); p = q * p * boost::math::conj(q);
x() = p.R_component_2(); getX() = p.R_component_2();
y() = p.R_component_3(); getY() = p.R_component_3();
z() = p.R_component_4(); getZ() = p.R_component_4();
} }
operator clipper::Coord_orth() const operator clipper::Coord_orth() const
{ {
return clipper::Coord_orth(x(), y(), z()); return clipper::Coord_orth(getX(), getY(), getZ());
} }
}; };
inline std::ostream& operator<<(std::ostream& os, const Point& pt)
inline std::ostream& operator<<(std::ostream& os, const point& pt)
{ {
os << '(' << pt.x() << ',' << pt.y() << ',' << pt.z() << ')'; os << '(' << pt.getX() << ',' << pt.getY() << ',' << pt.getZ() << ')';
return os; return os;
} }
inline point operator+(const point& lhs, const point& rhs) inline Point operator+(const Point& lhs, const Point& rhs)
{ {
return point(lhs.x() + rhs.x(), lhs.y() + rhs.y(), lhs.z() + rhs.z()); return Point(lhs.getX() + rhs.getX(), lhs.getY() + rhs.getY(), lhs.getZ() + rhs.getZ());
} }
inline point operator-(const point& lhs, const point& rhs) inline Point operator-(const Point& lhs, const Point& rhs)
{ {
return point(lhs.x() - rhs.x(), lhs.y() - rhs.y(), lhs.z() - rhs.z()); return Point(lhs.getX() - rhs.getX(), lhs.getY() - rhs.getY(), lhs.getZ() - rhs.getZ());
} }
inline point operator-(const point& pt) inline Point operator-(const Point& pt)
{ {
return point(-pt.x(), -pt.y(), -pt.z()); return Point(-pt.getX(), -pt.getY(), -pt.getZ());
} }
inline point operator*(const point& pt, float f) inline Point operator*(const Point& pt, float f)
{ {
return point(pt.x() * f, pt.y() * f, pt.z() * f); return Point(pt.getX() * f, pt.getY() * f, pt.getZ() * f);
} }
inline point operator/(const point& pt, float f) inline Point operator/(const Point& pt, float f)
{ {
return point(pt.x() / f, pt.y() / f, pt.z() / f); return Point(pt.getX() / f, pt.getY() / f, pt.getZ() / f);
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// several standard 3d operations // several standard 3d operations
inline double DistanceSquared(const point& a, const point& b) inline double DistanceSquared(const Point& a, const Point& b)
{ {
return return
(a.x() - b.x()) * (a.x() - b.x()) + (a.getX() - b.getX()) * (a.getX() - b.getX()) +
(a.y() - b.y()) * (a.y() - b.y()) + (a.getY() - b.getY()) * (a.getY() - b.getY()) +
(a.z() - b.z()) * (a.z() - b.z()); (a.getZ() - b.getZ()) * (a.getZ() - b.getZ());
} }
inline double Distance(const point& a, const point& b) inline double Distance(const Point& a, const Point& b)
{ {
return sqrt( return sqrt(
(a.x() - b.x()) * (a.x() - b.x()) + (a.getX() - b.getX()) * (a.getX() - b.getX()) +
(a.y() - b.y()) * (a.y() - b.y()) + (a.getY() - b.getY()) * (a.getY() - b.getY()) +
(a.z() - b.z()) * (a.z() - b.z())); (a.getZ() - b.getZ()) * (a.getZ() - b.getZ()));
} }
inline float DotProduct(const point& a, const point& b) inline float DotProduct(const Point& a, const Point& b)
{ {
return a.x() * b.x() + a.y() * b.y() + a.z() * b.z(); return a.getX() * b.getX() + a.getY() * b.getY() + a.getZ() * b.getZ();
} }
inline point CrossProduct(const point& a, const point& b) inline Point CrossProduct(const Point& a, const Point& b)
{ {
return point(a.y() * b.z() - b.y() * a.z(), return Point(a.getY() * b.getZ() - b.getY() * a.getZ(),
a.z() * b.x() - b.z() * a.x(), a.getZ() * b.getX() - b.getZ() * a.getX(),
a.x() * b.y() - b.x() * a.y()); a.getX() * b.getY() - b.getX() * a.getY());
} }
float DihedralAngle(const point& p1, const point& p2, const point& p3, const point& p4); float DihedralAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4);
float CosinusAngle(const point& p1, const point& p2, const point& p3, const point& p4); float CosinusAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// We use quaternions to do rotations in 3d space // We use quaternions to do rotations in 3d space
quaternion Normalize(quaternion q); quaternion Normalize(quaternion q);
//std::tuple<double,point> QuaternionToAngleAxis(quaternion q); //std::tuple<double,Point> QuaternionToAngleAxis(quaternion q);
point Centroid(std::vector<point>& points); Point Centroid(std::vector<Point>& Points);
point CenterPoints(std::vector<point>& points); Point CenterPoints(std::vector<Point>& Points);
quaternion AlignPoints(const std::vector<point>& a, const std::vector<point>& b); quaternion AlignPoints(const std::vector<Point>& a, const std::vector<Point>& b);
double RMSd(const std::vector<point>& a, const std::vector<point>& b); double RMSd(const std::vector<Point>& a, const std::vector<Point>& b);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// Helper class to generate evenly divided points on a sphere // Helper class to generate evenly divided Points on a sphere
// we use a fibonacci sphere to calculate even distribution of the dots // we use a fibonacci sphere to calculate even distribution of the dots
template<int N> template<int N>
class spherical_dots class SphericalDots
{ {
public: public:
enum { P = 2 * N + 1 }; enum { P = 2 * N + 1 };
typedef typename std::array<point,P> array_type; typedef typename std::array<Point,P> array_type;
typedef typename array_type::const_iterator iterator; typedef typename array_type::const_iterator iterator;
static spherical_dots& instance() static SphericalDots& instance()
{ {
static spherical_dots s_instance; static SphericalDots sInstance;
return s_instance; return sInstance;
} }
size_t size() const { return m_points.size(); } size_t size() const { return mPoints.size(); }
const point operator[](uint32 inIx) const { return m_points[inIx]; } const Point operator[](uint32 inIx) const { return mPoints[inIx]; }
iterator begin() const { return m_points.begin(); } iterator begin() const { return mPoints.begin(); }
iterator end() const { return m_points.end(); } iterator end() const { return mPoints.end(); }
double weight() const { return m_weight; } double weight() const { return mWeight; }
spherical_dots() SphericalDots()
{ {
using namespace std; using namespace std;
const double const double
kGoldenRatio = (1 + std::sqrt(5.0)) / 2; kGoldenRatio = (1 + std::sqrt(5.0)) / 2;
m_weight = (4 * kPI) / P; mWeight = (4 * kPI) / P;
auto p = m_points.begin(); auto p = mPoints.begin();
for (int32 i = -N; i <= N; ++i) for (int32 i = -N; i <= N; ++i)
{ {
...@@ -242,11 +241,11 @@ class spherical_dots ...@@ -242,11 +241,11 @@ class spherical_dots
private: private:
array_type m_points; array_type mPoints;
double m_weight; double mWeight;
}; };
typedef spherical_dots<50> spherical_dots_50; typedef SphericalDots<50> SphericalDots_50;
} }
// Lib for working with structures as contained in mmCIF and PDB files
#pragma once
#include <boost/filesystem/operations.hpp>
#include <boost/math/quaternion.hpp>
#include <boost/any.hpp>
#include "cif++/AtomType.h"
#include "cif++/Point.h"
#include "cif++/Compound.h"
/*
To modify a structure, you will have to use actions.
The currently supported actions are:
// - Move atom to new location
- Remove atom
// - Add new atom that was formerly missing
// - Add alternate Residue
-
Other important design principles:
- all objects here are references to the actual data. Not models of
the data itself. That means that if you copy an atom, you copy the
reference to an atom in the structure. You're not creating a new
atom. This may sound obvious, but it is not if you are used to
copy semantics in the C++ world.
*/
// forward declaration
namespace cif
{
class Datablock;
};
namespace libcif
{
class Atom;
class Residue;
class Monomer;
class Polymer;
class Structure;
class File;
// --------------------------------------------------------------------
// We do not want to introduce a dependency on cif++ here, we might want
// to change the backend storage in the future.
// So, in order to access the data we use properties based on boost::any
// Eventually this should be moved to std::variant, but that's only when
// c++17 is acceptable.
struct Property
{
Property() {}
Property(const std::string& name, const boost::any& value)
: name(name), value(value) {}
std::string name;
boost::any value;
};
typedef std::vector<Property> PropertyList;
// --------------------------------------------------------------------
class Atom
{
public:
// Atom(const structure& s, const std::string& id);
Atom(struct AtomImpl* impl);
Atom(const File& f, const std::string& id);
Atom(const Atom& rhs);
~Atom();
Atom& operator=(const Atom& rhs);
std::string id() const;
AtomType type() const;
point location() const;
const compound& comp() const;
const entity& ent() const;
bool isWater() const;
int charge() const;
boost::any property(const std::string& name) const;
void property(const std::string& name, const boost::any& value);
// specifications
std::string labelAtomId() const;
std::string labelCompId() const;
std::string labelAsymId() const;
int labelSeqId() const;
std::string labelAltId() const;
std::string authAtomId() const;
std::string authCompId() const;
std::string authAsymId() const;
int authSeqId() const;
std::string pdbxAuthInsCode() const;
std::string authAltId() const;
bool operator==(const Atom& rhs) const;
const File& getFile() const;
private:
struct AtomImpl* mImpl;
};
typedef std::vector<Atom> AtomView;
// --------------------------------------------------------------------
class Residue : public std::enable_shared_from_this<Residue>
{
public:
Residue(const compound& cmp) : mCompound(cmp) {}
const compound& comp() const { return mCompound; }
virtual AtomView atoms();
private:
const compound& mCompound;
};
//// --------------------------------------------------------------------
//// a monomer models a single Residue in a protein chain
//
//class monomer : public Residue
//{
// public:
// monomer(polymer& polymer, size_t seqId, const std::string& compId,
// const std::string& altId);
//
// int num() const { return mNum; }
//// polymer& getPolymer();
//
//// std::vector<monomer_ptr> alternates();
//
// private:
// polymer_ptr mPolymer;
// int mNum;
//};
//
//// --------------------------------------------------------------------
//
//class polymer : public std::enable_shared_from_this<polymer>
//{
// public:
// polymer(const polymerEntity& pe, const std::string& asymId);
//
// struct iterator : public std::iterator<std::random_access_iterator_tag, monomer>
// {
// typedef std::iterator<std::bidirectional_iterator_tag, monomer> base_type;
// typedef base_type::reference reference;
// typedef base_type::pointer pointer;
//
// iterator(polymer& list, uint32 index);
// iterator(iterator&& rhs);
// iterator(const iterator& rhs);
// iterator& operator=(const iterator& rhs);
// iterator& operator=(iterator&& rhs);
//
// reference operator*();
// pointer operator->();
//
// iterator& operator++();
// iterator operator++(int);
//
// iterator& operator--();
// iterator operator--(int);
//
// bool operator==(const iterator& rhs) const;
// bool operator!=(const iterator& rhs) const;
// };
//
// iterator begin();
// iterator end();
//
// private:
// polymer_entity mEntity;
// std::string mAsymId;
// std::vector<Residue_ptr> mMonomers;
//};
// --------------------------------------------------------------------
// file is a reference to the data stored in e.g. the cif file.
// This object is not copyable.
class File : public std::enable_shared_from_this<File>
{
public:
File();
File(boost::filesystem::path p);
~File();
File(const File&) = delete;
File& operator=(const File&) = delete;
void load(boost::filesystem::path p);
void save(boost::filesystem::path p);
structure* model(size_t nr = 1);
struct FileImpl& impl() const { return *mImpl; }
std::vector<const entity*> entities();
cif::datablock& data();
private:
struct FileImpl* mImpl;
};
// --------------------------------------------------------------------
class structure
{
public:
structure(File& p, uint32 modelNr = 1);
structure(const structure&);
structure& operator=(const structure&);
~structure();
File& getFile() const;
AtomView atoms() const;
AtomView waters() const;
Atom getAtomById(std::string id) const;
Atom getAtomByLocation(point pt, float maxDistance) const;
Atom getAtomForLabel(const std::string& atomId, const std::string& asymId,
const std::string& compId, int seqId, const std::string& altId = "");
Atom getAtomForAuth(const std::string& atomId, const std::string& asymId,
const std::string& compId, int seqId, const std::string& altId = "",
const std::string& pdbxAuthInsCode = "");
// map between auth and label locations
std::tuple<std::string,int,std::string> MapAuthToLabel(const std::string& asymId,
const std::string& seqId, const std::string& compId, const std::string& insCode = "");
std::tuple<std::string,std::string,std::string,std::string> MapLabelToAuth(
const std::string& asymId, int seqId, const std::string& compId);
// returns chain, seqnr
std::tuple<std::string,std::string> MapLabelToAuth(
const std::string& asymId, int seqId);
// returns chain,seqnr,comp,iCode
std::tuple<std::string,int,std::string,std::string> MapLabelToPDB(
const std::string& asymId, int seqId, const std::string& compId);
std::tuple<std::string,int,std::string,std::string> MapPDBToLabel(
const std::string& asymId, int seqId, const std::string& compId, const std::string& iCode);
// Actions
void removeAtom(Atom& a);
private:
struct StructureImpl* mImpl;
};
}
// cif parsing library
#include "libcif/cif++.h"
#include <boost/filesystem/path.hpp>
// the std regex of gcc is crashing....
#include <boost/regex.hpp>
#include <set>
namespace cif
{
struct validate_category;
// --------------------------------------------------------------------
class validation_error : public std::exception
{
public:
validation_error(const std::string& msg) : m_msg(msg) {}
const char* what() const noexcept { return m_msg.c_str(); }
std::string m_msg;
};
// --------------------------------------------------------------------
enum DDL_PrimitiveType
{
ptChar, ptUChar, ptNumb
};
DDL_PrimitiveType map_to_primitive_type(const std::string& s);
struct validate_type
{
std::string m_name;
DDL_PrimitiveType m_primitive_type;
boost::regex m_rx;
bool operator<(const validate_type& rhs) const
{
return icompare(m_name, rhs.m_name) < 0;
}
// compare values based on type
// int compare(const std::string& a, const std::string& b) const
// {
// return compare(a.c_str(), b.c_str());
// }
int compare(const char* a, const char* b) const;
};
struct validate_item
{
std::string m_tag;
bool m_mandatory;
const validate_type* m_type;
cif::iset m_enums;
validate_item* m_parent = nullptr;
std::set<validate_item*>
m_children;
validate_category* m_category = nullptr;
std::set<validate_item*>
m_foreign_keys;
void set_parent(validate_item* parent);
bool operator<(const validate_item& rhs) const
{
return icompare(m_tag, rhs.m_tag) < 0;
}
bool operator==(const validate_item& rhs) const
{
return iequals(m_tag, rhs.m_tag);
}
void operator()(std::string value) const;
};
struct validate_category
{
std::string m_name;
std::vector<string> m_keys;
cif::iset m_groups;
cif::iset m_mandatory_fields;
std::set<validate_item> m_item_validators;
bool operator<(const validate_category& rhs) const
{
return icompare(m_name, rhs.m_name) < 0;
}
void add_item_validator(validate_item&& v);
const validate_item* get_validator_for_item(std::string tag) const;
const std::set<validate_item>& item_validators() const
{
return m_item_validators;
}
};
// --------------------------------------------------------------------
class validator
{
public:
friend class dict_parser;
validator();
~validator();
validator(const validator& rhs) = delete;
validator& operator=(const validator& rhs) = delete;
validator(validator&& rhs);
validator& operator=(validator&& rhs);
void add_type_validator(validate_type&& v);
const validate_type* get_validator_for_type(std::string type_code) const;
void add_category_validator(validate_category&& v);
const validate_category* get_validator_for_category(std::string category) const;
void report_error(const std::string& msg);
std::string dict_name() const { return m_name; }
void dict_name(const std::string& name) { m_name = name; }
std::string dict_version() const { return m_version; }
void dict_version(const std::string& version) { m_version = version; }
private:
// name is fully qualified here:
validate_item* get_validator_for_item(std::string name) const;
std::string m_name;
std::string m_version;
bool m_strict = false;
// std::set<uint32> m_sub_categories;
std::set<validate_type> m_type_validators;
std::set<validate_category> m_category_validators;
};
}
#pragma once
#include "cif++.h"
void WritePDBFile(std::ostream& pdbFile, cif::file& cifFile);
// Lib for working with structures as contained in mmCIF and PDB files
#pragma once
#include <set>
#include <tuple>
#include <vector>
#include <map>
#include "libcif/atom_type.h"
namespace libcif
{
// --------------------------------------------------------------------
// The chemical composition of the structure in an mmCIF file is
// defined in the class composition. A compositon consists of
// entities. Each entity can be either a polymer, a non-polymer
// a macrolide or a water molecule.
// Entities themselves are made up of compounds. And compounds
// contain comp_atom records for each atom.
class composition;
class entity;
class compound;
struct comp_atom;
// --------------------------------------------------------------------
// struct containing information about an atom in a chemical compound
// This information comes from the CCP4 monomer library.
struct comp_atom
{
std::string id;
atom_type type_symbol;
std::string type_energy;
float partial_charge;
};
// --------------------------------------------------------------------
// a class that contains information about a chemical compound.
// This information is derived from the ccp4 monomer library by default.
// To create compounds, you'd best use the factory method.
class compound
{
public:
compound(const std::string& id, const std::string& name,
const std::string& group, std::vector<comp_atom>&& atoms,
std::map<std::tuple<std::string,std::string>,float>&& bonds)
: m_id(id), m_name(name), m_group(group)
, m_atoms(std::move(atoms)), m_bonds(std::move(bonds))
{
}
~compound();
// factory method, create a compound based on the three letter code
// (for amino acids) or the one-letter code (for bases) or the
// code as it is known in the CCP4 monomer library.
static const compound* create(const std::string& id);
// this second factory method can create a compound even if it is not
// recorded in the library. It will take the values from the CCP4 lib
// unless the value passed to this function is not empty.
static const compound* create(const std::string& id, const std::string& name,
const std::string& type, const std::string& formula);
// add an additional path to the monomer library.
static void add_monomer_library_path(const std::string& dir);
// accessors
std::string id() const { return m_id; }
std::string name() const { return m_name; }
std::string type() const;
// std::string group() const { return m_group; }
std::vector<comp_atom> atoms() const { return m_atoms; }
comp_atom get_atom_by_id(const std::string& atom_id) const;
bool atoms_bonded(const std::string& atom_id_1, const std::string& atom_id_2) const;
float atom_bond_value(const std::string& atom_id_1, const std::string& atom_id_2) const;
std::string formula() const;
float formula_weight() const;
int charge() const;
bool is_water() const;
private:
// entity& m_entity;
std::string m_id;
std::string m_name;
std::string m_group;
std::vector<comp_atom> m_atoms;
std::map<std::tuple<std::string,std::string>,float> m_bonds;
};
// --------------------------------------------------------------------
// an entity. This is a base class for polymer_entity and non_poly_entity
// The latter can be either a regular non-polymer (residue), a macrolide or
// water.
class entity
{
public:
entity(const std::string& id, const std::string& type, const std::string& description);
virtual ~entity();
std::string id() const;
std::string type() const;
std::string description() const;
virtual float formula_weight() const = 0;
private:
std::string m_id;
std::string m_type;
std::string m_description;
};
// --------------------------------------------------------------------
// A polymer entity
class polymer_entity : public entity
{
public:
polymer_entity(const std::string& id, const std::string& description);
~polymer_entity();
std::string seq_one_letter_code(bool cannonical) const;
std::string pdbx_strand_id() const;
virtual float formula_weight() const;
class monomer
{
public:
friend class polymer_entity;
size_t num() const; // sequence number
bool hetero() const; // whether this position contains alternate compounds
const compound& comp(size_t alt_nr) const; // the chemical compound of this monomer
private:
monomer* m_next;
monomer* m_alt;
size_t m_num;
compound* m_comp;
};
class iterator : public std::iterator<std::forward_iterator_tag, const monomer>
{
public:
typedef std::iterator<std::forward_iterator_tag, const monomer> base_type;
typedef base_type::reference reference;
typedef base_type::pointer pointer;
iterator(monomer* monomer = nullptr)
: m_cursor(monomer) {}
iterator(const iterator& rhs)
: m_cursor(rhs.m_cursor)
{
}
iterator& operator=(const iterator& rhs)
{
m_cursor = rhs.m_cursor;
return *this;
}
reference operator*() { return *m_cursor; }
pointer operator->() { return m_cursor; }
iterator& operator++() { m_cursor = m_cursor->m_next; return *this; }
iterator operator++(int)
{
iterator tmp(*this);
operator++();
return tmp;
}
bool operator==(const iterator& rhs) const { return m_cursor == rhs.m_cursor; }
bool operator!=(const iterator& rhs) const { return m_cursor != rhs.m_cursor; }
private:
monomer* m_cursor;
};
iterator begin() const { return iterator(m_seq); }
iterator end() const { return iterator(); }
const monomer& operator[](size_t index) const;
private:
entity& m_entity;
monomer* m_seq;
};
// --------------------------------------------------------------------
// non_poly entity
class non_poly_entity : public entity
{
public:
non_poly_entity(const std::string& id, const std::string& type, const std::string& description);
~non_poly_entity();
compound& comp() const;
virtual float formula_weight() const;
private:
compound* m_compound;
};
}
// Copyright Maarten L. Hekkelman 2006-2010
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef MRSRC_H
#define MRSRC_H
#include <string>
#include <list>
#include <exception>
/*
Resources are data sources for the application.
They are retrieved by name.
Basic usage:
mrsrc::rsrc rsrc("dialogs/my-cool-dialog.glade");
if (rsrc)
{
GladeXML* glade = glade_xml_new_from_buffer(rsrc.data(), rsrc.size(), NULL, "japi");
...
}
*/
namespace mrsrc {
struct rsrc_imp
{
unsigned int m_next;
unsigned int m_child;
unsigned int m_name;
unsigned int m_size;
unsigned int m_data;
};
}
// The following three variables are generated by the japi resource compiler:
extern const mrsrc::rsrc_imp gResourceIndex[];
extern const char gResourceData[];
extern const char gResourceName[];
namespace mrsrc
{
class rsrc_not_found_exception : public std::exception
{
public:
virtual const char* what() const throw() { return "resource not found"; }
};
class rsrc;
typedef std::list<rsrc> rsrc_list;
class rsrc
{
public:
rsrc() : m_impl(gResourceIndex) {}
rsrc(const rsrc& other)
: m_impl(other.m_impl) {}
rsrc& operator=(const rsrc& other)
{
m_impl = other.m_impl;
return *this;
}
rsrc(const std::string& path);
std::string name() const { return gResourceName + m_impl->m_name; }
const char* data() const { return gResourceData + m_impl->m_data; }
unsigned long size() const { return m_impl->m_size; }
operator bool () const { return m_impl != NULL and m_impl->m_size > 0; }
rsrc_list children() const;
private:
rsrc(const rsrc_imp* imp)
: m_impl(imp) {}
const rsrc_imp* m_impl;
};
inline
rsrc_list rsrc::children() const
{
rsrc_list result;
if (m_impl->m_child)
{
const rsrc_imp* impl = gResourceIndex + m_impl->m_child;
result.push_back(rsrc(impl));
while (impl->m_next)
{
impl = gResourceIndex + impl->m_next;
result.push_back(rsrc(impl));
}
}
return result;
}
inline
rsrc::rsrc(const std::string& path)
{
// static_assert(sizeof(m_impl->m_next) == 4, "invalid size for unsigned int");
m_impl = gResourceIndex;
std::string p(path);
// would love to use boost functions here, but then the dependancies
// should be minimal of course.
while (not p.empty())
{
if (m_impl->m_child == 0) // no children, this is an error
throw rsrc_not_found_exception();
m_impl = gResourceIndex + m_impl->m_child;
std::string::size_type s = p.find('/');
std::string name;
if (s != std::string::npos)
{
name = p.substr(0, s);
p.erase(0, s + 1);
}
else
std::swap(name, p);
while (name != gResourceName + m_impl->m_name)
{
if (m_impl->m_next == 0)
throw rsrc_not_found_exception();
m_impl = gResourceIndex + m_impl->m_next;
}
}
}
}
#endif
#pragma once
#include "cif++.h"
// --------------------------------------------------------------------
struct PDBRecord
{
PDBRecord* m_next;
uint32 m_line_nr;
char m_name[11];
size_t m_vlen;
char m_value[0];
PDBRecord(uint32 line_nr, const std::string& name, const std::string& value);
~PDBRecord();
void* operator new(size_t);
void* operator new(size_t size, size_t v_len);
void operator delete(void* p);
bool is(const char* name) const;
char v_c(size_t column);
std::string v_s(size_t column_first, size_t column_last = std::numeric_limits<size_t>::max());
int v_i(int column_first, int column_last);
std::string v_f(size_t column_first, size_t column_last);
};
// --------------------------------------------------------------------
void ReadPDBFile(std::istream& pdbFile, cif::file& cifFile);
// Lib for working with structures as contained in mmCIF and PDB files
#pragma once
#include <boost/filesystem/operations.hpp>
#include <boost/math/quaternion.hpp>
#include <boost/any.hpp>
#include "libcif/atom_type.h"
#include "libcif/point.h"
#include "libcif/compound.h"
/*
To modify a structure, you will have to use actions.
The currently supported actions are:
// - Move atom to new location
- Remove atom
// - Add new atom that was formerly missing
// - Add alternate residue
-
Other important design principles:
- all objects here are references to the actual data. Not models of
the data itself. That means that if you copy an atom, you copy the
reference to an atom in the structure. You're not creating a new
atom. This may sound obvious, but it is not if you are used to
copy semantics in the C++ world.
*/
// forward declaration
namespace cif
{
class datablock;
};
namespace libcif
{
class atom;
class residue;
class monomer;
class polymer;
class structure;
class file;
// --------------------------------------------------------------------
// We do not want to introduce a dependency on cif++ here, we might want
// to change the backend storage in the future.
// So, in order to access the data we use properties based on boost::any
// Eventually this should be moved to std::variant, but that's only when
// c++17 is acceptable.
struct property
{
property() {}
property(const std::string& name, const boost::any& value)
: name(name), value(value) {}
std::string name;
boost::any value;
};
typedef std::vector<property> property_list;
// --------------------------------------------------------------------
class atom
{
public:
// atom(const structure& s, const std::string& id);
atom(struct atom_impl* impl);
atom(const file& f, const std::string& id);
atom(const atom& rhs);
~atom();
atom& operator=(const atom& rhs);
std::string id() const;
atom_type type() const;
point location() const;
const compound& comp() const;
const entity& ent() const;
bool is_water() const;
int charge() const;
boost::any property(const std::string& name) const;
void property(const std::string& name, const boost::any& value);
// specifications
std::string label_atom_id() const;
std::string label_comp_id() const;
std::string label_asym_id() const;
int label_seq_id() const;
std::string label_alt_id() const;
std::string auth_atom_id() const;
std::string auth_comp_id() const;
std::string auth_asym_id() const;
int auth_seq_id() const;
std::string pdbx_auth_ins_code() const;
std::string auth_alt_id() const;
bool operator==(const atom& rhs) const;
const file& get_file() const;
private:
struct atom_impl* m_impl;
};
typedef std::vector<atom> atom_view;
// --------------------------------------------------------------------
class residue : public std::enable_shared_from_this<residue>
{
public:
residue(const compound& cmp) : m_compound(cmp) {}
const compound& comp() const { return m_compound; }
virtual atom_view atoms();
private:
const compound& m_compound;
};
//// --------------------------------------------------------------------
//// a monomer models a single residue in a protein chain
//
//class monomer : public residue
//{
// public:
// monomer(polymer& polymer, size_t seq_id, const std::string& comp_id,
// const std::string& alt_id);
//
// int num() const { return m_num; }
//// polymer& get_polymer();
//
//// std::vector<monomer_ptr> alternates();
//
// private:
// polymer_ptr m_polymer;
// int m_num;
//};
//
//// --------------------------------------------------------------------
//
//class polymer : public std::enable_shared_from_this<polymer>
//{
// public:
// polymer(const polymer_entity& pe, const std::string& asym_id);
//
// struct iterator : public std::iterator<std::random_access_iterator_tag, monomer>
// {
// typedef std::iterator<std::bidirectional_iterator_tag, monomer> base_type;
// typedef base_type::reference reference;
// typedef base_type::pointer pointer;
//
// iterator(polymer& list, uint32 index);
// iterator(iterator&& rhs);
// iterator(const iterator& rhs);
// iterator& operator=(const iterator& rhs);
// iterator& operator=(iterator&& rhs);
//
// reference operator*();
// pointer operator->();
//
// iterator& operator++();
// iterator operator++(int);
//
// iterator& operator--();
// iterator operator--(int);
//
// bool operator==(const iterator& rhs) const;
// bool operator!=(const iterator& rhs) const;
// };
//
// iterator begin();
// iterator end();
//
// private:
// polymer_entity m_entity;
// std::string m_asym_id;
// std::vector<residue_ptr> m_monomers;
//};
// --------------------------------------------------------------------
// file is a reference to the data stored in e.g. the cif file.
// This object is not copyable.
class file : public std::enable_shared_from_this<file>
{
public:
file();
file(boost::filesystem::path p);
~file();
file(const file&) = delete;
file& operator=(const file&) = delete;
void load(boost::filesystem::path p);
void save(boost::filesystem::path p);
structure* model(size_t nr = 1);
struct file_impl& impl() const { return *m_impl; }
std::vector<const entity*> entities();
cif::datablock& data();
private:
struct file_impl* m_impl;
};
// --------------------------------------------------------------------
class structure
{
public:
structure(file& p, uint32 model_nr = 1);
structure(const structure&);
structure& operator=(const structure&);
~structure();
file& get_file() const;
atom_view atoms() const;
atom_view waters() const;
atom get_atom_by_id(std::string id) const;
atom get_atom_by_location(point pt, float max_distance) const;
atom get_atom_for_label(const std::string& atom_id, const std::string& asym_id,
const std::string& comp_id, int seq_id, const std::string& alt_id = "");
atom get_atom_for_auth(const std::string& atom_id, const std::string& asym_id,
const std::string& comp_id, int seq_id, const std::string& alt_id = "",
const std::string& pdbx_auth_ins_code = "");
// map between auth and label locations
std::tuple<std::string,int,std::string> MapAuthToLabel(const std::string& asym_id,
const std::string& seq_id, const std::string& comp_id, const std::string& ins_code = "");
std::tuple<std::string,std::string,std::string,std::string> MapLabelToAuth(
const std::string& asym_id, int seq_id, const std::string& comp_id);
// returns chain, seqnr
std::tuple<std::string,std::string> MapLabelToAuth(
const std::string& asym_id, int seq_id);
// returns chain,seqnr,comp,iCode
std::tuple<std::string,int,std::string,std::string> MapLabelToPDB(
const std::string& asym_id, int seq_id, const std::string& comp_id);
std::tuple<std::string,int,std::string,std::string> MapPDBToLabel(
const std::string& asym_id, int seq_id, const std::string& comp_id, const std::string& iCode);
// Actions
void remove_atom(atom& a);
private:
friend class action;
struct structure_impl* m_impl;
};
}
// Lib for working with structures as contained in mmCIF and PDB files // Lib for working with structures as contained in mmCIF and PDB files
#include "libcif/atom_type.h" #include "cif++/AtomType.h"
#include "libcif/cif++.h" #include "cif++/Cif++.h"
using namespace std; using namespace std;
...@@ -10,7 +10,7 @@ namespace libcif ...@@ -10,7 +10,7 @@ namespace libcif
const float kNA = nan("1"); const float kNA = nan("1");
const atom_type_info kKnownAtoms[] = const AtomTypeInfo kKnownAtoms[] =
{ {
{ Nn, "Unknown", "Nn", 0, false, { kNA, kNA, kNA, kNA, kNA, kNA, kNA } }, { Nn, "Unknown", "Nn", 0, false, { kNA, kNA, kNA, kNA, kNA, kNA, kNA } },
{ H, "Hydro­gen", "H", 1.008, false, { 53, 25, 37, 32, kNA, kNA, 120 } }, { H, "Hydro­gen", "H", 1.008, false, { 53, 25, 37, 32, kNA, kNA, 120 } },
...@@ -133,35 +133,35 @@ const atom_type_info kKnownAtoms[] = ...@@ -133,35 +133,35 @@ const atom_type_info kKnownAtoms[] =
{ Lr, "Lawren­cium", "Lr", 266, true, { kNA, kNA, kNA, 161, 141, kNA, kNA } } { Lr, "Lawren­cium", "Lr", 266, true, { kNA, kNA, kNA, 161, 141, kNA, kNA } }
}; };
uint32 kKnownAtomsCount = sizeof(kKnownAtoms) / sizeof(atom_type_info); uint32 kKnownAtomsCount = sizeof(kKnownAtoms) / sizeof(AtomTypeInfo);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// atom_type_traits // AtomTypeTraits
atom_type_traits::atom_type_traits(const string& symbol) AtomTypeTraits::AtomTypeTraits(const string& symbol)
: m_info(nullptr) : mInfo(nullptr)
{ {
for (auto& i: kKnownAtoms) for (auto& i: kKnownAtoms)
{ {
if (cif::iequals(i.symbol, symbol)) if (cif::iequals(i.symbol, symbol))
{ {
m_info = &i; mInfo = &i;
break; break;
} }
} }
if (m_info == nullptr) if (mInfo == nullptr)
throw invalid_argument("Not a known element: " + symbol); throw invalid_argument("Not a known element: " + symbol);
} }
atom_type_traits::atom_type_traits(atom_type t) AtomTypeTraits::AtomTypeTraits(AtomType t)
{ {
if (t < H or t > Lr) if (t < H or t > Lr)
throw invalid_argument("atom_type out of range"); throw invalid_argument("atomType out of range");
m_info = &kKnownAtoms[t]; mInfo = &kKnownAtoms[t];
} }
bool atom_type_traits::is_element(const string& symbol) bool AtomTypeTraits::isElement(const string& symbol)
{ {
bool result = false; bool result = false;
...@@ -177,7 +177,7 @@ bool atom_type_traits::is_element(const string& symbol) ...@@ -177,7 +177,7 @@ bool atom_type_traits::is_element(const string& symbol)
return result; return result;
} }
bool atom_type_traits::is_metal(const std::string& symbol) bool AtomTypeTraits::isMetal(const string& symbol)
{ {
bool result = false; bool result = false;
......
This diff is collapsed. Click to expand it.
// CIF parser // CIF parser
#include "libcif/config.h" #include "cif++/Config.h"
#include <tuple> #include <tuple>
#include <iostream> #include <iostream>
#include <boost/algorithm/string.hpp> #include <boost/algorithm/string.hpp>
#include "libcif/cif-utils.h" #include "cif++/CifUtils.h"
using namespace std; using namespace std;
namespace ba = boost::algorithm; namespace ba = boost::algorithm;
......
...@@ -5,9 +5,9 @@ ...@@ -5,9 +5,9 @@
// since gcc's regex is crashing.... // since gcc's regex is crashing....
#include <boost/regex.hpp> #include <boost/regex.hpp>
#include "libcif/cif++.h" #include "cif++/Cif++.h"
#include "libcif/cif-parser.h" #include "cif++/CifParser.h"
#include "libcif/cif-validator.h" #include "cif++/CifValidator.h"
using namespace std; using namespace std;
namespace ba = boost::algorithm; namespace ba = boost::algorithm;
...@@ -17,7 +17,7 @@ extern int VERBOSE; ...@@ -17,7 +17,7 @@ extern int VERBOSE;
namespace cif namespace cif
{ {
DDL_PrimitiveType map_to_primitive_type(const string& s) DDL_PrimitiveType mapToPrimitiveType(const string& s)
{ {
DDL_PrimitiveType result; DDL_PrimitiveType result;
if (iequals(s, "char")) if (iequals(s, "char"))
...@@ -27,13 +27,13 @@ DDL_PrimitiveType map_to_primitive_type(const string& s) ...@@ -27,13 +27,13 @@ DDL_PrimitiveType map_to_primitive_type(const string& s)
else if (iequals(s, "numb")) else if (iequals(s, "numb"))
result = ptNumb; result = ptNumb;
else else
throw validation_error("Not a known primitive type"); throw ValidationError("Not a known primitive type");
return result; return result;
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
int validate_type::compare(const char* a, const char* b) const int ValidateType::compare(const char* a, const char* b) const
{ {
int result = 0; int result = 0;
...@@ -45,7 +45,7 @@ int validate_type::compare(const char* a, const char* b) const ...@@ -45,7 +45,7 @@ int validate_type::compare(const char* a, const char* b) const
{ {
try try
{ {
switch (m_primitive_type) switch (mPrimitiveType)
{ {
case ptNumb: case ptNumb:
{ {
...@@ -119,59 +119,59 @@ int validate_type::compare(const char* a, const char* b) const ...@@ -119,59 +119,59 @@ int validate_type::compare(const char* a, const char* b) const
// -------------------------------------------------------------------- // --------------------------------------------------------------------
void validate_item::set_parent(validate_item* parent) void ValidateItem::setParent(ValidateItem* parent)
{ {
m_parent = parent; mParent = parent;
if (m_type == nullptr and m_parent != nullptr) if (mType == nullptr and mParent != nullptr)
m_type = m_parent->m_type; mType = mParent->mType;
if (m_parent != nullptr) if (mParent != nullptr)
{ {
m_parent->m_children.insert(this); mParent->mChildren.insert(this);
if (m_category->m_keys == vector<string>{m_tag}) if (mCategory->mKeys == vector<string>{mTag})
m_parent->m_foreign_keys.insert(this); mParent->mForeignKeys.insert(this);
} }
} }
void validate_item::operator()(string value) const void ValidateItem::operator()(string value) const
{ {
if (VERBOSE >= 4) if (VERBOSE >= 4)
cout << "validating '" << value << "' for '" << m_tag << "'" << endl; cout << "validating '" << value << "' for '" << mTag << "'" << endl;
if (not value.empty() and value != "?" and value != ".") if (not value.empty() and value != "?" and value != ".")
{ {
if (m_type != nullptr and not boost::regex_match(value, m_type->m_rx)) if (mType != nullptr and not boost::regex_match(value, mType->mRx))
throw validation_error("Value '" + value + "' does not match type expression for type " + m_type->m_name + " in item " + m_tag); throw ValidationError("Value '" + value + "' does not match type expression for type " + mType->mName + " in item " + mTag);
if (not m_enums.empty()) if (not mEnums.empty())
{ {
if (m_enums.count(value) == 0) if (mEnums.count(value) == 0)
throw validation_error("Value '" + value + "' is not in the list of allowed values for item " + m_tag); throw ValidationError("Value '" + value + "' is not in the list of allowed values for item " + mTag);
} }
} }
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
void validate_category::add_item_validator(validate_item&& v) void ValidateCategory::addItemValidator(ValidateItem&& v)
{ {
if (v.m_mandatory) if (v.mMandatory)
m_mandatory_fields.insert(v.m_tag); mMandatoryFields.insert(v.mTag);
v.m_category = this; v.mCategory = this;
auto r = m_item_validators.insert(move(v)); auto r = mItemValidators.insert(move(v));
if (not r.second and VERBOSE >= 4) if (not r.second and VERBOSE >= 4)
cout << "Could not add validator for item " << v.m_tag << " to category " << m_name << endl; cout << "Could not add validator for item " << v.mTag << " to category " << mName << endl;
} }
const validate_item* validate_category::get_validator_for_item(string tag) const const ValidateItem* ValidateCategory::getValidatorForItem(string tag) const
{ {
const validate_item* result = nullptr; const ValidateItem* result = nullptr;
auto i = m_item_validators.find(validate_item{tag}); auto i = mItemValidators.find(ValidateItem{tag});
if (i != m_item_validators.end()) if (i != mItemValidators.end())
result = &*i; result = &*i;
else if (VERBOSE > 4) else if (VERBOSE > 4)
cout << "No validator for tag " << tag << endl; cout << "No validator for tag " << tag << endl;
...@@ -180,61 +180,61 @@ const validate_item* validate_category::get_validator_for_item(string tag) const ...@@ -180,61 +180,61 @@ const validate_item* validate_category::get_validator_for_item(string tag) const
// -------------------------------------------------------------------- // --------------------------------------------------------------------
validator::validator() Validator::Validator()
{ {
} }
validator::~validator() Validator::~Validator()
{ {
} }
void validator::add_type_validator(validate_type&& v) void Validator::addTypeValidator(ValidateType&& v)
{ {
auto r = m_type_validators.insert(move(v)); auto r = mTypeValidators.insert(move(v));
if (not r.second and VERBOSE > 4) if (not r.second and VERBOSE > 4)
cout << "Could not add validator for type " << v.m_name << endl; cout << "Could not add validator for type " << v.mName << endl;
} }
const validate_type* validator::get_validator_for_type(string type_code) const const ValidateType* Validator::getValidatorForType(string typeCode) const
{ {
const validate_type* result = nullptr; const ValidateType* result = nullptr;
auto i = m_type_validators.find(validate_type{ type_code, ptChar, boost::regex() }); auto i = mTypeValidators.find(ValidateType{ typeCode, ptChar, boost::regex() });
if (i != m_type_validators.end()) if (i != mTypeValidators.end())
result = &*i; result = &*i;
else if (VERBOSE > 4) else if (VERBOSE > 4)
cout << "No validator for type " << type_code << endl; cout << "No validator for type " << typeCode << endl;
return result; return result;
} }
void validator::add_category_validator(validate_category&& v) void Validator::addCategoryValidator(ValidateCategory&& v)
{ {
auto r = m_category_validators.insert(move(v)); auto r = mCategoryValidators.insert(move(v));
if (not r.second and VERBOSE > 4) if (not r.second and VERBOSE > 4)
cout << "Could not add validator for category " << v.m_name << endl; cout << "Could not add validator for category " << v.mName << endl;
} }
const validate_category* validator::get_validator_for_category(string category) const const ValidateCategory* Validator::getValidatorForCategory(string category) const
{ {
const validate_category* result = nullptr; const ValidateCategory* result = nullptr;
auto i = m_category_validators.find(validate_category{category}); auto i = mCategoryValidators.find(ValidateCategory{category});
if (i != m_category_validators.end()) if (i != mCategoryValidators.end())
result = &*i; result = &*i;
else if (VERBOSE > 4) else if (VERBOSE > 4)
cout << "No validator for category " << category << endl; cout << "No validator for category " << category << endl;
return result; return result;
} }
validate_item* validator::get_validator_for_item(string tag) const ValidateItem* Validator::getValidatorForItem(string tag) const
{ {
validate_item* result = nullptr; ValidateItem* result = nullptr;
string cat, item; string cat, item;
std::tie(cat, item) = split_tag_name(tag); std::tie(cat, item) = splitTagName(tag);
auto* cv = get_validator_for_category(cat); auto* cv = getValidatorForCategory(cat);
if (cv != nullptr) if (cv != nullptr)
result = const_cast<validate_item*>(cv->get_validator_for_item(item)); result = const_cast<ValidateItem*>(cv->getValidatorForItem(item));
if (result == nullptr and VERBOSE > 4) if (result == nullptr and VERBOSE > 4)
cout << "No validator for item " << tag << endl; cout << "No validator for item " << tag << endl;
...@@ -242,10 +242,10 @@ validate_item* validator::get_validator_for_item(string tag) const ...@@ -242,10 +242,10 @@ validate_item* validator::get_validator_for_item(string tag) const
return result; return result;
} }
void validator::report_error(const string& msg) void Validator::reportError(const string& msg)
{ {
if (m_strict) if (mStrict)
throw validation_error(msg); throw ValidationError(msg);
else if (VERBOSE) else if (VERBOSE)
cerr << msg << endl; cerr << msg << endl;
} }
......
#include "libpr.h"
#include <set>
#include <map>
#include <unordered_set>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/fstream.hpp>
#include <boost/algorithm/string.hpp>
#include "cif++.h"
#include "peptidedb.h"
using namespace std;
namespace fs = boost::filesystem;
namespace ba = boost::algorithm;
const map<string,char> kAAMap{
{ "ALA", 'A' },
{ "ARG", 'R' },
{ "ASN", 'N' },
{ "ASP", 'D' },
{ "CYS", 'C' },
{ "GLN", 'Q' },
{ "GLU", 'E' },
{ "GLY", 'G' },
{ "HIS", 'H' },
{ "ILE", 'I' },
{ "LEU", 'L' },
{ "LYS", 'K' },
{ "MET", 'M' },
{ "PHE", 'F' },
{ "PRO", 'P' },
{ "SER", 'S' },
{ "THR", 'T' },
{ "TRP", 'W' },
{ "TYR", 'Y' },
{ "VAL", 'V' },
{ "GLX", 'Z' },
{ "ASX", 'B' }
};
const map<string,char> kBaseMap{
{ "A", 'A' },
{ "C", 'C' },
{ "G", 'G' },
{ "T", 'T' },
{ "U", 'U' },
{ "DA", 'A' },
{ "DC", 'C' },
{ "DG", 'G' },
{ "DT", 'T' }
};
// --------------------------------------------------------------------
struct PeptideDBImpl
{
PeptideDBImpl(istream& data, PeptideDBImpl* next);
~PeptideDBImpl()
{
delete m_next;
}
/*unordered_*/set<string> m_known_peptides;
set<string> m_known_bases;
cif::file m_file;
cif::category& m_chem_comp;
PeptideDBImpl* m_next;
string name_for(const string& res_name) const
{
string result;
for (auto& chem_comp: m_chem_comp)
{
if (ba::iequals(chem_comp["three_letter_code"].as<string>(), res_name) == false)
continue;
result = chem_comp["name"].as<string>();
ba::trim(result);
break;
}
if (result.empty() and m_next)
result = m_next->name_for(res_name);
return result;
}
string formula_for(string res_name) const;
string unalias(const string& res_name) const
{
string result = res_name;
auto& e = const_cast<cif::file&>(m_file)["comp_synonym_list"];
for (auto& synonym: e["chem_comp_synonyms"])
{
if (ba::iequals(synonym["comp_alternative_id"].as<string>(), res_name) == false)
continue;
result = synonym["comp_id"].as<string>();
ba::trim(result);
break;
}
if (result.empty() and m_next)
result = m_next->unalias(res_name);
return result;
}
};
PeptideDBImpl::PeptideDBImpl(istream& data, PeptideDBImpl* next)
: m_file(data), m_chem_comp(m_file.first_datablock()["chem_comp"]), m_next(next)
{
for (auto& chem_comp: m_chem_comp)
{
string group, three_letter_code;
cif::tie(group, three_letter_code) = chem_comp.get("group", "three_letter_code");
if (group == "peptide" or group == "M-peptide" or group == "P-peptide")
m_known_peptides.insert(three_letter_code);
else if (group == "DNA" or group == "RNA")
m_known_bases.insert(three_letter_code);
}
}
string PeptideDBImpl::formula_for(string res) const
{
string result;
ba::to_upper(res);
for (auto& db: m_file)
{
if (db.name() != "comp_" + res)
continue;
auto& cat = db["chem_comp_atom"];
map<string,uint32> atoms;
for (auto r: cat)
atoms[r["type_symbol"].as<string>()] += 1;
for (auto a: atoms)
{
if (not result.empty())
result += ' ';
result += a.first;
if (a.second > 1)
result += to_string(a.second);
}
}
if (result.empty())
{
if (m_next != nullptr)
result = m_next->formula_for(res);
else
{
const char* clibd_mon = getenv("CLIBD_MON");
if (clibd_mon == nullptr)
throw runtime_error("Cannot locate peptide list, please souce the CCP4 environment");
fs::path resFile = fs::path(clibd_mon) / ba::to_lower_copy(res.substr(0, 1)) / (res + ".cif");
if (fs::exists(resFile))
{
fs::ifstream file(resFile);
if (file.is_open())
{
try
{
cif::file cf(file);
auto& cat = cf["comp_" + res]["chem_comp_atom"];
map<string,uint32> atoms;
for (auto r: cat)
atoms[r["type_symbol"].as<string>()] += 1;
for (auto a: atoms)
{
if (not result.empty())
result += ' ';
result += a.first;
if (a.second > 1)
result += to_string(a.second);
}
}
catch (exception& ex)
{
if (VERBOSE)
cerr << ex.what();
result.clear();
}
}
}
}
}
return result;
}
// --------------------------------------------------------------------
PeptideDB* PeptideDB::sInstance;
PeptideDB& PeptideDB::Instance()
{
if (sInstance == nullptr)
sInstance = new PeptideDB();
return *sInstance;
}
PeptideDB::PeptideDB()
{
const char* clibd_mon = getenv("CLIBD_MON");
if (clibd_mon == nullptr)
throw runtime_error("Cannot locate peptide list, please souce the CCP4 environment");
fs::path db = fs::path(clibd_mon) / "list" / "mon_lib_list.cif";
PushDictionary(db);
sInstance = this;
}
void PeptideDB::PushDictionary(boost::filesystem::path dict)
{
if (not fs::exists(dict))
throw runtime_error("file not found: " + dict.string());
fs::ifstream file(dict);
if (not file.is_open())
throw runtime_error("Could not open peptide list " + dict.string());
mImpl = new PeptideDBImpl(file, mImpl);
}
void PeptideDB::PopDictionary()
{
if (mImpl != nullptr)
{
auto i = mImpl;
mImpl = i->m_next;
i->m_next = nullptr;
delete i;
}
}
PeptideDB::~PeptideDB()
{
delete mImpl;
}
bool PeptideDB::IsKnownPeptide(const string& res_name) const
{
return mImpl->m_known_peptides.count(res_name) > 0;
}
bool PeptideDB::IsKnownBase(const string& res_name) const
{
return mImpl->m_known_bases.count(res_name) > 0;
}
string PeptideDB::GetNameForResidue(const string& res_name) const
{
return mImpl->name_for(res_name);
}
string PeptideDB::GetFormulaForResidue(const string& res_name) const
{
return mImpl->formula_for(res_name);
}
string PeptideDB::Unalias(const string& res_name) const
{
return mImpl->unalias(res_name);
}
// Lib for working with structures as contained in mmCIF and PDB files // Lib for working with structures as contained in mmCIF and PDB files
#include "libcif/point.h" #include "cif++/Point.h"
using namespace std; using namespace std;
...@@ -32,16 +32,16 @@ quaternion Normalize(quaternion q) ...@@ -32,16 +32,16 @@ quaternion Normalize(quaternion q)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
float DihedralAngle(const point& p1, const point& p2, const point& p3, const point& p4) float DihedralAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4)
{ {
point v12 = p1 - p2; // vector from p2 to p1 Point v12 = p1 - p2; // vector from p2 to p1
point v43 = p4 - p3; // vector from p3 to p4 Point v43 = p4 - p3; // vector from p3 to p4
point z = p2 - p3; // vector from p3 to p2 Point z = p2 - p3; // vector from p3 to p2
point p = CrossProduct(z, v12); Point p = CrossProduct(z, v12);
point x = CrossProduct(z, v43); Point x = CrossProduct(z, v43);
point y = CrossProduct(z, x); Point y = CrossProduct(z, x);
double u = DotProduct(x, x); double u = DotProduct(x, x);
double v = DotProduct(y, y); double v = DotProduct(y, y);
...@@ -58,10 +58,10 @@ float DihedralAngle(const point& p1, const point& p2, const point& p3, const poi ...@@ -58,10 +58,10 @@ float DihedralAngle(const point& p1, const point& p2, const point& p3, const poi
return result; return result;
} }
float CosinusAngle(const point& p1, const point& p2, const point& p3, const point& p4) float CosinusAngle(const Point& p1, const Point& p2, const Point& p3, const Point& p4)
{ {
point v12 = p1 - p2; Point v12 = p1 - p2;
point v34 = p3 - p4; Point v34 = p3 - p4;
double result = 0; double result = 0;
...@@ -74,7 +74,7 @@ float CosinusAngle(const point& p1, const point& p2, const point& p3, const poin ...@@ -74,7 +74,7 @@ float CosinusAngle(const point& p1, const point& p2, const point& p3, const poin
// -------------------------------------------------------------------- // --------------------------------------------------------------------
tuple<double,point> QuaternionToAngleAxis(quaternion q) tuple<double,Point> QuaternionToAngleAxis(quaternion q)
{ {
if (q.R_component_1() > 1) if (q.R_component_1() > 1)
q = Normalize(q); q = Normalize(q);
...@@ -88,58 +88,58 @@ tuple<double,point> QuaternionToAngleAxis(quaternion q) ...@@ -88,58 +88,58 @@ tuple<double,point> QuaternionToAngleAxis(quaternion q)
if (s < 0.001) if (s < 0.001)
s = 1; s = 1;
point axis(q.R_component_2() / s, q.R_component_3() / s, q.R_component_4() / s); Point axis(q.R_component_2() / s, q.R_component_3() / s, q.R_component_4() / s);
return make_tuple(angle, axis); return make_tuple(angle, axis);
} }
point CenterPoints(vector<point>& points) Point CenterPoints(vector<Point>& Points)
{ {
point t; Point t;
for (point& pt : points) for (Point& pt : Points)
{ {
t.x() += pt.x(); t.getX() += pt.getX();
t.y() += pt.y(); t.getY() += pt.getY();
t.z() += pt.z(); t.getZ() += pt.getZ();
} }
t.x() /= points.size(); t.getX() /= Points.size();
t.y() /= points.size(); t.getY() /= Points.size();
t.z() /= points.size(); t.getZ() /= Points.size();
for (point& pt : points) for (Point& pt : Points)
{ {
pt.x() -= t.x(); pt.getX() -= t.getX();
pt.y() -= t.y(); pt.getY() -= t.getY();
pt.z() -= t.z(); pt.getZ() -= t.getZ();
} }
return t; return t;
} }
point Centroid(vector<point>& points) Point Centroid(vector<Point>& Points)
{ {
point result; Point result;
for (point& pt : points) for (Point& pt : Points)
result += pt; result += pt;
result /= points.size(); result /= Points.size();
return result; return result;
} }
double RMSd(const vector<point>& a, const vector<point>& b) double RMSd(const vector<Point>& a, const vector<Point>& b)
{ {
double sum = 0; double sum = 0;
for (uint32 i = 0; i < a.size(); ++i) for (uint32 i = 0; i < a.size(); ++i)
{ {
valarray<double> d(3); valarray<double> d(3);
d[0] = b[i].x() - a[i].x(); d[0] = b[i].getX() - a[i].getX();
d[1] = b[i].y() - a[i].y(); d[1] = b[i].getY() - a[i].getY();
d[2] = b[i].z() - a[i].z(); d[2] = b[i].getZ() - a[i].getZ();
d *= d; d *= d;
...@@ -188,19 +188,19 @@ double LargestDepressedQuarticSolution(double a, double b, double c) ...@@ -188,19 +188,19 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
return t.max(); return t.max();
} }
//quaternion AlignPoints(const vector<point>& pa, const vector<point>& pb) //quaternion AlignPoints(const vector<Point>& pa, const vector<Point>& pb)
//{ //{
// // First calculate M, a 3x3 matrix containing the sums of products of the coordinates of A and B // // First calculate M, a 3x3 matrix containing the sums of products of the coordinates of A and B
// matrix<double> M(3, 3, 0); // matrix<double> M(3, 3, 0);
// //
// for (uint32 i = 0; i < pa.size(); ++i) // for (uint32 i = 0; i < pa.size(); ++i)
// { // {
// const point& a = pa[i]; // const Point& a = pa[i];
// const point& b = pb[i]; // const Point& b = pb[i];
// //
// M(0, 0) += a.x() * b.x(); M(0, 1) += a.x() * b.y(); M(0, 2) += a.x() * b.z(); // M(0, 0) += a.getX() * b.getX(); M(0, 1) += a.getX() * b.getY(); M(0, 2) += a.getX() * b.getZ();
// M(1, 0) += a.y() * b.x(); M(1, 1) += a.y() * b.y(); M(1, 2) += a.y() * b.z(); // M(1, 0) += a.getY() * b.getX(); M(1, 1) += a.getY() * b.getY(); M(1, 2) += a.getY() * b.getZ();
// M(2, 0) += a.z() * b.x(); M(2, 1) += a.z() * b.y(); M(2, 2) += a.z() * b.z(); // M(2, 0) += a.getZ() * b.getX(); M(2, 1) += a.getZ() * b.getY(); M(2, 2) += a.getZ() * b.getZ();
// } // }
// //
// // Now calculate N, a symmetric 4x4 matrix // // Now calculate N, a symmetric 4x4 matrix
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment