Skip to content

to_mmcif, to_pdb: keep link atoms present only as altlocs #443

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

Open
CV-GPhL wants to merge 1 commit into
base: master
Choose a base branch
from
Open

to_mmcif, to_pdb: keep link atoms present only as altlocs #443

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions include/gemmi/to_pdb.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

#include "model.hpp"
#include <ostream>
#include <utility> // for pair

namespace gemmi {

Expand Down Expand Up @@ -86,6 +87,21 @@ struct PdbWriteOptions {
/// @return true if HETATM should be used, false for ATOM.
GEMMI_DLL bool use_hetatm(const Residue& res);

/// @brief Find the atoms of a Connection in the first model, for writing.
///
/// Like Model::find_cra() for both partners, except for a partner without
/// altloc whose atom exists only in alternative conformations (not found by
/// the strict altloc matching): then the closest pair of atoms of the same
/// conformer is returned, so that the atom names, symmetry and distance can
/// still be written. The altloc of such an atom is not part of the Connection.
/// Atoms that are genuinely absent stay null.
///
/// @param st The Structure (must have at least one model).
/// @param con The Connection.
/// @return The CRAs of partner1 and partner2.
GEMMI_DLL std::pair<const_CRA, const_CRA> find_connection_atoms(const Structure& st,
const Connection& con);

/// @brief Write a Structure to PDB format on an output stream.
/// @param st The Structure to serialize.
/// @param os The output stream to write to.
Expand Down
12 changes: 8 additions & 4 deletions src/to_mmcif.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -410,8 +410,8 @@ void write_struct_conn(const Structure& st, cif::Block& block) {
if (use_ccp4_link_id)
conn_loop.tags.push_back("_struct_conn.ccp4_link_id");
for (const Connection& con : st.connections) {
const_CRA cra1 = st.models[0].find_cra(con.partner1, true);
const_CRA cra2 = st.models[0].find_cra(con.partner2, true);
const auto cras = find_connection_atoms(st, con);
const const_CRA& cra1 = cras.first, &cra2 = cras.second;
if (!cra1.residue || !cra2.residue)
continue;
const Atom* at1 = cra1.atom;
Expand All @@ -422,14 +422,18 @@ void write_struct_conn(const Structure& st, cif::Block& block) {
im_pdb_symbol = im.symmetry_code(true);
im_dist_str = to_str_prec<4>(im.dist());
}
// altloc as given in the Connection; an atom found for a partner
// without altloc (see find_connection_atoms) does not add one
char alt1 = at1 && con.partner1.altloc ? at1->altloc_or('?') : '?';
char alt2 = at2 && con.partner2.altloc ? at2->altloc_or('?') : '?';
auto& v = conn_loop.values;
v.emplace_back(string_or_qmark(con.name)); // id
v.emplace_back(connection_type_to_string(con.type)); // conn_type_id
v.emplace_back(subchain_or_dot(*cra1.residue)); // ptnr1_label_asym_id
v.emplace_back(cra1.residue->name); // ptnr1_label_comp_id
v.emplace_back(cra1.residue->label_seq.str('.')); // ptnr1_label_seq_id
v.emplace_back(at1 ? cif::quote(at1->name) : "?"); // ptnr1_label_atom_id
v.emplace_back(1, at1 ? at1->altloc_or('?') : '?'); // pdbx_ptnr1_label_alt_id
v.emplace_back(1, alt1); // pdbx_ptnr1_label_alt_id
v.emplace_back(qchain(con.partner1.chain_name)); // ptnr1_auth_asym_id
v.emplace_back(cra1.residue->seqid.num.str()); // ptnr1_auth_seq_id
v.emplace_back(pdbx_icode(con.partner1.res_id)); // ptnr1_PDB_ins_code
Expand All @@ -438,7 +442,7 @@ void write_struct_conn(const Structure& st, cif::Block& block) {
v.emplace_back(cra2.residue->name); // ptnr2_label_comp_id
v.emplace_back(cra2.residue->label_seq.str('.')); // ptnr2_label_seq_id
v.emplace_back(at2 ? cif::quote(at2->name) : "?"); // ptnr2_label_atom_id
v.emplace_back(1, at2 ? at2->altloc_or('?') : '?'); // pdbx_ptnr2_label_alt_id
v.emplace_back(1, alt2); // pdbx_ptnr2_label_alt_id
v.emplace_back(qchain(con.partner2.chain_name)); // ptnr2_auth_asym_id
v.emplace_back(cra2.residue->seqid.num.str()); // ptnr2_auth_seq_id
v.emplace_back(pdbx_icode(con.partner2.res_id)); // ptnr2_PDB_ins_code
Expand Down
52 changes: 46 additions & 6 deletions src/to_pdb.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@

#include <cassert>
#include <cctype> // for isdigit
#include <cmath> // for INFINITY
#include <cstring> // for memset, memcpy
#include <algorithm>
#include <array>
Expand Down Expand Up @@ -50,6 +51,43 @@ bool use_hetatm(const Residue& res) {
return !find_tabulated_residue(res.name).is_standard();
}

std::pair<const_CRA, const_CRA> find_connection_atoms(const Structure& st,
const Connection& con) {
const Model& model = st.models.at(0);
std::pair<const_CRA, const_CRA> cras{model.find_cra(con.partner1, true),
model.find_cra(con.partner2, true)};
const_CRA& cra1 = cras.first;
const_CRA& cra2 = cras.second;
if (!cra1.residue || !cra2.residue || (cra1.atom && cra2.atom))
return cras;
// A partner without altloc can name an atom that is present only as A, B...
// Pick the closest pair of the same conformer, as complete_ssbond() does.
auto candidates = [](const const_CRA& cra, const AtomAddress& ad) {
std::vector<const Atom*> atoms;
if (cra.atom)
atoms.push_back(cra.atom);
else if (ad.altloc == '\0' && !ad.atom_name.empty())
for (const Atom& a : cra.residue->atoms)
if (a.name == ad.atom_name)
atoms.push_back(&a);
return atoms;
};
const std::vector<const Atom*> atoms1 = candidates(cra1, con.partner1);
const std::vector<const Atom*> atoms2 = candidates(cra2, con.partner2);
double min_dist_sq = INFINITY;
for (const Atom* a1 : atoms1)
for (const Atom* a2 : atoms2)
if (a1->same_conformer(*a2)) {
double dist_sq = st.cell.find_nearest_image(a1->pos, a2->pos, con.asu).dist_sq;
if (dist_sq < min_dist_sq) {
min_dist_sq = dist_sq;
cra1.atom = a1;
cra2.atom = a2;
}
}
return cras;
}

namespace {

// works for non-negative values only
Expand Down Expand Up @@ -622,8 +660,8 @@ void write_pdb(const Structure& st, std::ostream& os, PdbWriteOptions opt) {
int counter = 0;
for (const Connection& con : st.connections)
if (con.type == Connection::Disulf) {
const_CRA cra1 = st.models[0].find_cra(con.partner1, true);
const_CRA cra2 = st.models[0].find_cra(con.partner2, true);
const auto cras = find_connection_atoms(st, con);
const const_CRA& cra1 = cras.first, &cra2 = cras.second;
if (!cra1.atom || !cra2.atom)
continue;
NearestImage im = st.cell.find_nearest_image(cra1.atom->pos,
Expand All @@ -645,8 +683,8 @@ void write_pdb(const Structure& st, std::ostream& os, PdbWriteOptions opt) {
for (const Connection& con : st.connections)
if (con.type == Connection::Covale || con.type == Connection::MetalC ||
con.type == Connection::Unknown) {
const_CRA cra1 = st.models[0].find_cra(con.partner1, true);
const_CRA cra2 = st.models[0].find_cra(con.partner2, true);
const auto cras = find_connection_atoms(st, con);
const const_CRA& cra1 = cras.first, &cra2 = cras.second;
// In special cases (LINKR gap) atoms are not there.
if (!cra1.residue || !cra2.residue)
continue;
Expand All @@ -667,12 +705,14 @@ void write_pdb(const Structure& st, std::ostream& os, PdbWriteOptions opt) {
snprintf_z(buf, 82, "LINK %-4s%c%3.3s%2s%5s "
" %-4s%c%3.3s%2s%5s %6s %6s %5s \n",
cra1.atom ? cra1.atom->padded_name().c_str() : "",
cra1.atom && cra1.atom->altloc ? std::toupper(cra1.atom->altloc) : ' ',
cra1.atom && cra1.atom->altloc && con.partner1.altloc
? std::toupper(cra1.atom->altloc) : ' ',
cra1.residue->name.c_str(),
con.partner1.chain_name.c_str(),
write_seq_id(cra1.residue->seqid).data(),
cra2.atom ? cra2.atom->padded_name().c_str() : "",
cra2.atom && cra2.atom->altloc ? std::toupper(cra2.atom->altloc) : ' ',
cra2.atom && cra2.atom->altloc && con.partner2.altloc
? std::toupper(cra2.atom->altloc) : ' ',
cra2.residue->name.c_str(),
con.partner2.chain_name.c_str(),
write_seq_id(cra2.residue->seqid).data(),
Expand Down
43 changes: 43 additions & 0 deletions tests/test_mol.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,23 @@ def is_written_to_pdb(line, via_cif):
ATOM 69 SG CYS A 11 36.989 25.994 19.570 1.00 14.23 16 S
"""

# links without altloc to atoms that are present only in conformers A and B
LINK_TO_ALTLOCS = """\
LINK ZN ZN A 301 NE2 HIS A 64 1555 1555 2.10
CRYST1 30.000 30.000 30.000 90.00 90.00 90.00 P 1
ATOM 1 NE2AHIS A 64 12.100 10.000 10.000 0.60 20.00 N
ATOM 2 NE2BHIS A 64 12.600 10.000 10.000 0.40 20.00 N
HETATM 3 ZN ZN A 301 10.000 10.000 10.000 1.00 20.00 ZN
"""

SSBOND_TO_ALTLOCS = """\
SSBOND 1 CYS A 5 CYS A 20 1555 1555 2.04
CRYST1 30.000 30.000 30.000 90.00 90.00 90.00 P 1
ATOM 1 SG ACYS A 5 10.000 10.000 10.000 0.70 20.00 S
ATOM 2 SG BCYS A 5 10.000 11.600 11.600 0.30 20.00 S
ATOM 3 SG CYS A 20 12.040 10.000 10.000 1.00 20.00 S
"""

# fragment of $CCP4/examples/toxd/toxd_mod_p1.pdb without chain names
BLANK_CHAIN_FRAGMENT = """\
CRYST1 50.000 50.00 50.00 90.00 90.00 90.00
Expand Down Expand Up @@ -780,6 +797,32 @@ def test_short_ssbond(self):
SHORT_SSBOND.splitlines()[0]
+ " 1555 1555 2.06 ")

def test_link_to_altloc_atoms(self):
# LINK has no altlocs, NE2 is present only as A and B
st = gemmi.read_pdb_string(LINK_TO_ALTLOCS)
self.assertEqual(st.connections[0].partner2.altloc, '\0')
out = [line.rstrip() for line in st.make_pdb_headers().splitlines()]
self.assertIn(LINK_TO_ALTLOCS.splitlines()[0], out)
doc = st.make_mmcif_document()
tags = ['ptnr1_label_atom_id', 'pdbx_ptnr1_label_alt_id',
'ptnr2_label_atom_id', 'pdbx_ptnr2_label_alt_id',
'ptnr2_symmetry', 'pdbx_dist_value']
row = doc[0].find('_struct_conn.', tags)[0]
# names written, no altloc added, distance from the closest conformer
self.assertEqual(list(row), ['ZN', '?', 'NE2', '?', '1_555', '2.1000'])

def test_ssbond_to_altloc_atoms(self):
st = gemmi.read_pdb_string(SSBOND_TO_ALTLOCS)
# as read from mmCIF with pdbx_ptnr1_label_alt_id '?'
st.connections[0].partner1.altloc = '\0'
out = st.make_pdb_headers().splitlines()
self.assertEqual(out[0].rstrip(), SSBOND_TO_ALTLOCS.splitlines()[0])
doc = st.make_mmcif_document()
row = doc[0].find('_struct_conn.', ['ptnr1_label_atom_id',
'pdbx_ptnr1_label_alt_id',
'pdbx_dist_value'])[0]
self.assertEqual(list(row), ['SG', '?', '2.0400'])

def test_add_remove(self):
st = gemmi.read_pdb_string(SSBOND_FRAGMENT)
st.add_model(st[0])
Expand Down
Loading