chiron started

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@215 a1961a4f-ab94-4bcc-80e8-33b5a54de466

include/cif++/Cif++.h

@@ -597,6 +597,24 @@ struct Condition
 namespace detail
 {
 
+struct KeyIsEmptyConditionImpl : public ConditionImpl
+{
+	KeyIsEmptyConditionImpl(const string& ItemTag)
+		: mItemTag(ItemTag) {}
+	
+	virtual bool test(const Category& c, const Row& r) const
+	{
+		return r[mItemTag].empty();
+	}
+	
+	virtual std::string str() const
+	{
+		return mItemTag + " == <empty>";
+	}
+	
+	string mItemTag;
+};
+
 template<typename T>
 struct KeyIsConditionImpl : public ConditionImpl
 {
@@ -782,7 +800,9 @@ inline Condition operator||(Condition&& a, Condition&& b)
 {
 	return Condition(new detail::orConditionImpl(std::move(a), std::move(b)));
 }
-	
+
+struct Empty {};
+
 struct Key
 {
 	Key(const string& ItemTag) : mItemTag(ItemTag) {}
@@ -800,6 +820,11 @@ struct Key
 		return Condition(new detail::KeyIsConditionImpl<std::string>(mItemTag, value));
 	}
 	
+	Condition operator==(Empty&) const
+	{
+		return Condition(new detail::KeyIsEmptyConditionImpl(mItemTag));
+	}
+	
 	template<typename T>
 	Condition operator!=(const T& v) const
 	{

include/cif++/Compound.h

@@ -42,15 +42,21 @@ struct CompoundAtom
 // This information is derived from the ccp4 monomer library by default.
 // To create compounds, you'd best use the factory method.
 
+enum ChiralVolumeSign { negativ, positiv, both };
+
 class Compound
 {
   public:
 
+	struct ChiralCentre;
+
 	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)
+		std::map<std::tuple<std::string,std::string>,float>&& bonds,
+		std::vector<ChiralCentre>&& chiralCentres)
 		: mId(id), mName(name), mGroup(group)
 		, mAtoms(std::move(atoms)), mBonds(std::move(bonds))
+		, mChiralCentres(std::move(chiralCentres))
 	{
 	}
 
@@ -88,13 +94,24 @@ class Compound
 	int charge() const;
 	bool isWater() const;
 
+	struct ChiralCentre
+	{
+		std::string			mID;
+		std::string			mAtomIDCentre;
+		std::string			mAtomID[3];
+		ChiralVolumeSign	mVolumeSign;
+	};
+
+	std::vector<ChiralCentre> chiralCentres() const			{ return mChiralCentres; }
+
   private:
-//	Entity&					mEntity;
-	std::string				mId;
-	std::string				mName;
-	std::string				mGroup;
+//	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;
+	std::vector<ChiralCentre>	mChiralCentres;
 };
 
 // --------------------------------------------------------------------

include/cif++/Structure.h

@@ -12,6 +12,7 @@
 #include "cif++/AtomType.h"
 #include "cif++/Point.h"
 #include "cif++/Compound.h"
+#include "cif++/Cif++.h"
 
 /*
 	To modify a structure, you will have to use actions.
@@ -35,12 +36,14 @@
 	
 */
 
-// forward declaration
-namespace cif
-{
-	class Datablock;
-	class File;
-};
+//// forward declaration
+//namespace cif
+//{
+//	class Category;
+//	class Datablock;
+//	class File;
+//	class RowSet;
+//};
 
 namespace libcif
 {
@@ -137,77 +140,121 @@ typedef std::vector<Atom> AtomView;
 
 // --------------------------------------------------------------------
 
-class Residue : public std::enable_shared_from_this<Residue>
+class Residue
 {
   public:
-	Residue(const Compound& cmp) : mCompound(cmp) {}
+	Residue(const Residue& rhs);
+	Residue& operator=(const Residue& rhs);
 
-	const Compound&		comp() const		{ return mCompound; }
-	virtual AtomView	atoms();
+	Residue(const Structure& structure)
+		: mStructure(&structure) {}
+	Residue(const Structure& structure, const std::string& compoundID,
+		const std::string& asymID, const std::string& seqID = "",
+		const std::string& altID = "")
+		: mStructure(&structure), mCompoundID(compoundID)
+		, mAsymID(asymID), mSeqID(seqID), mAltID(altID) {}
+	
+	const Compound&		comp() const;
+	AtomView			atoms() const;
+	
+	Atom				atomByID(const std::string& atomID) const;
 
-  private:
-	const Compound&		mCompound;
+	const std::string&	compoundID() const	{ return mCompoundID; }
+	const std::string&	asymID() const		{ return mAsymID; }
+	const std::string&	seqID() const		{ return mSeqID; }
+	const std::string&	altID() const		{ return mAltID; }
+
+  protected:
+
+	const Structure* mStructure;
+	std::string	mCompoundID, mAsymID, mSeqID, mAltID;
+	mutable AtomView mAtoms;
 };
 
-//// --------------------------------------------------------------------
-//// 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);
+// --------------------------------------------------------------------
+// a monomer models a single Residue in a protein chain 
+
+class Monomer : public Residue
+{
+  public:
+	Monomer(const Monomer& rhs);
+	Monomer& operator=(const Monomer& rhs);
+
+	Monomer(Polymer& polymer);
+	Monomer(Polymer& polymer, uint32 seqID,
+		const std::string& compoundID, const std::string& altID);
+
+  private:
+	Polymer*	mPolymer;
+};
+
+// --------------------------------------------------------------------
+
+class Polymer
+{
+  public:
+	Polymer(const Structure& s, const std::string& entityID, const std::string& asymID);
+	Polymer(const Polymer&);
+	Polymer& operator=(const Polymer&);
+	
+	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& p, 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;
-//};
+		
+		reference	operator*()											{ return mCurrent; }
+		pointer		operator->()										{ return &mCurrent; }
+		
+		iterator&	operator++();
+		iterator	operator++(int)
+		{
+			iterator result(*this);
+			operator++();
+			return result;
+		}
+
+		iterator&	operator--();
+		iterator	operator--(int)
+		{
+			iterator result(*this);
+			operator--();
+			return result;
+		}
+
+		bool		operator==(const iterator& rhs) const				{ return mPolymer == rhs.mPolymer and mIndex == rhs.mIndex; }
+		bool		operator!=(const iterator& rhs) const				{ return mPolymer != rhs.mPolymer or mIndex != rhs.mIndex; }
+
+	  private:
+		Polymer*	mPolymer;
+		uint32		mIndex;
+		Monomer		mCurrent;
+	};
+
+	iterator begin();
+	iterator end();
+
+	Structure* structure() const	{ return mStructure; }
+	
+	std::string asymID() const		{ return mAsymID; }
+	std::string entityID() const	{ return mEntityID; }
+
+  private:
+
+	friend struct iterator;
+
+	Structure*					mStructure;
+	std::string					mEntityID;
+	std::string					mAsymID;
+	cif::RowSet					mPolySeq;
+};
 
 // --------------------------------------------------------------------
 // file is a reference to the data stored in e.g. the cif file.
@@ -254,6 +301,9 @@ class Structure
 
 	AtomView atoms() const;
 	AtomView waters() const;
+	
+	std::vector<Polymer> polymers() const;
+	std::vector<Residue> nonPolymers() const;
 
 	Atom getAtomById(std::string id) const;
 	Atom getAtomByLocation(Point pt, float maxDistance) const;
@@ -286,8 +336,13 @@ class Structure
 	
 	// Actions
 	void removeAtom(Atom& a);
-
+	
   private:
+	friend Polymer;
+	friend Residue;
+
+	cif::Category& category(const char* name) const;
+
 	struct StructureImpl*	mImpl;
 };
 

