--- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -116,11 +116,9 @@ find_or_fetch_package(mcfp VERSION 1.3.5 GIT_REPOSITORY https://github.com/mhekkel/libmcfp GIT_TAG b6c62a3) find_or_fetch_package(cifpp VERSION 9 GIT_REPOSITORY https://github.com/PDB-REDO/libcifpp.git GIT_TAG v9.0.5 VARIABLES "CIFPP_SHARE_DIR") -if (NOT TARGET dssp) - add_subdirectory(dssp EXCLUDE_FROM_ALL) -endif() - add_executable(tortoize + ${PROJECT_SOURCE_DIR}/src/dssp.cpp + ${PROJECT_SOURCE_DIR}/src/dssp-io.cpp ${PROJECT_SOURCE_DIR}/src/tortoize.cpp ${PROJECT_SOURCE_DIR}/src/tortoize-main.cpp) @@ -137,7 +135,7 @@ target_include_directories(tortoize PRIVATE ${PROJECT_BINARY_DIR}) -target_link_libraries(tortoize dssp::dssp mcfp::mcfp cifpp::cifpp) +target_link_libraries(tortoize mcfp::mcfp cifpp::cifpp) install(TARGETS tortoize RUNTIME DESTINATION ${BIN_INSTALL_DIR} @@ -160,7 +158,7 @@ ${PROJECT_SOURCE_DIR}/src/tortoize-server-main.cpp ) - target_link_libraries(tortoized dssp::dssp mcfp::mcfp zeep::zeep) + target_link_libraries(tortoized mcfp::mcfp zeep::zeep) mrc_target_resources(tortoized ${RESOURCES} ${CMAKE_CURRENT_SOURCE_DIR}/docroot/) --- /dev/null +++ b/src/dssp.cpp @@ -0,0 +1,2298 @@ +/*- + * 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 + +#include "dssp.hpp" + +#include "dssp-io.hpp" + +#include +#include +#include +#include + +#ifdef near +# undef near +#endif + +using residue = dssp::residue; +using statistics = dssp::statistics; +using structure_type = dssp::structure_type; +using helix_type = dssp::helix_type; +using helix_position_type = dssp::helix_position_type; +using chain_break_type = dssp::chain_break_type; + +// -------------------------------------------------------------------- + +const double + kPI = 3.141592653589793238462643383279502884; + +struct point +{ + float mX, mY, mZ; +}; + +inline point operator-(const point &lhs, const point &rhs) +{ + return { lhs.mX - rhs.mX, lhs.mY - rhs.mY, lhs.mZ - rhs.mZ }; +} + +inline point operator*(const point &lhs, float rhs) +{ + return { lhs.mX * rhs, lhs.mY * rhs, lhs.mZ * rhs }; +} + +inline float distance_sq(const point &a, const point &b) +{ + return (a.mX - b.mX) * (a.mX - b.mX) + + (a.mY - b.mY) * (a.mY - b.mY) + + (a.mZ - b.mZ) * (a.mZ - b.mZ); +} + +inline float distance(const point &a, const point &b) +{ + return std::sqrt(distance_sq(a, b)); +} + +inline float dot_product(const point &a, const point &b) +{ + return a.mX * b.mX + a.mY * b.mY + a.mZ * b.mZ; +} + +inline point cross_product(const point &a, const point &b) +{ + return { + a.mY * b.mZ - b.mY * a.mZ, + a.mZ * b.mX - b.mZ * a.mX, + a.mX * b.mY - b.mX * a.mY + }; +} + +float dihedral_angle(const point &p1, const point &p2, const point &p3, const point &p4) +{ + point v12 = p1 - p2; // vector from p2 to p1 + point v43 = p4 - p3; // vector from p3 to p4 + + point z = p2 - p3; // vector from p3 to p2 + + point p = cross_product(z, v12); + point x = cross_product(z, v43); + point y = cross_product(z, x); + + float u = dot_product(x, x); + float v = dot_product(y, y); + + float result = 360; + if (u > 0 and v > 0) + { + u = dot_product(p, x) / std::sqrt(u); + v = dot_product(p, y) / std::sqrt(v); + if (u != 0 or v != 0) + result = std::atan2(v, u) * 180.f / static_cast(kPI); + } + + return result; +} + +float cosinus_angle(const point &p1, const point &p2, const point &p3, const point &p4) +{ + point v12 = p1 - p2; + point v34 = p3 - p4; + + float result = 0; + + float x = dot_product(v12, v12) * dot_product(v34, v34); + if (x > 0) + result = dot_product(v12, v34) / std::sqrt(x); + + return result; +} + +enum residue_type : char +{ + kUnknownResidue = 'X', + + // + kAlanine = 'A', // ala + kArginine = 'R', // arg + kAsparagine = 'N', // asn + kAsparticAcid = 'D', // asp + kCysteine = 'C', // cys + kGlutamicAcid = 'E', // glu + kGlutamine = 'Q', // gln + kGlycine = 'G', // gly + kHistidine = 'H', // his + kIsoleucine = 'I', // ile + kLeucine = 'L', // leu + kLysine = 'K', // lys + kMethionine = 'M', // met + kPhenylalanine = 'F', // phe + kProline = 'P', // pro + kSerine = 'S', // ser + kThreonine = 'T', // thr + kTryptophan = 'W', // trp + kTyrosine = 'Y', // tyr + kValine = 'V', // val +}; + +struct +{ + residue_type type; + char name[4]; +} const kResidueInfo[] = { + { kUnknownResidue, "UNK" }, + { kAlanine, "ALA" }, + { kArginine, "ARG" }, + { kAsparagine, "ASN" }, + { kAsparticAcid, "ASP" }, + { kCysteine, "CYS" }, + { kGlutamicAcid, "GLU" }, + { kGlutamine, "GLN" }, + { kGlycine, "GLY" }, + { kHistidine, "HIS" }, + { kIsoleucine, "ILE" }, + { kLeucine, "LEU" }, + { kLysine, "LYS" }, + { kMethionine, "MET" }, + { kPhenylalanine, "PHE" }, + { kProline, "PRO" }, + { kSerine, "SER" }, + { kThreonine, "THR" }, + { kTryptophan, "TRP" }, + { kTyrosine, "TYR" }, + { kValine, "VAL" } +}; + +constexpr residue_type MapResidue(std::string_view inName) +{ + residue_type result = kUnknownResidue; + + for (auto &ri : kResidueInfo) + { + if (inName == ri.name) + { + result = ri.type; + break; + } + } + + return result; +} + +struct HBond +{ + residue *res; + double energy; +}; + +enum class bridge_type +{ + None, + Parallel, + AntiParallel +}; + +struct bridge +{ + bridge_type type; + uint32_t sheet, ladder; + std::set link; + std::deque i, j; + std::string chainI, chainJ; + + bool operator<(const bridge &b) const { return chainI < b.chainI or (chainI == b.chainI and i.front() < b.i.front()); } +}; + +struct bridge_partner +{ + residue *m_residue; + uint32_t ladder; + bool parallel; +}; + +// -------------------------------------------------------------------- + +const float + // kSSBridgeDistance = 3.0f, + kMinimalDistance = 0.5f, + kMinimalCADistance = 9.0f, + kMinHBondEnergy = -9.9f, + kMaxHBondEnergy = -0.5f, + kCouplingConstant = -27.888f, // = -332 * 0.42 * 0.2 + kMaxPeptideBondLength = 2.5f; + +const float + kRadiusN = 1.65f, + kRadiusCA = 1.87f, + kRadiusC = 1.76f, + kRadiusO = 1.4f, + kRadiusSideAtom = 1.8f, + kRadiusWater = 1.4f; + +struct dssp::residue +{ + residue(residue &&) = default; + residue &operator=(residue &&) = default; + + residue(int model_nr, + std::string_view pdb_strand_id, int pdb_seq_num, std::string_view pdb_ins_code) + : mPDBStrandID(pdb_strand_id) + , mPDBSeqNum(pdb_seq_num) + , mPDBInsCode(pdb_ins_code) + , mChainBreak(chain_break_type::None) + , m_model_nr(model_nr) + { + // update the box containing all atoms + mBox[0].mX = mBox[0].mY = mBox[0].mZ = std::numeric_limits::max(); + mBox[1].mX = mBox[1].mY = mBox[1].mZ = -std::numeric_limits::max(); + + mH = point{ std::numeric_limits::max(), std::numeric_limits::max(), std::numeric_limits::max() }; + + for (auto &p : m_chiralAtoms) + p = {}; + } + + void addAtom(cif::row_handle atom) + { + std::string asymID, compID, atomID, type, authAsymID; + std::optional altID; + int seqID, authSeqID; + std::optional model; + float x, y, z; + + cif::tie(asymID, compID, atomID, altID, type, seqID, model, x, y, z, authAsymID, authSeqID) = + atom.get("label_asym_id", "label_comp_id", "label_atom_id", "label_alt_id", "type_symbol", "label_seq_id", + "pdbx_PDB_model_num", "Cartn_x", "Cartn_y", "Cartn_z", + "auth_asym_id", "auth_seq_id"); + + if (model and model != m_model_nr) + return; + + if (m_seen == 0) + { + mAsymID = asymID; + mCompoundID = compID; + mSeqID = seqID; + + mAuthSeqID = authSeqID; + mAuthAsymID = authAsymID; + + mType = MapResidue(mCompoundID); + + if (altID) + mAltID = *altID; + } + + if (atomID == "CA") + { + m_seen |= 1; + mCAlpha = { x, y, z }; + ExtendBox(mCAlpha, kRadiusCA + 2 * kRadiusWater); + + if (mType == kValine) + m_chiralAtoms[1] = mCAlpha; + } + else if (atomID == "C") + { + m_seen |= 2; + mC = { x, y, z }; + ExtendBox(mC, kRadiusC + 2 * kRadiusWater); + } + else if (atomID == "N") + { + m_seen |= 4; + mH = mN = { x, y, z }; + ExtendBox(mN, kRadiusN + 2 * kRadiusWater); + } + else if (atomID == "O") + { + m_seen |= 8; + mO = { x, y, z }; + ExtendBox(mO, kRadiusO + 2 * kRadiusWater); + } + else if (type != "H") + { + m_seen |= 16; + mSideChain.emplace_back(atomID, point{ x, y, z }); + ExtendBox(point{ x, y, z }, kRadiusSideAtom + 2 * kRadiusWater); + + if (mType == kLeucine) + { + if (atomID == "CG") + m_chiralAtoms[0] = { x, y, z }; + else if (atomID == "CB") + m_chiralAtoms[1] = { x, y, z }; + else if (atomID == "CD1") + m_chiralAtoms[2] = { x, y, z }; + else if (atomID == "CD2") + m_chiralAtoms[3] = { x, y, z }; + } + else if (mType == kValine) + { + if (atomID == "CB") + m_chiralAtoms[0] = { x, y, z }; + else if (atomID == "CG1") + m_chiralAtoms[2] = { x, y, z }; + else if (atomID == "CG2") + m_chiralAtoms[3] = { x, y, z }; + } + } + } + + void finish() + { + const int kSeenAll = (1 bitor 2 bitor 4 bitor 8); + mComplete = (m_seen bitand kSeenAll) == kSeenAll; + + if (mType == kValine or mType == kLeucine) + mChiralVolume = dot_product(m_chiralAtoms[1] - m_chiralAtoms[0], + cross_product(m_chiralAtoms[2] - m_chiralAtoms[0], m_chiralAtoms[3] - m_chiralAtoms[0])); + + mRadius = mBox[1].mX - mBox[0].mX; + if (mRadius < mBox[1].mY - mBox[0].mY) + mRadius = mBox[1].mY - mBox[0].mY; + if (mRadius < mBox[1].mZ - mBox[0].mZ) + mRadius = mBox[1].mZ - mBox[0].mZ; + + mCenter.mX = (mBox[0].mX + mBox[1].mX) / 2; + mCenter.mY = (mBox[0].mY + mBox[1].mY) / 2; + mCenter.mZ = (mBox[0].mZ + mBox[1].mZ) / 2; + } + + void assignHydrogen() + { + // assign the Hydrogen + mH = mN; + + if (mType != kProline and mPrev != nullptr) + { + auto pc = mPrev->mC; + auto po = mPrev->mO; + + float CODistance = static_cast(distance(pc, po)); + + mH.mX += (pc.mX - po.mX) / CODistance; + mH.mY += (pc.mY - po.mY) / CODistance; + mH.mZ += (pc.mZ - po.mZ) / CODistance; + } + } + + void SetSecondaryStructure(structure_type inSS) { mSecondaryStructure = inSS; } + structure_type GetSecondaryStructure() const { return mSecondaryStructure; } + + void SetBetaPartner(uint32_t n, residue &inResidue, uint32_t inLadder, bool inParallel) + { + assert(n == 0 or n == 1); + + mBetaPartner[n].m_residue = &inResidue; + mBetaPartner[n].ladder = inLadder; + mBetaPartner[n].parallel = inParallel; + } + + bridge_partner GetBetaPartner(uint32_t n) const + { + assert(n == 0 or n == 1); + return mBetaPartner[n]; + } + + void SetSheet(uint32_t inSheet) { mSheet = inSheet; } + uint32_t GetSheet() const { return mSheet; } + + void SetStrand(uint32_t inStrand) { mStrand = inStrand; } + uint32_t GetStrand() const { return mStrand; } + + bool IsBend() const { return mBend; } + void SetBend(bool inBend) { mBend = inBend; } + + helix_position_type GetHelixFlag(helix_type helixType) const + { + size_t stride = static_cast(helixType); + assert(stride < 4); + return mHelixFlags[stride]; + } + + bool IsHelixStart(helix_type helixType) const + { + size_t stride = static_cast(helixType); + assert(stride < 4); + return mHelixFlags[stride] == helix_position_type::Start or mHelixFlags[stride] == helix_position_type::StartAndEnd; + } + + void SetHelixFlag(helix_type helixType, helix_position_type inHelixFlag) + { + size_t stride = static_cast(helixType); + assert(stride < 4); + mHelixFlags[stride] = inHelixFlag; + } + + void SetSSBridgeNr(uint8_t inBridgeNr) + { + if (mType != kCysteine) + throw std::runtime_error("Only cysteine residues can form sulphur bridges"); + mSSBridgeNr = inBridgeNr; + } + + uint8_t GetSSBridgeNr() const + { + if (mType != kCysteine) + throw std::runtime_error("Only cysteine residues can form sulphur bridges"); + return mSSBridgeNr; + } + + float CalculateSurface(const std::vector &inResidues); + float CalculateSurface(const point &inAtom, float inRadius, const std::vector &inNeighbours); + + bool AtomIntersectsBox(const point &atom, float inRadius) const + { + return atom.mX + inRadius >= mBox[0].mX and + atom.mX - inRadius <= mBox[1].mX and + atom.mY + inRadius >= mBox[0].mY and + atom.mY - inRadius <= mBox[1].mY and + atom.mZ + inRadius >= mBox[0].mZ and + atom.mZ - inRadius <= mBox[1].mZ; + } + + void ExtendBox(const point &atom, float inRadius) + { + if (mBox[0].mX > atom.mX - inRadius) + mBox[0].mX = atom.mX - inRadius; + if (mBox[0].mY > atom.mY - inRadius) + mBox[0].mY = atom.mY - inRadius; + if (mBox[0].mZ > atom.mZ - inRadius) + mBox[0].mZ = atom.mZ - inRadius; + if (mBox[1].mX < atom.mX + inRadius) + mBox[1].mX = atom.mX + inRadius; + if (mBox[1].mY < atom.mY + inRadius) + mBox[1].mY = atom.mY + inRadius; + if (mBox[1].mZ < atom.mZ + inRadius) + mBox[1].mZ = atom.mZ + inRadius; + } + + point get_atom(std::string_view name) + { + if (name == "CA") + return mCAlpha; + else if (name == "C") + return mC; + else if (name == "N") + return mN; + else if (name == "O") + return mO; + else if (name == "H") + return mH; + else + { + for (const auto &[n, p] : mSideChain) + { + if (n == name) + return p; + } + } + + return {}; + } + + residue *mNext = nullptr; + residue *mPrev = nullptr; + + // const Monomer &mM; + + std::string mAsymID; + int mSeqID; + std::string mCompoundID; + std::string mAltID; + int mNumber; + bool mComplete = false; + + std::string mAuthAsymID; + int mAuthSeqID; + + std::string mPDBStrandID; + int mPDBSeqNum; + std::string mPDBInsCode; + + point mCAlpha, mC, mN, mO, mH; + point mBox[2] = {}; + float mRadius; + point mCenter; + std::vector> mSideChain; + float mAccessibility = 0; + float mChiralVolume = 0; + + // float mAlpha = 360, mKappa = 360, mPhi = 360, mPsi = 360, mTCO = 0, mOmega = 360; + std::optional mAlpha, mKappa, mPhi, mPsi, mTCO, mOmega; + + residue_type mType; + uint8_t mSSBridgeNr = 0; + structure_type mSecondaryStructure = structure_type::Loop; + HBond mHBondDonor[2] = {}, mHBondAcceptor[2] = {}; + bridge_partner mBetaPartner[2] = {}; + uint32_t mSheet = 0; + uint32_t mStrand = 0; // Added to ease the writing of mmCIF's struct_sheet and friends + helix_position_type mHelixFlags[4] = { helix_position_type::None, helix_position_type::None, helix_position_type::None, helix_position_type::None }; // + bool mBend = false; + chain_break_type mChainBreak = chain_break_type::None; + + int m_seen = 0, m_model_nr = 0; + std::array m_chiralAtoms; +}; + +// -------------------------------------------------------------------- + +class accumulator +{ + public: + struct candidate + { + point location; + double radius; + double distance; + + bool operator<(const candidate &rhs) const + { + return distance < rhs.distance; + } + }; + + void operator()(const point &a, const point &b, double d, double r) + { + double distance = distance_sq(a, b); + + d += kRadiusWater; + r += kRadiusWater; + + double test = d + r; + test *= test; + + if (distance < test and distance > 0.0001) + { + candidate c = { b - a, r * r, distance }; + + m_x.push_back(c); + push_heap(m_x.begin(), m_x.end()); + } + } + + void sort() + { + sort_heap(m_x.begin(), m_x.end()); + } + + std::vector m_x; +}; + +// we use a fibonacci sphere to calculate the even distribution of the dots +class MSurfaceDots +{ + public: + static MSurfaceDots &Instance(); + + size_t size() const { return mPoints.size(); } + const point &operator[](size_t inIx) const { return mPoints[inIx]; } + double weight() const { return mWeight; } + + private: + MSurfaceDots(int32_t inN); + + std::vector mPoints; + double mWeight; +}; + +MSurfaceDots &MSurfaceDots::Instance() +{ + const int32_t kN = 200; + + static MSurfaceDots sInstance(kN); + return sInstance; +} + +MSurfaceDots::MSurfaceDots(int32_t N) +{ + auto P = 2 * N + 1; + + const float kGoldenRatio = (1 + std::sqrt(5.0f)) / 2; + + mWeight = (4 * kPI) / P; + + for (auto i = -N; i <= N; ++i) + { + float lat = std::asin((2.0f * i) / P); + float lon = static_cast(std::fmod(i, kGoldenRatio) * 2 * kPI / kGoldenRatio); + + mPoints.emplace_back(point{ std::sin(lon) * std::cos(lat), std::cos(lon) * std::cos(lat), std::sin(lat) }); + } +} + +float residue::CalculateSurface(const point &inAtom, float inRadius, const std::vector &inNeighbours) +{ + accumulator accumulate; + + for (auto r : inNeighbours) + { + if (r->AtomIntersectsBox(inAtom, inRadius)) + { + accumulate(inAtom, r->mN, inRadius, kRadiusN); + accumulate(inAtom, r->mCAlpha, inRadius, kRadiusCA); + accumulate(inAtom, r->mC, inRadius, kRadiusC); + accumulate(inAtom, r->mO, inRadius, kRadiusO); + + for (const auto &[name, atom] : r->mSideChain) + accumulate(inAtom, atom, inRadius, kRadiusSideAtom); + } + } + + accumulate.sort(); + + float radius = inRadius + kRadiusWater; + float surface = 0; + + MSurfaceDots &surfaceDots = MSurfaceDots::Instance(); + + for (size_t i = 0; i < surfaceDots.size(); ++i) + { + point xx = surfaceDots[i] * radius; + + bool free = true; + for (size_t k = 0; free and k < accumulate.m_x.size(); ++k) + free = accumulate.m_x[k].radius < distance_sq(xx, accumulate.m_x[k].location); + + if (free) + surface += static_cast(surfaceDots.weight()); + } + + return surface * radius * radius; +} + +float residue::CalculateSurface(const std::vector &inResidues) +{ + std::vector neighbours; + + for (auto &r : inResidues) + { + point center = r.mCenter; + float radius = r.mRadius; + + if (distance_sq(mCenter, center) < (mRadius + radius) * (mRadius + radius)) + neighbours.push_back(const_cast(&r)); + } + + mAccessibility = CalculateSurface(mN, kRadiusN, neighbours) + + CalculateSurface(mCAlpha, kRadiusCA, neighbours) + + CalculateSurface(mC, kRadiusC, neighbours) + + CalculateSurface(mO, kRadiusO, neighbours); + + for (const auto &[name, atom] : mSideChain) + mAccessibility += CalculateSurface(atom, kRadiusSideAtom, neighbours); + + return mAccessibility; +} + +void CalculateAccessibilities(std::vector &inResidues, statistics &stats) +{ + stats.accessible_surface = 0; + for (auto &residue : inResidues) + stats.accessible_surface += residue.CalculateSurface(inResidues); +} + +// -------------------------------------------------------------------- +// TODO: use the angle to improve bond energy calculation. + +double CalculateHBondEnergy(residue &inDonor, residue &inAcceptor) +{ + double result = 0; + + if (inDonor.mType != kProline) + { + double distanceHO = distance(inDonor.mH, inAcceptor.mO); + double distanceHC = distance(inDonor.mH, inAcceptor.mC); + double distanceNC = distance(inDonor.mN, inAcceptor.mC); + double distanceNO = distance(inDonor.mN, inAcceptor.mO); + + if (distanceHO < kMinimalDistance or distanceHC < kMinimalDistance or distanceNC < kMinimalDistance or distanceNO < kMinimalDistance) + result = kMinHBondEnergy; + else + result = kCouplingConstant / distanceHO - kCouplingConstant / distanceHC + kCouplingConstant / distanceNC - kCouplingConstant / distanceNO; + + // DSSP compatibility mode: + result = std::round(result * 1000) / 1000; + + if (result < kMinHBondEnergy) + result = kMinHBondEnergy; + } + + // update donor + if (result < inDonor.mHBondAcceptor[0].energy) + { + inDonor.mHBondAcceptor[1] = inDonor.mHBondAcceptor[0]; + inDonor.mHBondAcceptor[0].res = &inAcceptor; + inDonor.mHBondAcceptor[0].energy = result; + } + else if (result < inDonor.mHBondAcceptor[1].energy) + { + inDonor.mHBondAcceptor[1].res = &inAcceptor; + inDonor.mHBondAcceptor[1].energy = result; + } + + // and acceptor + if (result < inAcceptor.mHBondDonor[0].energy) + { + inAcceptor.mHBondDonor[1] = inAcceptor.mHBondDonor[0]; + inAcceptor.mHBondDonor[0].res = &inDonor; + inAcceptor.mHBondDonor[0].energy = result; + } + else if (result < inAcceptor.mHBondDonor[1].energy) + { + inAcceptor.mHBondDonor[1].res = &inDonor; + inAcceptor.mHBondDonor[1].energy = result; + } + + return result; +} + +// -------------------------------------------------------------------- + +void CalculateHBondEnergies(std::vector &inResidues, std::vector> &q) +{ + std::unique_ptr progress; + if (cif::VERBOSE == 0 or cif::VERBOSE == 1) + progress.reset(new cif::progress_bar(q.size(), "calculate hbond energies")); + + for (const auto &[i, j] : q) + { + auto &ri = inResidues[i]; + auto &rj = inResidues[j]; + + CalculateHBondEnergy(ri, rj); + if (j != i + 1) + CalculateHBondEnergy(rj, ri); + + if (progress) + progress->consumed(1); + } +} + +// -------------------------------------------------------------------- + +bool NoChainBreak(const residue *a, const residue *b) +{ + bool result = a->mAsymID == b->mAsymID; + for (auto r = a; result and r != b; r = r->mNext) + { + auto next = r->mNext; + if (next == nullptr) + result = false; + else + result = next->mNumber == r->mNumber + 1; + } + return result; +} + +bool NoChainBreak(const residue &a, const residue &b) +{ + return NoChainBreak(&a, &b); +} + +// -------------------------------------------------------------------- + +bool TestBond(const residue *a, const residue *b) +{ + return (a->mHBondAcceptor[0].res == b and a->mHBondAcceptor[0].energy < kMaxHBondEnergy) or + (a->mHBondAcceptor[1].res == b and a->mHBondAcceptor[1].energy < kMaxHBondEnergy); +} + +bool test_bond(dssp::residue_info const &a, dssp::residue_info const &b) +{ + return a and b and TestBond(a.m_impl, b.m_impl); +} + +// -------------------------------------------------------------------- + +bridge_type TestBridge(const residue &r1, const residue &r2) +{ // I. a d II. a d parallel + auto a = r1.mPrev; // \ / + auto b = &r1; // b e b e + auto c = r1.mNext; // / \ .. + auto d = r2.mPrev; // c f c f + auto e = &r2; // + auto f = r2.mNext; // III. a <- f IV. a f antiparallel + // + bridge_type result = bridge_type::None; // b e b <-> e + // + // c -> d c d + + if (a and c and NoChainBreak(a, c) and d and f and NoChainBreak(d, f)) + { + if ((TestBond(c, e) and TestBond(e, a)) or (TestBond(f, b) and TestBond(b, d))) + result = bridge_type::Parallel; + else if ((TestBond(c, d) and TestBond(f, a)) or (TestBond(e, b) and TestBond(b, e))) + result = bridge_type::AntiParallel; + } + + return result; +} + +// -------------------------------------------------------------------- +// return true if any of the residues in bridge a is identical to any of the residues in bridge b +bool Linked(const bridge &a, const bridge &b) +{ + return find_first_of(a.i.begin(), a.i.end(), b.i.begin(), b.i.end()) != a.i.end() or + find_first_of(a.i.begin(), a.i.end(), b.j.begin(), b.j.end()) != a.i.end() or + find_first_of(a.j.begin(), a.j.end(), b.i.begin(), b.i.end()) != a.j.end() or + find_first_of(a.j.begin(), a.j.end(), b.j.begin(), b.j.end()) != a.j.end(); +} + +// -------------------------------------------------------------------- + +void CalculateBetaSheets(std::vector &inResidues, statistics &stats, std::vector> &q) +{ + // if (cif::VERBOSE) + // std::cerr << "calculating beta sheets" << std::endl; + + std::unique_ptr progress; + if (cif::VERBOSE == 0 or cif::VERBOSE == 1) + progress.reset(new cif::progress_bar(q.size(), "calculate beta sheets")); + + // Calculate Bridges + std::vector bridges; + + for (const auto &[i, j] : q) + { + if (progress) + progress->consumed(1); + + auto &ri = inResidues[i]; + auto &rj = inResidues[j]; + + bridge_type type = TestBridge(ri, rj); + if (type == bridge_type::None) + continue; + + bool found = false; + for (bridge &bridge : bridges) + { + if (type != bridge.type or i != bridge.i.back() + 1) + continue; + + if (type == bridge_type::Parallel and bridge.j.back() + 1 == j) + { + bridge.i.push_back(i); + bridge.j.push_back(j); + found = true; + break; + } + + if (type == bridge_type::AntiParallel and bridge.j.front() - 1 == j) + { + bridge.i.push_back(i); + bridge.j.push_front(j); + found = true; + break; + } + } + + if (not found) + { + bridge bridge = {}; + + bridge.type = type; + bridge.i.push_back(i); + bridge.chainI = ri.mAsymID; + bridge.j.push_back(j); + bridge.chainJ = rj.mAsymID; + + bridges.push_back(bridge); + } + } + + // extend ladders + std::sort(bridges.begin(), bridges.end()); + + for (uint32_t i = 0; i < bridges.size(); ++i) + { + for (uint32_t j = i + 1; j < bridges.size(); ++j) + { + uint32_t ibi = bridges[i].i.front(); + uint32_t iei = bridges[i].i.back(); + uint32_t jbi = bridges[i].j.front(); + uint32_t jei = bridges[i].j.back(); + uint32_t ibj = bridges[j].i.front(); + uint32_t iej = bridges[j].i.back(); + uint32_t jbj = bridges[j].j.front(); + uint32_t jej = bridges[j].j.back(); + + if (bridges[i].type != bridges[j].type or + NoChainBreak(inResidues[std::min(ibi, ibj)], inResidues[std::max(iei, iej)]) == false or + NoChainBreak(inResidues[std::min(jbi, jbj)], inResidues[std::max(jei, jej)]) == false or + ibj - iei >= 6 or + (iei >= ibj and ibi <= iej)) + { + continue; + } + + bool bulge; + if (bridges[i].type == bridge_type::Parallel) + bulge = ((jbj - jei < 6 and ibj - iei < 3) or (jbj - jei < 3)); + else + bulge = ((jbi - jej < 6 and ibj - iei < 3) or (jbi - jej < 3)); + + if (bulge) + { + bridges[i].i.insert(bridges[i].i.end(), bridges[j].i.begin(), bridges[j].i.end()); + if (bridges[i].type == bridge_type::Parallel) + bridges[i].j.insert(bridges[i].j.end(), bridges[j].j.begin(), bridges[j].j.end()); + else + bridges[i].j.insert(bridges[i].j.begin(), bridges[j].j.begin(), bridges[j].j.end()); + bridges.erase(bridges.begin() + j); + --j; + } + } + } + + // Sheet + std::set ladderset; + for (bridge &bridge : bridges) + { + ladderset.insert(&bridge); + + size_t n = bridge.i.size(); + if (n > dssp::kHistogramSize) + n = dssp::kHistogramSize; + + if (bridge.type == bridge_type::Parallel) + stats.histogram.parallel_bridges_per_ladder[n - 1] += 1; + else + stats.histogram.antiparallel_bridges_per_ladder[n - 1] += 1; + } + + uint32_t sheet = 1, ladder = 0; + while (not ladderset.empty()) + { + std::set sheetset; + sheetset.insert(*ladderset.begin()); + ladderset.erase(ladderset.begin()); + + bool done = false; + while (not done) + { + done = true; + for (bridge *a : sheetset) + { + for (bridge *b : ladderset) + { + if (Linked(*a, *b)) + { + sheetset.insert(b); + ladderset.erase(b); + done = false; + break; + } + } + if (not done) + break; + } + } + + for (bridge *bridge : sheetset) + { + bridge->ladder = ladder; + bridge->sheet = sheet; + bridge->link = sheetset; + + ++ladder; + } + + size_t nrOfLaddersPerSheet = sheetset.size(); + if (nrOfLaddersPerSheet > dssp::kHistogramSize) + nrOfLaddersPerSheet = dssp::kHistogramSize; + if (nrOfLaddersPerSheet == 1 and (*sheetset.begin())->i.size() > 1) + stats.histogram.ladders_per_sheet[0] += 1; + else if (nrOfLaddersPerSheet > 1) + stats.histogram.ladders_per_sheet[nrOfLaddersPerSheet - 1] += 1; + + ++sheet; + } + + for (bridge &bridge : bridges) + { + // find out if any of the i and j set members already have + // a bridge assigned, if so, we're assigning bridge 2 + + uint32_t betai = 0, betaj = 0; + + for (uint32_t l : bridge.i) + { + if (inResidues[l].GetBetaPartner(0).m_residue != nullptr) + { + betai = 1; + break; + } + } + + for (uint32_t l : bridge.j) + { + if (inResidues[l].GetBetaPartner(0).m_residue != nullptr) + { + betaj = 1; + break; + } + } + + structure_type ss = structure_type::Betabridge; + if (bridge.i.size() > 1) + ss = structure_type::Strand; + + if (bridge.type == bridge_type::Parallel) + { + stats.count.H_bonds_in_parallel_bridges += bridge.i.back() - bridge.i.front() + 2; + + std::deque::iterator j = bridge.j.begin(); + for (uint32_t i : bridge.i) + inResidues[i].SetBetaPartner(betai, inResidues[*j++], bridge.ladder, true); + + j = bridge.i.begin(); + for (uint32_t i : bridge.j) + inResidues[i].SetBetaPartner(betaj, inResidues[*j++], bridge.ladder, true); + } + else + { + stats.count.H_bonds_in_antiparallel_bridges += bridge.i.back() - bridge.i.front() + 2; + + std::deque::reverse_iterator j = bridge.j.rbegin(); + for (uint32_t i : bridge.i) + inResidues[i].SetBetaPartner(betai, inResidues[*j++], bridge.ladder, false); + + j = bridge.i.rbegin(); + for (uint32_t i : bridge.j) + inResidues[i].SetBetaPartner(betaj, inResidues[*j++], bridge.ladder, false); + } + + for (uint32_t i = bridge.i.front(); i <= bridge.i.back(); ++i) + { + if (inResidues[i].GetSecondaryStructure() != structure_type::Strand) + inResidues[i].SetSecondaryStructure(ss); + inResidues[i].SetSheet(bridge.sheet); + } + + for (uint32_t i = bridge.j.front(); i <= bridge.j.back(); ++i) + { + if (inResidues[i].GetSecondaryStructure() != structure_type::Strand) + inResidues[i].SetSecondaryStructure(ss); + inResidues[i].SetSheet(bridge.sheet); + } + } + + // Create 'strands'. A strand is a range of residues without a gap in between + // that belong to the same sheet. + + int strand = 0; + for (uint32_t iSheet = 1; iSheet < sheet; ++iSheet) + { + int lastNr = -1; + for (auto &res : inResidues) + { + if (res.mSheet != iSheet) + continue; + + if (lastNr + 1 < res.mNumber) + ++strand; + + res.mStrand = strand; + lastNr = res.mNumber; + } + } +} + +// -------------------------------------------------------------------- + +void CalculateAlphaHelices(std::vector &inResidues, statistics &stats, bool inPreferPiHelices = true) +{ + if (cif::VERBOSE) + std::cerr << "calculating alpha helices" << std::endl; + + // Helix and Turn + for (helix_type helixType : { helix_type::_3_10, helix_type::alpha, helix_type::pi }) + { + uint32_t stride = static_cast(helixType) + 3; + + for (uint32_t i = 0; i + stride < inResidues.size(); ++i) + { + if (NoChainBreak(inResidues[i], inResidues[i + stride]) and TestBond(&inResidues[i + stride], &inResidues[i])) + { + inResidues[i + stride].SetHelixFlag(helixType, helix_position_type::End); + for (uint32_t j = i + 1; j < i + stride; ++j) + { + if (inResidues[j].GetHelixFlag(helixType) == helix_position_type::None) + inResidues[j].SetHelixFlag(helixType, helix_position_type::Middle); + } + + if (inResidues[i].GetHelixFlag(helixType) == helix_position_type::End) + inResidues[i].SetHelixFlag(helixType, helix_position_type::StartAndEnd); + else + inResidues[i].SetHelixFlag(helixType, helix_position_type::Start); + } + } + } + + for (auto &r : inResidues) + { + if (r.mKappa.has_value()) + r.SetBend(*r.mKappa > 70); + } + + for (uint32_t i = 1; i + 4 < inResidues.size(); ++i) + { + if (inResidues[i].IsHelixStart(helix_type::alpha) and inResidues[i - 1].IsHelixStart(helix_type::alpha)) + { + for (uint32_t j = i; j <= i + 3; ++j) + inResidues[j].SetSecondaryStructure(structure_type::Alphahelix); + } + } + + for (uint32_t i = 1; i + 3 < inResidues.size(); ++i) + { + if (inResidues[i].IsHelixStart(helix_type::_3_10) and inResidues[i - 1].IsHelixStart(helix_type::_3_10)) + { + bool empty = true; + for (uint32_t j = i; empty and j <= i + 2; ++j) + empty = inResidues[j].GetSecondaryStructure() == structure_type::Loop or inResidues[j].GetSecondaryStructure() == structure_type::Helix_3; + if (empty) + { + for (uint32_t j = i; j <= i + 2; ++j) + inResidues[j].SetSecondaryStructure(structure_type::Helix_3); + } + } + } + + for (uint32_t i = 1; i + 5 < inResidues.size(); ++i) + { + if (inResidues[i].IsHelixStart(helix_type::pi) and inResidues[i - 1].IsHelixStart(helix_type::pi)) + { + bool empty = true; + for (uint32_t j = i; empty and j <= i + 4; ++j) + empty = inResidues[j].GetSecondaryStructure() == structure_type::Loop or inResidues[j].GetSecondaryStructure() == structure_type::Helix_5 or + (inPreferPiHelices and inResidues[j].GetSecondaryStructure() == structure_type::Alphahelix); + if (empty) + { + for (uint32_t j = i; j <= i + 4; ++j) + inResidues[j].SetSecondaryStructure(structure_type::Helix_5); + } + } + } + + for (uint32_t i = 1; i + 1 < inResidues.size(); ++i) + { + if (inResidues[i].GetSecondaryStructure() == structure_type::Loop) + { + bool isTurn = false; + for (helix_type helixType : { helix_type::_3_10, helix_type::alpha, helix_type::pi }) + { + uint32_t stride = 3 + static_cast(helixType); + for (uint32_t k = 1; k < stride and not isTurn; ++k) + isTurn = (i >= k) and inResidues[i - k].IsHelixStart(helixType); + } + + if (isTurn) + inResidues[i].SetSecondaryStructure(structure_type::Turn); + else if (inResidues[i].IsBend()) + inResidues[i].SetSecondaryStructure(structure_type::Bend); + } + } + + std::string asym; + size_t helixLength = 0; + for (auto &r : inResidues) + { + if (r.mAsymID != asym) + { + helixLength = 0; + asym = r.mAsymID; + } + + if (r.GetSecondaryStructure() == structure_type::Alphahelix) + ++helixLength; + else if (helixLength > 0) + { + if (helixLength > dssp::kHistogramSize) + helixLength = dssp::kHistogramSize; + + stats.histogram.residues_per_alpha_helix[helixLength - 1] += 1; + helixLength = 0; + } + } +} + +// -------------------------------------------------------------------- + +void CalculatePPHelices(std::vector &inResidues, statistics &stats, int stretch_length) +{ + if (cif::VERBOSE) + std::cerr << "calculating pp helices" << std::endl; + + size_t N = inResidues.size(); + + const float epsilon = 29; + const float phi_min = -75 - epsilon; + const float phi_max = -75 + epsilon; + const float psi_min = 145 - epsilon; + const float psi_max = 145 + epsilon; + + std::vector phi(N), psi(N); + + for (uint32_t i = 1; i + 1 < inResidues.size(); ++i) + { + phi[i] = static_cast(inResidues[i].mPhi.value_or(360)); + psi[i] = static_cast(inResidues[i].mPsi.value_or(360)); + } + + for (uint32_t i = 1; i + 3 < inResidues.size(); ++i) + { + switch (stretch_length) + { + case 2: + { + if (phi_min > phi[i + 0] or phi[i + 0] > phi_max or + phi_min > phi[i + 1] or phi[i + 1] > phi_max) + continue; + + if (psi_min > psi[i + 0] or psi[i + 0] > psi_max or + psi_min > psi[i + 1] or psi[i + 1] > psi_max) + continue; + + // auto phi_avg = (phi[i + 0] + phi[i + 1]) / 2; + // auto phi_sq = (phi[i + 0] - phi_avg) * (phi[i + 0] - phi_avg) + + // (phi[i + 1] - phi_avg) * (phi[i + 1] - phi_avg); + + // if (phi_sq >= 200) + // continue; + + // auto psi_avg = (psi[i + 0] + psi[i + 1]) / 2; + // auto psi_sq = (psi[i + 0] - psi_avg) * (psi[i + 0] - psi_avg) + + // (psi[i + 1] - psi_avg) * (psi[i + 1] - psi_avg); + + // if (psi_sq >= 200) + // continue; + + switch (inResidues[i].GetHelixFlag(helix_type::pp)) + { + case helix_position_type::None: + inResidues[i].SetHelixFlag(helix_type::pp, helix_position_type::Start); + break; + + case helix_position_type::End: + inResidues[i].SetHelixFlag(helix_type::pp, helix_position_type::Middle); + break; + + default: + break; + } + + inResidues[i + 1].SetHelixFlag(helix_type::pp, helix_position_type::End); + + if (inResidues[i].GetSecondaryStructure() == structure_type::Loop) + inResidues[i].SetSecondaryStructure(structure_type::Helix_PPII); + if (inResidues[i + 1].GetSecondaryStructure() == structure_type::Loop) + inResidues[i + 1].SetSecondaryStructure(structure_type::Helix_PPII); + } + break; + + case 3: + { + if (phi_min > phi[i + 0] or phi[i + 0] > phi_max or + phi_min > phi[i + 1] or phi[i + 1] > phi_max or + phi_min > phi[i + 2] or phi[i + 2] > phi_max) + continue; + + if (psi_min > psi[i + 0] or psi[i + 0] > psi_max or + psi_min > psi[i + 1] or psi[i + 1] > psi_max or + psi_min > psi[i + 2] or psi[i + 2] > psi_max) + continue; + + // auto phi_avg = (phi[i + 0] + phi[i + 1] + phi[i + 2]) / 3; + // auto phi_sq = (phi[i + 0] - phi_avg) * (phi[i + 0] - phi_avg) + + // (phi[i + 1] - phi_avg) * (phi[i + 1] - phi_avg) + + // (phi[i + 2] - phi_avg) * (phi[i + 2] - phi_avg); + + // if (phi_sq >= 300) + // continue; + + // auto psi_avg = (psi[i + 0] + psi[i + 1] + psi[i + 2]) / 3; + // auto psi_sq = (psi[i + 0] - psi_avg) * (psi[i + 0] - psi_avg) + + // (psi[i + 1] - psi_avg) * (psi[i + 1] - psi_avg) + + // (psi[i + 2] - psi_avg) * (psi[i + 2] - psi_avg); + + // if (psi_sq >= 300) + // continue; + + switch (inResidues[i].GetHelixFlag(helix_type::pp)) + { + case helix_position_type::None: + inResidues[i].SetHelixFlag(helix_type::pp, helix_position_type::Start); + break; + + case helix_position_type::End: + inResidues[i].SetHelixFlag(helix_type::pp, helix_position_type::StartAndEnd); + break; + + default: + break; + } + + inResidues[i + 1].SetHelixFlag(helix_type::pp, helix_position_type::Middle); + inResidues[i + 2].SetHelixFlag(helix_type::pp, helix_position_type::End); + + if (inResidues[i + 0].GetSecondaryStructure() == structure_type::Loop) + inResidues[i + 0].SetSecondaryStructure(structure_type::Helix_PPII); + + if (inResidues[i + 1].GetSecondaryStructure() == structure_type::Loop) + inResidues[i + 1].SetSecondaryStructure(structure_type::Helix_PPII); + + if (inResidues[i + 2].GetSecondaryStructure() == structure_type::Loop) + inResidues[i + 2].SetSecondaryStructure(structure_type::Helix_PPII); + + break; + } + + default: + throw std::runtime_error("Unsupported stretch length"); + } + } +} + +// -------------------------------------------------------------------- + +struct DSSP_impl +{ + DSSP_impl(const cif::datablock &db, int model_nr, int min_poly_proline_stretch_length); + + auto findRes(const std::string &asymID, int seqID) + { + return std::find_if(mResidues.begin(), mResidues.end(), [&](auto &r) + { return r.mAsymID == asymID and r.mSeqID == seqID; }); + } + + void calculateSurface(); + void calculateSecondaryStructure(); + + std::string GetPDBHEADERLine(); + std::string GetPDBCOMPNDLine(); + std::string GetPDBSOURCELine(); + std::string GetPDBAUTHORLine(); + + const cif::datablock &mDB; + std::vector mResidues; + std::vector> mSSBonds; + int m_min_poly_proline_stretch_length; + statistics mStats = {}; +}; + +// -------------------------------------------------------------------- + +DSSP_impl::DSSP_impl(const cif::datablock &db, int model_nr, int min_poly_proline_stretch_length) + : mDB(db) + , m_min_poly_proline_stretch_length(min_poly_proline_stretch_length) +{ + using namespace cif::literals; + + if (cif::VERBOSE) + std::cerr << "loading residues" << std::endl; + + int resNumber = 0; + + auto &pdbx_poly_seq_scheme = mDB["pdbx_poly_seq_scheme"]; + auto &atom_site = mDB["atom_site"]; + + using key_type = std::tuple; + using index_type = std::map; + + index_type index; + + mResidues.reserve(pdbx_poly_seq_scheme.size()); + + for (const auto &[asym_id, seq_id, pdb_strand_id, pdb_seq_num, pdb_ins_code] : pdbx_poly_seq_scheme.rows("asym_id", "seq_id", "pdb_strand_id", "pdb_seq_num", "pdb_ins_code")) + { + index[{ asym_id, seq_id }] = mResidues.size(); + mResidues.emplace_back(model_nr, pdb_strand_id, pdb_seq_num, pdb_ins_code); + } + + for (auto atom : atom_site) + { + std::string asym_id; + int seq_id; + + cif::tie(asym_id, seq_id) = atom.get("label_asym_id", "label_seq_id"); + auto i = index.find({ asym_id, seq_id }); + if (i == index.end()) + continue; + + mResidues[i->second].addAtom(atom); + } + + for (auto &residue : mResidues) + residue.finish(); + + mResidues.erase(std::remove_if(mResidues.begin(), mResidues.end(), [](const dssp::residue &r) + { return not r.mComplete; }), + mResidues.end()); + mStats.count.chains = 1; + + chain_break_type brk = chain_break_type::NewChain; + for (size_t i = 0; i < mResidues.size(); ++i) + { + auto &residue = mResidues[i]; + ++resNumber; + + if (i > 0) + { + if (distance(mResidues[i - 1].mC, mResidues[i].mN) > kMaxPeptideBondLength) + { + ++mStats.count.chains; + if (mResidues[i - 1].mAsymID == mResidues[i].mAsymID) + brk = chain_break_type::Gap; + else + brk = chain_break_type::NewChain; + + ++resNumber; + } + } + + residue.mChainBreak = brk; + residue.mNumber = resNumber; + + brk = chain_break_type::None; + } + + mStats.count.residues = static_cast(mResidues.size()); + + for (size_t i = 0; i + 1 < mResidues.size(); ++i) + { + auto &cur = mResidues[i]; + + auto &next = mResidues[i + 1]; + next.mPrev = &cur; + cur.mNext = &next; + } + + for (size_t i = 0; i < mResidues.size(); ++i) + { + auto &cur = mResidues[i]; + + if (i >= 2 and i + 2 < mResidues.size()) + { + auto &prevPrev = mResidues[i - 2]; + auto &nextNext = mResidues[i + 2]; + + if (NoChainBreak(prevPrev, nextNext) and prevPrev.mSeqID + 4 == nextNext.mSeqID) + { + float ckap = cosinus_angle( + cur.mCAlpha, + prevPrev.mCAlpha, + nextNext.mCAlpha, + cur.mCAlpha); + float skap = std::sqrt(1 - ckap * ckap); + + auto kappa = std::atan2(skap, ckap) * static_cast(180 / kPI); + if (not std::isnan(kappa)) + cur.mKappa = kappa; + } + } + + if (i + 1 < mResidues.size()) + { + auto &next = mResidues[i + 1]; + next.assignHydrogen(); + if (NoChainBreak(cur, next)) + { + cur.mPsi = dihedral_angle(cur.mN, cur.mCAlpha, cur.mC, next.mN); + cur.mOmega = dihedral_angle(cur.mCAlpha, cur.mC, next.mN, next.mCAlpha); + } + } + + if (i > 0) + { + auto &prev = mResidues[i - 1]; + if (NoChainBreak(prev, cur)) + { + cur.mTCO = cosinus_angle(cur.mC, cur.mO, prev.mC, prev.mO); + cur.mPhi = dihedral_angle(prev.mC, cur.mN, cur.mCAlpha, cur.mC); + } + } + + if (i >= 1 and i + 2 < mResidues.size()) + { + auto &prev = mResidues[i - 1]; + auto &next = mResidues[i + 1]; + auto &nextNext = mResidues[i + 2]; + + if (NoChainBreak(prev, nextNext)) + cur.mAlpha = dihedral_angle(prev.mCAlpha, cur.mCAlpha, next.mCAlpha, nextNext.mCAlpha); + } + } +} + +void DSSP_impl::calculateSecondaryStructure() +{ + if (cif::VERBOSE) + std::cerr << "calculating secondary structure" << std::endl; + + using namespace cif::literals; + + for (auto [asym1, seq1, asym2, seq2] : mDB["struct_conn"].find("conn_type_id"_key == "disulf", + "ptnr1_label_asym_id", "ptnr1_label_seq_id", "ptnr2_label_asym_id", "ptnr2_label_seq_id")) + { + auto r1 = findRes(asym1, seq1); + if (r1 == mResidues.end()) + { + if (cif::VERBOSE > 0) + std::cerr << "Missing (incomplete?) residue for SS bond when trying to find " << asym1 << '/' << seq1 << std::endl; + continue; + // throw std::runtime_error("Invalid file, missing residue for SS bond"); + } + + auto r2 = findRes(asym2, seq2); + if (r2 == mResidues.end()) + { + if (cif::VERBOSE > 0) + std::cerr << "Missing (incomplete?) residue for SS bond when trying to find " << asym2 << '/' << seq2 << std::endl; + continue; + // throw std::runtime_error("Invalid file, missing residue for SS bond"); + } + + mSSBonds.emplace_back(&*r1, &*r2); + } + + // Prefetch the c-alpha positions. No, really, that might be the trick + + std::vector cAlphas; + cAlphas.reserve(mResidues.size()); + for (auto &r : mResidues) + cAlphas.emplace_back(r.mCAlpha); + + std::unique_ptr progress; + if (cif::VERBOSE == 0 or cif::VERBOSE == 1) + progress.reset(new cif::progress_bar((mResidues.size() * (mResidues.size() - 1)) / 2, "calculate distances")); + + // Calculate the HBond energies + std::vector> near; + + for (uint32_t i = 0; i + 1 < mResidues.size(); ++i) + { + auto cai = cAlphas[i]; + + for (uint32_t j = i + 1; j < mResidues.size(); ++j) + { + auto caj = cAlphas[j]; + + if (distance_sq(cai, caj) > (kMinimalCADistance * kMinimalCADistance)) + continue; + + near.emplace_back(i, j); + } + + if (progress) + progress->consumed(mResidues.size() - i - 1); + } + + if (cif::VERBOSE > 0) + std::cerr << "Considering " << near.size() << " pairs of residues" << std::endl; + + progress.reset(nullptr); + + CalculateHBondEnergies(mResidues, near); + CalculateBetaSheets(mResidues, mStats, near); + CalculateAlphaHelices(mResidues, mStats); + CalculatePPHelices(mResidues, mStats, m_min_poly_proline_stretch_length); + + if (cif::VERBOSE > 1) + { + for (auto &r : mResidues) + { + char helix[5] = {}; + for (helix_type helixType : { helix_type::_3_10, helix_type::alpha, helix_type::pi, helix_type::pp }) + { + switch (r.GetHelixFlag(helixType)) + { + case helix_position_type::Start: helix[static_cast(helixType)] = '>'; break; + case helix_position_type::Middle: helix[static_cast(helixType)] = helixType == helix_type::pp ? 'P' : '3' + static_cast(helixType); break; + case helix_position_type::StartAndEnd: helix[static_cast(helixType)] = 'X'; break; + case helix_position_type::End: helix[static_cast(helixType)] = '<'; break; + case helix_position_type::None: helix[static_cast(helixType)] = ' '; break; + } + } + + auto id = r.mAsymID + ':' + std::to_string(r.mSeqID) + '/' + r.mCompoundID; + + std::cerr << id << std::string(12 - id.length(), ' ') + << char(r.mSecondaryStructure) << ' ' + << helix + << std::endl; + } + } + + // finish statistics + mStats.count.SS_bridges = static_cast(mSSBonds.size()); + + mStats.count.intra_chain_SS_bridges = 0; + uint8_t ssBondNr = 0; + for (const auto &[a, b] : mSSBonds) + { + if (a == b) + { + if (cif::VERBOSE > 0) + std::cerr << "In the SS bonds list, the residue " << a->mAsymID << ':' << a->mSeqID << " is bonded to itself" << std::endl; + continue; + } + + if (a->mAsymID == b->mAsymID and NoChainBreak(a, b)) + ++mStats.count.intra_chain_SS_bridges; + + a->mSSBridgeNr = b->mSSBridgeNr = ++ssBondNr; + } + + mStats.count.H_bonds = 0; + for (auto &r : mResidues) + { + auto donor = r.mHBondDonor; + + for (int i = 0; i < 2; ++i) + { + if (donor[i].res != nullptr and donor[i].energy < kMaxHBondEnergy) + { + ++mStats.count.H_bonds; + auto k = donor[i].res->mNumber - r.mNumber; + if (k >= -5 and k <= 5) + mStats.count.H_Bonds_per_distance[k + 5] += 1; + } + } + } +} + +void DSSP_impl::calculateSurface() +{ + CalculateAccessibilities(mResidues, mStats); +} + +// -------------------------------------------------------------------- + +// Truncate lines in pseudo PDB format to this length +const int kTruncateAt = 127; + +std::string FixStringLength(std::string s, std::string::size_type l = kTruncateAt) +{ + if (s.length() > l) + s = s.substr(0, l - 4) + "... "; + else if (s.length() < l) + s.append(l - s.length(), ' '); + + return s; +} + +std::string cif2pdbDate(const std::string &d) +{ + const std::regex rx(R"((\d{4})-(\d{2})(?:-(\d{2}))?)"); + const char *kMonths[12] = { + "JAN", "FEB", "MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP", "OCT", "NOV", "DEC" + }; + + std::smatch m; + std::ostringstream os; + + if (std::regex_match(d, m, rx)) + { + int year = std::stoi(m[1].str()); + int month = std::stoi(m[2].str()); + + if (m[3].matched) + os << std::setw(2) << std::setfill('0') << stoi(m[3].str()) << '-'; + os << kMonths[month - 1] << '-' << std::setw(2) << std::setfill('0') << (year % 100); + } + + return os.str(); +} + +std::string cif2pdbAuth(std::string name) +{ + const std::regex rx(R"(([^,]+), (\S+))"); + + std::smatch m; + if (std::regex_match(name, m, rx)) + name = m[2].str() + m[1].str(); + + return name; +} + +std::string DSSP_impl::GetPDBHEADERLine() +{ + std::string keywords; + auto &cat1 = mDB["struct_keywords"]; + + for (auto r : cat1) + { + keywords = FixStringLength(r["pdbx_keywords"].as(), 40); + break; + } + + std::string date; + for (auto r : mDB["pdbx_database_status"]) + { + date = r["recvd_initial_deposition_date"].as(); + if (date.empty()) + continue; + date = cif2pdbDate(date); + break; + } + + if (date.empty()) + { + for (auto r : mDB["database_PDB_rev"]) + { + date = r["date_original"].as(); + if (date.empty()) + continue; + date = cif2pdbDate(date); + break; + } + } + + date = FixStringLength(date, 9); + + // 0 1 2 3 4 5 6 7 8 + // HEADER xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxDDDDDDDDD IIII + char header[] = + "HEADER xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxDDDDDDDDD IIII"; + + std::copy(keywords.begin(), keywords.end(), header + 10); + std::copy(date.begin(), date.end(), header + 50); + + std::string id = mDB.name(); + if (id.length() < 4) + id.insert(id.end(), 4 - id.length(), ' '); + else if (id.length() > 4) + id.erase(id.begin() + 4, id.end()); + + std::copy(id.begin(), id.end(), header + 62); + + return FixStringLength(header); +} + +std::string DSSP_impl::GetPDBCOMPNDLine() +{ + // COMPND + using namespace std::placeholders; + using namespace cif::literals; + + int molID = 0; + std::vector cmpnd; + + for (auto r : mDB["entity"].find("type"_key == "polymer")) + { + std::string entityID = r["id"].as(); + + ++molID; + cmpnd.push_back("MOL_ID: " + std::to_string(molID)); + + std::string molecule = r["pdbx_description"].as(); + cmpnd.push_back("MOLECULE: " + molecule); + + auto poly = mDB["entity_poly"].find("entity_id"_key == entityID); + if (not poly.empty()) + { + std::string chains = poly.front()["pdbx_strand_id"].as(); + cif::replace_all(chains, ",", ", "); + cmpnd.push_back("CHAIN: " + chains); + } + + std::string fragment = r["pdbx_fragment"].as(); + if (not fragment.empty()) + cmpnd.push_back("FRAGMENT: " + fragment); + + for (auto sr : mDB["entity_name_com"].find("entity_id"_key == entityID)) + { + std::string syn = sr["name"].as(); + if (not syn.empty()) + cmpnd.push_back("SYNONYM: " + syn); + } + + std::string mutation = r["pdbx_mutation"].as(); + if (not mutation.empty()) + cmpnd.push_back("MUTATION: " + mutation); + + std::string ec = r["pdbx_ec"].as(); + if (not ec.empty()) + cmpnd.push_back("EC: " + ec); + + if (r["src_method"] == "man" or r["src_method"] == "syn") + cmpnd.push_back("ENGINEERED: YES"); + + std::string details = r["details"].as(); + if (not details.empty()) + cmpnd.push_back("OTHER_DETAILS: " + details); + } + + return FixStringLength("COMPND " + cif::join(cmpnd, "; "), kTruncateAt); +} + +std::string DSSP_impl::GetPDBSOURCELine() +{ + // SOURCE + + using namespace cif::literals; + + int molID = 0; + std::vector source; + + for (auto r : mDB["entity"]) + { + if (r["type"] != "polymer") + continue; + + std::string entityID = r["id"].as(); + + ++molID; + source.push_back("MOL_ID: " + std::to_string(molID)); + + if (r["src_method"] == "syn") + source.push_back("SYNTHETIC: YES"); + + auto &gen = mDB["entity_src_gen"]; + const std::pair kGenSourceMapping[] = { + { "gene_src_common_name", "ORGANISM_COMMON" }, + { "pdbx_gene_src_gene", "GENE" }, + { "gene_src_strain", "STRAIN" }, + { "pdbx_gene_src_cell_line", "CELL_LINE" }, + { "pdbx_gene_src_organelle", "ORGANELLE" }, + { "pdbx_gene_src_cellular_location", "CELLULAR_LOCATION" }, + { "pdbx_gene_src_scientific_name", "ORGANISM_SCIENTIFIC" }, + { "pdbx_gene_src_ncbi_taxonomy_id", "ORGANISM_TAXID" }, + { "pdbx_host_org_scientific_name", "EXPRESSION_SYSTEM" }, + { "pdbx_host_org_ncbi_taxonomy_id", "EXPRESSION_SYSTEM_TAXID" }, + { "pdbx_host_org_strain", "EXPRESSION_SYSTEM_STRAIN" }, + { "pdbx_host_org_variant", "EXPRESSION_SYSTEM_VARIANT" }, + { "pdbx_host_org_cellular_location", "EXPRESSION_SYSTEM_CELLULAR_LOCATION" }, + { "pdbx_host_org_vector_type", "EXPRESSION_SYSTEM_VECTOR_TYPE" }, + { "pdbx_host_org_vector", "EXPRESSION_SYSTEM_VECTOR" }, + { "pdbx_host_org_gene", "EXPRESSION_SYSTEM_GENE" }, + { "plasmid_name", "EXPRESSION_SYSTEM_PLASMID" } + }; + + for (auto gr : gen.find("entity_id"_key == entityID)) + { + for (auto m : kGenSourceMapping) + { + std::string cname, sname; + tie(cname, sname) = m; + + std::string s = gr[cname].as(); + if (not s.empty()) + source.push_back(sname + ": " + s); + } + } + + auto &nat = mDB["entity_src_nat"]; + const std::pair kNatSourceMapping[] = { + { "common_name", "ORGANISM_COMMON" }, + { "strain", "STRAIN" }, + { "pdbx_organism_scientific", "ORGANISM_SCIENTIFIC" }, + { "pdbx_ncbi_taxonomy_id", "ORGANISM_TAXID" }, + { "pdbx_cellular_location", "CELLULAR_LOCATION" }, + { "pdbx_plasmid_name", "PLASMID" }, + { "pdbx_organ", "ORGAN" }, + { "details", "OTHER_DETAILS" } + }; + + for (auto nr : nat.find("entity_id"_key == entityID)) + { + for (auto m : kNatSourceMapping) + { + std::string cname, sname; + tie(cname, sname) = m; + + std::string s = nr[cname].as(); + if (not s.empty()) + source.push_back(sname + ": " + s); + } + } + } + + return FixStringLength("SOURCE " + cif::join(source, "; "), kTruncateAt); +} + +std::string DSSP_impl::GetPDBAUTHORLine() +{ + // AUTHOR + std::vector author; + for (auto r : mDB["audit_author"]) + author.push_back(cif2pdbAuth(r["name"].as())); + + return FixStringLength("AUTHOR " + cif::join(author, "; "), kTruncateAt); +} + +// -------------------------------------------------------------------- + +std::string dssp::residue_info::asym_id() const +{ + return m_impl->mAsymID; +} + +std::string dssp::residue_info::compound_id() const +{ + return m_impl->mCompoundID; +} + +char dssp::residue_info::compound_letter() const +{ + return MapResidue(compound_id()); +} + +int dssp::residue_info::seq_id() const +{ + return m_impl->mSeqID; +} + +std::string dssp::residue_info::alt_id() const +{ + return m_impl->mAltID; +} + +std::string dssp::residue_info::auth_asym_id() const +{ + return m_impl->mAuthAsymID; +} + +int dssp::residue_info::auth_seq_id() const +{ + return m_impl->mAuthSeqID; +} + +std::string dssp::residue_info::pdb_strand_id() const +{ + return m_impl->mPDBStrandID; +} + +int dssp::residue_info::pdb_seq_num() const +{ + return m_impl->mPDBSeqNum; +} + +std::string dssp::residue_info::pdb_ins_code() const +{ + return m_impl->mPDBInsCode; +} + +std::optional dssp::residue_info::alpha() const +{ + return m_impl->mAlpha; +} + +std::optional dssp::residue_info::kappa() const +{ + return m_impl->mKappa; +} + +std::optional dssp::residue_info::omega() const +{ + return m_impl->mOmega; +} + +std::optional dssp::residue_info::phi() const +{ + return m_impl->mPhi; +} + +std::optional dssp::residue_info::psi() const +{ + return m_impl->mPsi; +} + +std::optional dssp::residue_info::tco() const +{ + return m_impl->mTCO; +} + +bool dssp::residue_info::is_pre_pro() const +{ + return m_impl->mType != kProline and m_impl->mNext != nullptr and m_impl->mNext->mType == kProline; +} + +float dssp::residue_info::chiral_volume() const +{ + return m_impl->mChiralVolume; +} + +const std::map> kChiAtomsMap = { + { MapResidue("ASP"), { "CG", "OD1" } }, + { MapResidue("ASN"), { "CG", "OD1" } }, + { MapResidue("ARG"), { "CG", "CD", "NE", "CZ" } }, + { MapResidue("HIS"), { "CG", "ND1" } }, + { MapResidue("GLN"), { "CG", "CD", "OE1" } }, + { MapResidue("GLU"), { "CG", "CD", "OE1" } }, + { MapResidue("SER"), { "OG" } }, + { MapResidue("THR"), { "OG1" } }, + { MapResidue("LYS"), { "CG", "CD", "CE", "NZ" } }, + { MapResidue("TYR"), { "CG", "CD1" } }, + { MapResidue("PHE"), { "CG", "CD1" } }, + { MapResidue("LEU"), { "CG", "CD1" } }, + { MapResidue("TRP"), { "CG", "CD1" } }, + { MapResidue("CYS"), { "SG" } }, + { MapResidue("ILE"), { "CG1", "CD1" } }, + { MapResidue("MET"), { "CG", "SD", "CE" } }, + { MapResidue("MSE"), { "CG", "SE", "CE" } }, + { MapResidue("PRO"), { "CG", "CD" } }, + { MapResidue("VAL"), { "CG1" } } +}; + +std::size_t dssp::residue_info::nr_of_chis() const +{ + auto i = kChiAtomsMap.find(m_impl->mType); + + return i != kChiAtomsMap.end() ? i->second.size() : 0; +} + +float dssp::residue_info::chi(std::size_t index) const +{ + float result = 0; + + auto type = m_impl->mType; + + auto i = kChiAtomsMap.find(type); + if (i != kChiAtomsMap.end() and index < i->second.size()) + { + std::vector atoms{ "N", "CA", "CB" }; + + atoms.insert(atoms.end(), i->second.begin(), i->second.end()); + + // in case we have a positive chiral volume we need to swap atoms + if (m_impl->mChiralVolume > 0) + { + if (type == kLeucine) + atoms.back() = "CD2"; + if (type == kValine) + atoms.back() = "CG2"; + } + + result = static_cast(dihedral_angle( + m_impl->get_atom(atoms[index + 0]), + m_impl->get_atom(atoms[index + 1]), + m_impl->get_atom(atoms[index + 2]), + m_impl->get_atom(atoms[index + 3]))); + } + + return result; +} + +std::tuple dssp::residue_info::ca_location() const +{ + return { m_impl->mCAlpha.mX, m_impl->mCAlpha.mY, m_impl->mCAlpha.mZ }; +} + +chain_break_type dssp::residue_info::chain_break() const +{ + return m_impl->mChainBreak; +} + +int dssp::residue_info::nr() const +{ + return m_impl->mNumber; +} + +structure_type dssp::residue_info::type() const +{ + return m_impl->mSecondaryStructure; +} + +int dssp::residue_info::ssBridgeNr() const +{ + return m_impl->mSSBridgeNr; +} + +helix_position_type dssp::residue_info::helix(helix_type helixType) const +{ + return m_impl->GetHelixFlag(helixType); +} + +bool dssp::residue_info::is_alpha_helix_end_before_start() const +{ + bool result = false; + + if (m_impl->mNext != nullptr) + result = m_impl->GetHelixFlag(helix_type::alpha) == helix_position_type::End and m_impl->mNext->GetHelixFlag(helix_type::alpha) == helix_position_type::Start; + + return result; +} + +bool dssp::residue_info::bend() const +{ + return m_impl->IsBend(); +} + +double dssp::residue_info::accessibility() const +{ + return m_impl->mAccessibility; +} + +std::tuple dssp::residue_info::bridge_partner(int i) const +{ + auto bp = m_impl->GetBetaPartner(i); + + residue_info ri(bp.m_residue); + + return std::make_tuple(std::move(ri), bp.ladder, bp.parallel); +} + +int dssp::residue_info::sheet() const +{ + return m_impl->GetSheet(); +} + +int dssp::residue_info::strand() const +{ + return m_impl->GetStrand(); +} + +std::tuple dssp::residue_info::acceptor(int i) const +{ + auto &a = m_impl->mHBondAcceptor[i]; + return { residue_info(a.res), a.energy }; +} + +std::tuple dssp::residue_info::donor(int i) const +{ + auto &d = m_impl->mHBondDonor[i]; + return { residue_info(d.res), d.energy }; +} + +dssp::residue_info dssp::residue_info::next() const +{ + return residue_info(m_impl ? m_impl->mNext : nullptr); +} + +// -------------------------------------------------------------------- + +dssp::iterator::iterator(residue *res) + : m_current(res) +{ +} + +dssp::iterator &dssp::iterator::operator++() +{ + ++m_current.m_impl; + return *this; +} + +dssp::iterator &dssp::iterator::operator--() +{ + --m_current.m_impl; + return *this; +} + +// -------------------------------------------------------------------- + +dssp::dssp(const cif::mm::structure &s, int min_poly_proline_stretch_length, bool calculateSurfaceAccessibility) + : dssp(s.get_datablock(), static_cast(s.get_model_nr()), min_poly_proline_stretch_length, calculateSurfaceAccessibility) +{ +} + +dssp::dssp(const cif::datablock &db, int model_nr, int min_poly_proline_stretch, bool calculateSurfaceAccessibility) + : m_impl(new DSSP_impl(db, model_nr, min_poly_proline_stretch)) +{ + if (calculateSurfaceAccessibility) + { + std::thread t(std::bind(&DSSP_impl::calculateSurface, m_impl)); + m_impl->calculateSecondaryStructure(); + t.join(); + } + else + m_impl->calculateSecondaryStructure(); +} + +dssp::~dssp() +{ + delete m_impl; +} + +dssp::iterator dssp::begin() const +{ + return iterator(m_impl->mResidues.empty() ? nullptr : m_impl->mResidues.data()); +} + +dssp::iterator dssp::end() const +{ + // careful now, MSVC is picky when it comes to dereferencing iterators that are at the end. + residue *res = nullptr; + if (not m_impl->mResidues.empty()) + { + res = m_impl->mResidues.data(); + res += m_impl->mResidues.size(); + } + + return iterator(res); +} + +dssp::residue_info dssp::operator[](const key_type &key) const +{ + auto i = std::find_if(begin(), end(), + [key](const residue_info &res) + { return res.asym_id() == std::get<0>(key) and res.seq_id() == std::get<1>(key); }); + + if (i == end()) + throw std::out_of_range("Could not find residue with supplied key"); + + return *i; +} + +dssp::statistics dssp::get_statistics() const +{ + return m_impl->mStats; +} + +std::string dssp::get_pdb_header_line(pdb_record_type pdb_record) const +{ + switch (pdb_record) + { + case pdb_record_type::HEADER: + return m_impl->GetPDBHEADERLine(); + case pdb_record_type::COMPND: + return m_impl->GetPDBCOMPNDLine(); + case pdb_record_type::SOURCE: + return m_impl->GetPDBSOURCELine(); + case pdb_record_type::AUTHOR: + return m_impl->GetPDBAUTHORLine(); + default: + return {}; + } +} + +// -------------------------------------------------------------------- + +void dssp::write_legacy_output(std::ostream &os) const +{ + writeDSSP(*this, os); +} + +void dssp::annotate(cif::datablock &db, bool writeOther, bool writeDSSPCategories) const +{ + annotateDSSP(db, *this, writeOther, writeDSSPCategories); +} --- /dev/null +++ b/src/dssp.hpp @@ -0,0 +1,283 @@ +/*- + * 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. + */ + +#pragma once + +/// \file DSSP.hpp +/// Calculate DSSP-like secondary structure information. + +#include + +#include + +class dssp +{ + public: + 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 + }; + + static constexpr 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 + }; + + dssp(const cif::datablock &db, int model_nr, int min_poly_proline_stretch_length, bool calculateSurfaceAccessibility); + dssp(const cif::mm::structure &s, 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::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; + char compound_letter() const; // Single letter for residue compound type, or 'X' in case it is not known + + 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; + + std::optional alpha() const; + std::optional kappa() const; + std::optional phi() const; + std::optional psi() const; + std::optional tco() const; + std::optional omega() const; + + bool is_pre_pro() const; + bool is_cis() const { return std::abs(omega().value_or(360)) < 30.0f; } + + float chiral_volume() const; + + std::size_t nr_of_chis() const; + float chi(std::size_t index) const; + + std::vector chis() const + { + std::vector result; + for (size_t i = 0; i < nr_of_chis(); ++i) + result.push_back(chi(i)); + return result; + } + + std::tuple 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 bridge_partner(int i) const; + + int sheet() const; + int strand() const; + + /// \brief return resinfo and the energy of the bond + std::tuple acceptor(int i) const; + std::tuple 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); + + residue_info next() const; + + 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; + + iterator begin() const; + iterator end() const; + + residue_info operator[](const key_type &key) const; + + bool empty() const { return begin() == end(); } + + // -------------------------------------------------------------------- + // Writing out the data, either in legacy format... + + void write_legacy_output(std::ostream &os) const; + + // ... or as annotation in the cif::datablock + void annotate(cif::datablock &db, bool writeOther, bool writeDSSPCategories) const; + + // 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; +}; --- a/src/tortoize.cpp +++ b/src/tortoize.cpp @@ -27,7 +27,7 @@ #include "tortoize.hpp" #include "revision.hpp" -#include +#include "dssp.hpp" #include #include @@ -1092,4 +1092,4 @@ } return data; -} \ No newline at end of file +} --- /dev/null +++ b/src/dssp-io.cpp @@ -0,0 +1,897 @@ +/*- + * 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. + */ + +#include "dssp-io.hpp" +#include "dssp-revision.hpp" + +#include +#include + +#include +#include +#include +#include + +// -------------------------------------------------------------------- + +std::string ResidueToDSSPLine(const dssp::residue_info &info) +{ + /* + This is the header line for the residue lines in a DSSP file: + + # RESIDUE AA STRUCTURE BP1 BP2 ACC N-H-->O O-->H-N N-H-->O O-->H-N TCO KAPPA ALPHA PHI PSI X-CA Y-CA Z-CA + */ + + // auto& residue = info.residue(); + auto &residue = info; + + if (residue.pdb_strand_id().length() > 1) + throw std::runtime_error("This file contains data that won't fit in the original DSSP format"); + + char code = residue.compound_letter(); + + if (code == 'C') // a cysteine + { + auto ssbridgenr = info.ssBridgeNr(); + if (ssbridgenr) + code = 'a' + ((ssbridgenr - 1) % 26); + } + + char ss = ' '; + switch (info.type()) + { + case dssp::structure_type::Alphahelix: ss = 'H'; break; + case dssp::structure_type::Betabridge: ss = 'B'; break; + case dssp::structure_type::Strand: ss = 'E'; break; + case dssp::structure_type::Helix_3: ss = 'G'; break; + case dssp::structure_type::Helix_5: ss = 'I'; break; + case dssp::structure_type::Helix_PPII: ss = 'P'; break; + case dssp::structure_type::Turn: ss = 'T'; break; + case dssp::structure_type::Bend: ss = 'S'; break; + case dssp::structure_type::Loop: ss = ' '; break; + } + + char helix[4] = { ' ', ' ', ' ', ' ' }; + for (dssp::helix_type helixType : { dssp::helix_type::_3_10, dssp::helix_type::alpha, dssp::helix_type::pi, dssp::helix_type::pp }) + { + switch (info.helix(helixType)) + { + case dssp::helix_position_type::None: helix[static_cast(helixType)] = ' '; break; + case dssp::helix_position_type::Start: helix[static_cast(helixType)] = '>'; break; + case dssp::helix_position_type::End: helix[static_cast(helixType)] = '<'; break; + case dssp::helix_position_type::StartAndEnd: helix[static_cast(helixType)] = 'X'; break; + case dssp::helix_position_type::Middle: helix[static_cast(helixType)] = (helixType == dssp::helix_type::pp ? 'P' : static_cast('3' + static_cast(helixType))); break; + } + } + + char bend = ' '; + if (info.bend()) + bend = 'S'; + + double alpha = residue.alpha().value_or(360); + char chirality = alpha == 360 ? ' ' : (alpha < 0 ? '-' : '+'); + + uint32_t bp[2] = {}; + char bridgelabel[2] = { ' ', ' ' }; + for (uint32_t i : { 0, 1 }) + { + const auto &[p, ladder, parallel] = info.bridge_partner(i); + if (not p) + continue; + + bp[i] = p.nr() % 10000; // won't fit otherwise... + bridgelabel[i] = (parallel ? 'a' : 'A') + ladder % 26; + } + + char sheet = ' '; + if (info.sheet() != 0) + sheet = 'A' + (info.sheet() - 1) % 26; + + std::string NHO[2], ONH[2]; + for (int i : { 0, 1 }) + { + const auto &[donor, donorE] = info.donor(i); + const auto &[acceptor, acceptorE] = info.acceptor(i); + + NHO[i] = ONH[i] = "0, 0.0"; + + if (acceptor) + { + auto d = acceptor.nr() - info.nr(); + NHO[i] = cif::format("{:d},{:3.1f}", d, acceptorE); + } + + if (donor) + { + auto d = donor.nr() - info.nr(); + ONH[i] = cif::format("{:d},{:3.1f}", d, donorE); + } + } + + // auto ca = residue.atomByID("CA"); + auto const &[cax, cay, caz] = residue.ca_location(); + + return cif::format("{:5d}{:5d}{:1.1s}{:1.1s} {:1c} {:1c}{:1c}{:1c}{:1c}{:1c}{:1c}{:1c}{:1c}{:1c}{:4d}{:4d}{:1c}{:4.0f} {:>11s}{:>11s}{:>11s}{:>11s} {:6.3f}{:6.1f}{:6.1f}{:6.1f}{:6.1f} {:6.1f} {:6.1f} {:6.1f}", + info.nr(), residue.pdb_seq_num(), residue.pdb_ins_code(), residue.pdb_strand_id(), code, + ss, helix[3], helix[0], helix[1], helix[2], bend, chirality, bridgelabel[0], bridgelabel[1], + bp[0], bp[1], sheet, floor(info.accessibility() + 0.5), + NHO[0], ONH[0], NHO[1], ONH[1], + residue.tco().value_or(0), + residue.kappa().value_or(360), + residue.alpha().value_or(360), + residue.phi().value_or(360), + residue.psi().value_or(360), + cax, cay, caz); +} + +void writeDSSP(const dssp &dssp, std::ostream &os) +{ + using namespace std::chrono; + + auto stats = dssp.get_statistics(); + + std::time_t today = system_clock::to_time_t(system_clock::now()); + std::tm *tm = std::gmtime(&today); + + std::string version = klibdsspVersionNumber; + if (version.length() < 10) + version.insert(version.end(), 10 - version.length(), ' '); + + os << "==== Secondary Structure Definition by the program DSSP, NKI version " << version << " ==== DATE=" << std::put_time(tm, "%F") << " ." << std::endl + << "REFERENCE M.L. HEKKELMAN ET AL, PROTEIN SCIENCE 34.8 (2025) e70208; W. KABSCH AND C.SANDER, BIOPOLYMERS 22 (1983) 2577-2637    ." << std::endl + << dssp.get_pdb_header_line(dssp::pdb_record_type::HEADER) << '.' << std::endl + << dssp.get_pdb_header_line(dssp::pdb_record_type::COMPND) << '.' << std::endl + << dssp.get_pdb_header_line(dssp::pdb_record_type::SOURCE) << '.' << std::endl + << dssp.get_pdb_header_line(dssp::pdb_record_type::AUTHOR) << '.' << std::endl; + + os << cif::format("{:5d}{:3d}{:3d}{:3d}{:3d} TOTAL NUMBER OF RESIDUES, NUMBER OF CHAINS, NUMBER OF SS-BRIDGES(TOTAL,INTRACHAIN,INTERCHAIN) .", + stats.count.residues, stats.count.chains, stats.count.SS_bridges, stats.count.intra_chain_SS_bridges, (stats.count.SS_bridges - stats.count.intra_chain_SS_bridges)) + << std::endl; + + os << cif::format("{:8.1f} ACCESSIBLE SURFACE OF PROTEIN (ANGSTROM**2) .", stats.accessible_surface) << std::endl; + + // hydrogenbond summary + + os << cif::format("{:5d}{:5.1f} TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(J) , SAME NUMBER PER 100 RESIDUES .", stats.count.H_bonds, (stats.count.H_bonds * 100.0 / stats.count.residues)) << std::endl; + + os << cif::format("{:5d}{:5.1f} TOTAL NUMBER OF HYDROGEN BONDS IN PARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES .", stats.count.H_bonds_in_parallel_bridges, (stats.count.H_bonds_in_parallel_bridges * 100.0 / stats.count.residues)) << std::endl; + + os << cif::format("{:5d}{:5.1f} TOTAL NUMBER OF HYDROGEN BONDS IN ANTIPARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES .", stats.count.H_bonds_in_antiparallel_bridges, (stats.count.H_bonds_in_antiparallel_bridges * 100.0 / stats.count.residues)) << std::endl; + + for (int k = 0; k < 11; ++k) + os << cif::format("{:5d}{:5.1f} TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(I{:1c}{:1d}), SAME NUMBER PER 100 RESIDUES .", stats.count.H_Bonds_per_distance[k], (stats.count.H_Bonds_per_distance[k] * 100.0 / stats.count.residues), (k - 5 < 0 ? '-' : '+'), abs(k - 5)) << std::endl; + + // histograms... + os << " 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 *** HISTOGRAMS OF *** ." << std::endl; + + for (auto hi : stats.histogram.residues_per_alpha_helix) + os << cif::format("{:3d}", hi); + os << " RESIDUES PER ALPHA HELIX ." << std::endl; + + for (auto hi : stats.histogram.parallel_bridges_per_ladder) + os << cif::format("{:3d}", hi); + os << " PARALLEL BRIDGES PER LADDER ." << std::endl; + + for (auto hi : stats.histogram.antiparallel_bridges_per_ladder) + os << cif::format("{:3d}", hi); + os << " ANTIPARALLEL BRIDGES PER LADDER ." << std::endl; + + for (auto hi : stats.histogram.ladders_per_sheet) + os << cif::format("{:3d}", hi); + os << " LADDERS PER SHEET ." << std::endl; + + // per residue information + + os << " # RESIDUE AA STRUCTURE BP1 BP2 ACC N-H-->O O-->H-N N-H-->O O-->H-N TCO KAPPA ALPHA PHI PSI X-CA Y-CA Z-CA" << std::endl; + + int last = 0; + for (auto ri : dssp) + { + // insert a break line whenever we detect missing residues + // can be the transition to a different chain, or missing residues in the current chain + + if (ri.nr() != last + 1) + os << cif::format("{:5d} !{:1c} 0 0 0 0, 0.0 0, 0.0 0, 0.0 0, 0.0 0.000 360.0 360.0 360.0 360.0 0.0 0.0 0.0", + (last + 1), (ri.chain_break() == dssp::chain_break_type::NewChain ? '*' : ' ')) + << std::endl; + + os << ResidueToDSSPLine(ri) << std::endl; + last = ri.nr(); + } +} + +// -------------------------------------------------------------------- + +void writeBridgePairs(cif::datablock &db, const dssp &dssp) +{ + auto &hb = db["dssp_struct_bridge_pairs"]; + + hb.add_item("id"); + hb.add_item("label_comp_id"); + hb.add_item("label_seq_id"); + hb.add_item("label_asym_id"); + hb.add_item("auth_seq_id"); + hb.add_item("auth_asym_id"); + hb.add_item("pdbx_PDB_ins_code"); + + // force right order + for (std::string da : { "acceptor_", "donor_" }) + { + for (std::string i : { "1_", "2_" }) + { + for (std::string n : { "label_comp_id", "label_seq_id", "label_asym_id", "auth_seq_id", "auth_asym_id", "pdbx_PDB_ins_code", "energy" }) + hb.add_item(da + i + n); + } + } + + for (auto &res : dssp) + { + cif::row_initializer data({ { "id", hb.get_unique_id("") }, + { "label_comp_id", res.compound_id() }, + { "label_seq_id", res.seq_id() }, + { "label_asym_id", res.asym_id() }, + // { "auth_comp_id", res.compound_id() }, + { "auth_seq_id", res.auth_seq_id() }, + { "auth_asym_id", res.auth_asym_id() }, + { "pdbx_PDB_ins_code", res.pdb_ins_code() } }); + + for (int i : { 0, 1 }) + { + const auto &&[acceptor, acceptorEnergy] = res.acceptor(i); + const auto &&[donor, donorEnergy] = res.donor(i); + + if (acceptor) + { + if (i == 0) + { + data.emplace_back("acceptor_1_label_comp_id", acceptor.compound_id()); + data.emplace_back("acceptor_1_label_seq_id", acceptor.seq_id()); + data.emplace_back("acceptor_1_label_asym_id", acceptor.asym_id()); + // data.emplace_back("acceptor_1_auth_comp_id", acceptor.compound_id()); + data.emplace_back("acceptor_1_auth_seq_id", acceptor.auth_seq_id()); + data.emplace_back("acceptor_1_auth_asym_id", acceptor.auth_asym_id()); + data.emplace_back("acceptor_1_pdbx_PDB_ins_code", acceptor.pdb_ins_code()); + data.emplace_back("acceptor_1_energy", acceptorEnergy, 1); + } + else + { + data.emplace_back("acceptor_2_label_comp_id", acceptor.compound_id()); + data.emplace_back("acceptor_2_label_seq_id", acceptor.seq_id()); + data.emplace_back("acceptor_2_label_asym_id", acceptor.asym_id()); + // data.emplace_back("acceptor_2_auth_comp_id", acceptor.compound_id()); + data.emplace_back("acceptor_2_auth_seq_id", acceptor.auth_seq_id()); + data.emplace_back("acceptor_2_auth_asym_id", acceptor.auth_asym_id()); + data.emplace_back("acceptor_2_pdbx_PDB_ins_code", acceptor.pdb_ins_code()); + data.emplace_back("acceptor_2_energy", acceptorEnergy, 1); + } + } + + if (donor) + { + if (i == 0) + { + data.emplace_back("donor_1_label_comp_id", donor.compound_id()); + data.emplace_back("donor_1_label_seq_id", donor.seq_id()); + data.emplace_back("donor_1_label_asym_id", donor.asym_id()); + // data.emplace_back("donor_1_auth_comp_id", donor.compound_id()); + data.emplace_back("donor_1_auth_seq_id", donor.auth_seq_id()); + data.emplace_back("donor_1_auth_asym_id", donor.auth_asym_id()); + data.emplace_back("donor_1_pdbx_PDB_ins_code", donor.pdb_ins_code()); + data.emplace_back("donor_1_energy", donorEnergy, 1); + } + else + { + data.emplace_back("donor_2_label_comp_id", donor.compound_id()); + data.emplace_back("donor_2_label_seq_id", donor.seq_id()); + data.emplace_back("donor_2_label_asym_id", donor.asym_id()); + // data.emplace_back("donor_2_auth_comp_id", donor.compound_id()); + data.emplace_back("donor_2_auth_seq_id", donor.auth_seq_id()); + data.emplace_back("donor_2_auth_asym_id", donor.auth_asym_id()); + data.emplace_back("donor_2_pdbx_PDB_ins_code", donor.pdb_ins_code()); + data.emplace_back("donor_2_energy", donorEnergy, 1); + } + } + } + + hb.emplace(std::move(data)); + } +} + +void writeSheets(cif::datablock &db, const dssp &dssp) +{ + using namespace cif::literals; + + using res_list = std::vector; + + // clean up old info first + + for (auto sheet_cat : { "struct_sheet", "struct_sheet_order", "struct_sheet_range", "struct_sheet_hbond", "pdbx_struct_sheet_hbond" }) + db.erase(remove_if(db.begin(), db.end(), [sheet_cat](const cif::category &cat) + { return cat.name() == sheet_cat; }), + db.end()); + + // create a list of strands, based on the SS info in DSSP. Store sheet number along with the strand. + + std::map, res_list> strands; + std::set sheetNrs; + + for (auto &res : dssp) + { + if (res.type() != dssp::structure_type::Strand and res.type() != dssp::structure_type::Betabridge) + continue; + + strands[{ res.sheet(), res.strand() }].emplace_back(res); + sheetNrs.insert(res.sheet()); + } + + // -------------------------------------------------------------------- + + auto &struct_sheet = db["struct_sheet"]; + auto &struct_sheet_range = db["struct_sheet_range"]; + + for (int sheetNr : sheetNrs) + { + auto sheetID = cif::cif_id_for_number(sheetNr - 1); + + struct_sheet.emplace({ { "id", sheetID }, + { "number_strands", + std::count_if(strands.begin(), strands.end(), [nr = sheetNr](std::tuple, res_list> const &s) + { + const auto &[strandID, strand] = s; + return strand.front().sheet() == nr; }) } }); + + for (auto &&[strandTuple, strand] : strands) + { + if (strand.front().sheet() != sheetNr) + continue; + + std::string strandID = cif::cif_id_for_number(strand.front().strand() - 1); + + std::sort(strand.begin(), strand.end(), [](dssp::residue_info const &a, dssp::residue_info const &b) + { return a.nr() < b.nr(); }); + + auto &beg = strand.front(); + auto &end = strand.back(); + + struct_sheet_range.emplace({ { "sheet_id", sheetID }, + { "id", strandID }, + { "beg_label_comp_id", beg.compound_id() }, + { "beg_label_asym_id", beg.asym_id() }, + { "beg_label_seq_id", beg.seq_id() }, + { "pdbx_beg_PDB_ins_code", beg.pdb_ins_code() }, + { "end_label_comp_id", end.compound_id() }, + { "end_label_asym_id", end.asym_id() }, + { "end_label_seq_id", end.seq_id() }, + { "pdbx_end_PDB_ins_code", end.pdb_ins_code() }, + { "beg_auth_comp_id", beg.compound_id() }, + { "beg_auth_asym_id", beg.auth_asym_id() }, + { "beg_auth_seq_id", beg.auth_seq_id() }, + { "end_auth_comp_id", end.compound_id() }, + { "end_auth_asym_id", end.auth_asym_id() }, + { "end_auth_seq_id", end.auth_seq_id() } }); + } + } +} + +void writeLadders(cif::datablock &db, const dssp &dssp) +{ + // Write out the DSSP ladders + struct ladder_info + { + ladder_info(int label, int sheet, bool parallel, const dssp::residue_info &a, const dssp::residue_info &b) + : ladder(label) + , sheet(sheet) + , parallel(parallel) + , pairs({ { a, b } }) + { + } + + bool operator<(const ladder_info &rhs) const + { + return ladder < rhs.ladder; + } + + int ladder; + int sheet; + bool parallel; + std::vector> pairs; + }; + + std::vector ladders; + + for (auto res : dssp) + { + for (int i : { 0, 1 }) + { + const auto &[p, ladder, parallel] = res.bridge_partner(i); + + if (not p) + continue; + + bool is_new = true; + for (auto &l : ladders) + { + if (l.ladder != ladder) + continue; + + assert(l.parallel == parallel); + + if (find_if(l.pairs.begin(), l.pairs.end(), [na = p.nr(), nb = res.nr()](const auto &p) + { return p.first.nr() == na and p.second.nr() == nb; }) != l.pairs.end()) + { + is_new = false; + break; + } + + l.pairs.emplace_back(res, p); + is_new = false; + break; + } + + if (not is_new) + continue; + + ladders.emplace_back(ladder, res.sheet() - 1, parallel, res, p); + } + } + + std::sort(ladders.begin(), ladders.end()); + + auto &dssp_struct_ladder = db["dssp_struct_ladder"]; + + for (const auto &l : ladders) + { + const auto &[beg1, beg2] = l.pairs.front(); + const auto &[end1, end2] = l.pairs.back(); + + dssp_struct_ladder.emplace({ { "id", cif::cif_id_for_number(l.ladder) }, + { "sheet_id", cif::cif_id_for_number(l.sheet) }, + { "range_id_1", cif::cif_id_for_number(beg1.strand() - 1) }, + { "range_id_2", cif::cif_id_for_number(beg2.strand() - 1) }, + { "type", l.parallel ? "parallel" : "anti-parallel" }, + + { "beg_1_label_comp_id", beg1.compound_id() }, + { "beg_1_label_asym_id", beg1.asym_id() }, + { "beg_1_label_seq_id", beg1.seq_id() }, + { "pdbx_beg_1_PDB_ins_code", beg1.pdb_ins_code() }, + { "end_1_label_comp_id", end1.compound_id() }, + { "end_1_label_asym_id", end1.asym_id() }, + { "end_1_label_seq_id", end1.seq_id() }, + { "pdbx_end_1_PDB_ins_code", end1.pdb_ins_code() }, + { "beg_1_auth_comp_id", beg1.compound_id() }, + { "beg_1_auth_asym_id", beg1.auth_asym_id() }, + { "beg_1_auth_seq_id", beg1.auth_seq_id() }, + { "end_1_auth_comp_id", end1.compound_id() }, + { "end_1_auth_asym_id", end1.auth_asym_id() }, + { "end_1_auth_seq_id", end1.auth_seq_id() }, + + { "beg_2_label_comp_id", beg2.compound_id() }, + { "beg_2_label_asym_id", beg2.asym_id() }, + { "beg_2_label_seq_id", beg2.seq_id() }, + { "pdbx_beg_2_PDB_ins_code", beg2.pdb_ins_code() }, + { "end_2_label_comp_id", end2.compound_id() }, + { "end_2_label_asym_id", end2.asym_id() }, + { "end_2_label_seq_id", end2.seq_id() }, + { "pdbx_end_2_PDB_ins_code", end2.pdb_ins_code() }, + { "beg_2_auth_comp_id", beg2.compound_id() }, + { "beg_2_auth_asym_id", beg2.auth_asym_id() }, + { "beg_2_auth_seq_id", beg2.auth_seq_id() }, + { "end_2_auth_comp_id", end2.compound_id() }, + { "end_2_auth_asym_id", end2.auth_asym_id() }, + { "end_2_auth_seq_id", end2.auth_seq_id() } }); + } +} + +void writeStatistics(cif::datablock &db, const dssp &dssp) +{ + using namespace std::literals; + + auto stats = dssp.get_statistics(); + + auto &dssp_statistics = db["dssp_statistics"]; + + std::optional surface_accessibility; + if (stats.accessible_surface > 0) + surface_accessibility = stats.accessible_surface; + + auto stats_i = dssp_statistics.emplace({ // + { "entry_id", db.name() }, + { "nr_of_residues", stats.count.residues }, + { "nr_of_chains", stats.count.chains }, + { "nr_of_ss_bridges_total", stats.count.SS_bridges }, + { "nr_of_ss_bridges_intra_chain", stats.count.intra_chain_SS_bridges }, + { "nr_of_ss_bridges_inter_chain", stats.count.SS_bridges - stats.count.intra_chain_SS_bridges }, + { "accessible_surface_of_protein", surface_accessibility } }); + + auto &dssp_struct_hbonds = db["dssp_statistics_hbond"]; + + dssp_struct_hbonds.emplace({ { "entry_id", db.name() }, + { "type", "O(I)-->H-N(J)" }, + { "count", stats.count.H_bonds }, + { "count_per_100", stats.count.H_bonds * 100.0 / stats.count.residues, 1 } }); + + dssp_struct_hbonds.emplace({ { "entry_id", db.name() }, + { "type", "PARALLEL BRIDGES" }, + { "count", stats.count.H_bonds_in_parallel_bridges }, + { "count_per_100", stats.count.H_bonds_in_parallel_bridges * 100.0 / stats.count.residues, 1 } }); + + dssp_struct_hbonds.emplace({ { "entry_id", db.name() }, + { "type", "ANTIPARALLEL BRIDGES" }, + { "count", stats.count.H_bonds_in_antiparallel_bridges }, + { "count_per_100", stats.count.H_bonds_in_antiparallel_bridges * 100.0 / stats.count.residues, 1 } }); + + for (int k = 0; k < 11; ++k) + dssp_struct_hbonds.emplace({ { "entry_id", db.name() }, + { "type", "O(I)-->H-N(I"s + char(k - 5 < 0 ? '-' : '+') + std::to_string(abs(k - 5)) + ")" }, + { "count", stats.count.H_Bonds_per_distance[k] }, + { "count_per_100", stats.count.H_Bonds_per_distance[k] * 100.0 / stats.count.residues, 1 } }); + + auto &dssp_statistics_histogram = db["dssp_statistics_histogram"]; + + using histogram_data_type = std::tuple; + for (const auto &[type, values] : std::vector{ + { "residues_per_alpha_helix", stats.histogram.residues_per_alpha_helix }, + { "parallel_bridges_per_ladder", stats.histogram.parallel_bridges_per_ladder }, + { "antiparallel_bridges_per_ladder", stats.histogram.antiparallel_bridges_per_ladder }, + { "ladders_per_sheet", stats.histogram.ladders_per_sheet } }) + { + auto hi = dssp_statistics_histogram.emplace({ + { "entry_id", db.name() }, + { "type", type }, + { "1", values[0] }, + { "2", values[1] }, + { "3", values[2] }, + { "4", values[3] }, + { "5", values[4] }, + { "6", values[5] }, + { "7", values[6] }, + { "8", values[7] }, + { "9", values[8] }, + { "10", values[9] }, + { "11", values[10] }, + { "12", values[11] }, + { "13", values[12] }, + { "14", values[13] }, + { "15", values[14] }, + { "16", values[15] }, + { "17", values[16] }, + { "18", values[17] }, + { "19", values[18] }, + { "20", values[19] }, + { "21", values[20] }, + { "22", values[21] }, + { "23", values[22] }, + { "24", values[23] }, + { "25", values[24] }, + { "26", values[25] }, + { "27", values[26] }, + { "28", values[27] }, + { "29", values[28] }, + { "30", values[29] }, + }); + } +} + +void writeSummary(cif::datablock &db, const dssp &dssp) +{ + bool writeAccessibility = dssp.get_statistics().accessible_surface > 0; + + // A approximation of the old format + + auto &dssp_struct_summary = db["dssp_struct_summary"]; + + // prime the category with the field labels we need, this is to ensure proper order in writing out the data. + + for (auto label : { "entry_id", "label_comp_id", "label_asym_id", "label_seq_id", "secondary_structure", + "ss_bridge", "helix_3_10", "helix_alpha", "helix_pi", "helix_pp", "bend", "chirality", "sheet", + "strand", "ladder_1", "ladder_2", "accessibility", "TCO", "kappa", "alpha", "phi", "psi", + "x_ca", "y_ca", "z_ca" }) + dssp_struct_summary.add_item(label); + + for (auto res : dssp) + { + + /* + This is the header line for the residue lines in a DSSP file: + + # RESIDUE AA STRUCTURE BP1 BP2 ACC N-H-->O O-->H-N N-H-->O O-->H-N TCO KAPPA ALPHA PHI PSI X-CA Y-CA Z-CA + */ + + std::string ss_bridge = "."; + if (res.ssBridgeNr() != 0) + ss_bridge = std::to_string(res.ssBridgeNr()); + + std::string ss = { res.type() == dssp::structure_type::Loop ? '.' : static_cast(res.type()) }; + + std::string helix[4] = { ".", ".", ".", "." }; + for (dssp::helix_type helixType : { dssp::helix_type::_3_10, dssp::helix_type::alpha, dssp::helix_type::pi, dssp::helix_type::pp }) + { + switch (res.helix(helixType)) + { + // case dssp::helix_position_type::None: helix[static_cast(helixType)] = ' '; break; + case dssp::helix_position_type::Start: + helix[static_cast(helixType)] = { '>' }; + break; + + case dssp::helix_position_type::End: + helix[static_cast(helixType)] = { '<' }; + break; + + case dssp::helix_position_type::StartAndEnd: + helix[static_cast(helixType)] = { 'X' }; + break; + + case dssp::helix_position_type::Middle: + if (helixType == dssp::helix_type::pp) + helix[static_cast(helixType)] = { 'P' }; + else + helix[static_cast(helixType)] = { static_cast('3' + static_cast(helixType)) }; + break; + + default: + break; + } + } + + std::string bend = "."; + if (res.bend()) + bend = "S"; + + std::string chirality = "."; + if (res.alpha().has_value()) + chirality = *res.alpha() < 0 ? "-" : "+"; + + std::string ladders[2] = { ".", "." }; + + for (uint32_t i : { 0, 1 }) + { + const auto &[p, ladder, parallel] = res.bridge_partner(i); + if (not p) + continue; + + ladders[i] = cif::cif_id_for_number(ladder); + } + + auto const &[cax, cay, caz] = res.ca_location(); + + cif::row_initializer data{ + { "entry_id", db.name() }, + { "label_comp_id", res.compound_id() }, + { "label_asym_id", res.asym_id() }, + { "label_seq_id", res.seq_id() }, + + { "secondary_structure", ss }, + + { "ss_bridge", ss_bridge }, + + { "helix_3_10", helix[0] }, + { "helix_alpha", helix[1] }, + { "helix_pi", helix[2] }, + { "helix_pp", helix[3] }, + + { "bend", bend }, + { "chirality", chirality }, + + { "sheet", res.sheet() ? cif::cif_id_for_number(res.sheet() - 1) : "." }, + { "strand", res.strand() ? cif::cif_id_for_number(res.strand() - 1) : "." }, + { "ladder_1", ladders[0] }, + { "ladder_2", ladders[1] }, + + { "x_ca", cax, 1 }, + { "y_ca", cay, 1 }, + { "z_ca", caz, 1 }, + }; + + if (writeAccessibility) + data.emplace_back("accessibility", res.accessibility(), 1); + + if (res.tco().has_value()) + data.emplace_back("TCO", *res.tco(), 3); + else + data.emplace_back("TCO", "."); + + if (res.kappa().has_value()) + data.emplace_back("kappa", *res.kappa(), 1); + else + data.emplace_back("kappa", "."); + + if (res.alpha().has_value()) + data.emplace_back("alpha", *res.alpha(), 1); + else + data.emplace_back("alpha", "."); + + if (res.phi().has_value()) + data.emplace_back("phi", *res.phi(), 1); + else + data.emplace_back("phi", "."); + + if (res.psi().has_value()) + data.emplace_back("psi", *res.psi(), 1); + else + data.emplace_back("psi", "."); + + dssp_struct_summary.emplace(std::move(data)); + } +} + +void annotateDSSP(cif::datablock &db, const dssp &dssp, bool writeOther, bool writeNewFormat) +{ + using namespace std::literals; + + auto &audit_conform = db["audit_conform"]; + + if (audit_conform.empty()) + { + auto &cf = cif::validator_factory::instance(); + cf.get("mmcif_pdbx.dic").fill_audit_conform(audit_conform); + } + + audit_conform.erase(cif::key("dict_name") == "dssp-extension.dic"); + audit_conform.emplace({ // + { "dict_name", "dssp-extension.dic" }, + { "dict_version", "1.1.1" }, + { "dict_location", "https://pdb-redo.eu/dssp/dssp-extension.dic" } }); + + // Re-load the dictionary + db.load_dictionary(); + + if (dssp.empty()) + { + if (cif::VERBOSE > 0) + std::cout << "No secondary structure information found" << std::endl; + } + else + { + for (auto cat : std::initializer_list{ + "dssp_struct_bridge_pairs", + "dssp_struct_ladder", + "dssp_statistics", + "dssp_statistics_hbond", + "dssp_statistics_histogram", + "dssp_struct_summary" + }) + { + db.erase(std::remove_if(db.begin(), db.end(), [cat] (cif::category &c) { return c.name() == cat; }), db.end()); + } + + if (writeNewFormat) + { + writeBridgePairs(db, dssp); + writeSheets(db, dssp); + writeLadders(db, dssp); + writeStatistics(db, dssp); + writeSummary(db, dssp); + } + + // replace all struct_conf and struct_conf_type records + auto &structConfType = db["struct_conf_type"]; + structConfType.clear(); + + auto &structConf = db["struct_conf"]; + structConf.clear(); + + std::map foundTypes; + + auto st = dssp.begin(), lt = st; + auto lastSS = st->type(); + + for (auto t = dssp.begin();; lt = t, ++t) + { + bool stop = t == dssp.end(); + + bool flush = (stop or t->type() != lastSS); + + if (flush and (writeOther or lastSS != dssp::structure_type::Loop)) + { + auto &rb = *st; + auto &re = *lt; + + std::string id; + switch (lastSS) + { + case dssp::structure_type::Helix_3: + id = "HELX_RH_3T_P"; + break; + + case dssp::structure_type::Alphahelix: + id = "HELX_RH_AL_P"; + break; + + case dssp::structure_type::Helix_5: + id = "HELX_RH_PI_P"; + break; + + case dssp::structure_type::Helix_PPII: + id = "HELX_LH_PP_P"; + break; + + case dssp::structure_type::Turn: + id = "TURN_TY1_P"; + break; + + case dssp::structure_type::Bend: + id = "BEND"; + break; + + case dssp::structure_type::Betabridge: + case dssp::structure_type::Strand: + id = "STRN"; + break; + + case dssp::structure_type::Loop: + id = "OTHER"; + break; + } + + if (foundTypes.count(id) == 0) + { + structConfType.emplace({ { "id", id }, + { "criteria", "DSSP" } }); + foundTypes[id] = 1; + } + + structConf.emplace({ // + { "conf_type_id", id }, + { "id", id + std::to_string(foundTypes[id]++) }, + // { "pdbx_PDB_helix_id", vS(12, 14) }, + { "beg_label_comp_id", rb.compound_id() }, + { "beg_label_asym_id", rb.asym_id() }, + { "beg_label_seq_id", rb.seq_id() }, + { "pdbx_beg_PDB_ins_code", rb.pdb_ins_code() }, + { "end_label_comp_id", re.compound_id() }, + { "end_label_asym_id", re.asym_id() }, + { "end_label_seq_id", re.seq_id() }, + { "pdbx_end_PDB_ins_code", re.pdb_ins_code() }, + + { "beg_auth_comp_id", rb.compound_id() }, + { "beg_auth_asym_id", rb.auth_asym_id() }, + { "beg_auth_seq_id", rb.auth_seq_id() }, + { "end_auth_comp_id", re.compound_id() }, + { "end_auth_asym_id", re.auth_asym_id() }, + { "end_auth_seq_id", re.auth_seq_id() } }); + + st = t; + } + + if (stop) + break; + + if (lastSS != t->type()) + { + st = t; + lastSS = t->type(); + } + } + } + + auto &software = db["software"]; + software.emplace({ // + { "pdbx_ordinal", software.get_unique_id("") }, + { "name", "dssp" }, + { "version", klibdsspVersionNumber }, + { "date", klibdsspRevisionDate }, + { "classification", "model annotation" } }); +} --- /dev/null +++ b/src/dssp-io.hpp @@ -0,0 +1,32 @@ +/*- + * 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. + */ + +#pragma once + +#include "dssp.hpp" + +void writeDSSP(const dssp& dssp, std::ostream& os); +void annotateDSSP(cif::datablock &db, const dssp& dssp, bool writeOther, bool writeNewFormat); --- /dev/null +++ b/src/dssp-revision.hpp @@ -0,0 +1,121 @@ +// This file was generated by VersionString.cmake + +#pragma once + +#include + +constexpr const char klibdsspProjectName[] = "libdssp"; +constexpr const char klibdsspVersionNumber[] = "4.4.10"; +constexpr int klibdsspBuildNumber = 306; +constexpr const char klibdsspRevisionGitTag[] = "7b87dac"; +constexpr const char klibdsspRevisionDate[] = "2024-10-28T08:00:57Z"; + +#ifndef VERSION_INFO_DEFINED +#define VERSION_INFO_DEFINED 1 + +namespace version_info_v1_1 +{ + +class version_info_base +{ + public: + static void write_version_string(std::ostream &os, bool verbose) + { + auto s_main = registered_main(); + if (s_main != nullptr) + s_main->write(os, verbose); + + if (verbose) + { + for (auto lib = registered_libraries(); lib != nullptr; lib = lib->m_next) + { + os << "-\n"; + lib->write(os, verbose); + } + } + } + + protected: + version_info_base(const char *name, const char *version, int build_number, const char *git_tag, const char *revision_date, bool is_main) + : m_name(name) + , m_version(version) + , m_build_number(build_number) + , m_git_tag(git_tag) + , m_revision_date(revision_date) + { + if (is_main) + registered_main() = this; + else + { + auto &s_head = registered_libraries(); + m_next = s_head; + s_head = this; + } + } + + void write(std::ostream &os, bool verbose) + { + os << m_name << " version " << m_version << '\n'; + + if (verbose) + { + if (m_build_number != 0) + { + os << "build: " << m_build_number << ' ' << m_revision_date << '\n'; + if (m_git_tag[0] != 0) + os << "git tag: " << m_git_tag << '\n'; + } + } + } + + using version_info_ptr = version_info_base *; + + static version_info_ptr ®istered_main() + { + static version_info_ptr s_main = nullptr; + return s_main; + } + + static version_info_ptr ®istered_libraries() + { + static version_info_ptr s_head = nullptr; + return s_head; + } + + const char *m_name; + const char *m_version; + int m_build_number; + const char *m_git_tag; + const char *m_revision_date; + version_info_base *m_next = nullptr; +}; + +template +class version_info : public version_info_base +{ + public: + using implementation_type = T; + + version_info(const char *name, const char *version, int build_number, const char *git_tag, const char *revision_date, bool is_main) + : version_info_base(name, version, build_number, git_tag, revision_date, is_main) + { + } +}; + +} // namespace version_info_v1_1 + +inline void write_version_string(std::ostream &os, bool verbose) +{ + version_info_v1_1::version_info_base::write_version_string(os, verbose); +} + +#endif + +const class version_info_libdssp_impl : public version_info_v1_1::version_info +{ + public: + version_info_libdssp_impl() + : version_info(klibdsspProjectName, klibdsspVersionNumber, klibdsspBuildNumber, klibdsspRevisionGitTag, klibdsspRevisionDate, false) + { + } +} s_version_info_libdssp_instance; --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -2,13 +2,15 @@ add_executable(tortoize-unit-test + ${PROJECT_SOURCE_DIR}/src/dssp.cpp + ${PROJECT_SOURCE_DIR}/src/dssp-io.cpp ${PROJECT_SOURCE_DIR}/test/tortoize-unit-test.cpp ${PROJECT_SOURCE_DIR}/src/tortoize.cpp) target_compile_definitions(tortoize-unit-test PUBLIC NOMINMAX=1) target_include_directories(tortoize-unit-test PRIVATE ${PROJECT_SOURCE_DIR}/src ${PROJECT_BINARY_DIR}) -target_link_libraries(tortoize-unit-test dssp::dssp cifpp::cifpp Catch2::Catch2) +target_link_libraries(tortoize-unit-test cifpp::cifpp Catch2::Catch2) add_test(NAME tortoize-unit-test COMMAND $ --data-dir ${PROJECT_SOURCE_DIR}/test)