Skip to content

L
libcifpp
Commit 99671072 by maarten
Options

water numbering scheme... 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@367 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent b5f159c3
Showing with 107 additions and 71 deletions
include/cif++/Statistics.h
......@@ -22,8 +22,7 @@ struct ResidueStatistics
std::string asymID;
int seqID;
std::string compID;
std::string pdbID; // comp_chain_seqnr''icode
std::string authSeqID;
double RSR, SRSR, RSCCS, EDIAm, OPIA;
int ngrid;
......
include/cif++/Structure.h
......@@ -108,7 +108,7 @@ class Atom
std::string authAtomId() const;
std::string authCompId() const;
std::string authAsymId() const;
int authSeqId() const;
std::string authSeqId() const;
std::string pdbxAuthInsCode() const;
std::string authAltId() const;
......@@ -169,7 +169,8 @@ class Residue
public:
// constructors should be private, but that's not possible for now (needed in emplace)
// constructor for waters
Residue(const Structure& structure, const std::string& asymID, Atom water, int ndbSeqID);
Residue(const Structure& structure, const std::string& compoundID,
const std::string& asymID, const std::string& authSeqID);
Residue(const Structure& structure, const std::string& compoundID,
const std::string& asymID, int seqID = 0);
......@@ -191,7 +192,7 @@ class Residue
const std::string& asymID() const { return mAsymID; }
int seqID() const { return mSeqID; }
int authSeqID() const;
std::string authSeqID() const;
std::string authInsCode() const;
// return a human readable PDB-like auth id (chain+seqnr+iCode)
......@@ -222,7 +223,8 @@ class Residue
const Structure* mStructure = nullptr;
std::string mCompoundID, mAsymID;
int mSeqID = 0, mNDBSeqID = 0;
int mSeqID = 0;
std::string mAuthSeqID;
AtomView mAtoms;
};
......@@ -377,7 +379,8 @@ class Structure
// returns chain,seqnr,comp,iCode
std::tuple<std::string,int,std::string,std::string> MapLabelToPDB(
const std::string& asymId, int seqId, const std::string& compId) const;
const std::string& asymId, int seqId, const std::string& compId,
const std::string& authSeqID) const;
std::tuple<std::string,int,std::string> MapPDBToLabel(
const std::string& asymId, int seqId, const std::string& compId, const std::string& iCode) const;
......
src/Cif++.cpp
......@@ -1162,6 +1162,16 @@ size_t Category::getColumnIndex(const string& name) const
break;
}
if (result == mColumns.size() and mCatValidator != nullptr) // validate the name, if it is known at all (since it was not found)
{
auto iv = mCatValidator->getValidatorForItem(name);
if (iv == nullptr)
{
assert(false);
throw logic_error("Invalid name used? " + name + " is not a known column in " + mName);
}
}
return result;
}
......
src/Statistics.cpp
......@@ -25,8 +25,12 @@ using clipper::Xmap;
ostream& operator<<(ostream& os, const ResidueStatistics& st)
{
os << st.asymID << '_' << st.seqID << '_' << st.compID << '\t'
<< st.RSR << '\t'
if (st.compID == "HOH")
os << st.asymID << '_' << st.authSeqID << '_' << st.compID << '\t';
else
os << st.asymID << '_' << st.seqID << '_' << st.compID << '\t';
os << st.RSR << '\t'
<< st.SRSR << '\t'
<< st.RSCCS << '\t'
<< st.ngrid << '\t'
......@@ -577,10 +581,13 @@ vector<ResidueStatistics> StatsCollector::collect() const
for (auto atom: mStructure.atoms())
{
auto k = make_tuple(atom.labelAsymId(), atom.labelSeqId(), atom.labelCompId(), atom.pdbID());
if (atom.isWater())
continue;
auto k = make_tuple(atom.labelAsymId(), atom.labelSeqId(), atom.labelCompId(), atom.authSeqId());
// auto k = make_tuple(atom.authAsymId(), atom.property<string>("auth_seq_id"), atom.authCompId());
if (not atom.isWater() and (residues.empty() or residues.back() != k))
if (residues.empty() or residues.back() != k)
{
residues.emplace_back(move(k));
atoms.emplace_back(move(atom));
......@@ -598,9 +605,14 @@ vector<ResidueStatistics> StatsCollector::collect(const string& asymID, int resF
for (auto atom: mStructure.atoms())
{
if (atom.isWater())
continue;
if (authNameSpace)
{
if (atom.authAsymId() != asymID or atom.authSeqId() < resFirst or atom.authSeqId() > resLast)
int authSeqID = stoi(atom.authSeqId());
if (atom.authAsymId() != asymID or authSeqID < resFirst or authSeqID > resLast)
continue;
}
else
......@@ -609,10 +621,10 @@ vector<ResidueStatistics> StatsCollector::collect(const string& asymID, int resF
continue;
}
auto k = make_tuple(atom.labelAsymId(), atom.labelSeqId(), atom.labelCompId(), atom.pdbID());
auto k = make_tuple(atom.labelAsymId(), atom.labelSeqId(), atom.labelCompId(), atom.authSeqId());
// auto k = make_tuple(atom.authAsymId(), atom.property<string>("auth_seq_id"), atom.authCompId());
if (not atom.isWater() and (residues.empty() or residues.back() != k))
if (residues.empty() or residues.back() != k)
{
residues.emplace_back(move(k));
atoms.emplace_back(move(atom));
......@@ -661,8 +673,8 @@ vector<ResidueStatistics> StatsCollector::collect(const vector<tuple<string,int,
for (auto r: residues)
{
int seqID;
string asymID, compID, pdbID;
tie(asymID, seqID, compID, pdbID) = r;
string asymID, compID, authSeqID;
tie(asymID, seqID, compID, authSeqID) = r;
AtomDataSums sums;
double ediaSum = 0;
......@@ -718,7 +730,7 @@ vector<ResidueStatistics> StatsCollector::collect(const vector<tuple<string,int,
}
result.emplace_back(ResidueStatistics{asymID, seqID, compID,
pdbID,
authSeqID,
(sums.rfSums[0] / sums.rfSums[1]), // rsr
sums.srg(), // srsr
sums.cc(), // rsccs
......@@ -737,7 +749,7 @@ vector<ResidueStatistics> StatsCollector::collect(const vector<tuple<string,int,
continue;
result.emplace_back(ResidueStatistics{d.asymID, d.seqID, "HOH",
atom.pdbID(),
atom.authSeqId(),
(d.sums.rfSums[0] / d.sums.rfSums[1]), // rsr
d.sums.srg(), // srsr
d.sums.cc(), // rsccs
......
src/Structure.cpp
......@@ -515,9 +515,9 @@ string Atom::authCompId() const
return property<string>("auth_comp_id");
}
int Atom::authSeqId() const
string Atom::authSeqId() const
{
return property<int>("auth_seq_id");
return property<string>("auth_seq_id");
}
string Atom::labelID() const
......@@ -588,11 +588,35 @@ float Atom::radius() const
// --------------------------------------------------------------------
// residue
Residue::Residue(const Structure& structure, const std::string& compoundID,
const std::string& asymID, int seqID)
Residue::Residue(const Structure& structure, const string& compoundID,
const string& asymID, const string& authSeqID)
: mStructure(&structure), mCompoundID(compoundID)
, mAsymID(asymID), mAuthSeqID(authSeqID)
{
assert(mCompoundID == "HOH");
for (auto& a: mStructure->atoms())
{
if (a.labelAsymId() != mAsymID or
a.labelCompId() != mCompoundID)
continue;
if (not mAuthSeqID.empty() and a.authSeqId() != mAuthSeqID) // water!
continue;
mAtoms.push_back(a);
}
assert(not mAtoms.empty());
}
Residue::Residue(const Structure& structure, const string& compoundID,
const string& asymID, int seqID)
: mStructure(&structure), mCompoundID(compoundID)
, mAsymID(asymID), mSeqID(seqID)
{
assert(mCompoundID != "HOH");
for (auto& a: mStructure->atoms())
{
if (mSeqID > 0 and a.labelSeqId() != mSeqID)
......@@ -606,15 +630,9 @@ Residue::Residue(const Structure& structure, const std::string& compoundID,
}
}
Residue::Residue(const Structure& structure, const std::string& asymID, Atom water, int ndbSeqID)
: mStructure(&structure), mCompoundID("HOH")
, mAsymID(asymID), mNDBSeqID(ndbSeqID), mAtoms({water})
{
}
Residue::Residue(Residue&& rhs)
: mStructure(rhs.mStructure), mCompoundID(move(rhs.mCompoundID)), mAsymID(move(rhs.mAsymID))
, mSeqID(rhs.mSeqID), mNDBSeqID(rhs.mNDBSeqID), mAtoms(move(rhs.mAtoms))
, mSeqID(rhs.mSeqID), mAuthSeqID(rhs.mAuthSeqID), mAtoms(move(rhs.mAtoms))
{
//cerr << "move constructor residue" << endl;
rhs.mStructure = nullptr;
......@@ -627,7 +645,7 @@ Residue& Residue::operator=(Residue&& rhs)
mCompoundID = move(rhs.mCompoundID);
mAsymID = move(rhs.mAsymID);
mSeqID = rhs.mSeqID;
mNDBSeqID = rhs.mNDBSeqID;
mAuthSeqID = rhs.mAuthSeqID;
mAtoms = move(rhs.mAtoms);
return *this;
......@@ -657,11 +675,11 @@ string Residue::authInsCode() const
return result;
}
int Residue::authSeqID() const
string Residue::authSeqID() const
{
assert(mStructure);
int result = 0;
string result;
try
{
......@@ -740,7 +758,7 @@ string Residue::authID() const
string Residue::labelID() const
{
if (mCompoundID == "HOH")
return mAsymID + to_string(mNDBSeqID);
return mAsymID + mAuthSeqID;
else
return mAsymID + to_string(mSeqID);
}
......@@ -1276,37 +1294,17 @@ void Structure::loadData()
auto& nonPolyScheme = category("pdbx_nonpoly_scheme");
vector<string> waterIDs;
for (auto& a: mAtoms)
{
if (a.isWater())
waterIDs.push_back(a.id());
}
auto waterIDIter = waterIDs.begin();
for (auto& r: nonPolyScheme)
{
string asymID, monID;
int ndbSeqNum;
cif::tie(asymID, monID, ndbSeqNum) =
r.get("asym_id", "mon_id", "ndb_seq_num");
string asymID, monID, pdbSeqNum;
cif::tie(asymID, monID, pdbSeqNum) =
r.get("asym_id", "mon_id", "pdb_seq_num");
if (monID == "HOH")
{
if (waterIDIter == waterIDs.end())
throw runtime_error("pdbx_nonpoly_scheme contains more waters than water atoms are available in atom_site");
Atom water = getAtomById(*waterIDIter);
++waterIDIter;
mNonPolymers.emplace_back(*this, asymID, water, ndbSeqNum);
}
mNonPolymers.emplace_back(*this, monID, asymID, pdbSeqNum);
else if (mNonPolymers.empty() or mNonPolymers.back().asymID() != asymID)
mNonPolymers.emplace_back(*this, monID, asymID);
}
if (waterIDIter != waterIDs.end())
throw runtime_error("there are more waters in atom_site than referenced in pdbx_nonpoly_scheme");
}
void Structure::updateAtomIndex()
......@@ -1421,28 +1419,42 @@ tuple<char,int,char> Structure::MapLabelToAuth(
}
tuple<string,int,string,string> Structure::MapLabelToPDB(
const string& asymId, int seqId, const string& monId) const
const string& asymId, int seqId, const string& monId,
const string& authSeqID) const
{
auto& db = datablock();
tuple<string,int,string,string> result;
for (auto r: db["pdbx_poly_seq_scheme"].find(
cif::Key("asym_id") == asymId and
cif::Key("seq_id") == seqId and
cif::Key("mon_id") == monId))
if (monId == "HOH")
{
result = r.get("pdb_strand_id", "pdb_seq_num", "pdb_mon_id", "pdb_ins_code");
break;
for (auto r: db["pdbx_nonpoly_scheme"].find(
cif::Key("asym_id") == asymId and
cif::Key("pdb_seq_num") == authSeqID and
cif::Key("mon_id") == monId))
{
result = r.get("pdb_strand_id", "pdb_seq_num", "pdb_mon_id", "pdb_ins_code");
break;
}
}
for (auto r: db["pdbx_nonpoly_scheme"].find(
cif::Key("asym_id") == asymId and
cif::Key("seq_id") == seqId and
cif::Key("mon_id") == monId))
else
{
result = r.get("pdb_strand_id", "pdb_seq_num", "pdb_mon_id", "pdb_ins_code");
break;
for (auto r: db["pdbx_poly_seq_scheme"].find(
cif::Key("asym_id") == asymId and
cif::Key("seq_id") == seqId and
cif::Key("mon_id") == monId))
{
result = r.get("pdb_strand_id", "pdb_seq_num", "pdb_mon_id", "pdb_ins_code");
break;
}
for (auto r: db["pdbx_nonpoly_scheme"].find(
cif::Key("asym_id") == asymId and
cif::Key("mon_id") == monId))
{
result = r.get("pdb_strand_id", "pdb_seq_num", "pdb_mon_id", "pdb_ins_code");
break;
}
}
return result;
......
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