Skip to content

L
libcifpp
Commit 43eda65d by Maarten L. Hekkelman
Options

WIP on trunk: 7c5bf010 calculate surface only when needed, added deuterium

parents 7c5bf010 3a943847
Showing with 92 additions and 83 deletions
include/cif++/AtomType.hpp
......@@ -14,6 +14,7 @@ enum AtomType : uint8_t
Nn = 0, // Unknown
H = 1, // Hydro­gen
D = 129, // Deuterium
He = 2, // He­lium
Li = 3, // Lith­ium
......
include/cif++/Secondary.hpp
......@@ -91,7 +91,7 @@ enum class ChainBreak
class DSSP
{
public:
DSSP(const Structure& s, int min_poly_proline_stretch_length);
DSSP(const Structure& s, int min_poly_proline_stretch_length, bool calculateSurfaceAccessibility);
~DSSP();
DSSP(const DSSP&) = delete;
......
src/AtomType.cpp
......@@ -17,6 +17,7 @@ const AtomTypeInfo kKnownAtoms[] =
{
{ Nn, "Unknown", "Nn", 0, false, { kNA, kNA, kNA, kNA, kNA, kNA, kNA } }, // 0 Nn Unknown
{ H, "Hydrogen", "H", 1.008, false, { 53, 25, 37, 32, kNA, kNA, 120 } }, // 1 H Hydro­gen
{ D, "Deuterium", "D", 2.014, false, { 53, 25, 37, 32, kNA, kNA, 120 } }, // 1 D Deuterium
{ He, "Helium", "He", 4.0026, false, { 31, kNA, 32, 46, kNA, kNA, 140 } }, // 2 He He­lium
{ Li, "Lithium", "Li", 6.94, true, { 167, 145, 134, 133, 124, kNA, 182 } }, // 3 Li Lith­ium
{ Be, "Beryllium", "Be", 9.0122, true, { 112, 105, 90, 102, 90, 85, kNA } }, // 4 Be Beryl­lium
......
src/Secondary.cpp
......@@ -1150,6 +1150,9 @@ struct DSSPImpl
return std::find_if(mResidues.begin(), mResidues.end(), [&](auto& r) { return r.mM.asymID() == asymID and r.mM.seqID() == seqID; });
}
void calculateSurface();
void calculateSecondaryStructure();
DSSP_Statistics mStats = {};
};
......@@ -1208,116 +1211,112 @@ DSSPImpl::DSSPImpl(const Structure& s, int min_poly_proline_stretch_length)
mResidues[i + 1].assignHydrogen();
}
}
std::thread ta(std::bind(&CalculateAccessibilities, std::ref(mResidues), std::ref(mStats)));
try
void DSSPImpl::calculateSecondaryStructure()
{
auto& db = mStructure.getFile().data();
for (auto r: db["struct_conn"].find(cif::Key("conn_type_id") == "disulf"))
{
auto& db = s.getFile().data();
for (auto r: db["struct_conn"].find(cif::Key("conn_type_id") == "disulf"))
{
std::string asym1, asym2;
int seq1, seq2;
cif::tie(asym1, seq1, asym2, seq2) = r.get("ptnr1_label_asym_id", "ptnr1_label_seq_id", "ptnr2_label_asym_id", "ptnr2_label_seq_id");
std::string asym1, asym2;
int seq1, seq2;
cif::tie(asym1, seq1, asym2, seq2) = r.get("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())
throw std::runtime_error("Invalid file, missing residue for SS bond");
auto r1 = findRes(asym1, seq1);
if (r1 == mResidues.end())
throw std::runtime_error("Invalid file, missing residue for SS bond");
auto r2 = findRes(asym2, seq2);
if (r2 == mResidues.end())
throw std::runtime_error("Invalid file, missing residue for SS bond");
auto r2 = findRes(asym2, seq2);
if (r2 == mResidues.end())
throw std::runtime_error("Invalid file, missing residue for SS bond");
mSSBonds.emplace_back(&*r1, &*r2);
}
mSSBonds.emplace_back(&*r1, &*r2);
}
if (cif::VERBOSE) std::cerr << ".";
CalculateHBondEnergies(mResidues);
if (cif::VERBOSE) std::cerr << ".";
CalculateHBondEnergies(mResidues);
if (cif::VERBOSE) std::cerr << ".";
CalculateBetaSheets(mResidues, mStats);
if (cif::VERBOSE) std::cerr << ".";
CalculateBetaSheets(mResidues, mStats);
if (cif::VERBOSE) std::cerr << ".";
CalculateAlphaHelices(mResidues, mStats);
if (cif::VERBOSE) std::cerr << ".";
CalculateAlphaHelices(mResidues, mStats);
if (cif::VERBOSE) std::cerr << ".";
CalculatePPHelices(mResidues, mStats, m_min_poly_proline_stretch_length);
if (cif::VERBOSE) std::cerr << ".";
CalculatePPHelices(mResidues, mStats, m_min_poly_proline_stretch_length);
if (cif::VERBOSE) std::cerr << std::endl;
if (cif::VERBOSE) std::cerr << std::endl;
if (cif::VERBOSE > 1)
if (cif::VERBOSE > 1)
{
for (auto& r: mResidues)
{
for (auto& r: mResidues)
auto& m = r.mM;
char helix[5] = { };
for (HelixType helixType: { HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi, HelixType::rh_pp })
{
auto& m = r.mM;
char helix[5] = { };
for (HelixType helixType: { HelixType::rh_3_10, HelixType::rh_alpha, HelixType::rh_pi, HelixType::rh_pp })
switch (r.GetHelixFlag(helixType))
{
switch (r.GetHelixFlag(helixType))
{
case Helix::Start: helix[static_cast<int>(helixType)] = '>'; break;
case Helix::Middle: helix[static_cast<int>(helixType)] = helixType == HelixType::rh_pp ? 'P' : '3' + static_cast<int>(helixType); break;
case Helix::StartAndEnd: helix[static_cast<int>(helixType)] = 'X'; break;
case Helix::End: helix[static_cast<int>(helixType)] = '<'; break;
case Helix::None: helix[static_cast<int>(helixType)] = ' '; break;
}
case Helix::Start: helix[static_cast<int>(helixType)] = '>'; break;
case Helix::Middle: helix[static_cast<int>(helixType)] = helixType == HelixType::rh_pp ? 'P' : '3' + static_cast<int>(helixType); break;
case Helix::StartAndEnd: helix[static_cast<int>(helixType)] = 'X'; break;
case Helix::End: helix[static_cast<int>(helixType)] = '<'; break;
case Helix::None: helix[static_cast<int>(helixType)] = ' '; break;
}
auto id = m.asymID() + ':' + std::to_string(m.seqID()) + '/' + m.compoundID();
std::cerr << id << std::string(12 - id.length(), ' ')
<< char(r.mSecondaryStructure) << ' '
<< helix
<< std::endl;
}
auto id = m.asymID() + ':' + std::to_string(m.seqID()) + '/' + m.compoundID();
std::cerr << id << std::string(12 - id.length(), ' ')
<< char(r.mSecondaryStructure) << ' '
<< helix
<< std::endl;
}
}
// finish statistics
mStats.nrOfSSBridges = mSSBonds.size();
// finish statistics
mStats.nrOfSSBridges = mSSBonds.size();
mStats.nrOfIntraChainSSBridges = 0;
int ssBondNr = 0;
for (const auto& [a, b]: mSSBonds)
mStats.nrOfIntraChainSSBridges = 0;
int ssBondNr = 0;
for (const auto& [a, b]: mSSBonds)
{
if (a == b)
{
if (a == b)
{
if (cif::VERBOSE)
std::cerr << "In the SS bonds list, the residue " << a->mM << " is bonded to itself" << std::endl;
continue;
}
// throw std::runtime_error("first and second residue are the same");
if (cif::VERBOSE)
std::cerr << "In the SS bonds list, the residue " << a->mM << " is bonded to itself" << std::endl;
continue;
}
// throw std::runtime_error("first and second residue are the same");
if (a->mM.asymID() == b->mM.asymID() and NoChainBreak(a, b))
++mStats.nrOfIntraChainSSBridges;
if (a->mM.asymID() == b->mM.asymID() and NoChainBreak(a, b))
++mStats.nrOfIntraChainSSBridges;
a->mSSBridgeNr = b->mSSBridgeNr = ++ssBondNr;
}
a->mSSBridgeNr = b->mSSBridgeNr = ++ssBondNr;
}
mStats.nrOfHBonds = 0;
for (auto& r: mResidues)
{
auto donor = r.mHBondDonor;
mStats.nrOfHBonds = 0;
for (auto& r: mResidues)
{
auto donor = r.mHBondDonor;
for (int i = 0; i < 2; ++i)
for (int i = 0; i < 2; ++i)
{
if (donor[i].residue != nullptr and donor[i].energy < kMaxHBondEnergy)
{
if (donor[i].residue != nullptr and donor[i].energy < kMaxHBondEnergy)
{
++mStats.nrOfHBonds;
auto k = donor[i].residue->mNumber - r.mNumber;
if (k >= -5 and k <= 5)
mStats.nrOfHBondsPerDistance[k + 5] += 1;
}
++mStats.nrOfHBonds;
auto k = donor[i].residue->mNumber - r.mNumber;
if (k >= -5 and k <= 5)
mStats.nrOfHBondsPerDistance[k + 5] += 1;
}
}
}
catch (const std::exception& ex)
{
ta.join();
throw;
}
}
ta.join();
void DSSPImpl::calculateSurface()
{
CalculateAccessibilities(mResidues, mStats);
}
// --------------------------------------------------------------------
......@@ -1424,9 +1423,17 @@ DSSP::iterator& DSSP::iterator::operator++()
// --------------------------------------------------------------------
DSSP::DSSP(const Structure& s, int min_poly_proline_stretch)
DSSP::DSSP(const Structure& s, int min_poly_proline_stretch, bool calculateSurfaceAccessibility)
: mImpl(new DSSPImpl(s, min_poly_proline_stretch))
{
if (calculateSurfaceAccessibility)
{
std::thread t(std::bind(&DSSPImpl::calculateSurface, mImpl));
mImpl->calculateSecondaryStructure();
t.join();
}
else
mImpl->calculateSecondaryStructure();
}
DSSP::~DSSP()
......
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