src/Compound.cpp

@@ -271,8 +271,38 @@ const Compound* CompoundFactory::create(std::string id)
 				
 				bonds[make_tuple(atomId_1, atomId_2)] = value;
 			}
+
+			auto& compChir = cf["comp_" + id]["chem_comp_chir"];
 			
-			result = new Compound(id, name, group, move(atoms), move(bonds));
+			vector<Compound::ChiralCentre> chiralCentres;
+			for (auto row: compChir)
+			{
+				Compound::ChiralCentre cc;
+				string volumeSign;
+				
+				cif::tie(cc.mID, cc.mAtomIDCentre, cc.mAtomID[0],
+					cc.mAtomID[1], cc.mAtomID[2], volumeSign) = 
+					row.get("id", "atom_id_centre", "atom_id_1",
+						"atom_id_2", "atom_id_3", "volume_sign");
+				
+				if (volumeSign == "negativ")
+					cc.mVolumeSign = negativ;
+				else if (volumeSign == "positiv")
+					cc.mVolumeSign = positiv;
+				else if (volumeSign == "both")
+					cc.mVolumeSign = both;
+				else
+				{
+					if (VERBOSE)
+						cerr << "Unimplemented chem_comp_chir.volume_sign " << volumeSign << " in file " << resFile << endl;
+					continue;
+				}
+				
+				chiralCentres.push_back(cc);
+			}
+			
+			result = new Compound(id, name, group, move(atoms), move(bonds),
+				move(chiralCentres));
 			mCompounds.push_back(result);
 		}
 	}

