Refactored dssp to be standalone

CMakeLists.txt

@@ -212,10 +212,11 @@ set(project_sources
 	# ${PROJECT_SOURCE_DIR}/src/pdb/PDB2Cif.cpp
 	# ${PROJECT_SOURCE_DIR}/src/pdb/PDB2CifRemark3.cpp
 
+	${PROJECT_SOURCE_DIR}/src/dssp/DSSP.cpp
+
 	${PROJECT_SOURCE_DIR}/src/structure/AtomType.cpp
 	# ${PROJECT_SOURCE_DIR}/src/structure/BondMap.cpp
 	${PROJECT_SOURCE_DIR}/src/structure/Compound.cpp
-	${PROJECT_SOURCE_DIR}/src/structure/DSSP.cpp
 	# ${PROJECT_SOURCE_DIR}/src/structure/Structure.cpp
 	${PROJECT_SOURCE_DIR}/src/structure/Symmetry.cpp
 	${PROJECT_SOURCE_DIR}/src/structure/TlsParser.cpp
@@ -242,12 +243,12 @@ set(project_headers
 	${PROJECT_SOURCE_DIR}/include/cif++/cif/condition.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/cif/category.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/cif/row.hpp
+	${PROJECT_SOURCE_DIR}/include/cif++/dssp/DSSP.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/structure/AtomType.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/structure/BondMap.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/structure/TlsParser.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/structure/Symmetry.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/structure/Structure.hpp
-	${PROJECT_SOURCE_DIR}/include/cif++/structure/DSSP.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/structure/Compound.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/pdb/PDB2Cif.hpp
 	${PROJECT_SOURCE_DIR}/include/cif++/pdb/PDB2CifRemark3.hpp

include/cif++/cif/item.hpp

@@ -227,6 +227,13 @@ struct item_handle
 		return item_value_as<T>::compare(*this, value, icase);
 	}
 
+	template <typename T>
+	bool operator==(const T &value) const
+	{
+		// TODO: icase or not icase?
+		return item_value_as<T>::compare(*this, value, true) == 0;
+	}
+
 	// empty means either null or unknown
 	bool empty() const
 	{

include/cif++/dssp/DSSP.hpp

+/*-
+ * SPDX-License-Identifier: BSD-2-Clause
+ *
+ * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright notice, this
+ *    list of conditions and the following disclaimer
+ * 2. Redistributions in binary form must reproduce the above copyright notice,
+ *    this list of conditions and the following disclaimer in the documentation
+ *    and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+ * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/// \file DSSP.hpp
+/// Calculate DSSP-like secondary structure information.
+
+#pragma once
+
+#include <cif++/cif.hpp>
+
+namespace dssp
+{
+
+struct residue;
+
+enum class structure_type : char
+{
+	Loop = ' ',
+	Alphahelix = 'H',
+	Betabridge = 'B',
+	Strand = 'E',
+	Helix_3 = 'G',
+	Helix_5 = 'I',
+	Helix_PPII = 'P',
+	Turn = 'T',
+	Bend = 'S'
+};
+
+enum class helix_type
+{
+	_3_10,
+	alpha,
+	pi,
+	pp
+};
+
+enum class helix_position_type
+{
+	None,
+	Start,
+	End,
+	StartAndEnd,
+	Middle
+};
+
+const size_t
+	kHistogramSize = 30;
+
+struct statistics
+{
+	struct
+	{
+		uint32_t residues, chains, SS_bridges, intra_chain_SS_bridges, H_bonds;
+		uint32_t H_bonds_in_antiparallel_bridges, H_bonds_in_parallel_bridges;
+		uint32_t H_Bonds_per_distance[11];
+	} count;
+
+	double accessible_surface;
+
+	struct
+	{
+		uint32_t residues_per_alpha_helix[kHistogramSize];
+		uint32_t parallel_bridges_per_ladder[kHistogramSize];
+		uint32_t antiparallel_bridges_per_ladder[kHistogramSize];
+		uint32_t ladders_per_sheet[kHistogramSize];
+	} histogram;
+};
+
+enum class chain_break_type
+{
+	None,
+	NewChain,
+	Gap
+};
+
+class DSSP
+{
+  public:
+	DSSP(const cif::datablock &db, int model_nr, int min_poly_proline_stretch_length, bool calculateSurfaceAccessibility);
+	~DSSP();
+
+	DSSP(const DSSP &) = delete;
+	DSSP &operator=(const DSSP &) = delete;
+
+	statistics get_statistics() const;
+
+	class iterator;
+	using res_iterator = typename std::vector<residue>::iterator;
+
+	class residue_info
+	{
+	  public:
+		friend class iterator;
+
+		residue_info() = default;
+		residue_info(const residue_info &rhs) = default;
+		residue_info &operator=(const residue_info &rhs) = default;
+
+		explicit operator bool() const { return not empty(); }
+		bool empty() const { return m_impl == nullptr; }
+
+		std::string asym_id() const;
+		int seq_id() const;
+		std::string alt_id() const;
+		std::string compound_id() const;
+
+		std::string auth_asym_id() const;
+		int auth_seq_id() const;
+
+		std::string pdb_strand_id() const;
+		int pdb_seq_num() const;
+		std::string pdb_ins_code() const;
+
+		float alpha() const;
+		float kappa() const;
+		float phi() const;
+		float psi() const;
+		float tco() const;
+		std::tuple<float, float, float> ca_location() const;
+
+		chain_break_type chain_break() const;
+
+		/// \brief the internal number in DSSP
+		int nr() const;
+
+		structure_type type() const;
+
+		int ssBridgeNr() const;
+
+		helix_position_type helix(helix_type helixType) const;
+
+		bool is_alpha_helix_end_before_start() const;
+
+		bool bend() const;
+
+		double accessibility() const;
+
+		/// \brief returns resinfo, ladder and parallel
+		std::tuple<residue_info, int, bool> bridge_partner(int i) const;
+
+		int sheet() const;
+
+		/// \brief return resinfo and the energy of the bond
+		std::tuple<residue_info, double> acceptor(int i) const;
+		std::tuple<residue_info, double> donor(int i) const;
+
+		/// \brief Simple compare equals
+		bool operator==(const residue_info &rhs) const
+		{
+			return m_impl == rhs.m_impl;
+		}
+
+		/// \brief Returns \result true if there is a bond between two residues
+		friend bool test_bond(residue_info const &a, residue_info const &b);
+
+	  private:
+		residue_info(residue *res)
+			: m_impl(res)
+		{
+		}
+
+		residue *m_impl = nullptr;
+	};
+
+	class iterator
+	{
+	  public:
+		using iterator_category = std::bidirectional_iterator_tag;
+		using value_type = residue_info;
+		using difference_type = std::ptrdiff_t;
+		using pointer = value_type *;
+		using reference = value_type &;
+
+		iterator(const iterator &i) = default;
+		iterator(residue *res);
+		iterator &operator=(const iterator &i) = default;
+
+		reference operator*() { return m_current; }
+		pointer operator->() { return &m_current; }
+
+		iterator &operator++();
+		iterator operator++(int)
+		{
+			auto tmp(*this);
+			this->operator++();
+			return tmp;
+		}
+
+		iterator &operator--();
+		iterator operator--(int)
+		{
+			auto tmp(*this);
+			this->operator--();
+			return tmp;
+		}
+
+		bool operator==(const iterator &rhs) const { return m_current.m_impl == rhs.m_current.m_impl; }
+		bool operator!=(const iterator &rhs) const { return m_current.m_impl != rhs.m_current.m_impl; }
+
+	  private:
+		residue_info m_current;
+	};
+
+	using value_type = residue_info;
+
+	// To access residue info by key, i.e. LabelAsymID and LabelSeqID
+	using key_type = std::tuple<std::string,int>;
+
+	iterator begin() const;
+	iterator end() const;
+
+	residue_info operator[](const key_type &key) const;
+
+	bool empty() const { return begin() == end(); }
+
+	// convenience method, when creating old style DSSP files
+
+	enum class pdb_record_type { HEADER, COMPND, SOURCE, AUTHOR };
+
+	std::string get_pdb_header_line(pdb_record_type pdb_record) const;
+
+  private:
+	struct DSSP_impl *m_impl;
+};
+
+} // namespace dssp

include/cif++/structure/DSSP.hpp

-/*-
- * SPDX-License-Identifier: BSD-2-Clause
- *
- * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- *
- * 1. Redistributions of source code must retain the above copyright notice, this
- *    list of conditions and the following disclaimer
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- *    this list of conditions and the following disclaimer in the documentation
- *    and/or other materials provided with the distribution.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
- * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
- * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
- * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
- * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
- * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
- * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
- * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
- * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
-// Calculate DSSP-like secondary structure information
-
-#pragma once
-
-#include <cif++/cif.hpp>
-
-namespace mmcif
-{
-
-struct Res;
-
-extern const float
-	kCouplingConstant,
-	kMinHBondEnergy, kMaxHBondEnergy;
-
-enum SecondaryStructureType : char
-{
-	ssLoop = ' ',
-	ssAlphahelix = 'H',
-	ssBetabridge = 'B',
-	ssStrand = 'E',
-	ssHelix_3 = 'G',
-	ssHelix_5 = 'I',
-	ssHelix_PPII = 'P',
-	ssTurn = 'T',
-	ssBend = 'S'
-};
-
-enum class HelixType
-{
-	rh_3_10,
-	rh_alpha,
-	rh_pi,
-	rh_pp
-};
-
-enum class Helix
-{
-	None,
-	Start,
-	End,
-	StartAndEnd,
-	Middle
-};
-
-// struct HBond
-//{
-//	std::string 				labelAsymID;
-//	int							labelSeqID;
-//	double						energy;
-// };
-//
-// struct BridgePartner
-//{
-//	std::string					labelAsymID;
-//	int							labelSeqID;
-//	int							ladder;
-//	bool						parallel;
-// };
-
-struct SecondaryStructure
-{
-	SecondaryStructureType type;
-	//	HBond						donor[2], acceptor[2];
-	//	BridgePartner				beta[2];
-	//	int							sheet;
-	//	bool						bend;
-};
-
-// void CalculateSecondaryStructure(Structure& s);
-
-const size_t
-	kHistogramSize = 30;
-
-struct DSSP_Statistics
-{
-	uint32_t nrOfResidues, nrOfChains, nrOfSSBridges, nrOfIntraChainSSBridges, nrOfHBonds;
-	uint32_t nrOfHBondsInAntiparallelBridges, nrOfHBondsInParallelBridges;
-	uint32_t nrOfHBondsPerDistance[11] = {};
-	double accessibleSurface = 0;
-
-	uint32_t residuesPerAlphaHelixHistogram[kHistogramSize] = {};
-	uint32_t parallelBridgesPerLadderHistogram[kHistogramSize] = {};
-	uint32_t antiparallelBridgesPerLadderHistogram[kHistogramSize] = {};
-	uint32_t laddersPerSheetHistogram[kHistogramSize] = {};
-};
-
-enum class ChainBreak
-{
-	None,
-	NewChain,
-	Gap
-};
-
-class DSSP
-{
-  public:
-	DSSP(const cif::datablock &db, int model_nr, int min_poly_proline_stretch_length, bool calculateSurfaceAccessibility);
-	~DSSP();
-
-	DSSP(const DSSP &) = delete;
-	DSSP &operator=(const DSSP &) = delete;
-
-	SecondaryStructureType operator()(const std::string &inAsymID, int inSeqID) const;
-	double accessibility(const std::string &inAsymID, int inSeqID) const;
-	bool isAlphaHelixEndBeforeStart(const std::string &inAsymID, int inSeqID) const;
-
-	DSSP_Statistics GetStatistics() const;
-
-	class iterator;
-	using res_iterator = typename std::vector<Res>::iterator;
-
-	class ResidueInfo
-	{
-	  public:
-		friend class iterator;
-
-		ResidueInfo()
-			: mImpl(nullptr)
-		{
-		}
-
-		ResidueInfo(const ResidueInfo &rhs)
-			: mImpl(rhs.mImpl)
-		{
-		}
-
-		ResidueInfo &operator=(const ResidueInfo &rhs)
-		{
-			mImpl = rhs.mImpl;
-			return *this;
-		}
-
-		explicit operator bool() const { return not empty(); }
-		bool empty() const { return mImpl == nullptr; }
-
-		std::string asym_id() const;
-		int seq_id() const;
-		std::string alt_id() const;
-		std::string compound_id() const;
-
-		std::string auth_asym_id() const;
-		int auth_seq_id() const;
-
-		std::string pdb_strand_id() const;
-		int pdb_seq_num() const;
-		std::string pdb_ins_code() const;
-
-		float alpha() const;
-		float kappa() const;
-		float phi() const;
-		float psi() const;
-		float tco() const;
-		std::tuple<float,float,float> ca_location() const;
-
-		/// \brief return 0 if not a break, ' ' in case of a new chain and '*' in case of a broken chain
-		ChainBreak chainBreak() const;
-
-		/// \brief the internal number in DSSP
-		int nr() const;
-
-		SecondaryStructureType ss() const;
-
-		int ssBridgeNr() const;
-
-		Helix helix(HelixType helixType) const;
-
-		bool bend() const;
-
-		double accessibility() const;
-
-		/// \brief returns resinfo, ladder and parallel
-		std::tuple<ResidueInfo, int, bool> bridgePartner(int i) const;
-
-		int sheet() const;
-
-		/// \brief return resinfo and the energy of the bond
-		std::tuple<ResidueInfo, double> acceptor(int i) const;
-		std::tuple<ResidueInfo, double> donor(int i) const;
-
-		/// \brief Simple compare equals
-		bool operator==(const ResidueInfo &rhs) const
-		{
-			return mImpl == rhs.mImpl;
-		}
-
-		/// \brief Returns \result true if there is a bond between two residues
-		friend bool TestBond(ResidueInfo const &a, ResidueInfo const &b);
-
-	  private:
-		ResidueInfo(Res *res)
-			: mImpl(res)
-		{
-		}
-
-		Res *mImpl;
-	};
-
-	class iterator
-	{
-	  public:
-		using iterator_category = std::bidirectional_iterator_tag;
-		using value_type = ResidueInfo;
-		using difference_type = std::ptrdiff_t;
-		using pointer = value_type *;
-		using reference = value_type &;
-
-		iterator(const iterator &i);
-		iterator(Res *res);
-		iterator &operator=(const iterator &i);
-
-		reference operator*() { return mCurrent; }
-		pointer operator->() { return &mCurrent; }
-
-		iterator &operator++();
-		iterator operator++(int)
-		{
-			auto tmp(*this);
-			this->operator++();
-			return tmp;
-		}
-
-		iterator &operator--();
-		iterator operator--(int)
-		{
-			auto tmp(*this);
-			this->operator--();
-			return tmp;
-		}
-
-		bool operator==(const iterator &rhs) const { return mCurrent.mImpl == rhs.mCurrent.mImpl; }
-		bool operator!=(const iterator &rhs) const { return mCurrent.mImpl != rhs.mCurrent.mImpl; }
-
-	  private:
-		ResidueInfo mCurrent;
-	};
-
-	iterator begin() const;
-	iterator end() const;
-
-	bool empty() const { return begin() == end(); }
-
-  private:
-	struct DSSPImpl *mImpl;
-};
-
-} // namespace mmcif

test/dssp-test.cpp

@@ -31,7 +31,7 @@
 #include <stdexcept>
 
 #include <cif++/cif.hpp>
-#include <cif++/structure/DSSP.hpp>
+#include <cif++/dssp/DSSP.hpp>
 
 namespace tt = boost::test_tools;
 
@@ -85,7 +85,7 @@ BOOST_AUTO_TEST_CASE(dssp_1)
 
 	std::ifstream t(gTestDir / "1cbs-dssp-test.tsv");
 
-	mmcif::DSSP dssp(f.front(), 1, 3, true);
+	dssp::DSSP dssp(f.front(), 1, 3, true);
 
 	for (auto residue : dssp)
 	{
@@ -109,6 +109,6 @@ BOOST_AUTO_TEST_CASE(dssp_1)
 
 		BOOST_CHECK_EQUAL(residue.asym_id(), asymID);
 		BOOST_CHECK_EQUAL(residue.seq_id(), seqID);
-		BOOST_CHECK_EQUAL((char)residue.ss(), secstr.front());
+		BOOST_CHECK_EQUAL((char)residue.type(), secstr.front());
 	}
 }
\ No newline at end of file