Skip to content

D
dssp
Commit 1bd675e1 by Maarten L. Hekkelman
Options

Merge branch 'develop' into trunk

parents 49306a57 722c616a
Showing with 176 additions and 568 deletions
CMakeLists.txt
......@@ -25,7 +25,7 @@
cmake_minimum_required(VERSION 3.15)
# set the project name
project(mkdssp VERSION 4.3.2 LANGUAGES CXX)
project(mkdssp VERSION 4.4.0 LANGUAGES CXX)
list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
......@@ -211,6 +211,10 @@ install(TARGETS ${PROJECT_NAME}
RUNTIME DESTINATION ${BIN_INSTALL_DIR}
)
set(CIFPP_DATA_DIR "${CMAKE_INSTALL_FULL_DATADIR}/libcifpp")
install(FILES "${CMAKE_CURRENT_SOURCE_DIR}/mmcif_pdbx/dssp-extension.dic"
DESTINATION ${CIFPP_DATA_DIR})
# manual
find_program(PANDOC pandoc)
......@@ -265,7 +269,9 @@ if(ENABLE_TESTING)
add_executable(unit-test-dssp ${PROJECT_SOURCE_DIR}/test/unit-test-dssp.cpp ${PROJECT_SOURCE_DIR}/src/dssp-io.cpp)
if(USE_RSRC)
mrc_target_resources(unit-test-dssp ${CIFPP_SHARE_DIR}/mmcif_pdbx.dic ${CIFPP_SHARE_DIR}/mmcif_ddl.dic)
mrc_target_resources(unit-test-dssp
${CIFPP_SHARE_DIR}/mmcif_pdbx.dic ${CIFPP_SHARE_DIR}/mmcif_ddl.dic
${CMAKE_CURRENT_SOURCE_DIR}/mmcif_pdbx/dssp-extension.dic)
endif()
target_include_directories(unit-test-dssp PRIVATE
......
changelog
Version 4.3.2
Version 4.4.0
- Packaging (split out the web server again)
- Added markdown doc file, all doc files are now created from mkdssp.1
using pandoc, if that is installed.
......
include/dssp.hpp
......@@ -172,6 +172,7 @@ class dssp
std::tuple<residue_info, int, bool> bridge_partner(int i) const;
int sheet() const;
int strand() const;
/// \brief return resinfo and the energy of the bond
std::tuple<residue_info, double> acceptor(int i) const;
......
mmcif_pdbx/dssp-extension.dic
......@@ -477,6 +477,9 @@ save_dssp_struct_ladder
;
;
_dssp_struct_ladder.id A
_dssp_struct_ladder.sheet_id A
_dssp_struct_ladder.range_id_1 A
_dssp_struct_ladder.range_id_2 B
_dssp_struct_ladder.type anti-parallel
_dssp_struct_ladder.beg_1_label_comp_id GLY
_dssp_struct_ladder.beg_1_label_asym_id A
......@@ -520,6 +523,39 @@ save__dssp_struct_ladder.id
_item_type.code code
save_
save__dssp_struct_ladder.sheet_id
_item.description
; This data item is a pointer to _struct_sheet.id in the
STRUCT_SHEET category.
;
_item.name '_dssp_struct_ladder.sheet_id'
_item.category_id dssp_struct_ladder
_item.mandatory_code yes
_item_type.code code
save_
save__dssp_struct_ladder.range_id_1
_item.description
; This data item is a pointer to _struct_sheet_range.id in
the STRUCT_SHEET_RANGE category.
;
_item.name '_dssp_struct_ladder.range_id_1'
_item.category_id dssp_struct_ladder
_item.mandatory_code yes
_item_type.code code
save_
save__dssp_struct_ladder.range_id_2
_item.description
; This data item is a pointer to _struct_sheet_range.id in
the STRUCT_SHEET_RANGE category.
;
_item.name '_dssp_struct_ladder.range_id_2'
_item.category_id dssp_struct_ladder
_item.mandatory_code yes
_item_type.code code
save_
save__dssp_struct_ladder.type
_item.description
; The type of the ladder, be it parallel or anti-parallel
......@@ -1466,9 +1502,10 @@ save_dssp_struct_summary
_dssp_struct_summary.helix_pp .
_dssp_struct_summary.bend S
_dssp_struct_summary.chirality +
_dssp_struct_summary.sheet .
_dssp_struct_summary.strand .
_dssp_struct_summary.ladder_1 .
_dssp_struct_summary.ladder_2 .
_dssp_struct_summary.sheet .
_dssp_struct_summary.accessibility 23.1
_dssp_struct_summary.TCO 0.940
_dssp_struct_summary.kappa 73.9
......@@ -1655,31 +1692,43 @@ save__dssp_struct_summary.chirality
'-'
save_
save__dssp_struct_summary.ladder_1
save__dssp_struct_summary.sheet
_item.description
; Label for the first beta bridge of which this residue is part.
; This data item is a pointer to _struct_sheet.id in
the STRUCT_SHEET category.
;
_item.name '_dssp_struct_summary.ladder_1'
_item.name '_dssp_struct_summary.sheet'
_item.category dssp_struct_summary
_item.mandatory_code no
_item_type.code code
save_
save__dssp_struct_summary.ladder_2
save__dssp_struct_summary.strand
_item.description
; Label for the second beta bridge of which this residue is part.
; This data item is a pointer to _struct_sheet_range.id in
the STRUCT_SHEET_RANGE category.
;
_item.name '_dssp_struct_summary.ladder_2'
_item.name '_dssp_struct_summary.strand'
_item.category dssp_struct_summary
_item.mandatory_code no
_item_type.code code
save_
save__dssp_struct_summary.sheet
save__dssp_struct_summary.ladder_1
_item.description
; Label for the sheet of which this residue is part.
; Label for the first beta bridge of which this residue is part.
;
_item.name '_dssp_struct_summary.sheet'
_item.name '_dssp_struct_summary.ladder_1'
_item.category dssp_struct_summary
_item.mandatory_code no
_item_type.code code
save_
save__dssp_struct_summary.ladder_2
_item.description
; Label for the second beta bridge of which this residue is part.
;
_item.name '_dssp_struct_summary.ladder_2'
_item.category dssp_struct_summary
_item.mandatory_code no
_item_type.code code
......
src/dssp-io.cpp
......@@ -158,7 +158,11 @@ void writeDSSP(const dssp &dssp, std::ostream &os)
std::time_t today = system_clock::to_time_t(system_clock::now());
std::tm *tm = std::gmtime(&today);
os << "==== Secondary Structure Definition by the program DSSP, NKI version 4.3 ==== DATE=" << std::put_time(tm, "%F") << " ." << std::endl
std::string version = kVersionNumber;
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 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
......@@ -322,8 +326,9 @@ void writeBridgePairs(cif::datablock &db, const dssp &dssp)
void writeSheets(cif::datablock &db, const dssp &dssp)
{
using namespace cif::literals;
using res_list = std::vector<dssp::residue_info>;
using ss_type = dssp::structure_type;
// clean up old info first
......@@ -334,423 +339,44 @@ void writeSheets(cif::datablock &db, const dssp &dssp)
// create a list of strands, based on the SS info in DSSP. Store sheet number along with the strand.
std::vector<std::tuple<int, res_list>> strands;
std::map<std::tuple<std::string, int>, int> sheetMap; // mapping from bridge number (=info.sheet()) in DSSP to sheet ID
ss_type ss = ss_type::Loop;
std::map<std::tuple<int,int>, res_list> strands;
std::set<int> sheetNrs;
for (auto &res : dssp)
{
std::tuple<std::string, int> sheetID{ res.asym_id(), res.sheet() };
ss_type iss = res.type();
if (iss == ss and ss == ss_type::Strand and
sheetMap.count(sheetID) > 0 and sheetMap[sheetID] == std::get<0>(strands.back()))
{
std::get<1>(strands.back()).emplace_back(res);
if (res.type() != dssp::structure_type::Strand and res.type() != dssp::structure_type::Betabridge)
continue;
}
ss = iss;
if (ss != ss_type::Strand)
continue;
if (not sheetMap.count(sheetID))
sheetMap[sheetID] = static_cast<int>(sheetMap.size());
strands.emplace_back(sheetMap[sheetID], res_list{ res });
strands[{res.sheet(), res.strand()}].emplace_back(res);
sheetNrs.insert(res.sheet());
}
// sort the strands vector
std::sort(strands.begin(), strands.end(), [](auto &a, auto &b)
{
const auto &[sheetA, strandsA] = a;
const auto &[sheetB, strandsB] = b;
int d = sheetA - sheetB;
if (d == 0)
d = strandsA.front().nr() - strandsB.front().nr();
return d < 0; });
// --------------------------------------------------------------------
// write out the struct_sheet, since all info is available now
auto &struct_sheet = db["struct_sheet"];
auto &struct_sheet_range = db["struct_sheet_range"];
int lastSheet = -1;
for (const auto &[sheetNr, strand] : strands)
{
if (sheetNr != lastSheet)
{
struct_sheet.emplace({
{ "id", cif::cif_id_for_number(sheetNr) },
{ "number_strands",
std::count_if(strands.begin(), strands.end(), [nr = sheetNr](std::tuple<int, res_list> const &s)
{ return std::get<0>(s) == nr; })
}
});
lastSheet = sheetNr;
}
}
// Each residue resides in a single strand which is part of a single sheet
// this function returns the sequence number inside the sheet for the strand
// containing res
std::map<int,int> strandMap;
int strandNr = 0;
for (auto &&[sheet, strand] : strands)
{
for (auto &r : strand)
{
if (r.type() != ss_type::Strand)
continue;
strandMap[r.nr()] = strandNr;
}
++strandNr;
}
// This map is used to record the sense of a ladder, and can be used
// to detect ladders already seen.
std::map<std::tuple<int, int, int>, std::tuple<int, bool>> ladderSense;
for (auto &res : dssp)
{
if (res.type() != ss_type::Strand)
continue;
int s1 = strandMap[res.nr()];
for (int i : { 0, 1 })
{
const auto &[p, ladder, parallel] = res.bridge_partner(i);
if (not p or p.asym_id() != res.asym_id() or p.sheet() != res.sheet() or p.type() != ss_type::Strand)
continue;
int s2 = strandMap[p.nr()];
// assert(s1 != s2);
if (s2 == s1)
continue;
std::tuple<std::string, int> sheetID{ res.asym_id(), res.sheet() };
int sheet = sheetMap[sheetID];
auto k = s1 > s2 ? std::make_tuple(sheet, s2, s1) : std::make_tuple(sheet, s1, s2);
if (ladderSense.count(k))
{
// assert(ladderSense[k] == std::make_tuple(ladder, parallel));
assert(std::get<1>(ladderSense[k]) == parallel);
continue;
}
ladderSense.emplace(k, std::make_tuple(ladder, parallel));
}
}
auto &struct_sheet_order = db["struct_sheet_order"];
auto &struct_sheet_hbond = db["struct_sheet_hbond"];
auto &pdbx_struct_sheet_hbond = db["pdbx_struct_sheet_hbond"];
auto &pdbx_poly_seq_scheme = db["pdbx_poly_seq_scheme"];
for (const auto &[key, value] : ladderSense)
for (int sheetNr : sheetNrs)
{
const auto &[sheet, s1, s2] = key;
const auto &[ladder, parallel] = value;
struct_sheet_order.emplace({ { "sheet_id", cif::cif_id_for_number(sheet) },
// { "dssp_struct_ladder_id", cif::cifIdForNumber(ladder) },
{ "range_id_1", s1 + 1 },
{ "range_id_2", s2 + 1 },
{ "sense", parallel ? "parallel" : "anti-parallel" } });
res_list strand1, strand2;
for (int strandIx = 0; static_cast<size_t>(strandIx) < strands.size(); ++strandIx)
{
const auto &[sSheet, strand] = strands[strandIx];
if (sSheet != sheet)
continue;
if (strandIx == s1)
strand1 = strand;
else if (strandIx == s2)
{
strand2 = strand;
break;
}
}
assert(not(strand1.empty() or strand2.empty()));
int beg1SeqID = 0, beg2SeqID = 0, end1SeqID = 0, end2SeqID = 0;
std::string beg1AtomID, beg2AtomID, end1AtomID, end2AtomID;
if (parallel)
{
// I. a d II. a d parallel
// \ /
// b e b e <= the residues forming the bridge
// / \ ..
// c f c f
for (auto b : strand1)
{
for (int i : { 0, 1 })
{
const auto &[e, ladder1, parallel1] = b.bridge_partner(i);
auto esi = std::find(strand2.begin(), strand2.end(), e);
if (esi == strand2.end())
continue;
auto bi = std::find(dssp.begin(), dssp.end(), b);
assert(bi != dssp.end() and bi != dssp.begin());
auto a = *std::prev(bi);
auto c = *std::next(bi);
auto ei = std::find(dssp.begin(), dssp.end(), e);
assert(ei != dssp.end() and ei != dssp.begin());
auto d = *std::prev(ei);
auto f = *std::next(ei);
if (test_bond(e, a) and test_bond(c, e)) // case I.
{
beg1SeqID = a.seq_id();
beg2SeqID = e.seq_id();
beg1AtomID = "O";
beg2AtomID = "N";
}
else if (test_bond(b, d) and test_bond(f, b)) // case II.
{
beg1SeqID = b.seq_id();
beg2SeqID = d.seq_id();
beg1AtomID = "N";
beg2AtomID = "O";
}
break;
}
if (beg1SeqID)
break;
}
auto sheetID = cif::cif_id_for_number(sheetNr - 1);
std::reverse(strand1.begin(), strand1.end());
std::reverse(strand2.begin(), strand2.end());
for (auto b : strand1)
{
for (int i : { 0, 1 })
struct_sheet.emplace({
{ "id", sheetID },
{ "number_strands",
std::count_if(strands.begin(), strands.end(), [nr = sheetNr](std::tuple<std::tuple<int,int>, res_list> const &s)
{
const auto &[e, ladder1, parallel1] = b.bridge_partner(i);
auto esi = std::find(strand2.begin(), strand2.end(), e);
if (esi == strand2.end())
continue;
auto bi = std::find(dssp.begin(), dssp.end(), b);
assert(bi != dssp.end() and bi != dssp.begin());
auto a = *std::next(bi);
auto c = *std::prev(bi);
auto ei = std::find(dssp.begin(), dssp.end(), e);
assert(ei != dssp.end() and ei != dssp.begin());
auto d = *std::next(ei);
auto f = *std::prev(ei);
if (test_bond(a, e) and test_bond(e, c)) // case I.
{
end1SeqID = a.seq_id();
end2SeqID = e.seq_id();
end1AtomID = "N";
end2AtomID = "O";
}
else if (test_bond(d, b) and test_bond(b, f)) // case II.
{
end1SeqID = b.seq_id();
end2SeqID = d.seq_id();
end1AtomID = "O";
end2AtomID = "N";
}
break;
}
if (end1SeqID)
break;
const auto &[strandID, strand] = s;
return strand.front().sheet() == nr;
})
}
}
else
{
// III. a <- f IV. a f antiparallel
//
// b e b <-> e <= the residues forming the bridge
//
// c -> d c d
std::reverse(strand2.begin(), strand2.end());
for (auto b : strand1)
{
for (int i : { 0, 1 })
{
const auto &[e, ladder1, parallel1] = b.bridge_partner(i);
auto esi = std::find(strand2.begin(), strand2.end(), e);
if (esi == strand2.end())
continue;
auto bi = std::find(dssp.begin(), dssp.end(), b);
assert(bi != dssp.end() and bi != dssp.begin());
auto a = *std::prev(bi);
auto c = *std::next(bi);
auto ei = std::find(dssp.begin(), dssp.end(), e);
assert(ei != dssp.end() and ei != dssp.begin());
auto d = *std::prev(ei);
auto f = *std::next(ei);
if (test_bond(f, a) and test_bond(c, d)) // case III.
{
beg1SeqID = a.seq_id();
beg2SeqID = f.seq_id();
beg1AtomID = "O";
beg2AtomID = "N";
}
else if (test_bond(b, e) and test_bond(e, b)) // case IV.
{
beg1SeqID = b.seq_id();
beg2SeqID = e.seq_id();
beg1AtomID = "N";
beg2AtomID = "O";
}
break;
}
if (beg1SeqID)
break;
}
std::reverse(strand1.begin(), strand1.end());
std::reverse(strand2.begin(), strand2.end());
for (auto b : strand1)
{
for (int i : { 0, 1 })
{
const auto &[e, ladder1, parallel1] = b.bridge_partner(i);
auto esi = std::find(strand2.begin(), strand2.end(), e);
if (esi == strand2.end())
continue;
auto bi = std::find(dssp.begin(), dssp.end(), b);
assert(bi != dssp.end() and bi != dssp.begin());
auto a = *std::next(bi);
auto c = *std::prev(bi);
auto ei = std::find(dssp.begin(), dssp.end(), e);
assert(ei != dssp.end() and ei != dssp.begin());
auto d = *std::next(ei);
auto f = *std::prev(ei);
if (test_bond(a, f) and test_bond(d, c)) // case III.
{
end1SeqID = a.seq_id();
end2SeqID = f.seq_id();
end1AtomID = "N";
end2AtomID = "O";
}
else if (test_bond(b, e) and test_bond(e, b)) // case IV.
{
end1SeqID = b.seq_id();
end2SeqID = e.seq_id();
end1AtomID = "N";
end2AtomID = "O";
}
break;
}
if (end1SeqID)
break;
}
}
struct_sheet_hbond.emplace({ { "sheet_id", cif::cif_id_for_number(sheet) },
{ "range_id_1", s1 + 1 },
{ "range_id_2", s2 + 1 },
{ "range_1_beg_label_seq_id", beg1SeqID },
{ "range_1_beg_label_atom_id", beg1AtomID },
{ "range_2_beg_label_seq_id", beg2SeqID },
{ "range_2_beg_label_atom_id", beg2AtomID },
{ "range_1_end_label_seq_id", end1SeqID },
{ "range_1_end_label_atom_id", end1AtomID },
{ "range_2_end_label_seq_id", end2SeqID },
{ "range_2_end_label_atom_id", end2AtomID } });
using namespace cif::literals;
auto b1 = pdbx_poly_seq_scheme.find_first("asym_id"_key == strand1.front().asym_id() and "seq_id"_key == beg1SeqID);
auto e1 = pdbx_poly_seq_scheme.find_first("asym_id"_key == strand1.front().asym_id() and "seq_id"_key == end1SeqID);
auto b2 = pdbx_poly_seq_scheme.find_first("asym_id"_key == strand2.front().asym_id() and "seq_id"_key == beg2SeqID);
auto e2 = pdbx_poly_seq_scheme.find_first("asym_id"_key == strand2.front().asym_id() and "seq_id"_key == end2SeqID);
if (not(b1 and e1 and b2 and e2))
{
if (cif::VERBOSE > 0)
std::cerr << "error looking up strand ends" << std::endl;
continue;
}
pdbx_struct_sheet_hbond.emplace({
{ "sheet_id", cif::cif_id_for_number(sheet) },
{ "range_id_1", s1 + 1 },
{ "range_id_2", s2 + 1 },
{ "range_1_label_atom_id", beg1AtomID },
{ "range_1_label_comp_id", b1["mon_id"].as<std::string>() },
{ "range_1_label_asym_id", b1["asym_id"].as<std::string>() },
{ "range_1_label_seq_id", b1["seq_id"].as<std::string>() },
{ "range_1_PDB_ins_code", b1["pdb_ins_code"].as<std::string>() },
{ "range_1_auth_atom_id", beg1AtomID },
{ "range_1_auth_comp_id", b1["auth_mon_id"].as<std::string>() },
{ "range_1_auth_asym_id", b1["pdb_strand_id"].as<std::string>() },
{ "range_1_auth_seq_id", b1["auth_seq_num"].as<std::string>() },
{ "range_2_label_atom_id", beg2AtomID },
{ "range_2_label_comp_id", b2["mon_id"].as<std::string>() },
{ "range_2_label_asym_id", b2["asym_id"].as<std::string>() },
{ "range_2_label_seq_id", b2["seq_id"].as<std::string>() },
{ "range_2_PDB_ins_code", b2["pdb_ins_code"].as<std::string>() },
{ "range_2_auth_atom_id", beg2AtomID },
{ "range_2_auth_comp_id", b2["auth_mon_id"].as<std::string>() },
{ "range_2_auth_asym_id", b2["pdb_strand_id"].as<std::string>() },
{ "range_2_auth_seq_id", b2["auth_seq_num"].as<std::string>() },
});
}
auto &struct_sheet_range = db["struct_sheet_range"];
for (const auto &[key, iSheet] : sheetMap)
{
for (auto &&[sheet, strand] : strands)
for (auto &&[strandTuple, strand] : strands)
{
if (sheet != iSheet)
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(); });
......@@ -759,8 +385,8 @@ void writeSheets(cif::datablock &db, const dssp &dssp)
auto &end = strand.back();
struct_sheet_range.emplace({
{ "sheet_id", cif::cif_id_for_number(sheet) },
{ "id", strandMap[strand.front().nr()] + 1 },
{ "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() },
......@@ -784,14 +410,21 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
// Write out the DSSP ladders
struct ladder_info
{
ladder_info(const std::string &label, bool parallel, const dssp::residue_info &a, const dssp::residue_info &b)
: label(label)
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 } })
{
}
std::string label;
bool operator<(const ladder_info &rhs) const
{
return ladder < rhs.ladder;
}
int ladder;
int sheet;
bool parallel;
std::vector<std::pair<dssp::residue_info, dssp::residue_info>> pairs;
};
......@@ -807,12 +440,10 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
if (not p)
continue;
auto label = cif::cif_id_for_number(ladder);
bool is_new = true;
for (auto &l : ladders)
{
if (l.label != label)
if (l.ladder != ladder)
continue;
assert(l.parallel == parallel);
......@@ -832,10 +463,12 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
if (not is_new)
continue;
ladders.emplace_back(label, parallel, res, p);
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)
......@@ -843,7 +476,11 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
const auto &[beg1, beg2] = l.pairs.front();
const auto &[end1, end2] = l.pairs.back();
dssp_struct_ladder.emplace({ { "id", l.label },
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() },
......@@ -975,7 +612,10 @@ void writeSummary(cif::datablock &db, const dssp &dssp)
// 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", "ladder_1", "ladder_2", "sheet", "accessibility", "TCO", "kappa", "alpha", "phi", "psi", "x_ca", "y_ca", "z_ca"})
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_column(label);
for (auto res : dssp)
......@@ -1062,11 +702,11 @@ void writeSummary(cif::datablock &db, const dssp &dssp)
{ "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] },
{ "sheet", res.sheet() ? cif::cif_id_for_number(res.sheet() - 1) : "." },
{ "x_ca", cax, 1 },
{ "y_ca", cay, 1 },
{ "z_ca", caz, 1 },
......@@ -1218,10 +858,6 @@ void annotateDSSP(cif::datablock &db, const dssp &dssp, bool writeOther, bool wr
{ "end_auth_comp_id", re.compound_id() },
{ "end_auth_asym_id", re.auth_asym_id() },
{ "end_auth_seq_id", re.auth_seq_id() }
// { "pdbx_PDB_helix_class", vS(39, 40) },
// { "details", vS(41, 70) },
// { "pdbx_PDB_helix_length", vI(72, 76) }
});
st = t;
......
src/dssp.cpp
......@@ -27,7 +27,6 @@
// Calculate DSSP-like secondary structure information
#include "dssp.hpp"
#include "queue.hpp"
#include "dssp-io.hpp"
......@@ -43,8 +42,6 @@ using helix_type = dssp::helix_type;
using helix_position_type = dssp::helix_position_type;
using chain_break_type = dssp::chain_break_type;
using queue_type = blocking_queue<std::tuple<uint32_t,uint32_t>>;
// --------------------------------------------------------------------
const double
......@@ -421,6 +418,9 @@ struct dssp::residue
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; }
......@@ -548,6 +548,7 @@ struct dssp::residue
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;
......@@ -1099,6 +1100,26 @@ void CalculateBetaSheets(std::vector<residue> &inResidues, statistics &stats, st
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;
}
}
}
// --------------------------------------------------------------------
......@@ -1487,7 +1508,10 @@ DSSP_impl::DSSP_impl(const cif::datablock &db, int model_nr, int min_poly_prolin
nextNext.mCAlpha,
cur.mCAlpha);
float skap = std::sqrt(1 - ckap * ckap);
cur.mKappa = std::atan2(skap, ckap) * static_cast<float>(180 / kPI);
auto kappa = std::atan2(skap, ckap) * static_cast<float>(180 / kPI);
if (not std::isnan(kappa))
cur.mKappa = kappa;
}
}
......@@ -2138,6 +2162,11 @@ 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, double> dssp::residue_info::acceptor(int i) const
{
auto &a = m_impl->mHBondAcceptor[i];
......
src/queue.hpp deleted 100644 → 0
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
* Copyright (c) 2021 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 <condition_variable>
#include <mutex>
#include <queue>
template <typename T, size_t N = 100>
class blocking_queue
{
public:
void push(T const &value)
{
std::unique_lock<std::mutex> lock(m_guard);
while (m_queue.size() >= N)
m_full_signal.wait(lock);
m_queue.push(value);
m_empty_signal.notify_one();
}
T pop()
{
std::unique_lock<std::mutex> lock(m_guard);
while (m_queue.empty())
m_empty_signal.wait(lock);
auto value = m_queue.front();
m_queue.pop();
m_full_signal.notify_one();
return value;
}
template<class Rep, class Period>
std::tuple<bool,T> pop(const std::chrono::duration<Rep, Period>& wait_for)
{
std::unique_lock<std::mutex> lock(m_guard);
while (m_queue.empty())
{
auto now = std::chrono::system_clock::now();
if (m_empty_signal.wait_until(lock, now + wait_for) == std::cv_status::timeout)
return { true , T{} };
}
auto value = m_queue.front();
m_queue.pop();
m_full_signal.notify_one();
return { false, value };
}
bool is_full() const
{
std::unique_lock<std::mutex> lock(m_guard);
return m_queue.size() >= N;
}
private:
std::queue<T> m_queue;
mutable std::mutex m_guard;
std::condition_variable m_empty_signal, m_full_signal;
};
template <typename T, size_t N = 10>
class non_blocking_queue
{
public:
bool push(T const &value)
{
bool result = false;
std::unique_lock<std::mutex> lock(m_guard);
if (m_queue.size() < N)
{
m_queue.push(value);
m_empty_signal.notify_one();
result = true;
}
return result;
}
T pop()
{
std::unique_lock<std::mutex> lock(m_guard);
while (m_queue.empty())
m_empty_signal.wait(lock);
auto value = m_queue.front();
m_queue.pop();
return value;
}
private:
std::queue<T> m_queue;
mutable std::mutex m_guard;
std::condition_variable m_empty_signal;
};
\ No newline at end of file
test/1cbs.dssp
==== Secondary Structure Definition by the program DSSP, NKI version 4.0 ==== DATE=2021-08-25 .
==== Secondary Structure Definition by the program DSSP, NKI version 4.4.0 ==== DATE=2021-08-25 .
REFERENCE W. KABSCH AND C.SANDER, BIOPOLYMERS 22 (1983) 2577-2637 .
HEADER RETINOIC-ACID TRANSPORT 28-SEP-94 1CBS .
COMPND MOL_ID: 1; MOLECULE: CELLULAR RETINOIC ACID BINDING PROTEIN TYPE II; CHAIN: A; ENGINEERED: YES .
......
test/unit-test-dssp.cpp
......@@ -32,6 +32,8 @@
#include "dssp.hpp"
#include "dssp-io.hpp"
#include <cif++/dictionary_parser.hpp>
namespace fs = std::filesystem;
// --------------------------------------------------------------------
......@@ -97,7 +99,7 @@ BOOST_AUTO_TEST_CASE(ut_dssp)
std::string line_t, line_r;
BOOST_CHECK(std::getline(test, line_t) and std::getline(reference, line_r));
const char *kHeaderLineStart = "==== Secondary Structure Definition by the program DSSP, NKI version 4.3 ====";
const char *kHeaderLineStart = "==== Secondary Structure Definition by the program DSSP, NKI version 4.4.0 ====";
BOOST_CHECK(line_t.compare(0, std::strlen(kHeaderLineStart), kHeaderLineStart) == 0);
// BOOST_CHECK(line_r.compare(0, std::strlen(kHeaderLineStart), kHeaderLineStart) == 0);
......@@ -225,4 +227,19 @@ BOOST_AUTO_TEST_CASE(dssp_2)
BOOST_CHECK_EQUAL(ri.seq_id(), seqID);
BOOST_CHECK_EQUAL((char)ri.type(), secstr.front());
}
}
// --------------------------------------------------------------------
BOOST_AUTO_TEST_CASE(dssp_3)
{
cif::file f(gTestDir / "1cbs.cif.gz");
BOOST_ASSERT(f.is_valid());
dssp dssp(f.front(), 1, 3, true);
dssp.annotate(f.front(), true, true);
BOOST_TEST(f.is_valid());
}
\ No newline at end of file
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