Skip to content

L
libcifpp
Commit 38a852ca by maarten
Options

backup 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@277 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent a5b64ec8
Showing with 493 additions and 46 deletions
include/cif++/Compound.h
...@@ -184,6 +184,8 @@ struct LinkAtom ...@@ -184,6 +184,8 @@ struct LinkAtom
{ {
int compID; int compID;
std::string atomID; std::string atomID;
bool operator==(const LinkAtom& rhs) const { return compID == rhs.compID and atomID == rhs.atomID; }
}; };
struct LinkBond struct LinkBond
...@@ -261,6 +263,10 @@ class Link ...@@ -261,6 +263,10 @@ class Link
std::vector<LinkPlane> planes() const { return mPlanes; } std::vector<LinkPlane> planes() const { return mPlanes; }
std::vector<LinkTorsion> torsions() const { return mTorsions; } std::vector<LinkTorsion> torsions() const { return mTorsions; }
float atomBondValue(const LinkAtom& atomId_1, const LinkAtom& atomId_2) const;
float bondAngle(const LinkAtom& atomId_1, const LinkAtom& atomId_2, const LinkAtom& atomId_3) const;
float chiralVolume(const std::string& id) const;
private: private:
~Link(); ~Link();
......
include/cif++/MapMaker.h
...@@ -93,6 +93,12 @@ class MapMaker ...@@ -93,6 +93,12 @@ class MapMaker
std::initializer_list<std::string> foLabels = { "FP", "SIGFP" }, std::initializer_list<std::string> foLabels = { "FP", "SIGFP" },
std::initializer_list<std::string> freeLabels = { "FREE" }); std::initializer_list<std::string> freeLabels = { "FREE" });
void recalc(const Structure& structure,
bool noBulk = false, AnisoScalingFlag anisoScaling = as_None,
float samplingRate = 4.5, bool electronScattering = false);
void printStats();
void writeMTZ(const boost::filesystem::path& file, void writeMTZ(const boost::filesystem::path& file,
const std::string& project, const std::string& crystal); const std::string& project, const std::string& crystal);
...@@ -109,6 +115,7 @@ class MapMaker ...@@ -109,6 +115,7 @@ class MapMaker
const Cell& cell() const { return mHKLInfo.cell(); } const Cell& cell() const { return mHKLInfo.cell(); }
const Grid_sampling& gridSampling() const { return mGrid; } const Grid_sampling& gridSampling() const { return mGrid; }
private: private:
void loadFoFreeFromReflectionsFile(const boost::filesystem::path& hklin); void loadFoFreeFromReflectionsFile(const boost::filesystem::path& hklin);
...@@ -116,12 +123,7 @@ class MapMaker ...@@ -116,12 +123,7 @@ class MapMaker
std::initializer_list<std::string> foLabels, std::initializer_list<std::string> foLabels,
std::initializer_list<std::string> freeLabels); std::initializer_list<std::string> freeLabels);
void recalc(const Structure& structure,
bool noBulk, AnisoScalingFlag anisoScaling,
float samplingRate = 4.5, bool electronScattering = false);
void fixMTZ(); void fixMTZ();
void printStats();
MapType mFb, mFd; MapType mFb, mFd;
Grid_sampling mGrid; Grid_sampling mGrid;
......
include/cif++/Structure.h
...@@ -142,6 +142,11 @@ class Atom ...@@ -142,6 +142,11 @@ class Atom
struct AtomImpl* mImpl; struct AtomImpl* mImpl;
}; };
inline double Distance(const Atom& a, const Atom& b)
{
return Distance(a.location(), b.location());
}
typedef std::vector<Atom> AtomView; typedef std::vector<Atom> AtomView;
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -176,6 +181,8 @@ class Residue ...@@ -176,6 +181,8 @@ class Residue
const Structure& structure() const { return *mStructure; } const Structure& structure() const { return *mStructure; }
protected: protected:
const Structure* mStructure; const Structure* mStructure;
...@@ -185,6 +192,45 @@ class Residue ...@@ -185,6 +192,45 @@ class Residue
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// We have DSSP built-in
struct HBond
{
Monomer* residue;
double energy;
};
enum BridgeType
{
btNoBridge, btParallel, btAntiParallel
};
struct BridgeParner
{
Monomer* residue;
uint32 ladder;
bool parallel;
};
enum HelixFlag
{
helixNone, helixStart, helixEnd, helixStartAndEnd, helixMiddle
};
enum SecondaryStructure : char
{
undefined = 0,
loop = ' ',
alphahelix = 'H',
betabridge = 'B',
strand = 'E',
helix_3 = 'G',
helix_5 = 'I',
turn = 'T',
bend = 'S'
};
// --------------------------------------------------------------------
// a monomer models a single Residue in a protein chain // a monomer models a single Residue in a protein chain
class Monomer : public Residue class Monomer : public Residue
...@@ -211,9 +257,23 @@ class Monomer : public Residue ...@@ -211,9 +257,23 @@ class Monomer : public Residue
Atom O() const { return atomByID("O"); } Atom O() const { return atomByID("O"); }
Atom H() const { return atomByID("H"); } Atom H() const { return atomByID("H"); }
SecondaryStructure secondaryStructure() const { return mSS; }
static double CalculateHBondEnergy(Monomer& donor, Monomer& acceptor);
private: private:
Polymer* mPolymer; Polymer* mPolymer;
uint32 mIndex; uint32 mIndex;
// DSSP info
SecondaryStructure mSS;
uint8 mSSBridgeNr;
// double mAccessibility;
HBond mHBondDonor[2], mHBondAcceptor[2];
BridgeParner mBetaPartner[2];
uint32 mSheet;
HelixFlag mHelixFlags[3]; //
bool mBend;
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -372,6 +432,9 @@ class Structure ...@@ -372,6 +432,9 @@ class Structure
void changeResidue(Residue& res, const std::string& newCompound, void changeResidue(Residue& res, const std::string& newCompound,
const std::vector<std::tuple<std::string,std::string>>& remappedAtoms); const std::vector<std::tuple<std::string,std::string>>& remappedAtoms);
// DSSP
void CalculateSecondaryStructure(bool inPreferPiHelices = true);
private: private:
friend Polymer; friend Polymer;
friend Residue; friend Residue;
......
src/Compound.cpp
...@@ -750,39 +750,6 @@ float Compound::chiralVolume(const string& centreID) const ...@@ -750,39 +750,6 @@ float Compound::chiralVolume(const string& centreID) const
float cosc = cos(gamma * kPI / 180); float cosc = cos(gamma * kPI / 180);
result = (a * b * c * sqrt(1 + 2 * cosa * cosb * cosc - (cosa * cosa) - (cosb * cosb) - (cosc * cosc))) / 6; result = (a * b * c * sqrt(1 + 2 * cosa * cosb * cosc - (cosa * cosa) - (cosb * cosb) - (cosc * cosc))) / 6;
//
//
// // the bond angles for the top of the tetrahedron
// float b[3] =
// {
// bondAngle(cv.atomID[0], cv.atomIDCentre, cv.atomID[1]),
// bondAngle(cv.atomID[1], cv.atomIDCentre, cv.atomID[2]),
// bondAngle(cv.atomID[2], cv.atomIDCentre, cv.atomID[0])
// };
//
// // the edges
//
// float a[6] =
// {
// atomBondValue(cv.atomIDCentre, cv.atomID[0]),
// atomBondValue(cv.atomIDCentre, cv.atomID[1]),
// atomBondValue(cv.atomIDCentre, cv.atomID[2])
// };
//
// a[3] = calcC(a[0], a[1], b[0]);
// a[4] = calcC(a[1], a[2], b[1]);
// a[5] = calcC(a[2], a[0], b[2]);
//
// for (auto& i: a)
// i *= i;
//
// float v2 =
// a[0] * a[5] * (a[1] + a[2] + a[3] + a[5] - a[0] - a[4]) +
// a[1] * a[5] * (a[0] + a[2] + a[3] + a[4] - a[1] - a[5]) +
// a[2] * a[3] * (a[0] + a[1] + a[4] + a[5] - a[2] - a[3]) -
// a[0] * a[1] * a[3] - a[1] * a[2] * a[4] - a[0] * a[2] * a[5] - a[3] * a[4] * a[5];
//
// result = sqrt(v2 / 144);
if (cv.volumeSign == negativ) if (cv.volumeSign == negativ)
result = -result; result = -result;
...@@ -926,6 +893,74 @@ Link::~Link() ...@@ -926,6 +893,74 @@ Link::~Link()
{ {
} }
float Link::atomBondValue(const LinkAtom& atom1, const LinkAtom& atom2) const
{
auto i = find_if(mBonds.begin(), mBonds.end(),
[&](auto& b)
{
return (b.atom[0] == atom1 and b.atom[1] == atom2)
or (b.atom[0] == atom2 and b.atom[1] == atom1);
});
return i != mBonds.end() ? i->distance : 0;
}
float Link::bondAngle(const LinkAtom& atom1, const LinkAtom& atom2, const LinkAtom& atom3) const
{
float result = nanf("1");
for (auto& a: mAngles)
{
if (not (a.atom[1] == atom2 and
((a.atom[0] == atom1 and a.atom[2] == atom3) or
(a.atom[2] == atom1 and a.atom[0] == atom3))))
continue;
result = a.angle;
break;
}
return result;
}
float Link::chiralVolume(const string& centreID) const
{
float result = 0;
for (auto& cv: mChiralCentres)
{
if (cv.id != centreID)
continue;
// calculate the expected chiral volume
// the edges
float a = atomBondValue(cv.atomCentre, cv.atom[0]);
float b = atomBondValue(cv.atomCentre, cv.atom[1]);
float c = atomBondValue(cv.atomCentre, cv.atom[2]);
// the angles for the top of the tetrahedron
float alpha = bondAngle(cv.atom[0], cv.atomCentre, cv.atom[1]);
float beta = bondAngle(cv.atom[1], cv.atomCentre, cv.atom[2]);
float gamma = bondAngle(cv.atom[2], cv.atomCentre, cv.atom[0]);
float cosa = cos(alpha * kPI / 180);
float cosb = cos(beta * kPI / 180);
float cosc = cos(gamma * kPI / 180);
result = (a * b * c * sqrt(1 + 2 * cosa * cosb * cosc - (cosa * cosa) - (cosb * cosb) - (cosc * cosc))) / 6;
if (cv.volumeSign == negativ)
result = -result;
break;
}
return result;
}
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// a factory class to generate compounds // a factory class to generate compounds
......
src/DistanceMap.cpp
...@@ -71,7 +71,6 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -71,7 +71,6 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
dim = atoms.size(); dim = atoms.size();
vector<clipper::Coord_orth> locations(dim); vector<clipper::Coord_orth> locations(dim);
vector<bool> isWater(dim, false);
// bounding box // bounding box
Point pMin(numeric_limits<float>::max(), numeric_limits<float>::max(), numeric_limits<float>::max()), Point pMin(numeric_limits<float>::max(), numeric_limits<float>::max(), numeric_limits<float>::max()),
...@@ -82,8 +81,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -82,8 +81,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
size_t ix = index.size(); size_t ix = index.size();
index[atom.id()] = ix; index[atom.id()] = ix;
locations[ix] = atom.location(); locations[ix] = toCell(atom.location(), cell);
isWater[ix] = atom.isWater();
auto p = atom.location(); auto p = atom.location();
...@@ -125,19 +123,16 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro ...@@ -125,19 +123,16 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
if (i >= locations.size()) if (i >= locations.size())
break; break;
clipper::Coord_orth pi = toCell(locations[i], cell); clipper::Coord_orth pi = locations[i];
for (size_t j = i + 1; j < locations.size(); ++j) for (size_t j = i + 1; j < locations.size(); ++j)
{ {
// if (not (isWater[i] or isWater[j]))
// continue;
// find nearest location based on spacegroup/cell // find nearest location based on spacegroup/cell
double minR2 = numeric_limits<double>::max(); double minR2 = numeric_limits<double>::max();
for (auto rt: rtOrth) for (auto rt: rtOrth)
{ {
auto pj = toCell(locations[j], cell); auto pj = locations[j];
pj = pj.transform(rt); pj = pj.transform(rt);
double r2 = (pi - pj).lengthsq(); double r2 = (pi - pj).lengthsq();
......
src/Structure.cpp
...@@ -25,6 +25,18 @@ extern int VERBOSE; ...@@ -25,6 +25,18 @@ extern int VERBOSE;
namespace mmcif namespace mmcif
{ {
// DSSP constants
const double
kSSBridgeDistance = 3.0,
kMinimalDistance = 0.5,
kMinimalCADistance = 9.0,
kMinHBondEnergy = -9.9,
kMaxHBondEnergy = -0.5,
kCouplingConstant = -27.888, // = -332 * 0.42 * 0.2
kMaxPeptideBondLength = 2.5;
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// FileImpl // FileImpl
...@@ -942,6 +954,8 @@ struct StructureImpl ...@@ -942,6 +954,8 @@ struct StructureImpl
AtomView mAtoms; AtomView mAtoms;
}; };
// --------------------------------------------------------------------
void StructureImpl::insertCompound(const string& compoundID, bool isEntity) void StructureImpl::insertCompound(const string& compoundID, bool isEntity)
{ {
auto compound = Compound::create(compoundID); auto compound = Compound::create(compoundID);
...@@ -1084,6 +1098,8 @@ void StructureImpl::changeResidue(Residue& res, const string& newCompound, ...@@ -1084,6 +1098,8 @@ void StructureImpl::changeResidue(Residue& res, const string& newCompound,
} }
} }
// --------------------------------------------------------------------
Structure::Structure(File& f, uint32 modelNr) Structure::Structure(File& f, uint32 modelNr)
: mImpl(new StructureImpl(*this, f, modelNr)) : mImpl(new StructureImpl(*this, f, modelNr))
{ {
...@@ -1166,6 +1182,336 @@ vector<Residue> Structure::nonPolymers() const ...@@ -1166,6 +1182,336 @@ vector<Residue> Structure::nonPolymers() const
return result; return result;
} }
double Monomer::CalculateHBondEnergy(Monomer& inDonor, Monomer& inAcceptor)
{
double result = 0;
if (inDonor.compoundID() != "PRO")
{
double distanceHO = Distance(inDonor.H(), inAcceptor.O());
double distanceHC = Distance(inDonor.H(), inAcceptor.C());
double distanceNC = Distance(inDonor.N(), inAcceptor.C());
double distanceNO = Distance(inDonor.N(), inAcceptor.O());
if (distanceHO < kMinimalDistance or distanceHC < kMinimalDistance or distanceNC < kMinimalDistance or distanceNO < kMinimalDistance)
result = kMinHBondEnergy;
else
result = kCouplingConstant / distanceHO - kCouplingConstant / distanceHC + kCouplingConstant / distanceNC - kCouplingConstant / distanceNO;
// DSSP compatibility mode:
result = round(result * 1000) / 1000;
if (result < kMinHBondEnergy)
result = kMinHBondEnergy;
}
// update donor
if (result < inDonor.mHBondAcceptor[0].energy)
{
inDonor.mHBondAcceptor[1] = inDonor.mHBondAcceptor[0];
inDonor.mHBondAcceptor[0].residue = &inAcceptor;
inDonor.mHBondAcceptor[0].energy = result;
}
else if (result < inDonor.mHBondAcceptor[1].energy)
{
inDonor.mHBondAcceptor[1].residue = &inAcceptor;
inDonor.mHBondAcceptor[1].energy = result;
}
// and acceptor
if (result < inAcceptor.mHBondDonor[0].energy)
{
inAcceptor.mHBondDonor[1] = inAcceptor.mHBondDonor[0];
inAcceptor.mHBondDonor[0].residue = &inDonor;
inAcceptor.mHBondDonor[0].energy = result;
}
else if (result < inAcceptor.mHBondDonor[1].energy)
{
inAcceptor.mHBondDonor[1].residue = &inDonor;
inAcceptor.mHBondDonor[1].energy = result;
}
return result;
}
struct Bridge
{
BridgeType type;
uint32 sheet, ladder;
set<Bridge*> link;
deque<uint32> i, j;
string chainI, chainJ;
bool operator<(const Bridge& b) const { return chainI < b.chainI or (chainI == b.chainI and i.front() < b.i.front()); }
};
void Structure::CalculateSecondaryStructure(bool inPreferPiHelices)
{
auto polies = polymers();
// CalculateHBondEnergies
for (auto& poly: polies)
{
for (size_t i = 0; i + 1 < poly.size(); ++i)
{
auto ri = poly[i];
for (size_t j = i + 1; j < poly.size(); ++j)
{
auto rj = poly[j];
if (Distance(ri.CAlpha().location(), rj.CAlpha().location()) < kMinimalCADistance)
{
Monomer::CalculateHBondEnergy(ri, rj);
if (j != i + 1)
Monomer::CalculateHBondEnergy(rj, ri);
}
}
}
}
// CalculateBetaSheets
vector<Bridge> bridges;
if (inResidues.size() > 4)
{
for (uint32 i = 1; i + 4 < inResidues.size(); ++i)
{
MResidue* ri = inResidues[i];
for (uint32 j = i + 3; j + 1 < inResidues.size(); ++j)
{
MResidue* rj = inResidues[j];
BridgeType type = ri->TestBridge(rj);
if (type == btNoBridge)
continue;
bool found = false;
for (Bridge& bridge : bridges)
{
if (type != bridge.type or i != bridge.i.back() + 1)
continue;
if (type == btParallel and bridge.j.back() + 1 == j)
{
bridge.i.push_back(i);
bridge.j.push_back(j);
found = true;
break;
}
if (type == btAntiParallel and bridge.j.front() - 1 == j)
{
bridge.i.push_back(i);
bridge.j.push_front(j);
found = true;
break;
}
}
if (not found)
{
Bridge bridge = {};
bridge.type = type;
bridge.i.push_back(i);
bridge.chainI = ri->GetChainID();
bridge.j.push_back(j);
bridge.chainJ = rj->GetChainID();
bridges.push_back(bridge);
}
}
}
}
// extend ladders
sort(bridges.begin(), bridges.end());
for (uint32 i = 0; i < bridges.size(); ++i)
{
for (uint32 j = i + 1; j < bridges.size(); ++j)
{
uint32 ibi = bridges[i].i.front();
uint32 iei = bridges[i].i.back();
uint32 jbi = bridges[i].j.front();
uint32 jei = bridges[i].j.back();
uint32 ibj = bridges[j].i.front();
uint32 iej = bridges[j].i.back();
uint32 jbj = bridges[j].j.front();
uint32 jej = bridges[j].j.back();
if (bridges[i].type != bridges[j].type or
MResidue::NoChainBreak(inResidues[min(ibi, ibj)], inResidues[max(iei, iej)]) == false or
MResidue::NoChainBreak(inResidues[min(jbi, jbj)], inResidues[max(jei, jej)]) == false or
ibj - iei >= 6 or
(iei >= ibj and ibi <= iej))
{
continue;
}
bool bulge;
if (bridges[i].type == btParallel)
bulge = ((jbj - jei < 6 and ibj - iei < 3) or (jbj - jei < 3));
else
bulge = ((jbi - jej < 6 and ibj - iei < 3) or (jbi - jej < 3));
if (bulge)
{
bridges[i].i.insert(bridges[i].i.end(), bridges[j].i.begin(), bridges[j].i.end());
if (bridges[i].type == btParallel)
bridges[i].j.insert(bridges[i].j.end(), bridges[j].j.begin(), bridges[j].j.end());
else
bridges[i].j.insert(bridges[i].j.begin(), bridges[j].j.begin(), bridges[j].j.end());
bridges.erase(bridges.begin() + j);
--j;
}
}
}
// Sheet
set<Bridge*> ladderset;
for (Bridge& bridge : bridges)
{
ladderset.insert(&bridge);
uint32 n = bridge.i.size();
if (n > kHistogramSize)
n = kHistogramSize;
if (bridge.type == btParallel)
mParallelBridgesPerLadderHistogram[n - 1] += 1;
else
mAntiparallelBridgesPerLadderHistogram[n - 1] += 1;
}
uint32 sheet = 1, ladder = 0;
while (not ladderset.empty())
{
set<Bridge*> sheetset;
sheetset.insert(*ladderset.begin());
ladderset.erase(ladderset.begin());
bool done = false;
while (not done)
{
done = true;
for (Bridge* a : sheetset)
{
for (Bridge* b : ladderset)
{
if (Linked(*a, *b))
{
sheetset.insert(b);
ladderset.erase(b);
done = false;
break;
}
}
if (not done)
break;
}
}
for (Bridge* bridge : sheetset)
{
bridge->ladder = ladder;
bridge->sheet = sheet;
bridge->link = sheetset;
++ladder;
}
uint32 nrOfLaddersPerSheet = sheetset.size();
if (nrOfLaddersPerSheet > kHistogramSize)
nrOfLaddersPerSheet = kHistogramSize;
if (nrOfLaddersPerSheet == 1 and (*sheetset.begin())->i.size() > 1)
mLaddersPerSheetHistogram[0] += 1;
else if (nrOfLaddersPerSheet > 1)
mLaddersPerSheetHistogram[nrOfLaddersPerSheet - 1] += 1;
++sheet;
}
for (Bridge& bridge : bridges)
{
// find out if any of the i and j set members already have
// a bridge assigned, if so, we're assigning bridge 2
uint32 betai = 0, betaj = 0;
for (uint32 l : bridge.i)
{
if (inResidues[l]->GetBetaPartner(0).residue != nullptr)
{
betai = 1;
break;
}
}
for (uint32 l : bridge.j)
{
if (inResidues[l]->GetBetaPartner(0).residue != nullptr)
{
betaj = 1;
break;
}
}
SecondaryStructure ss = betabridge;
if (bridge.i.size() > 1)
ss = strand;
if (bridge.type == btParallel)
{
mNrOfHBondsInParallelBridges += bridge.i.back() - bridge.i.front() + 2;
deque<uint32>::iterator j = bridge.j.begin();
for (uint32 i : bridge.i)
inResidues[i]->SetBetaPartner(betai, inResidues[*j++], bridge.ladder, true);
j = bridge.i.begin();
for (uint32 i : bridge.j)
inResidues[i]->SetBetaPartner(betaj, inResidues[*j++], bridge.ladder, true);
}
else
{
mNrOfHBondsInAntiparallelBridges += bridge.i.back() - bridge.i.front() + 2;
deque<uint32>::reverse_iterator j = bridge.j.rbegin();
for (uint32 i : bridge.i)
inResidues[i]->SetBetaPartner(betai, inResidues[*j++], bridge.ladder, false);
j = bridge.i.rbegin();
for (uint32 i : bridge.j)
inResidues[i]->SetBetaPartner(betaj, inResidues[*j++], bridge.ladder, false);
}
for (uint32 i = bridge.i.front(); i <= bridge.i.back(); ++i)
{
if (inResidues[i]->GetSecondaryStructure() != strand)
inResidues[i]->SetSecondaryStructure(ss);
inResidues[i]->SetSheet(bridge.sheet);
}
for (uint32 i = bridge.j.front(); i <= bridge.j.back(); ++i)
{
if (inResidues[i]->GetSecondaryStructure() != strand)
inResidues[i]->SetSecondaryStructure(ss);
inResidues[i]->SetSheet(bridge.sheet);
}
}
// CalculateAlphaHelices
}
Atom Structure::getAtomById(string id) const Atom Structure::getAtomById(string id) const
{ {
for (auto& a: mImpl->mAtoms) for (auto& a: mImpl->mAtoms)
......
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