From 4867eda022846b6c0f773f307cf5ef9a1ad66a5c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 03:00:55 +0100 Subject: [PATCH] RESULT-CHANGING: RefinementNcs: copies of a molecule must also lie within 2.5 A rms of each other A chain is a copy of a molecule only where, besides the sequence rule, its C-alpha (P) atoms lie within 2.5 A rms of the molecule's first copy after superposition (REFMAC5's second criterion): the same sequence in a conformation further away is another molecule. 4NWV (a T=3 capsid expanded from 14 strict-NCS operators) has 45 protein chains of ONE sequence: A, B and C are its quasi-equivalent subunits (A starts 14 residues earlier; A is 3.3-3.6 A from B and C, B and C 1.2 A apart). All 45 were one molecule, paired in model order A-B-C, copies 2.9 A apart, and NCS was refused. Now A x 15 and B/C x 30: NCS kept (sigma 0.1 A), --mode refine R-free 0.2565 -> 0.2450. Every other NCS set (open arm and the MD starts) has all its copies within 2.5 A, so its groups and results are unchanged. Ordering each molecule's copies as a chain of nearest neighbours (strict copies beside each other) was tried too: 4NWV 0.2401, but one MD start of a four-copy set lost 0.039 and a twelve-copy set 0.001 on three starts. Not kept. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/ACKNOWLEDGEMENT.md | 3 +- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 12 ++++---- .../structure_refinement/RefinementNcs.cpp | 29 +++++++++++++++---- .../structure_refinement/RefinementNcs.h | 16 +++++----- tests/RefinementTest.cpp | 15 ++++++++-- 5 files changed, 55 insertions(+), 20 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 96c2b3c59..a7f8749da 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -631,7 +631,8 @@ what the data support rather than a rule by resolution: E. A. Merritt, "To B or resolution?" (2012), Acta Cryst. D68, 468-477 [doi:10.1107/S0907444911028320](https://doi.org/10.1107/S0907444911028320). **NCS restraints** — copies of a molecule in the model are found by sequence alignment with REFMAC5's criterion -(more than 15 residues aligned, 80% identity - here of the longer chain; Murshudov et al. 2011, above) and restrained locally, on pairs of +(more than 15 residues aligned, 80% identity - here of the longer chain - and 2.5 Å rms after superposition; +Murshudov et al. 2011, above) and restrained locally, on pairs of corresponding interatomic distances under the Geman-McClure robust loss, as REFMAC5's local NCS restraints are, with the distance range of R. A. Nicholls, F. Long and G. N. Murshudov, "Low-resolution refinement tools in REFMAC5" (2012), Acta Cryst. D68, 404-417 [doi:10.1107/S090744491105606X](https://doi.org/10.1107/S090744491105606X), and diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index 0e93175e8..01274b2b0 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -236,11 +236,13 @@ the ceiling of §15.1. **NCS restraints.** Where the model holds several copies of one molecule, they are found from the model alone: two polymer chains of one kind (protein, or nucleic acid) are copies where their sequences align with more than 15 -residues identical and identical in 80% of the longer chain or more - REFMAC5's criterion, taken on the longer chain. In the residues aligned to a residue of the same name, the -atoms of the same name correspond - not a hydrogen, not an alternative conformation, both refinable. Each copy is -paired with the next of its molecule (and the last with the first, where there are three or more), and each pair is -restrained **locally**, on its interatomic distances rather than on a superposition: every distance within a chain -shorter than 4.2 Å in either copy, except those a bond or an angle of the dictionary already fixes, against the same +residues identical and identical in 80% of the longer chain or more, and their C-alpha (P) atoms lie within 2.5 Å +rms of each other after superposition - REFMAC5's criterion, the identity taken on the longer chain. The same +sequence in a conformation further away than that, such as a quasi-equivalent subunit of a virus capsid, is another +molecule. In the residues aligned to a residue of the same name, the atoms of the same name correspond - not a +hydrogen, not an alternative conformation, both refinable. Each copy is paired with the next of its molecule (and the +last with the first, where there are three or more), and each pair is restrained **locally**, on its interatomic +distances rather than on a superposition: every distance within a chain shorter than 4.2 Å in either copy, except those a bond or an angle of the dictionary already fixes, against the same distance in the other copy, with the robust loss of Geman and McClure, $$ diff --git a/image_analysis/structure_refinement/RefinementNcs.cpp b/image_analysis/structure_refinement/RefinementNcs.cpp index 7590a7c08..51c566ad9 100644 --- a/image_analysis/structure_refinement/RefinementNcs.cpp +++ b/image_analysis/structure_refinement/RefinementNcs.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -26,6 +27,10 @@ constexpr float NCS_RELEASE = 0.5f; // A // Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 constexpr int NCS_MIN_ALIGNED = 15; constexpr double NCS_MIN_IDENTITY = 80.0; // % +// and their C-alpha (P) atoms no further apart than this after superposition: the same sequence further away (a +// quasi-equivalent subunit of a virus capsid) is another molecule. +// Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 +constexpr double NCS_MAX_RMSD = 2.5; // A // One polymer chain: its residues and, per residue, its atoms' indices by name. struct NcsChain { @@ -65,15 +70,30 @@ void Correspond(const NcsChain &a, const NcsChain &b, std::vector &pos) { + const double none = std::numeric_limits::infinity(); if (a.type != b.type || a.sequence.size() <= NCS_MIN_ALIGNED || b.sequence.size() <= NCS_MIN_ALIGNED) - return false; + return none; const gemmi::AlignmentResult al = gemmi::align_string_sequences(a.sequence, b.sequence, {}, gemmi::AlignmentScoring::simple()); // identity over the longer chain: over the shorter one, a short RNA aligned with gaps into a long one finds most of // its four letters matched const size_t longer = std::max(a.sequence.size(), b.sequence.size()); - return al.match_count > NCS_MIN_ALIGNED && 100.0 * al.match_count >= NCS_MIN_IDENTITY * double(longer); + if (al.match_count <= NCS_MIN_ALIGNED || 100.0 * al.match_count < NCS_MIN_IDENTITY * double(longer)) + return none; + std::vector> atoms, anchors; + std::vector residue_of; + Correspond(a, b, atoms, residue_of, anchors); + if (anchors.size() < 3) + return none; + std::vector pa, pb; + for (const std::array &p : anchors) { + pa.push_back(pos[p[0]]); + pb.push_back(pos[p[1]]); + } + return gemmi::superpose_positions(pa.data(), pb.data(), pa.size(), nullptr).rmsd; } // The pairs (i, j), i < j, of `atoms` closer than `reach` at `pos`, by a grid of cells of that size. @@ -146,7 +166,7 @@ NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &top for (size_t k = 0; k < chains.size(); k++) { bool placed = false; for (std::vector &g : groups) - if (AreCopies(chains[g[0]], chains[k])) { + if (CopyRmsd(chains[g[0]], chains[k], pos) <= NCS_MAX_RMSD) { g.push_back(static_cast(k)); placed = true; break; @@ -154,7 +174,6 @@ NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &top if (!placed) groups.push_back({static_cast(k)}); } - // The distances the dictionary fixes: bonded atoms and those bonded to a common atom. std::set> fixed; for (const RestraintRecord &r : topo.geometry.records) diff --git a/image_analysis/structure_refinement/RefinementNcs.h b/image_analysis/structure_refinement/RefinementNcs.h index 0d8592ca4..322592cfc 100644 --- a/image_analysis/structure_refinement/RefinementNcs.h +++ b/image_analysis/structure_refinement/RefinementNcs.h @@ -30,13 +30,15 @@ struct NcsCopies { }; // The copies of `model` (whose entities are set up, as BuildRefinementTopology leaves them): polymer chains whose -// sequences align with more than 15 residues identical, and identical in 80% of the longer chain's residues, are -// copies of one molecule (REFMAC's criterion, taken on the longer chain). Their residues of the same name aligned to -// each other correspond, and in them the atoms of the same name - neither a hydrogen nor an alternative conformation, -// all refinable (topo.refinable). Each copy is paired with the next of its molecule (the last with the first, where -// there are three or more), and each pair restrained on the distances within a chain shorter than NCS_DISTANCE_MAX in -// either copy, those the dictionary's bonds and angles fix excepted. Atoms are numbered in model order, as -// ModelPositions() lists them. +// sequences align with more than 15 residues identical, identical in 80% of the longer chain's residues, and whose +// C-alpha (P) atoms are within 2.5 A rms of each other after superposition are copies of one molecule (REFMAC's +// criterion, the identity taken on the longer chain); the same sequence in a conformation further away - a +// quasi-equivalent subunit of a virus capsid - is another molecule. Their residues of the same name aligned to each +// other correspond, and in them the atoms of the same name - neither a hydrogen nor an alternative conformation, all +// refinable (topo.refinable). Each copy is paired with the next of its molecule (the last with the first, where there +// are three or more), and each pair restrained on the distances within a chain shorter than NCS_DISTANCE_MAX in either +// copy, those the dictionary's bonds and angles fix excepted. Atoms are numbered in model order, as ModelPositions() +// lists them. NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &topo); // The restraints at `sigma` (A). diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp index 6b7932280..5ddd593cd 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -973,8 +973,8 @@ TEST_CASE("Refinement_ion_form_factors", "[Refinement]") { } } -// Two copies of one RNA chain in the asymmetric unit, the second turned and moved, and a fragment of it: the two found -// as copies of one molecule, the fragment not, held +// Two copies of one RNA chain in the asymmetric unit, the second turned and moved, a fragment of it and the chain in +// another conformation: the two found as copies of one molecule, the fragment and the other conformation not, held // to each other by local distance restraints whose gradient is right on the CPU and the GPU, and let go where the // copies really differ. A single copy has none. TEST_CASE("RefinementNcs_two_copies", "[Refinement]") { @@ -1005,6 +1005,17 @@ TEST_CASE("RefinementNcs_two_copies", "[Refinement]") { for (gemmi::Atom &at : r.atoms) at.pos += gemmi::Position(-30, 40, 20); st.models[0].chains.push_back(fragment); + // and the whole molecule in another conformation, each residue moved 3 A its own way: the same sequence, not a copy + gemmi::Chain other = st.models[0].chains.at(0); + other.name = "D"; + std::mt19937 shift_rng(11); + std::uniform_real_distribution shift(-1.0, 1.0); + for (gemmi::Residue &r : other.residues) { + const gemmi::Vec3 v = gemmi::Vec3(shift(shift_rng), shift(shift_rng), shift(shift_rng)).normalized() * 3.0; + for (gemmi::Atom &at : r.atoms) + at.pos += gemmi::Position(v) + gemmi::Position(30, -40, 25); + } + st.models[0].chains.push_back(other); st.connections.clear(); // the LINK records name chain A only const RefinementTopology topo = BuildRefinementTopology(st); const NcsCopies ncs = FindNcsCopies(st.models[0], topo);