src/Structure.cpp

@@ -451,17 +451,200 @@ float Atom::radius() const
 // --------------------------------------------------------------------
 // residue
 
-//AtomView residue::Atoms()
-//{
-//	assert(false);
-//}
+Residue::Residue(const Residue& rhs)
+	: mStructure(rhs.mStructure)
+	, mCompoundID(rhs.mCompoundID), mAsymID(rhs.mAsymID), mSeqID(rhs.mSeqID), mAltID(rhs.mAltID)
+{
+}
+
+Residue& Residue::operator=(const Residue& rhs)
+{
+	if (this != &rhs)
+	{
+		mStructure = rhs.mStructure;
+		mCompoundID = rhs.mCompoundID;
+		mAsymID = rhs.mAsymID;
+		mSeqID = rhs.mSeqID;
+		mAltID = rhs.mAltID;
+	}
+	
+	return *this;
+}
+
+const Compound& Residue::comp() const
+{
+	auto compound = Compound::create(mCompoundID);
+	if (compound == nullptr)
+		throw runtime_error("Failed to create compound " + mCompoundID);
+	return *compound;
+}
+
+AtomView Residue::atoms() const
+{
+	if (mStructure == nullptr)
+		throw runtime_error("Invalid Residue object");
+
+	if (mAtoms.empty())
+	{
+		for (auto& a: mStructure->atoms())
+		{
+			if (a.property<string>("label_asym_id") != mAsymID or
+				a.property<string>("label_comp_id") != mCompoundID)
+					continue;
+			
+			if (not mSeqID.empty() and a.property<string>("label_seq_id") != mSeqID)
+				continue;
+			
+			if (not mAltID.empty() and a.property<string>("label_alt_id") != mAltID)
+				continue;
+			
+			mAtoms.push_back(a);
+		}
+	}
+	
+	return mAtoms;
+	
+//	auto& atomSites = mStructure->category("atom_site");
+//	
+//	auto query = cif::Key("label_asym_id") == mAsymID and cif::Key("label_comp_id") == mCompoundID;
+//	
+//	if (not mSeqID.empty())
+//		query = move(query) and cif::Key("label_seq_id") == mSeqID;
+//	
+//	if (not mAltID.empty())
+//		query = move(query) and (cif::Key("label_alt_id") == cif::Empty() or cif::Key("label_alt_id") == mAltID);
+//
+//	return atomSites.find(move(query));
+}
+
+Atom Residue::atomByID(const string& atomID) const
+{
+	for (auto& a: atoms())
+	{
+		if (a.property<string>("label_atom_id") == atomID)
+			return a;
+	}
+	
+	throw runtime_error("Atom with atom_id " + atomID + " not found in residue " + mAsymID + ':' + mSeqID);
+}
 
 // --------------------------------------------------------------------
 // monomer
 
+Monomer::Monomer(const Monomer& rhs)
+	: Residue(rhs), mPolymer(rhs.mPolymer)
+{
+}
+
+Monomer& Monomer::operator=(const Monomer& rhs)
+{
+	if (this != &rhs)
+	{
+		Residue::operator=(rhs);
+		mPolymer = rhs.mPolymer;
+	}
+	
+	return *this;
+}
+
+Monomer::Monomer(Polymer& polymer)
+	: Residue(*polymer.structure())
+	, mPolymer(&polymer)
+{
+	
+}
+
+Monomer::Monomer(Polymer& polymer, uint32 seqID, const string& compoundID, const string& altID)
+	: Residue(*polymer.structure(), compoundID, polymer.asymID(), to_string(seqID), altID)
+	, mPolymer(&polymer)
+{
+}
+
 // --------------------------------------------------------------------
 // polymer
 
