Skip to content

D
dssp
Commit 0cedf7b1 by Maarten L. Hekkelman
Options

Alternate sheet writing, again

parent 7460007b
Showing with 116 additions and 84 deletions
src/dssp-io.cpp
...@@ -340,7 +340,7 @@ void writeSheets(cif::datablock &db, const dssp &dssp) ...@@ -340,7 +340,7 @@ void writeSheets(cif::datablock &db, const dssp &dssp)
for (auto &res : dssp) for (auto &res : dssp)
{ {
if (res.type() != dssp::structure_type::Strand) if (res.type() != dssp::structure_type::Strand and res.type() != dssp::structure_type::Betabridge)
continue; continue;
strands[{res.sheet(), res.strand()}].emplace_back(res); strands[{res.sheet(), res.strand()}].emplace_back(res);
...@@ -377,7 +377,7 @@ void writeSheets(cif::datablock &db, const dssp &dssp) ...@@ -377,7 +377,7 @@ void writeSheets(cif::datablock &db, const dssp &dssp)
if (strand.front().sheet() != sheetNr) if (strand.front().sheet() != sheetNr)
continue; continue;
std::string strandID = "strand_" + cif::cif_id_for_number(strand.front().strand() - 1); 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) std::sort(strand.begin(), strand.end(), [](dssp::residue_info const &a, dssp::residue_info const &b)
{ return a.nr() < b.nr(); }); { return a.nr() < b.nr(); });
...@@ -409,7 +409,8 @@ void writeSheets(cif::datablock &db, const dssp &dssp) ...@@ -409,7 +409,8 @@ void writeSheets(cif::datablock &db, const dssp &dssp)
for (int i : { 0, 1 }) for (int i : { 0, 1 })
{ {
const auto &[bp, ladder, parallel] = res.bridge_partner(i); const auto &[bp, ladder, parallel] = res.bridge_partner(i);
if (not bp or bp.type() != dssp::structure_type::Strand) // if (not bp or bp.type() != dssp::structure_type::Strand)
if (not bp)
continue; continue;
auto s1 = res.strand(); auto s1 = res.strand();
...@@ -431,13 +432,13 @@ void writeSheets(cif::datablock &db, const dssp &dssp) ...@@ -431,13 +432,13 @@ void writeSheets(cif::datablock &db, const dssp &dssp)
bool parallel = parallel_test > 0; bool parallel = parallel_test > 0;
std::string strandID_1 = "strand_" + cif::cif_id_for_number(s1 - 1); std::string strandID_1 = cif::cif_id_for_number(s1 - 1);
std::string strandID_2 = "strand_" + cif::cif_id_for_number(s2 - 1); std::string strandID_2 = cif::cif_id_for_number(s2 - 1);
struct_sheet_order.emplace({ struct_sheet_order.emplace({
{ "sheet_id", sheetID }, { "sheet_id", sheetID },
{ "range_id_1", "strand_" + strandID_1 }, { "range_id_1", strandID_1 },
{ "range_id_2", "strand_" + strandID_2 }, { "range_id_2", strandID_2 },
{ "sense", parallel > 0 ? "parallel" : "anti-parallel" } { "sense", parallel > 0 ? "parallel" : "anti-parallel" }
}); });
...@@ -748,14 +749,21 @@ void writeLadders(cif::datablock &db, const dssp &dssp) ...@@ -748,14 +749,21 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
// Write out the DSSP ladders // Write out the DSSP ladders
struct ladder_info struct ladder_info
{ {
ladder_info(const std::string &label, bool parallel, const dssp::residue_info &a, const dssp::residue_info &b) ladder_info(int label, int sheet, bool parallel, const dssp::residue_info &a, const dssp::residue_info &b)
: label(label) : ladder(label)
, sheet(sheet)
, parallel(parallel) , parallel(parallel)
, pairs({ { a, b } }) , pairs({ { a, b } })
{ {
} }
std::string label; bool operator<(const ladder_info &rhs) const
{
return ladder < rhs.ladder;
}
int ladder;
int sheet;
bool parallel; bool parallel;
std::vector<std::pair<dssp::residue_info, dssp::residue_info>> pairs; std::vector<std::pair<dssp::residue_info, dssp::residue_info>> pairs;
}; };
...@@ -771,12 +779,10 @@ void writeLadders(cif::datablock &db, const dssp &dssp) ...@@ -771,12 +779,10 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
if (not p) if (not p)
continue; continue;
auto label = cif::cif_id_for_number(ladder);
bool is_new = true; bool is_new = true;
for (auto &l : ladders) for (auto &l : ladders)
{ {
if (l.label != label) if (l.ladder != ladder)
continue; continue;
assert(l.parallel == parallel); assert(l.parallel == parallel);
...@@ -796,10 +802,12 @@ void writeLadders(cif::datablock &db, const dssp &dssp) ...@@ -796,10 +802,12 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
if (not is_new) if (not is_new)
continue; 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"]; auto &dssp_struct_ladder = db["dssp_struct_ladder"];
for (const auto &l : ladders) for (const auto &l : ladders)
...@@ -807,7 +815,11 @@ void writeLadders(cif::datablock &db, const dssp &dssp) ...@@ -807,7 +815,11 @@ void writeLadders(cif::datablock &db, const dssp &dssp)
const auto &[beg1, beg2] = l.pairs.front(); const auto &[beg1, beg2] = l.pairs.front();
const auto &[end1, end2] = l.pairs.back(); 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_1", cif::cif_id_for_number(beg1.strand() - 1) },
{ "range_2", cif::cif_id_for_number(beg2.strand() - 1) },
{ "type", l.parallel ? "parallel" : "anti-parallel" }, { "type", l.parallel ? "parallel" : "anti-parallel" },
{ "beg_1_label_comp_id", beg1.compound_id() }, { "beg_1_label_comp_id", beg1.compound_id() },
...@@ -939,7 +951,7 @@ void writeSummary(cif::datablock &db, const dssp &dssp) ...@@ -939,7 +951,7 @@ 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. // 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); dssp_struct_summary.add_column(label);
for (auto res : dssp) for (auto res : dssp)
...@@ -1026,11 +1038,11 @@ void writeSummary(cif::datablock &db, const dssp &dssp) ...@@ -1026,11 +1038,11 @@ void writeSummary(cif::datablock &db, const dssp &dssp)
{ "bend", bend }, { "bend", bend },
{ "chirality", chirality }, { "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_1", ladders[0] },
{ "ladder_2", ladders[1] }, { "ladder_2", ladders[1] },
{ "sheet", res.sheet() ? cif::cif_id_for_number(res.sheet() - 1) : "." },
{ "x_ca", cax, 1 }, { "x_ca", cax, 1 },
{ "y_ca", cay, 1 }, { "y_ca", cay, 1 },
{ "z_ca", caz, 1 }, { "z_ca", caz, 1 },
......
src/dssp.cpp
...@@ -1104,85 +1104,105 @@ void CalculateBetaSheets(std::vector<residue> &inResidues, statistics &stats, st ...@@ -1104,85 +1104,105 @@ void CalculateBetaSheets(std::vector<residue> &inResidues, statistics &stats, st
} }
} }
// Construct the 'strands' // Second attempt to create 'strands'. A strand is a range of residues
// For mmCIF output, this is needed and since we now have the information available // without a gap in between in the same chain that belong to the same sheet.
// it is best to do the calculation here.
// strands are ranges of residues of length > 1 that form beta bridges in a sheet.
int strand = 0;
for (uint32_t iSheet = 1; iSheet < sheet; ++iSheet) for (uint32_t iSheet = 1; iSheet < sheet; ++iSheet)
{ {
std::vector<std::tuple<uint32_t,uint32_t>> strands; int lastNr = -1;
for (auto &bridge : bridges) for (auto &res : inResidues)
{ {
if (bridge.sheet != iSheet) if (res.mSheet != iSheet)
continue; continue;
for (auto &range : { bridge.i, bridge.j}) if (lastNr + 1 < res.mNumber)
{ ++strand;
auto imin = range.front();
auto imax = range.back();
if (imin == imax)
continue;
if (imin > imax)
std::swap(imin, imax);
auto ii = find_if(strands.begin(), strands.end(), [a = imin, b = imax] (std::tuple<uint32_t,uint32_t> &t)
{
auto &&[start, end] = t;
bool result = false;
if (start <= b and end >= a)
{
result = true;
if (start > a)
start = a;
if (end < b)
end = b;
}
return result; res.mStrand = strand;
}); lastNr = res.mNumber;
if (ii == strands.end())
strands.emplace_back(imin, imax);
} }
} }
std::sort(strands.begin(), strands.end()); // // Construct the 'strands'
// // For mmCIF output, this is needed and since we now have the information available
// collapse ranges that overlap // // it is best to do the calculation here.
if (strands.size() > 1)
{
auto si = strands.begin();
while (std::next(si) != strands.end())
{
auto &&[afirst, alast] = *si;
auto &&[bfirst, blast] = *(std::next(si));
if (alast >= bfirst) // // strands are ranges of residues of length > 1 that form beta bridges in a sheet.
{
bfirst = afirst;
si = strands.erase(si);
continue;
}
++si; // for (uint32_t iSheet = 1; iSheet < sheet; ++iSheet)
} // {
} // std::vector<std::tuple<uint32_t,uint32_t>> strands;
// for (auto &bridge : bridges)
// {
// if (bridge.sheet != iSheet)
// continue;
for (size_t i = 0; i < strands.size(); ++i) // for (auto &range : { bridge.i, bridge.j})
{ // {
const auto &[first, last] = strands[i]; // auto imin = range.front();
for (auto nr = first; nr <= last; ++nr) // auto imax = range.back();
{
assert(inResidues[nr].mStrand == 0); // // if (imin == imax)
inResidues[nr].SetStrand(i + 1); // // continue;
}
} // if (imin > imax)
} // std::swap(imin, imax);
// auto ii = find_if(strands.begin(), strands.end(), [a = imin, b = imax] (std::tuple<uint32_t,uint32_t> &t)
// {
// auto &&[start, end] = t;
// bool result = false;
// if (start <= b and end >= a)
// {
// result = true;
// if (start > a)
// start = a;
// if (end < b)
// end = b;
// }
// return result;
// });
// if (ii == strands.end())
// strands.emplace_back(imin, imax);
// }
// }
// std::sort(strands.begin(), strands.end());
// // collapse ranges that overlap
// if (strands.size() > 1)
// {
// auto si = strands.begin();
// while (std::next(si) != strands.end())
// {
// auto &&[afirst, alast] = *si;
// auto &&[bfirst, blast] = *(std::next(si));
// if (alast >= bfirst)
// {
// bfirst = afirst;
// si = strands.erase(si);
// continue;
// }
// ++si;
// }
// }
// for (size_t i = 0; i < strands.size(); ++i)
// {
// const auto &[first, last] = strands[i];
// for (auto nr = first; nr <= last; ++nr)
// {
// assert(inResidues[nr].mStrand == 0);
// inResidues[nr].SetStrand(i + 1);
// }
// }
// }
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
......
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