Repository navigation
Write CONECT records for bonds between residues differing in insertion code - #946
Conversation
padix-key
left a comment
padix-key
left a comment
There was a problem hiding this comment.
Good catch! Thank you for the fix 👍
| np.testing.assert_array_equal(actual_bonds, expected_bonds) | ||
|
|
||
|
|
||
| def test_inter_residue_bond_with_ins_code(): |
There was a problem hiding this comment.
Would it be possible to generalize this test? For example we could check that PDBFile.set_structure() adds the same CONECT records as there are in the PDB file downloaded from RCSB. The test assets do contain structures with insertion codes, so this should catch it. If for some reason complete CONECT equivalence is not possible, we could at least compare the number of CONECT records.
Merging this PR will degrade performance by 16.1%2 warnings
Performance Changes
Comparing 
Footnotes
|
|
Good idea. I replaced the synthetic test with one over every PDB test asset. It checks that the written One thing I found while doing this: the RCSB comparison alone doesn't catch the bug. RCSB only writes The failing "databases" job is the |
| test_file.set_structure(atoms) | ||
| test_pairs = _conect_pairs(test_file) | ||
|
|
||
| assert _conect_pairs(ref_file) <= test_pairs |
There was a problem hiding this comment.
Shouldn't these two sets be equal?
| assert _conect_pairs(ref_file) <= test_pairs | |
| assert _conect_pairs(ref_file) = test_pairs |
There was a problem hiding this comment.
Sorry for the slow reply on this one.
I tried it with == and it fails on all 19 PDB test assets, but only one way. None of the original CONECT pairs are missing from the written file. The extra ones are all peptide C-N or O3'-P bonds (7279 pairs across the 19 files, e.g. LYS1.C - VAL2.N in 1aki). The original files have no records for backbone links, while set_structure() writes CONECT for every inter-residue bond. main does that too, so it's not from this PR.
So I made it an exact check instead:
assert test_pairs == _conect_pairs(ref_file) | inter_residue_pairsThis replaces both subset asserts. It passes on all 19 files and still fails on main only for 1igy (the 8 bonds like 82 - 82A).
Also rebased on main, the only conflict was the import block.
About CodSpeed: the 3 flagged benchmarks are the k-mer matching ones, which this PR doesn't touch, and benchmark_match[KmerTable-None] flips between ~294 and ~328 µs on unrelated PRs too, in both directions. Locally the residue start lookup adds ~30 µs to _remove_non_conect_bonds() on 1aki, which doesn't show up in set_structure() timings.
…n code _remove_non_conect_bonds() compared only res_id and chain_id to find inter-residue bonds, so a bond between e.g. residues 52 and 52A was treated as intra-residue and not written. Compare residue starts instead, like the PDBx writer does.
Compare the written CONECT records with the CONECT records of the original file and with the inter-residue bonds of the structure, matched by atom identity since serials are renumbered.
Keeping all altlocs and the original atom IDs makes the identity matching unnecessary.
The written records are exactly the original CONECT records plus the inter-residue bonds. The original files have no records for peptide and phosphodiester bonds, so plain equality with them cannot hold.
9777950 to
aad7e2e
Compare
Fixes #945
_remove_non_conect_bonds()found inter-residue bonds by comparingres_idandchain_idonly. Two residues that differ just inins_code(e.g.52and52A) looked like one residue, so a bond between them got noCONECTrecord and was lost when writing a PDB file.This now compares residue starts via
get_residue_starts_for(), the same way_filter_bonds()does for the PDBx writer.Added one small test with a disulfide between
52and52A. It fails onmainand passes here.tests/structure/io/test_pdb.py,test_pdbx.pyandtests/structure/test_bonds.pypass,ruff0.9.7 andpyrightare clean on the changed file.