From 20b60f5bfc5dccd28043b5b306d946198fed5659 Mon Sep 17 00:00:00 2001 From: Mike German Date: Wed, 29 Jul 2026 23:51:38 -0400 Subject: [PATCH 1/3] Fix dihedral_backbone() returning garbage angles across chain breaks dihedral_backbone() and nucleotide_dihedral_backbone() built their phi/psi/omega (and alpha/beta/gamma/delta/epsilon/zeta) coordinate arrays purely from positional adjacency between residues, with no check that consecutive residues were actually bonded. When two residues were adjacent in the AtomArray but not chemically connected (missing loop, concatenated chains), the dihedrals spanning that gap came out as finite garbage instead of NaN. Add a bond-distance check (C-N for the peptide backbone, O3'-P for the phosphate backbone) using the same 1.2-1.8 A tolerance already used by check_backbone_continuity()/filter_linear_bond_continuity(), and NaN out the angles that span a junction where the distance falls outside that range. Applied to both dihedral_backbone() (fixes #744) and the identical defect in nucleotide_dihedral_backbone(), which was not filed but has the same root cause. Fixes #744. --- src/biotite/structure/geometry.py | 48 ++++++++++++++++++++ tests/structure/test_geometry.py | 74 +++++++++++++++++++++++++++++++ 2 files changed, 122 insertions(+) diff --git a/src/biotite/structure/geometry.py b/src/biotite/structure/geometry.py index 32453deb8..1c89610ec 100644 --- a/src/biotite/structure/geometry.py +++ b/src/biotite/structure/geometry.py @@ -60,6 +60,14 @@ NDArray3, ) +# Bond length range (in Angstrom) used to decide whether two positionally +# adjacent residues are actually connected via a peptide (C-N) or +# phosphodiester (O3'-P) bond. +# Mirrors the defaults of `filter_linear_bond_continuity()` / +# `check_backbone_continuity()` in `biotite.structure.integrity`. +_BACKBONE_BOND_LENGTH_MIN = 1.2 +_BACKBONE_BOND_LENGTH_MAX = 1.8 + # The names of the atoms participating in chi angle _CHI_ATOMS = { "ARG": [ @@ -647,6 +655,9 @@ def dihedral_backbone( `phi` is not defined at the N-terminus, `psi` and `omega` are not defined at the C-terminus. In these places the arrays have *NaN* values. + The same is true if two consecutive residues are not actually connected, + e.g. due to a chain break or missing loop, as indicated by an implausible + C-N bond length. If an :class:`AtomArrayStack` is given, the output angles are 2-dimensional, the first dimension corresponds to the model number. """ @@ -694,6 +705,23 @@ def dihedral_backbone( psi = dihedral(*(coord_for_psi[..., i] for i in range(4))) omg = dihedral(*(coord_for_omg[..., i] for i in range(4))) + # Two residues that are merely positionally adjacent in the atom array + # (e.g. due to a missing loop or concatenated chains) are not + # necessarily bonded to each other. + # Hence, the C-N distance between them is checked to only compute + # dihedral angles for backbone atoms that are actually connected. + c_n_dist = np.linalg.norm(coord_c[..., :-1, :] - coord_n[..., 1:, :], axis=-1) + is_discontinuous = ~( + (c_n_dist >= _BACKBONE_BOND_LENGTH_MIN) + & (c_n_dist <= _BACKBONE_BOND_LENGTH_MAX) + ) + # `psi` and `omega` of residue `i` as well as `phi` of residue `i + 1` + # are defined using atoms from both sides of the junction between + # residue `i` and `i + 1` + psi[..., :-1][is_discontinuous] = np.nan + omg[..., :-1][is_discontinuous] = np.nan + phi[..., 1:][is_discontinuous] = np.nan + # `dihedral`'s union return includes a scalar `np.floating` branch # that only fires on rank-0 inputs; here the inputs are always at # least 2D so the scalar branch is unreachable. @@ -851,6 +879,9 @@ def nucleotide_dihedral_backbone( :math:`\alpha` is not defined at the 5'-terminus, :math:`\epsilon` and :math:`\zeta` are not defined at the 3'-terminus. In these places the arrays have *NaN* values. + The same is true if two consecutive residues are not actually connected, + e.g. due to a chain break or missing residues, as indicated by an + implausible O3'-P bond length. If an :class:`AtomArrayStack` is given, the output angles are 2-dimensional, the first dimension corresponds to the model number. @@ -932,6 +963,23 @@ def nucleotide_dihedral_backbone( epsilon = dihedral(*(coord_for_epsilon[..., i] for i in range(4))) zeta = dihedral(*(coord_for_zeta[..., i] for i in range(4))) + # Two residues that are merely positionally adjacent in the atom array + # (e.g. due to a missing loop or concatenated chains) are not + # necessarily bonded to each other. + # Hence, the O3'-P distance between them is checked to only compute + # dihedral angles for backbone atoms that are actually connected. + o3p_p_dist = np.linalg.norm(coord_o3p[..., :-1, :] - coord_p[..., 1:, :], axis=-1) + is_discontinuous = ~( + (o3p_p_dist >= _BACKBONE_BOND_LENGTH_MIN) + & (o3p_p_dist <= _BACKBONE_BOND_LENGTH_MAX) + ) + # `epsilon` and `zeta` of residue `i` as well as `alpha` of residue + # `i + 1` are defined using atoms from both sides of the junction + # between residue `i` and `i + 1` + epsilon[..., :-1][is_discontinuous] = np.nan + zeta[..., :-1][is_discontinuous] = np.nan + alpha[..., 1:][is_discontinuous] = np.nan + # See note in `dihedral_backbone` about the scalar branch of # `dihedral`'s return type being unreachable here. return alpha, beta, gamma, delta, epsilon, zeta # pyright: ignore[reportReturnType] diff --git a/tests/structure/test_geometry.py b/tests/structure/test_geometry.py index 46f47e4be..a315e54ed 100644 --- a/tests/structure/test_geometry.py +++ b/tests/structure/test_geometry.py @@ -98,6 +98,80 @@ def test_dihedral_backbone_consistency(multi_model): assert test_ome == pytest.approx(ref_omega, abs=1e-3, nan_ok=True) +@pytest.mark.parametrize("multi_model", [False, True]) +def test_dihedral_backbone_chain_break(multi_model): + """ + :func:`dihedral_backbone()` must not compute dihedral angles across a + chain break, i.e. two residues that are positionally adjacent in the + :class:`AtomArray` but not actually bonded. + This is simulated by taking two unrelated fragments of the same chain + and translating one of them far away, while keeping the residue IDs + perfectly continuous, to ensure the detection is based on the actual + bond distance and not on residue numbering. + """ + pdbx_file = pdbx.BinaryCIFFile.read(data_dir("structure") / "pdb" / "1l2y.bcif") + atoms = pdbx.get_structure(pdbx_file, model=1) + + fragment_1 = atoms[np.isin(atoms.res_id, [1, 2])].copy() + fragment_2 = atoms[np.isin(atoms.res_id, [15, 16])].copy() + # Renumber so residue IDs are contiguous with `fragment_1`, although + # the two fragments are not physically connected + fragment_2.res_id = fragment_2.res_id - 15 + 3 + fragment_2.coord = fragment_2.coord + np.array([1000, 0, 0], dtype=np.float32) + combined = fragment_1 + fragment_2 + if multi_model: + combined = struc.stack([combined] * 2) + + phi, psi, omega = struc.dihedral_backbone(combined) + + # The junction between residue index 1 (res_id 2) and residue index 2 + # (res_id 3) is not an actual bond -> angles spanning it must be NaN + assert np.all(np.isnan(psi[..., 1])) + assert np.all(np.isnan(omega[..., 1])) + assert np.all(np.isnan(phi[..., 2])) + # All other angles within a fragment are unaffected and must remain + # finite + assert np.all(np.isfinite(phi[..., 1])) + assert np.all(np.isfinite(psi[..., 0])) + assert np.all(np.isfinite(omega[..., 0])) + assert np.all(np.isfinite(psi[..., 2])) + assert np.all(np.isfinite(omega[..., 2])) + assert np.all(np.isfinite(phi[..., 3])) + + +@pytest.mark.parametrize("multi_model", [False, True]) +def test_nucleotide_dihedral_backbone_chain_break(multi_model): + """ + Analogous to `test_dihedral_backbone_chain_break()`, but for + :func:`nucleotide_dihedral_backbone()`. + """ + pdbx_file = pdbx.BinaryCIFFile.read(data_dir("structure") / "pdb" / "4p5j.bcif") + atoms = pdbx.get_structure(pdbx_file, model=1) + atoms = atoms[struc.filter_canonical_nucleotides(atoms)] + + fragment_1 = atoms[np.isin(atoms.res_id, [1, 2])].copy() + fragment_2 = atoms[np.isin(atoms.res_id, [15, 16])].copy() + fragment_2.res_id = fragment_2.res_id - 15 + 3 + fragment_2.coord = fragment_2.coord + np.array([1000, 0, 0], dtype=np.float32) + combined = fragment_1 + fragment_2 + if multi_model: + combined = struc.stack([combined] * 2) + + alpha, beta, gamma, delta, epsilon, zeta = struc.nucleotide_dihedral_backbone( + combined + ) + + assert np.all(np.isnan(epsilon[..., 1])) + assert np.all(np.isnan(zeta[..., 1])) + assert np.all(np.isnan(alpha[..., 2])) + assert np.all(np.isfinite(epsilon[..., 0])) + assert np.all(np.isfinite(zeta[..., 0])) + assert np.all(np.isfinite(alpha[..., 1])) + assert np.all(np.isfinite(epsilon[..., 2])) + assert np.all(np.isfinite(zeta[..., 2])) + assert np.all(np.isfinite(alpha[..., 3])) + + @pytest.mark.parametrize("multi_model", [False, True]) def test_dihedral_side_chain_consistency(multi_model): """ From 282704cfbdea16bf559ee3deabc20736e2838127 Mon Sep 17 00:00:00 2001 From: Mike German Date: Thu, 6 Aug 2026 19:17:18 -0400 Subject: [PATCH 2/3] Combine the chain break tests into one parametrized test The two chain break tests differed only in the input structure and the names of the returned angles, so parametrize a single test on dihedral_backbone/nucleotide_dihedral_backbone instead. Look the returned angles up by name so the assertions can be grouped per angle in ascending residue index order. --- tests/structure/test_geometry.py | 85 ++++++++++++++------------------ 1 file changed, 37 insertions(+), 48 deletions(-) diff --git a/tests/structure/test_geometry.py b/tests/structure/test_geometry.py index a315e54ed..05cc5419a 100644 --- a/tests/structure/test_geometry.py +++ b/tests/structure/test_geometry.py @@ -99,18 +99,39 @@ def test_dihedral_backbone_consistency(multi_model): @pytest.mark.parametrize("multi_model", [False, True]) -def test_dihedral_backbone_chain_break(multi_model): +@pytest.mark.parametrize( + "function_name", ["dihedral_backbone", "nucleotide_dihedral_backbone"] +) +def test_dihedral_backbone_chain_break(function_name, multi_model): """ - :func:`dihedral_backbone()` must not compute dihedral angles across a - chain break, i.e. two residues that are positionally adjacent in the - :class:`AtomArray` but not actually bonded. + :func:`dihedral_backbone()` and :func:`nucleotide_dihedral_backbone()` + must not compute dihedral angles across a chain break, i.e. two + residues that are positionally adjacent in the :class:`AtomArray` but + not actually bonded. This is simulated by taking two unrelated fragments of the same chain and translating one of them far away, while keeping the residue IDs perfectly continuous, to ensure the detection is based on the actual bond distance and not on residue numbering. """ - pdbx_file = pdbx.BinaryCIFFile.read(data_dir("structure") / "pdb" / "1l2y.bcif") + if function_name == "dihedral_backbone": + pdb_id = "1l2y" + angle_names = ["phi", "psi", "omega"] + # The angles reaching into the following residue and the angle + # reaching back into the preceding one + trailing_names = ["psi", "omega"] + leading_name = "phi" + else: + pdb_id = "4p5j" + angle_names = ["alpha", "beta", "gamma", "delta", "epsilon", "zeta"] + trailing_names = ["epsilon", "zeta"] + leading_name = "alpha" + + pdbx_file = pdbx.BinaryCIFFile.read( + data_dir("structure") / "pdb" / f"{pdb_id}.bcif" + ) atoms = pdbx.get_structure(pdbx_file, model=1) + if function_name == "nucleotide_dihedral_backbone": + atoms = atoms[struc.filter_canonical_nucleotides(atoms)] fragment_1 = atoms[np.isin(atoms.res_id, [1, 2])].copy() fragment_2 = atoms[np.isin(atoms.res_id, [15, 16])].copy() @@ -122,54 +143,22 @@ def test_dihedral_backbone_chain_break(multi_model): if multi_model: combined = struc.stack([combined] * 2) - phi, psi, omega = struc.dihedral_backbone(combined) + angles = dict( + zip(angle_names, getattr(struc, function_name)(combined), strict=True) + ) # The junction between residue index 1 (res_id 2) and residue index 2 # (res_id 3) is not an actual bond -> angles spanning it must be NaN - assert np.all(np.isnan(psi[..., 1])) - assert np.all(np.isnan(omega[..., 1])) - assert np.all(np.isnan(phi[..., 2])) + for name in trailing_names: + assert np.all(np.isnan(angles[name][..., 1])) + assert np.all(np.isnan(angles[leading_name][..., 2])) # All other angles within a fragment are unaffected and must remain # finite - assert np.all(np.isfinite(phi[..., 1])) - assert np.all(np.isfinite(psi[..., 0])) - assert np.all(np.isfinite(omega[..., 0])) - assert np.all(np.isfinite(psi[..., 2])) - assert np.all(np.isfinite(omega[..., 2])) - assert np.all(np.isfinite(phi[..., 3])) - - -@pytest.mark.parametrize("multi_model", [False, True]) -def test_nucleotide_dihedral_backbone_chain_break(multi_model): - """ - Analogous to `test_dihedral_backbone_chain_break()`, but for - :func:`nucleotide_dihedral_backbone()`. - """ - pdbx_file = pdbx.BinaryCIFFile.read(data_dir("structure") / "pdb" / "4p5j.bcif") - atoms = pdbx.get_structure(pdbx_file, model=1) - atoms = atoms[struc.filter_canonical_nucleotides(atoms)] - - fragment_1 = atoms[np.isin(atoms.res_id, [1, 2])].copy() - fragment_2 = atoms[np.isin(atoms.res_id, [15, 16])].copy() - fragment_2.res_id = fragment_2.res_id - 15 + 3 - fragment_2.coord = fragment_2.coord + np.array([1000, 0, 0], dtype=np.float32) - combined = fragment_1 + fragment_2 - if multi_model: - combined = struc.stack([combined] * 2) - - alpha, beta, gamma, delta, epsilon, zeta = struc.nucleotide_dihedral_backbone( - combined - ) - - assert np.all(np.isnan(epsilon[..., 1])) - assert np.all(np.isnan(zeta[..., 1])) - assert np.all(np.isnan(alpha[..., 2])) - assert np.all(np.isfinite(epsilon[..., 0])) - assert np.all(np.isfinite(zeta[..., 0])) - assert np.all(np.isfinite(alpha[..., 1])) - assert np.all(np.isfinite(epsilon[..., 2])) - assert np.all(np.isfinite(zeta[..., 2])) - assert np.all(np.isfinite(alpha[..., 3])) + for name in trailing_names: + assert np.all(np.isfinite(angles[name][..., 0])) + assert np.all(np.isfinite(angles[name][..., 2])) + assert np.all(np.isfinite(angles[leading_name][..., 1])) + assert np.all(np.isfinite(angles[leading_name][..., 3])) @pytest.mark.parametrize("multi_model", [False, True]) From 74200bb7eddc5df479d8e45b8fff4f2a8dacb991 Mon Sep 17 00:00:00 2001 From: Mike German Date: Mon, 10 Aug 2026 08:38:51 -0400 Subject: [PATCH 3/3] Narrow dihedral results to arrays so the chain-break masking type-checks --- src/biotite/structure/geometry.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/biotite/structure/geometry.py b/src/biotite/structure/geometry.py index 1c89610ec..900b13fdb 100644 --- a/src/biotite/structure/geometry.py +++ b/src/biotite/structure/geometry.py @@ -701,9 +701,12 @@ def dihedral_backbone( coord_for_omg[..., 0:-1, :, 3] = coord_ca[..., 1:, :] # fmt: on - phi = dihedral(*(coord_for_phi[..., i] for i in range(4))) - psi = dihedral(*(coord_for_psi[..., i] for i in range(4))) - omg = dihedral(*(coord_for_omg[..., i] for i in range(4))) + # `np.asarray` strips the unreachable scalar (`np.floating`) branch of + # `dihedral`'s return type (see the note on the return statement below), + # so the in-place NaN masking further down type-checks. + phi = np.asarray(dihedral(*(coord_for_phi[..., i] for i in range(4)))) + psi = np.asarray(dihedral(*(coord_for_psi[..., i] for i in range(4)))) + omg = np.asarray(dihedral(*(coord_for_omg[..., i] for i in range(4)))) # Two residues that are merely positionally adjacent in the atom array # (e.g. due to a missing loop or concatenated chains) are not @@ -956,12 +959,16 @@ def nucleotide_dihedral_backbone( coord_for_zeta[..., 0:-1, :, 3] = coord_o5p[..., 1:, :] # fmt: on - alpha = dihedral(*(coord_for_alpha[..., i] for i in range(4))) + # `alpha`, `epsilon` and `zeta` are wrapped in `np.asarray` to strip the + # unreachable scalar (`np.floating`) branch of `dihedral`'s return type + # (see the note on the return statement below), so the in-place NaN + # masking further down type-checks. + alpha = np.asarray(dihedral(*(coord_for_alpha[..., i] for i in range(4)))) beta = dihedral(*(coord_for_beta[..., i] for i in range(4))) gamma = dihedral(*(coord_for_gamma[..., i] for i in range(4))) delta = dihedral(*(coord_for_delta[..., i] for i in range(4))) - epsilon = dihedral(*(coord_for_epsilon[..., i] for i in range(4))) - zeta = dihedral(*(coord_for_zeta[..., i] for i in range(4))) + epsilon = np.asarray(dihedral(*(coord_for_epsilon[..., i] for i in range(4)))) + zeta = np.asarray(dihedral(*(coord_for_zeta[..., i] for i in range(4)))) # Two residues that are merely positionally adjacent in the atom array # (e.g. due to a missing loop or concatenated chains) are not