Skip to content

D
dssp
Commit 122d7b02 by Maarten L. Hekkelman
Options

towards v5 of libcifpp

parent fa880e3d
Showing with 151 additions and 161 deletions
CMakeLists.txt
...@@ -120,7 +120,7 @@ find_package(Threads) ...@@ -120,7 +120,7 @@ find_package(Threads)
# Note: use -DBoost_USE_STATIC_LIBS=ON to use boost static libraries # Note: use -DBoost_USE_STATIC_LIBS=ON to use boost static libraries
find_package(cifpp 4.0.1 REQUIRED) find_package(cifpp 5.0.0 REQUIRED)
find_package(Boost COMPONENTS date_time program_options) find_package(Boost COMPONENTS date_time program_options)
add_executable(mkdssp add_executable(mkdssp
......
src/dssp.cpp
/*- /*-
* SPDX-License-Identifier: BSD-2-Clause * SPDX-License-Identifier: BSD-2-Clause
* *
* Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute * Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
* *
* Redistribution and use in source and binary forms, with or without * Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met: * modification, are permitted provided that the following conditions are met:
* *
* 1. Redistributions of source code must retain the above copyright notice, this * 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer * list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice, * 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation * this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution. * and/or other materials provided with the distribution.
* *
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND * 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 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
...@@ -29,41 +29,42 @@ ...@@ -29,41 +29,42 @@
#endif #endif
#include <exception> #include <exception>
#include <iostream>
#include <filesystem> #include <filesystem>
#include <fstream> #include <fstream>
#include <iostream>
#include <boost/format.hpp>
#include <boost/date_time/gregorian/formatters.hpp> #include <boost/date_time/gregorian/formatters.hpp>
#include <boost/format.hpp>
#include <cif++/CifUtils.hpp> #include <cif++/structure/Compound.hpp>
#include <cif++/Cif2PDB.hpp> #include <cif++/utilities.hpp>
#include "dssp.hpp" #include "dssp.hpp"
#include "revision.hpp" #include "revision.hpp"
// -------------------------------------------------------------------- // --------------------------------------------------------------------
std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info) std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo &info)
{ {
/* /*
This is the header line for the residue lines in a DSSP file: 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 # 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
*/ */
boost::format kDSSPResidueLine( boost::format kDSSPResidueLine(
"%5.5d%5.5d%1.1s%1.1s %c %c%c%c%c%c%c%c%c%c%4.4d%4.4d%c%4.4d %11s%11s%11s%11s %6.3f%6.1f%6.1f%6.1f%6.1f %6.1f %6.1f %6.1f"); "%5.5d%5.5d%1.1s%1.1s %c %c%c%c%c%c%c%c%c%c%4.4d%4.4d%c%4.4d %11s%11s%11s%11s %6.3f%6.1f%6.1f%6.1f%6.1f %6.1f %6.1f %6.1f");
auto& residue = info.residue(); // auto& residue = info.residue();
auto &residue = info;
if (residue.asymID().length() > 1) if (residue.asym_id().length() > 1)
throw std::runtime_error("This file contains data that won't fit in the original DSSP format"); throw std::runtime_error("This file contains data that won't fit in the original DSSP format");
char code = 'X'; char code = 'X';
if (mmcif::kAAMap.find(residue.compoundID()) != mmcif::kAAMap.end()) if (mmcif::kAAMap.find(residue.compound_id()) != mmcif::kAAMap.end())
code = mmcif::kAAMap.at(residue.compoundID()); code = mmcif::kAAMap.at(residue.compound_id());
if (code == 'C') // a cysteine if (code == 'C') // a cysteine
{ {
auto ssbridgenr = info.ssBridgeNr(); auto ssbridgenr = info.ssBridgeNr();
if (ssbridgenr) if (ssbridgenr)
...@@ -73,27 +74,27 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info) ...@@ -73,27 +74,27 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info)
char ss; char ss;
switch (info.ss()) switch (info.ss())
{ {
case mmcif::ssAlphahelix: ss = 'H'; break; case mmcif::ssAlphahelix: ss = 'H'; break;
case mmcif::ssBetabridge: ss = 'B'; break; case mmcif::ssBetabridge: ss = 'B'; break;
case mmcif::ssStrand: ss = 'E'; break; case mmcif::ssStrand: ss = 'E'; break;
case mmcif::ssHelix_3: ss = 'G'; break; case mmcif::ssHelix_3: ss = 'G'; break;
case mmcif::ssHelix_5: ss = 'I'; break; case mmcif::ssHelix_5: ss = 'I'; break;
case mmcif::ssHelix_PPII: ss = 'P'; break; case mmcif::ssHelix_PPII: ss = 'P'; break;
case mmcif::ssTurn: ss = 'T'; break; case mmcif::ssTurn: ss = 'T'; break;
case mmcif::ssBend: ss = 'S'; break; case mmcif::ssBend: ss = 'S'; break;
case mmcif::ssLoop: ss = ' '; break; case mmcif::ssLoop: ss = ' '; break;
} }
char helix[4] = { ' ', ' ', ' ', ' ' }; char helix[4] = {' ', ' ', ' ', ' '};
for (mmcif::HelixType helixType: { mmcif::HelixType::rh_3_10, mmcif::HelixType::rh_alpha, mmcif::HelixType::rh_pi, mmcif::HelixType::rh_pp }) for (mmcif::HelixType helixType : {mmcif::HelixType::rh_3_10, mmcif::HelixType::rh_alpha, mmcif::HelixType::rh_pi, mmcif::HelixType::rh_pp})
{ {
switch (info.helix(helixType)) switch (info.helix(helixType))
{ {
case mmcif::Helix::None: helix[static_cast<int>(helixType)] = ' '; break; case mmcif::Helix::None: helix[static_cast<int>(helixType)] = ' '; break;
case mmcif::Helix::Start: helix[static_cast<int>(helixType)] = '>'; break; case mmcif::Helix::Start: helix[static_cast<int>(helixType)] = '>'; break;
case mmcif::Helix::End: helix[static_cast<int>(helixType)] = '<'; break; case mmcif::Helix::End: helix[static_cast<int>(helixType)] = '<'; break;
case mmcif::Helix::StartAndEnd: helix[static_cast<int>(helixType)] = 'X'; break; case mmcif::Helix::StartAndEnd: helix[static_cast<int>(helixType)] = 'X'; break;
case mmcif::Helix::Middle: helix[static_cast<int>(helixType)] = (helixType == mmcif::HelixType::rh_pp ? 'P' : static_cast<char>('3' + static_cast<int>(helixType))); break; case mmcif::Helix::Middle: helix[static_cast<int>(helixType)] = (helixType == mmcif::HelixType::rh_pp ? 'P' : static_cast<char>('3' + static_cast<int>(helixType))); break;
} }
} }
...@@ -103,16 +104,16 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info) ...@@ -103,16 +104,16 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info)
double alpha = residue.alpha(); double alpha = residue.alpha();
char chirality = alpha == 360 ? ' ' : (alpha < 0 ? '-' : '+'); char chirality = alpha == 360 ? ' ' : (alpha < 0 ? '-' : '+');
uint32_t bp[2] = {}; uint32_t bp[2] = {};
char bridgelabel[2] = { ' ', ' ' }; char bridgelabel[2] = {' ', ' '};
for (uint32_t i: { 0, 1 }) for (uint32_t i : {0, 1})
{ {
const auto& [p, ladder, parallel] = info.bridgePartner(i); const auto &[p, ladder, parallel] = info.bridgePartner(i);
if (not p) if (not p)
continue; continue;
bp[i] = p.nr() % 10000; // won't fit otherwise... bp[i] = p.nr() % 10000; // won't fit otherwise...
bridgelabel[i] = (parallel ? 'a' : 'A') + ladder % 26; bridgelabel[i] = (parallel ? 'a' : 'A') + ladder % 26;
} }
...@@ -121,11 +122,11 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info) ...@@ -121,11 +122,11 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info)
sheet = 'A' + (info.sheet() - 1) % 26; sheet = 'A' + (info.sheet() - 1) % 26;
std::string NHO[2], ONH[2]; std::string NHO[2], ONH[2];
for (int i: { 0, 1 }) for (int i : {0, 1})
{ {
const auto& [donor, donorE] = info.donor(i); const auto &[donor, donorE] = info.donor(i);
const auto& [acceptor, acceptorE] = info.acceptor(i); const auto &[acceptor, acceptorE] = info.acceptor(i);
NHO[i] = ONH[i] = "0, 0.0"; NHO[i] = ONH[i] = "0, 0.0";
if (acceptor) if (acceptor)
...@@ -141,21 +142,20 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info) ...@@ -141,21 +142,20 @@ std::string ResidueToDSSPLine(const mmcif::DSSP::ResidueInfo& info)
} }
} }
auto ca = residue.atomByID("CA"); // auto ca = residue.atomByID("CA");
auto const& [cax, cay, caz] = ca.location(); auto const &[cax, cay, caz] = residue.ca_location();
return (kDSSPResidueLine % info.nr() % ca.authSeqID() % ca.pdbxAuthInsCode() % ca.authAsymID() % code % return (kDSSPResidueLine % 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] % 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) % bp[0] % bp[1] % sheet % floor(info.accessibility() + 0.5) %
NHO[0] % ONH[0] % NHO[1] % ONH[1] % NHO[0] % ONH[0] % NHO[1] % ONH[1] %
residue.tco() % residue.kappa() % alpha % residue.phi() % residue.psi() % residue.tco() % residue.kappa() % alpha % residue.phi() % residue.psi() %
cax % cay % caz).str(); cax % cay % caz)
.str();
} }
void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std::ostream& os) void writeDSSP(const cif::datablock &db, const mmcif::DSSP &dssp, std::ostream &os)
{ {
auto &db = structure.datablock();
const std::string kFirstLine("==== Secondary Structure Definition by the program DSSP, NKI version 4.0 ==== "); const std::string kFirstLine("==== Secondary Structure Definition by the program DSSP, NKI version 4.0 ==== ");
boost::format kHeaderLine("%1% %|127t|%2%"); boost::format kHeaderLine("%1% %|127t|%2%");
...@@ -169,28 +169,23 @@ void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std:: ...@@ -169,28 +169,23 @@ void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std::
os << kHeaderLine % (kFirstLine + "DATE=" + to_iso_extended_string(today)) % '.' << std::endl os << kHeaderLine % (kFirstLine + "DATE=" + to_iso_extended_string(today)) % '.' << std::endl
<< kHeaderLine % "REFERENCE W. KABSCH AND C.SANDER, BIOPOLYMERS 22 (1983) 2577-2637" % '.' << std::endl << kHeaderLine % "REFERENCE W. KABSCH AND C.SANDER, BIOPOLYMERS 22 (1983) 2577-2637" % '.' << std::endl
<< GetPDBHEADERLine(db, 127) << '.' << std::endl << /*GetPDBHEADERLine(db, 127) << */ '.' << std::endl
<< GetPDBCOMPNDLine(db, 127) << '.' << std::endl << /*GetPDBCOMPNDLine(db, 127) << */ '.' << std::endl
<< GetPDBSOURCELine(db, 127) << '.' << std::endl << /*GetPDBSOURCELine(db, 127) << */ '.' << std::endl
<< GetPDBAUTHORLine(db, 127) << '.' << std::endl; << /*GetPDBAUTHORLine(db, 127) << */ '.' << std::endl;
os << boost::format("%5.5d%3.3d%3.3d%3.3d%3.3d TOTAL NUMBER OF RESIDUES, NUMBER OF CHAINS, NUMBER OF SS-BRIDGES(TOTAL,INTRACHAIN,INTERCHAIN) %|127t|%c") % os << boost::format("%5.5d%3.3d%3.3d%3.3d%3.3d TOTAL NUMBER OF RESIDUES, NUMBER OF CHAINS, NUMBER OF SS-BRIDGES(TOTAL,INTRACHAIN,INTERCHAIN) %|127t|%c") %
stats.nrOfResidues % stats.nrOfChains % stats.nrOfSSBridges % stats.nrOfIntraChainSSBridges % (stats.nrOfSSBridges - stats.nrOfIntraChainSSBridges) % '.' << std::endl; stats.nrOfResidues % stats.nrOfChains % stats.nrOfSSBridges % stats.nrOfIntraChainSSBridges % (stats.nrOfSSBridges - stats.nrOfIntraChainSSBridges) % '.'
os << kHeaderLine % (boost::format("%8.1f ACCESSIBLE SURFACE OF PROTEIN (ANGSTROM**2)") % stats.accessibleSurface) % '.' << std::endl; << std::endl;
os << kHeaderLine % (boost::format("%8.1f ACCESSIBLE SURFACE OF PROTEIN (ANGSTROM**2)") % stats.accessibleSurface) % '.' << std::endl;
// hydrogenbond summary // hydrogenbond summary
os << kHeaderLine % ( os << kHeaderLine % (boost::format("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(J) , SAME NUMBER PER 100 RESIDUES") % stats.nrOfHBonds % (stats.nrOfHBonds * 100.0 / stats.nrOfResidues)) % '.' << std::endl;
boost::format("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(J) , SAME NUMBER PER 100 RESIDUES")
% stats.nrOfHBonds % (stats.nrOfHBonds * 100.0 / stats.nrOfResidues)) % '.' << std::endl;
os << kHeaderLine % ( os << kHeaderLine % (boost::format("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS IN PARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES") % stats.nrOfHBondsInParallelBridges % (stats.nrOfHBondsInParallelBridges * 100.0 / stats.nrOfResidues)) % '.' << std::endl;
boost::format("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS IN PARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES")
% stats.nrOfHBondsInParallelBridges % (stats.nrOfHBondsInParallelBridges * 100.0 / stats.nrOfResidues)) % '.' << std::endl;
os << kHeaderLine % ( os << kHeaderLine % (boost::format("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS IN ANTIPARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES") % stats.nrOfHBondsInAntiparallelBridges % (stats.nrOfHBondsInAntiparallelBridges * 100.0 / stats.nrOfResidues)) % '.' << std::endl;
boost::format("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS IN ANTIPARALLEL BRIDGES, SAME NUMBER PER 100 RESIDUES")
% stats.nrOfHBondsInAntiparallelBridges % (stats.nrOfHBondsInAntiparallelBridges * 100.0 / stats.nrOfResidues)) % '.' << std::endl;
boost::format kHBondsLine("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(I%c%1.1d), SAME NUMBER PER 100 RESIDUES"); boost::format kHBondsLine("%5.5d%5.1f TOTAL NUMBER OF HYDROGEN BONDS OF TYPE O(I)-->H-N(I%c%1.1d), SAME NUMBER PER 100 RESIDUES");
for (int k = 0; k < 11; ++k) for (int k = 0; k < 11; ++k)
...@@ -199,19 +194,19 @@ void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std:: ...@@ -199,19 +194,19 @@ void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std::
// histograms... // 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; 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.residuesPerAlphaHelixHistogram) for (auto hi : stats.residuesPerAlphaHelixHistogram)
os << boost::format("%3.3d") % hi; os << boost::format("%3.3d") % hi;
os << " RESIDUES PER ALPHA HELIX ." << std::endl; os << " RESIDUES PER ALPHA HELIX ." << std::endl;
for (auto hi: stats.parallelBridgesPerLadderHistogram) for (auto hi : stats.parallelBridgesPerLadderHistogram)
os << boost::format("%3.3d") % hi; os << boost::format("%3.3d") % hi;
os << " PARALLEL BRIDGES PER LADDER ." << std::endl; os << " PARALLEL BRIDGES PER LADDER ." << std::endl;
for (auto hi: stats.antiparallelBridgesPerLadderHistogram) for (auto hi : stats.antiparallelBridgesPerLadderHistogram)
os << boost::format("%3.3d") % hi; os << boost::format("%3.3d") % hi;
os << " ANTIPARALLEL BRIDGES PER LADDER ." << std::endl; os << " ANTIPARALLEL BRIDGES PER LADDER ." << std::endl;
for (auto hi: stats.laddersPerSheetHistogram) for (auto hi : stats.laddersPerSheetHistogram)
os << boost::format("%3.3d") % hi; os << boost::format("%3.3d") % hi;
os << " LADDERS PER SHEET ." << std::endl; os << " LADDERS PER SHEET ." << std::endl;
...@@ -222,23 +217,21 @@ void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std:: ...@@ -222,23 +217,21 @@ void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std::
"%5.5d !%c 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"); "%5.5d !%c 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");
int last = 0; int last = 0;
for (auto ri: dssp) for (auto ri : dssp)
{ {
// insert a break line whenever we detect missing residues // insert a break line whenever we detect missing residues
// can be the transition to a different chain, or missing residues in the current chain // can be the transition to a different chain, or missing residues in the current chain
if (ri.nr() != last + 1) if (ri.nr() != last + 1)
os << (kDSSPResidueLine % (last + 1) % (ri.chainBreak() == mmcif::ChainBreak::NewChain ? '*' : ' ')) << std::endl; os << (kDSSPResidueLine % (last + 1) % (ri.chainBreak() == mmcif::ChainBreak::NewChain ? '*' : ' ')) << std::endl;
os << ResidueToDSSPLine(ri) << std::endl; os << ResidueToDSSPLine(ri) << std::endl;
last = ri.nr(); last = ri.nr();
} }
} }
void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool writeOther, std::ostream& os) void annotateDSSP(cif::datablock &db, const mmcif::DSSP &dssp, bool writeOther, std::ostream &os)
{ {
auto& db = structure.datablock();
if (dssp.empty()) if (dssp.empty())
{ {
if (cif::VERBOSE) if (cif::VERBOSE)
...@@ -247,18 +240,18 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri ...@@ -247,18 +240,18 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri
else else
{ {
// replace all struct_conf and struct_conf_type records // replace all struct_conf and struct_conf_type records
auto& structConfType = db["struct_conf_type"]; auto &structConfType = db["struct_conf_type"];
structConfType.clear(); structConfType.clear();
auto& structConf = db["struct_conf"]; auto &structConf = db["struct_conf"];
structConf.clear(); structConf.clear();
std::map<std::string,int> foundTypes; std::map<std::string, int> foundTypes;
auto st = dssp.begin(), lt = st; auto st = dssp.begin(), lt = st;
auto lastSS = st->ss(); auto lastSS = st->ss();
for (auto t = dssp.begin(); ; lt = t, ++t) for (auto t = dssp.begin();; lt = t, ++t)
{ {
bool stop = t == dssp.end(); bool stop = t == dssp.end();
...@@ -266,8 +259,8 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri ...@@ -266,8 +259,8 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri
if (flush and (writeOther or lastSS != mmcif::SecondaryStructureType::ssLoop)) if (flush and (writeOther or lastSS != mmcif::SecondaryStructureType::ssLoop))
{ {
auto& rb = st->residue(); auto &rb = *st;
auto& re = lt->residue(); auto &re = *lt;
std::string id; std::string id;
switch (lastSS) switch (lastSS)
...@@ -308,36 +301,34 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri ...@@ -308,36 +301,34 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri
if (foundTypes.count(id) == 0) if (foundTypes.count(id) == 0)
{ {
structConfType.emplace({ structConfType.emplace({{"id", id},
{ "id", id }, {"criteria", "DSSP"}});
{ "criteria", "DSSP" }
});
foundTypes[id] = 1; foundTypes[id] = 1;
} }
structConf.emplace({ structConf.emplace({
{ "conf_type_id", id }, {"conf_type_id", id},
{ "id", id + std::to_string(foundTypes[id]++) }, {"id", id + std::to_string(foundTypes[id]++)},
// { "pdbx_PDB_helix_id", vS(12, 14) }, // { "pdbx_PDB_helix_id", vS(12, 14) },
{ "beg_label_comp_id", rb.compoundID() }, {"beg_label_comp_id", rb.compound_id()},
{ "beg_label_asym_id", rb.asymID() }, {"beg_label_asym_id", rb.asym_id()},
{ "beg_label_seq_id", rb.seqID() }, {"beg_label_seq_id", rb.seq_id()},
{ "pdbx_beg_PDB_ins_code", rb.authInsCode() }, {"pdbx_beg_PDB_ins_code", rb.pdb_ins_code()},
{ "end_label_comp_id", re.compoundID() }, {"end_label_comp_id", re.compound_id()},
{ "end_label_asym_id", re.asymID() }, {"end_label_asym_id", re.asym_id()},
{ "end_label_seq_id", re.seqID() }, {"end_label_seq_id", re.seq_id()},
{ "pdbx_end_PDB_ins_code", re.authInsCode() }, {"pdbx_end_PDB_ins_code", re.pdb_ins_code()},
{ "beg_auth_comp_id", rb.compoundID() }, {"beg_auth_comp_id", rb.compound_id()},
{ "beg_auth_asym_id", rb.authAsymID() }, {"beg_auth_asym_id", rb.auth_asym_id()},
{ "beg_auth_seq_id", rb.authSeqID() }, {"beg_auth_seq_id", rb.auth_seq_id()},
{ "end_auth_comp_id", re.compoundID() }, {"end_auth_comp_id", re.compound_id()},
{ "end_auth_asym_id", re.authAsymID() }, {"end_auth_asym_id", re.auth_asym_id()},
{ "end_auth_seq_id", re.authSeqID() } {"end_auth_seq_id", re.auth_seq_id()}
// { "pdbx_PDB_helix_class", vS(39, 40) }, // { "pdbx_PDB_helix_class", vS(39, 40) },
// { "details", vS(41, 70) }, // { "details", vS(41, 70) },
// { "pdbx_PDB_helix_length", vI(72, 76) } // { "pdbx_PDB_helix_length", vI(72, 76) }
}); });
st = t; st = t;
...@@ -354,7 +345,10 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri ...@@ -354,7 +345,10 @@ void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool wri
} }
} }
db.add_software("dssp", "other", kVersionNumber, kBuildDate); // db.add_software("dssp", "other", kVersionNumber, kBuildDate);
// db["software"].emplace({
// });
db.write(os); db.write(os);
} }
src/dssp.hpp
...@@ -26,9 +26,8 @@ ...@@ -26,9 +26,8 @@
#pragma once #pragma once
#include <cif++/Structure.hpp> #include <cif++/structure/DSSP.hpp>
#include <cif++/Secondary.hpp>
void writeDSSP(const mmcif::Structure& structure, const mmcif::DSSP& dssp, std::ostream& os); void writeDSSP(const cif::datablock &db, const mmcif::DSSP& dssp, std::ostream& os);
void annotateDSSP(mmcif::Structure& structure, const mmcif::DSSP& dssp, bool writeOther, std::ostream& os); void annotateDSSP(cif::datablock &db, const mmcif::DSSP& dssp, bool writeOther, std::ostream& os);
src/mkdssp.cpp
...@@ -33,23 +33,20 @@ ...@@ -33,23 +33,20 @@
#include <filesystem> #include <filesystem>
#include <fstream> #include <fstream>
#include <gxrio.hpp>
#include <boost/format.hpp> #include <boost/format.hpp>
#include <boost/date_time/gregorian/formatters.hpp> #include <boost/date_time/gregorian/formatters.hpp>
#include <cif++/Structure.hpp> #include <cif++/structure/DSSP.hpp>
#include <cif++/Secondary.hpp> #include <cif++/structure/Compound.hpp>
#include <cif++/CifUtils.hpp>
#include <cif++/Cif2PDB.hpp>
#include <boost/program_options.hpp> #include <boost/program_options.hpp>
#include <boost/iostreams/filter/gzip.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include "dssp.hpp" #include "dssp.hpp"
#include "revision.hpp" #include "revision.hpp"
namespace fs = std::filesystem; namespace fs = std::filesystem;
namespace io = boost::iostreams;
namespace po = boost::program_options; namespace po = boost::program_options;
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -148,9 +145,9 @@ int d_main(int argc, const char* argv[]) ...@@ -148,9 +145,9 @@ int d_main(int argc, const char* argv[])
// Load extra CCD definitions, if any // Load extra CCD definitions, if any
if (vm.count("compounds")) if (vm.count("compounds"))
cif::addFileResource("components.cif", vm["compounds"].as<std::string>()); cif::add_file_resource("components.cif", vm["compounds"].as<std::string>());
else if (vm.count("components")) else if (vm.count("components"))
cif::addFileResource("components.cif", vm["components"].as<std::string>()); cif::add_file_resource("components.cif", vm["components"].as<std::string>());
if (vm.count("extra-compounds")) if (vm.count("extra-compounds"))
mmcif::CompoundFactory::instance().pushDictionary(vm["extra-compounds"].as<std::string>()); mmcif::CompoundFactory::instance().pushDictionary(vm["extra-compounds"].as<std::string>());
...@@ -158,7 +155,7 @@ int d_main(int argc, const char* argv[]) ...@@ -158,7 +155,7 @@ int d_main(int argc, const char* argv[])
// And perhaps a private mmcif_pdbx dictionary // And perhaps a private mmcif_pdbx dictionary
if (vm.count("mmcif-dictionary")) if (vm.count("mmcif-dictionary"))
cif::addFileResource("mmcif_pdbx_v50.dic", vm["mmcif-dictionary"].as<std::string>()); cif::add_file_resource("mmcif_pdbx_v50.dic", vm["mmcif-dictionary"].as<std::string>());
if (vm.count("dict")) if (vm.count("dict"))
{ {
...@@ -166,8 +163,17 @@ int d_main(int argc, const char* argv[]) ...@@ -166,8 +163,17 @@ int d_main(int argc, const char* argv[])
mmcif::CompoundFactory::instance().pushDictionary(dict); mmcif::CompoundFactory::instance().pushDictionary(dict);
} }
mmcif::File f(vm["xyzin"].as<std::string>()); cif::file f(vm["xyzin"].as<std::string>());
mmcif::Structure structure(f, 1, mmcif::StructureOpenOptions::SkipHydrogen);
f.load_dictionary("mmcif_pdbx_v50");
if (not f.is_valid())
{
std::cerr << "File not valid" << std::endl;
exit(1);
}
// mmcif::Structure structure(f, 1, mmcif::StructureOpenOptions::SkipHydrogen);
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -185,9 +191,7 @@ int d_main(int argc, const char* argv[]) ...@@ -185,9 +191,7 @@ int d_main(int argc, const char* argv[])
{ {
fs::path output = vm["output"].as<std::string>(); fs::path output = vm["output"].as<std::string>();
if (output.extension() == ".gz") if (output.extension() == ".gz" or output.extension() == ".xz")
output = output.stem();
else if (output.extension() == ".bz2")
output = output.stem(); output = output.stem();
if (output.extension() == ".dssp") if (output.extension() == ".dssp")
...@@ -197,37 +201,30 @@ int d_main(int argc, const char* argv[]) ...@@ -197,37 +201,30 @@ int d_main(int argc, const char* argv[])
} }
mmcif::DSSP dssp(structure, pp_stretch, fmt == "dssp"); mmcif::DSSP dssp(f.front(), 1, pp_stretch, fmt == "dssp");
if (vm.count("output")) if (vm.count("output"))
{ {
fs::path output = vm["output"].as<std::string>(); fs::path output = vm["output"].as<std::string>();
std::ofstream of(output, std::ios_base::out | std::ios_base::binary); gxrio::ofstream out(output);
if (not of.is_open()) if (not out.is_open())
{ {
std::cerr << "Could not open output file" << std::endl; std::cerr << "Could not open output file" << std::endl;
exit(1); exit(1);
} }
io::filtering_stream<io::output> out;
if (output.extension() == ".gz")
out.push(io::gzip_compressor());
out.push(of);
if (fmt == "dssp") if (fmt == "dssp")
writeDSSP(structure, dssp, out); writeDSSP(f.front(), dssp, out);
else else
annotateDSSP(structure, dssp, writeOther, out); annotateDSSP(f.front(), dssp, writeOther, out);
} }
else else
{ {
if (fmt == "dssp") if (fmt == "dssp")
writeDSSP(structure, dssp, std::cout); writeDSSP(f.front(), dssp, std::cout);
else else
annotateDSSP(structure, dssp, writeOther, std::cout); annotateDSSP(f.front(), dssp, writeOther, std::cout);
} }
return 0; return 0;
...@@ -242,7 +239,7 @@ int main(int argc, const char* argv[]) ...@@ -242,7 +239,7 @@ int main(int argc, const char* argv[])
try try
{ {
#if defined(DATA_DIR) #if defined(DATA_DIR)
cif::addDataDirectory(DATA_DIR); cif::add_data_directory(DATA_DIR);
#endif #endif
result = d_main(argc, argv); result = d_main(argc, argv);
} }
......
test/unit-test.cpp
...@@ -30,10 +30,7 @@ ...@@ -30,10 +30,7 @@
#include <stdexcept> #include <stdexcept>
#include <cif++/Structure.hpp> #include <cif++/structure/DSSP.hpp>
#include <cif++/Secondary.hpp>
#include <cif++/CifUtils.hpp>
#include <cif++/Cif2PDB.hpp>
#include "dssp.hpp" #include "dssp.hpp"
...@@ -45,14 +42,16 @@ BOOST_AUTO_TEST_CASE(ut_dssp) ...@@ -45,14 +42,16 @@ BOOST_AUTO_TEST_CASE(ut_dssp)
{ {
using namespace std::literals; using namespace std::literals;
mmcif::File f("1cbs.cif.gz"); cif::file f("1cbs.cif.gz");
mmcif::Structure structure(f, 1, mmcif::StructureOpenOptions::SkipHydrogen); f.load_dictionary("mmcif_pdbx_v50");
BOOST_CHECK(f.is_valid());
// mmcif::Structure structure(f, 1, mmcif::StructureOpenOptions::SkipHydrogen);
mmcif::DSSP dssp(structure, 3, true); mmcif::DSSP dssp(f.front(), 1, 3, true);
std::stringstream test; std::stringstream test;
writeDSSP(structure, dssp, test); writeDSSP(f.front(), dssp, test);
std::ifstream reference("1cbs.dssp"); std::ifstream reference("1cbs.dssp");
...@@ -91,25 +90,26 @@ BOOST_AUTO_TEST_CASE(ut_mmcif_2) ...@@ -91,25 +90,26 @@ BOOST_AUTO_TEST_CASE(ut_mmcif_2)
using namespace std::literals; using namespace std::literals;
using namespace cif::literals; using namespace cif::literals;
mmcif::File f("1cbs.cif.gz"); cif::file f("1cbs.cif.gz");
f.loadDictionary("mmcif_pdbx_v50"); f.load_dictionary("mmcif_pdbx_v50");
BOOST_CHECK(f.is_valid());
mmcif::Structure structure(f, 1, mmcif::StructureOpenOptions::SkipHydrogen); // mmcif::Structure structure(f, 1, mmcif::StructureOpenOptions::SkipHydrogen);
mmcif::DSSP dssp(structure, 3, true); mmcif::DSSP dssp(f.front(), 1, 3, true);
std::stringstream test; std::stringstream test;
annotateDSSP(structure, dssp, true, test); annotateDSSP(f.front(), dssp, true, test);
mmcif::File rf("1cbs-dssp.cif"); cif::file rf("1cbs-dssp.cif");
mmcif::Structure rs(rf, 1, mmcif::StructureOpenOptions::SkipHydrogen); // mmcif::Structure rs(rf, 1, mmcif::StructureOpenOptions::SkipHydrogen);
structure.datablock()["software"].erase("name"_key == "dssp"); // structure.datablock()["software"].erase("name"_key == "dssp");
rs.datablock()["software"].erase("name"_key == "dssp"); // rs.datablock()["software"].erase("name"_key == "dssp");
// generate some output on different files: // generate some output on different files:
cif::VERBOSE = 2; cif::VERBOSE = 2;
BOOST_CHECK(structure.datablock() == rs.datablock()); // BOOST_CHECK(f.front() == rf.front());
} }
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment