Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve PDB formatting with incomplete Monomer info #7286

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Improve PDB formatting with incomplete Monomer info #7286

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
257b3c0
Fix PDB formatting if residue name is set but atom name is not.
Mar 19, 2024
ad2ce4c
Merge branch 'master' of github.com:fwaibl/rdkit
Mar 21, 2024
03d1414
Add a test for missing residue name and short atom names
Mar 22, 2024
ab2fcb9
Fix atom name alignment in unit test
Mar 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
75 changes: 45 additions & 30 deletions Code/GraphMol/FileParsers/PDBWriter.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,13 +38,17 @@
// flavor & 32 : Write TER record

namespace RDKit {

std::array<char, 2> GetDefaultAtomNumber(const Atom* atom,
std::map<unsigned int, unsigned int> &elem);

std::string GetPDBAtomLine(const Atom *atom, const Conformer *conf,
std::map<unsigned int, unsigned int> &elem) {
PRECONDITION(atom, "bad atom");
std::stringstream ss;

std::string symb = atom->getSymbol();
char at1, at2, at3, at4;
char at1, at2;
switch (symb.length()) {
case 0:
at1 = ' ';
Expand All @@ -68,13 +72,19 @@ std::string GetPDBAtomLine(const Atom *atom, const Conformer *conf,
ss << (info->getIsHeteroAtom() ? "HETATM" : "ATOM ");
ss << std::setw(5) << atom->getIdx() + 1;
ss << ' ';
ss << info->getName(); // Always 4 characters?
const std::string& name = info->getName();
if (name.size() == 0) {
std::array<char, 2> atnum = GetDefaultAtomNumber(atom, elem);
ss << at1 << at2 << atnum[0] << atnum[1];
} else {
ss << std::setw(4) << name;
}
const char *ptr = info->getAltLoc().c_str();
if (*ptr == '\0') {
ptr = " ";
}
ss << *ptr;
ss << info->getResidueName(); // Always 3 characters?
ss << std::setw(3) << info->getResidueName();
ss << ' ';
ptr = info->getChainId().c_str();
if (*ptr == '\0') {
Expand All @@ -90,36 +100,11 @@ std::string GetPDBAtomLine(const Atom *atom, const Conformer *conf,
ss << " ";
} else {
info = (AtomPDBResidueInfo *)nullptr;
unsigned int atno = atom->getAtomicNum();
if (elem.find(atno) == elem.end()) {
elem[atno] = 1;
at3 = '1';
at4 = ' ';
} else {
unsigned int tmp = elem[atno] + 1;
elem[atno] = tmp;
if (tmp < 10) {
at3 = tmp + '0';
at4 = ' ';
} else if (tmp < 100) {
at3 = (tmp / 10) + '0';
at4 = (tmp % 10) + '0';
} else if (tmp < 360) {
at3 = ((tmp - 100) / 10) + 'A';
at4 = ((tmp - 100) % 10) + '0';
} else if (tmp < 1036) {
at3 = ((tmp - 360) / 26) + 'A';
at4 = ((tmp - 360) % 26) + 'A';
} else {
at3 = ' ';
at4 = ' ';
}
}

std::array<char, 2> atnum = GetDefaultAtomNumber(atom, elem);
ss << "HETATM";
ss << std::setw(5) << atom->getIdx() + 1;
ss << ' ';
ss << at1 << at2 << at3 << at4;
ss << at1 << at2 << atnum[0] << atnum[1];
ss << " UNL 1 ";
}

Expand Down Expand Up @@ -153,6 +138,36 @@ std::string GetPDBAtomLine(const Atom *atom, const Conformer *conf,
return ss.str();
}

std::array<char, 2> GetDefaultAtomNumber(const Atom* atom, std::map<unsigned int, unsigned int> &elem) {
std::array<char, 2> ret;
unsigned int atno = atom->getAtomicNum();
if (elem.find(atno) == elem.end()) {
elem[atno] = 1;
ret[0] = '1';
ret[1] = ' ';
} else {
unsigned int tmp = elem[atno] + 1;
elem[atno] = tmp;
if (tmp < 10) {
ret[0] = tmp + '0';
ret[1] = ' ';
} else if (tmp < 100) {
ret[0] = (tmp / 10) + '0';
ret[1] = (tmp % 10) + '0';
} else if (tmp < 360) {
ret[0] = ((tmp - 100) / 10) + 'A';
ret[1] = ((tmp - 100) % 10) + '0';
} else if (tmp < 1036) {
ret[0] = ((tmp - 360) / 26) + 'A';
ret[1] = ((tmp - 360) % 26) + 'A';
} else {
ret[0] = ' ';
ret[1] = ' ';
}
}
return ret;
}

std::string GetPDBBondLines(const Atom *atom, bool all, bool both, bool mult,
unsigned int &conect_count) {
PRECONDITION(atom, "bad atom");
Expand Down
42 changes: 42 additions & 0 deletions Code/GraphMol/FileParsers/file_parsers_catch.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
#include <GraphMol/FileParsers/PNGParser.h>
#include <GraphMol/FileParsers/MolFileStereochem.h>
#include <GraphMol/FileParsers/MolWriters.h>
#include <GraphMol/MonomerInfo.h>
#include <RDGeneral/FileParseException.h>
#include <boost/algorithm/string.hpp>

Expand Down Expand Up @@ -3019,6 +3020,47 @@ TEST_CASE("test bond flavors when writing PDBs", "[bug]") {
}
}

TEST_CASE("test output with incomplete monomer info", "[bug][writer]") {
SECTION("basics") {
{
std::unique_ptr<RWMol> m{SmilesToMol("Cl")};
REQUIRE(m);
std::string pdb = MolToPDBBlock(*m, -1);
CHECK(pdb.find("HETATM 1 CL1 UNL 1") != std::string::npos);
}
{
std::unique_ptr<RWMol> m{SmilesToMol("Cl")};
// will get deleted by ~Atom()
AtomPDBResidueInfo *info = new AtomPDBResidueInfo();
info->setResidueName("HCL");
m->getAtomWithIdx(0)->setMonomerInfo(info);
std::string pdb = MolToPDBBlock(*m, -1);
CHECK(pdb.find("ATOM 1 CL1 HCL 0") != std::string::npos);
}
{
std::unique_ptr<RWMol> m{SmilesToMol("Cl")};
// will get deleted by ~Atom()
AtomPDBResidueInfo *info = new AtomPDBResidueInfo();
info->setResidueName("HCL");
info->setName("Cl1");
m->getAtomWithIdx(0)->setMonomerInfo(info);
std::string pdb = MolToPDBBlock(*m, -1);
CHECK(pdb.find("ATOM 1 Cl1 HCL 0") != std::string::npos);
}
{
// 1. should add spaces for padding if the atom name is too short
// 2. should add spaces for missing residue name.
std::unique_ptr<RWMol> m{SmilesToMol("Cl")};
// will get deleted by ~Atom()
AtomPDBResidueInfo *info = new AtomPDBResidueInfo();
info->setName("CL");
m->getAtomWithIdx(0)->setMonomerInfo(info);
std::string pdb = MolToPDBBlock(*m, -1);
CHECK(pdb.find("ATOM 1 CL 0") != std::string::npos);
}
}
}

TEST_CASE(
"github #3599: Add explicit support for remaining CTAB query bond types",
"[feature]") {
Expand Down