+Polymer::iterator::iterator(Polymer& p, uint32 index)
+	: mPolymer(&p), mIndex(index), mCurrent(p)
+{
+	auto& polySeq = mPolymer->mPolySeq;
+
+	if (index < polySeq.size())
+	{
+		int seqID;
+		string asymID, monID;
+		cif::tie(asymID, seqID, monID) =
+			polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
+		
+		mCurrent = Monomer(*mPolymer, seqID, monID, "");
+	}
+}
+
+Polymer::iterator::iterator(const iterator& rhs)
+	: mPolymer(rhs.mPolymer), mIndex(rhs.mIndex), mCurrent(rhs.mCurrent)
+{
+}
+
+Polymer::iterator& Polymer::iterator::operator=(const iterator& rhs)
+{
+	if (this != &rhs)
+	{
+		mPolymer = rhs.mPolymer;
+		mIndex = rhs.mIndex;
+		mCurrent = rhs.mCurrent;
+	}
+	
+	return *this;
+}
+
+Polymer::iterator& Polymer::iterator::operator++()
+{
+	auto& polySeq = mPolymer->mPolySeq;
+
+	if (mIndex < polySeq.size())
+		++mIndex;
+
+	if (mIndex < polySeq.size())
+	{
+		int seqID;
+		string asymID, monID;
+		cif::tie(asymID, seqID, monID) =
+			polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
+		
+		mCurrent = Monomer(*mPolymer, seqID, monID, "");
+	}
+	
+	return *this;
+}
+
+Polymer::Polymer(const Polymer& rhs)
+	: mStructure(rhs.mStructure), mEntityID(rhs.mEntityID), mAsymID(rhs.mAsymID), mPolySeq(rhs.mPolySeq)
+{
+	
+}
+
+Polymer::Polymer(const Structure& s, const string& entityID, const string& asymID)
+	: mStructure(const_cast<Structure*>(&s)), mEntityID(entityID), mAsymID(asymID)
+	, mPolySeq(s.category("pdbx_poly_seq_scheme").find(cif::Key("asym_id") == mAsymID and cif::Key("entity_id") == mEntityID))
+{
+}
+
+Polymer::iterator Polymer::begin()
+{
+	return iterator(*this, 0);
+}
+
+Polymer::iterator Polymer::end()
+{
+	return iterator(*this, mPolySeq.size());
+}
+
+//cif::RowSet Polymer::polySeqRows() const
+//{
+//	auto& cat = mStructure->category("pdbx_poly_seq_scheme");
+//	
+//	return cat.find(cif::Key("asym_id") == mAsymID and cif::Key("entityID") == mEntityID);
+//}
+//
 // --------------------------------------------------------------------
 // File
 
@@ -655,6 +838,44 @@ AtomView Structure::waters() const
 	return result;
 }
 
+vector<Polymer> Structure::polymers() const
+{
+	vector<Polymer> result;
+	
+	auto& polySeqScheme = category("pdbx_poly_seq_scheme");
+	
+	for (auto& r: polySeqScheme)
+	{
+		string asymID, entityID, seqID, monID;
+		cif::tie(asymID, entityID, seqID, monID) =
+			r.get("asym_id", "entity_id", "seq_id", "mon_id");
+		
+		if (result.empty() or result.back().asymID() != asymID or result.back().entityID() != entityID)
+			result.emplace_back(*this, entityID, asymID);
+	}
+	
+	return result;
+}
+
+vector<Residue> Structure::nonPolymers() const
+{
+	vector<Residue> result;
+
+	auto& nonPolyScheme = category("pdbx_nonpoly_scheme");
+	
+	for (auto& r: nonPolyScheme)
+	{
+		string asymID, monID;
+		cif::tie(asymID, monID) =
+			r.get("asym_id", "mon_id");
+		
+		if (result.empty() or result.back().asymID() != asymID)
+			result.emplace_back(*this, monID, asymID);
+	}
+	
+	return result;
+}
+
 Atom Structure::getAtomById(string id) const
 {
 	for (auto& a: mImpl->mAtoms)
@@ -671,6 +892,12 @@ File& Structure::getFile() const
 	return *mImpl->mFile;
 }
 
+cif::Category& Structure::category(const char* name) const
+{
+	auto& db = *getFile().impl().mDb;
+	return db[name];
+}
+
 //tuple<string,string> Structure::MapLabelToAuth(
 //	const string& asymId, int seqId)
 //{