Skip to content

L
libcifpp
Commit a9671546 by Maarten L. Hekkelman
Options

PP helix assignment

parent 78dd9a3c
Showing with 82 additions and 164 deletions
include/cif++/Secondary.hpp
...@@ -91,7 +91,7 @@ enum class ChainBreak ...@@ -91,7 +91,7 @@ enum class ChainBreak
class DSSP class DSSP
{ {
public: public:
DSSP(const Structure& s); DSSP(const Structure& s, int min_poly_proline_stretch_length);
~DSSP(); ~DSSP();
DSSP(const DSSP&) = delete; DSSP(const DSSP&) = delete;
......
src/Secondary.cpp
...@@ -1000,9 +1000,8 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, ...@@ -1000,9 +1000,8 @@ void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats,
// -------------------------------------------------------------------- // --------------------------------------------------------------------
void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats) void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, int stretch_length)
{ {
// less strict variant: three residues instead of four
size_t N = inResidues.size(); size_t N = inResidues.size();
const float epsilon = 29; const float epsilon = 29;
...@@ -1021,213 +1020,131 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats) ...@@ -1021,213 +1020,131 @@ void CalculatePPHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats)
for (uint32_t i = 1; i + 3 < inResidues.size(); ++i) for (uint32_t i = 1; i + 3 < inResidues.size(); ++i)
{ {
if (phi_min > phi[i - 1] or phi[i - 1] > phi_max or switch (stretch_length)
phi_min > phi[i + 0] or phi[i + 0] > phi_max) {
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; continue;
if (psi_min > psi[i - 1] or psi[i - 1] > psi_max or if (psi_min > psi[i + 0] or psi[i + 0] > psi_max or
psi_min > psi[i + 0] or psi[i + 0] > psi_max) psi_min > psi[i + 1] or psi[i + 1] > psi_max)
continue; continue;
auto phi_avg = (phi[i - 1] + phi[i]) / 2; // auto phi_avg = (phi[i + 0] + phi[i + 1]) / 2;
auto phi_sq = (phi[i - 1] - phi_avg) * (phi[i - 1] - phi_avg) + // auto phi_sq = (phi[i + 0] - phi_avg) * (phi[i + 0] - phi_avg) +
(phi[i + 0] - phi_avg) * (phi[i + 0] - phi_avg); // (phi[i + 1] - phi_avg) * (phi[i + 1] - phi_avg);
if (phi_sq >= 200) // if (phi_sq >= 200)
continue; // continue;
auto psi_avg = (psi[i - 1] + psi[i]) / 2; // auto psi_avg = (psi[i + 0] + psi[i + 1]) / 2;
auto psi_sq = (psi[i - 1] - psi_avg) * (psi[i - 1] - psi_avg) + // auto psi_sq = (psi[i + 0] - psi_avg) * (psi[i + 0] - psi_avg) +
(psi[i + 0] - psi_avg) * (psi[i + 0] - psi_avg); // (psi[i + 1] - psi_avg) * (psi[i + 1] - psi_avg);
if (psi_sq >= 200) // if (psi_sq >= 200)
continue; // continue;
switch (inResidues[i - 1].GetHelixFlag(HelixType::rh_pp)) switch (inResidues[i].GetHelixFlag(HelixType::rh_pp))
{ {
case Helix::None: case Helix::None:
inResidues[i - 1].SetHelixFlag(HelixType::rh_pp, Helix::Start); inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Start);
break; break;
case Helix::End: case Helix::End:
inResidues[i - 1].SetHelixFlag(HelixType::rh_pp, Helix::StartAndEnd); inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Middle);
break; break;
default: default:
break; break;
} }
inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Middle);
inResidues[i + 1].SetHelixFlag(HelixType::rh_pp, Helix::End); inResidues[i + 1].SetHelixFlag(HelixType::rh_pp, Helix::End);
if (inResidues[i - 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop and if (inResidues[i].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
inResidues[i + 0].GetSecondaryStructure() == SecondaryStructureType::ssLoop and inResidues[i].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
inResidues[i + 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop) if (inResidues[i + 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
{
inResidues[i - 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
inResidues[i + 0].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
inResidues[i + 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII); inResidues[i + 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_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;
// // less strict variant: three residues instead of four if (psi_min > psi[i + 0] or psi[i + 0] > psi_max or
// size_t N = inResidues.size(); psi_min > psi[i + 1] or psi[i + 1] > psi_max or
psi_min > psi[i + 2] or psi[i + 2] > psi_max)
continue;
// std::vector<float> phi(N), psi(N); // 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);
// for (uint32_t i = 1; i + 1 < inResidues.size(); ++i) // if (phi_sq >= 300)
// { // continue;
// // if (not NoChainBreak(inResidues[i], inResidues[i]))
// // continue;
// phi[i] = inResidues[i].mM.phi(); // auto psi_avg = (psi[i + 0] + psi[i + 1] + psi[i + 2]) / 3;
// psi[i] = inResidues[i].mM.psi(); // 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);
// for (uint32_t i = 1; i + 3 < inResidues.size(); ++i) // if (psi_sq >= 300)
// {
// if (phi_min > phi[i - 1] or phi[i - 1] > phi_max or
// phi_min > phi[i + 0] or phi[i + 0] > phi_max or
// phi_min > phi[i + 1] or phi[i + 1] > phi_max)
// continue; // continue;
// if (psi_min > psi[i - 1] or psi[i - 1] > psi_max or switch (inResidues[i].GetHelixFlag(HelixType::rh_pp))
// psi_min > psi[i + 0] or psi[i + 0] > psi_max or {
// psi_min > psi[i + 1] or psi[i + 1] > psi_max) case Helix::None:
// continue; inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Start);
break;
// // auto phi_avg = (phi[i - 1] + phi[i] + phi[i + 1]) / 3; case Helix::End:
// // auto phi_sq = (phi[i - 1] - phi_avg) * (phi[i - 1] - phi_avg) + inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::StartAndEnd);
// // (phi[i + 0] - phi_avg) * (phi[i + 0] - phi_avg) + break;
// // (phi[i + 1] - phi_avg) * (phi[i + 1] - phi_avg);
// // if (phi_sq >= 300)
// // continue;
// // auto psi_avg = (psi[i - 1] + psi[i] + psi[i + 1]) / 3;
// // auto psi_sq = (psi[i - 1] - psi_avg) * (psi[i - 1] - psi_avg) +
// // (psi[i + 0] - psi_avg) * (psi[i + 0] - psi_avg) +
// // (psi[i + 1] - psi_avg) * (psi[i + 1] - psi_avg);
// // if (psi_sq >= 300)
// // continue;
// switch (inResidues[i - 1].GetHelixFlag(HelixType::rh_pp))
// {
// case Helix::None:
// inResidues[i - 1].SetHelixFlag(HelixType::rh_pp, Helix::Start);
// break;
// case Helix::End:
// inResidues[i - 1].SetHelixFlag(HelixType::rh_pp, Helix::StartAndEnd);
// break;
// default:
// break;
// }
// inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Middle);
// inResidues[i + 1].SetHelixFlag(HelixType::rh_pp, Helix::End);
// if (inResidues[i - 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop and
// inResidues[i + 0].GetSecondaryStructure() == SecondaryStructureType::ssLoop and
// inResidues[i + 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
// {
// inResidues[i - 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// inResidues[i + 0].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// inResidues[i + 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// }
// }
// size_t N = inResidues.size();
// std::vector<float> phi(N), psi(N);
// for (uint32_t i = 1; i + 1 < inResidues.size(); ++i)
// {
// // if (not NoChainBreak(inResidues[i], inResidues[i]))
// // continue;
// phi[i] = inResidues[i].mM.phi();
// psi[i] = inResidues[i].mM.psi();
// }
// for (uint32_t i = 1; i + 3 < inResidues.size(); ++i)
// {
// if (phi_min > phi[i - 1] or phi[i - 1] > phi_max or
// 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 - 1] or psi[i - 1] > psi_max or default:
// psi_min > psi[i + 0] or psi[i + 0] > psi_max or break;
// 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 - 1] + phi[i] + phi[i + 1] + phi[i + 2]) / 4; inResidues[i + 1].SetHelixFlag(HelixType::rh_pp, Helix::Middle);
// auto phi_sq = (phi[i - 1] - phi_avg) * (phi[i - 1] - phi_avg) + inResidues[i + 2].SetHelixFlag(HelixType::rh_pp, Helix::End);
// (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 >= 400) if (inResidues[i + 0].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
// continue; inResidues[i + 0].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// auto psi_avg = (psi[i - 1] + psi[i] + psi[i + 1] + psi[i + 2]) / 4; if (inResidues[i + 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
// auto psi_sq = (psi[i - 1] - psi_avg) * (psi[i - 1] - psi_avg) + inResidues[i + 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// (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 >= 400) if (inResidues[i + 2].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
// continue; inResidues[i + 2].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// switch (inResidues[i - 1].GetHelixFlag(HelixType::rh_pp)) break;
// { }
// case Helix::None:
// inResidues[i - 1].SetHelixFlag(HelixType::rh_pp, Helix::Start); default:
// break; throw std::runtime_error("Unsupported stretch length");
}
// case Helix::End: }
// inResidues[i - 1].SetHelixFlag(HelixType::rh_pp, Helix::StartAndEnd);
// break;
// default:
// break;
// }
// inResidues[i].SetHelixFlag(HelixType::rh_pp, Helix::Middle);
// inResidues[i + 1].SetHelixFlag(HelixType::rh_pp, Helix::Middle);
// inResidues[i + 2].SetHelixFlag(HelixType::rh_pp, Helix::End);
// if (inResidues[i - 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop and
// inResidues[i + 0].GetSecondaryStructure() == SecondaryStructureType::ssLoop and
// inResidues[i + 1].GetSecondaryStructure() == SecondaryStructureType::ssLoop and
// inResidues[i + 2].GetSecondaryStructure() == SecondaryStructureType::ssLoop)
// {
// inResidues[i - 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// inResidues[i + 0].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// inResidues[i + 1].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// inResidues[i + 2].SetSecondaryStructure(SecondaryStructureType::ssHelix_PPII);
// }
// }
} }
// -------------------------------------------------------------------- // --------------------------------------------------------------------
struct DSSPImpl struct DSSPImpl
{ {
DSSPImpl(const Structure& s); DSSPImpl(const Structure& s, int min_poly_proline_stretch_length);
const Structure& mStructure; const Structure& mStructure;
const std::list<Polymer>& mPolymers; const std::list<Polymer>& mPolymers;
std::vector<Res> mResidues; std::vector<Res> mResidues;
std::vector<std::pair<Res*,Res*>> mSSBonds; std::vector<std::pair<Res*,Res*>> mSSBonds;
int m_min_poly_proline_stretch_length;
auto findRes(const std::string& asymID, int seqID) auto findRes(const std::string& asymID, int seqID)
{ {
return std::find_if(mResidues.begin(), mResidues.end(), [&](auto& r) { return r.mM.asymID() == asymID and r.mM.seqID() == seqID; }); return std::find_if(mResidues.begin(), mResidues.end(), [&](auto& r) { return r.mM.asymID() == asymID and r.mM.seqID() == seqID; });
...@@ -1238,9 +1155,10 @@ struct DSSPImpl ...@@ -1238,9 +1155,10 @@ struct DSSPImpl
// -------------------------------------------------------------------- // --------------------------------------------------------------------
DSSPImpl::DSSPImpl(const Structure& s) DSSPImpl::DSSPImpl(const Structure& s, int min_poly_proline_stretch_length)
: mStructure(s) : mStructure(s)
, mPolymers(mStructure.polymers()) , mPolymers(mStructure.polymers())
, m_min_poly_proline_stretch_length(min_poly_proline_stretch_length)
{ {
if (cif::VERBOSE) if (cif::VERBOSE)
std::cerr << "Calculating DSSP "; std::cerr << "Calculating DSSP ";
...@@ -1321,7 +1239,7 @@ DSSPImpl::DSSPImpl(const Structure& s) ...@@ -1321,7 +1239,7 @@ DSSPImpl::DSSPImpl(const Structure& s)
CalculateAlphaHelices(mResidues, mStats); CalculateAlphaHelices(mResidues, mStats);
if (cif::VERBOSE) std::cerr << "."; if (cif::VERBOSE) std::cerr << ".";
CalculatePPHelices(mResidues, mStats); CalculatePPHelices(mResidues, mStats, m_min_poly_proline_stretch_length);
if (cif::VERBOSE) std::cerr << std::endl; if (cif::VERBOSE) std::cerr << std::endl;
...@@ -1493,8 +1411,8 @@ DSSP::iterator& DSSP::iterator::operator++() ...@@ -1493,8 +1411,8 @@ DSSP::iterator& DSSP::iterator::operator++()
// -------------------------------------------------------------------- // --------------------------------------------------------------------
DSSP::DSSP(const Structure& s) DSSP::DSSP(const Structure& s, int min_poly_proline_stretch)
: mImpl(new DSSPImpl(s)) : mImpl(new DSSPImpl(s, min_poly_proline_stretch))
{ {
} }
......
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