From af4d60a0ea37cf85acb51d1917be68effb40a3ce Mon Sep 17 00:00:00 2001 From: choijein0410 <246378271+choijein0410@users.noreply.github.com> Date: Sun, 4 Oct 2026 10:50:26 -0400 Subject: [PATCH] Prevent torsions and couplings across filtered residues --- barnaba/nucleic.py | 10 +- test/test_skipped_residues.py | 168 ++++++++++++++++++++++++++++++++++ 2 files changed, 173 insertions(+), 5 deletions(-) create mode 100644 test/test_skipped_residues.py diff --git a/barnaba/nucleic.py b/barnaba/nucleic.py index 5a1c482..1af4b9d 100644 --- a/barnaba/nucleic.py +++ b/barnaba/nucleic.py @@ -115,12 +115,13 @@ def get_bb_torsion_idx(self, residues=None): res = self.ok_residues[i] rr.append(self.rna_seq[i]) + # Filtering must not make nonadjacent topology residues neighbors. # alpha try: res_minus = self.ok_residues[i-1] #print res,res_minus, res_minus.chain.index,res.chain.index idxs_tmp = [res_minus.atom("O3'").index,res.atom("P").index,res.atom("O5'").index,res.atom("C5'").index] - if(i!=0 and (res_minus.chain.index == res.chain.index)): + if(i!=0 and res_minus.chain.index == res.chain.index and res_minus.index + 1 == res.index): #print "mm", idxs_tmp idxs[j,0] = idxs_tmp else: @@ -144,7 +145,7 @@ def get_bb_torsion_idx(self, residues=None): try: res_plus = self.ok_residues[i+1] idxs_tmp = [res.atom("C4'").index,res.atom("C3'").index,res.atom("O3'").index,res_plus.atom("P").index] - if(res_plus.chain.index == res.chain.index): + if(res_plus.chain.index == res.chain.index and res_plus.index == res.index + 1): idxs[j,4] = idxs_tmp else: pass @@ -155,7 +156,7 @@ def get_bb_torsion_idx(self, residues=None): try: res_plus = self.ok_residues[i+1] idxs_tmp = [res.atom("C3'").index,res.atom("O3'").index,res_plus.atom("P").index,res_plus.atom("O5'").index] - if(res_plus.chain.index == res.chain.index): + if(res_plus.chain.index == res.chain.index and res_plus.index == res.index + 1): idxs[j,5] = idxs_tmp else: pass @@ -278,7 +279,7 @@ def get_coupling_idx(self,residues=None): try: res_plus = self.ok_residues[i+1] idxs_tmp = [res.atom("C4'").index,res.atom("C3'").index,res.atom("O3'").index,res_plus.atom("P").index] - if(res_plus.chain.index == res.chain.index): + if(res_plus.chain.index == res.chain.index and res_plus.index == res.index + 1): idxs[j,5] = idxs_tmp else: pass @@ -296,4 +297,3 @@ def get_coupling_idx(self,residues=None): return idxs, rr - diff --git a/test/test_skipped_residues.py b/test/test_skipped_residues.py new file mode 100644 index 0000000..04f5cdf --- /dev/null +++ b/test/test_skipped_residues.py @@ -0,0 +1,168 @@ +"""Inter-residue torsions must not bridge a residue skipped by Nucleic.""" +from __future__ import absolute_import, division, print_function + +import os +import tempfile + +import barnaba as bb +import mdtraj as md +import numpy as np + + +DATA = os.path.join(os.path.dirname(__file__), "data", "UUCG.pdb") +EPSILON_COUPLINGS = ["C4Pe", "H3P"] + + +def _trajectory(): + return md.load(DATA) + + +def _skip_unknown_residue(): + traj = _trajectory() + traj.topology.residue(3).name = "UNK" + return traj + + +def _skip_incomplete_residue(): + traj = _trajectory() + missing = traj.topology.residue(3).atom("C2").index + return traj.atom_slice([atom.index for atom in traj.topology.atoms + if atom.index != missing]) + + +def _check_skipped_backbone(traj): + original, names = bb.backbone_angles_traj(_trajectory()) + actual, recognized = bb.backbone_angles_traj(traj) + kept = [i for i in range(len(names)) if i != 3] + assert recognized == [names[i] for i in kept] + expected = original[:, kept, :].copy() + # U3 cannot obtain epsilon/zeta from C5 across the skipped U4. + expected[:, 2, 4:6] = np.nan + # C5 cannot obtain alpha from U3 across the skipped U4. + expected[:, 3, 0] = np.nan + np.testing.assert_array_equal(actual, expected) + + +def _check_skipped_couplings(traj): + for raw in [False, True]: + original, names = bb.jcouplings_traj( + _trajectory(), couplings=EPSILON_COUPLINGS, raw=raw) + actual, recognized = bb.jcouplings_traj( + traj, couplings=EPSILON_COUPLINGS, raw=raw) + kept = [i for i in range(len(names)) if i != 3] + assert recognized == [names[i] for i in kept] + expected = original[:, kept, :].copy() + expected[:, 2, :] = np.nan + # Every unaffected value and pre-existing NaN must be unchanged. + np.testing.assert_array_equal(actual, expected) + + +def test_unknown_residue_does_not_bridge_backbone_torsions(): + _check_skipped_backbone(_skip_unknown_residue()) + + +def test_incomplete_residue_does_not_bridge_backbone_torsions(): + _check_skipped_backbone(_skip_incomplete_residue()) + + +def test_unknown_residue_does_not_bridge_epsilon_couplings(): + _check_skipped_couplings(_skip_unknown_residue()) + + +def test_incomplete_residue_does_not_bridge_epsilon_couplings(): + _check_skipped_couplings(_skip_incomplete_residue()) + + +def test_skipped_residue_subset_preserves_requested_order_and_missing_angles(): + traj = _skip_unknown_residue() + full, recognized = bb.backbone_angles_traj(traj) + requested = [recognized[3], recognized[2]] + selected, names = bb.backbone_angles_traj( + traj, residues=requested, angles=["zeta", "alpha", "eps"]) + assert names == requested + np.testing.assert_array_equal(selected, full[:, [3, 2], :][:, :, [5, 0, 4]]) + assert np.isnan(selected[0, 0, 1]) + assert np.isnan(selected[0, 1, [0, 2]]).all() + couplings, names = bb.jcouplings_traj( + traj, residues=requested, couplings=EPSILON_COUPLINGS) + assert names == requested + assert np.isfinite(couplings[0, 0]).all() + assert np.isnan(couplings[0, 1]).all() + + +def test_contiguous_fixture_has_only_existing_terminal_missing_angles(): + traj = _trajectory() + angles, names = bb.backbone_angles_traj( + traj, angles=["alpha", "eps", "zeta"]) + assert len(names) == 8 + expected_missing = np.zeros((1, 8, 3), dtype=bool) + expected_missing[:, 0, 0] = True + expected_missing[:, -1, 1:] = True + np.testing.assert_array_equal(np.isnan(angles), expected_missing) + couplings, _ = bb.jcouplings_traj(traj, couplings=EPSILON_COUPLINGS) + assert np.isfinite(couplings[:, :-1, :]).all() + assert np.isnan(couplings[:, -1, :]).all() + + +def test_residue_numbers_and_insertion_codes_do_not_break_topology_neighbors(): + original = _trajectory() + renamed = md.Trajectory(original.xyz.copy(), original.topology.copy()) + for residue, number in zip(renamed.topology.residues, + [10, 20, 40, 41, 75, 100, 150, 200]): + residue.resSeq = number + expected, _ = bb.backbone_angles_traj(original) + actual, _ = bb.backbone_angles_traj(renamed) + np.testing.assert_array_equal(actual, expected) + expected, _ = bb.jcouplings_traj(original, couplings=EPSILON_COUPLINGS) + actual, _ = bb.jcouplings_traj(renamed, couplings=EPSILON_COUPLINGS) + np.testing.assert_array_equal(actual, expected) + with open(DATA) as source: + lines = source.readlines() + renamed = [] + for line in lines: + if line.startswith(("ATOM ", "HETATM")): + index = int(line[22:26]) - 1 + number = 10 + index // 2 + insertion = "A" if index % 2 else " " + line = line[:22] + "%4d%s" % (number, insertion) + line[27:] + renamed.append(line) + with tempfile.TemporaryDirectory() as directory: + path = os.path.join(directory, "insertions.pdb") + with open(path, "w") as target: + target.writelines(renamed) + inserted = md.load(path) + original = _trajectory() + assert inserted.n_residues == original.n_residues + np.testing.assert_array_equal(inserted.xyz, original.xyz) + expected, _ = bb.backbone_angles_traj(original) + actual, _ = bb.backbone_angles_traj(inserted) + np.testing.assert_array_equal(actual, expected) + expected, _ = bb.jcouplings_traj(original, couplings=EPSILON_COUPLINGS) + actual, _ = bb.jcouplings_traj(inserted, couplings=EPSILON_COUPLINGS) + np.testing.assert_array_equal(actual, expected) + + +def test_chain_boundaries_still_leave_inter_residue_angles_missing(): + original = _trajectory() + topology = md.Topology() + chains = [topology.add_chain(), topology.add_chain()] + atoms = {} + for residue in original.topology.residues: + target = topology.add_residue(residue.name, chains[int(residue.index >= 4)], + resSeq=residue.resSeq) + for atom in residue.atoms: + atoms[atom.index] = topology.add_atom(atom.name, atom.element, target) + for a, b in original.topology.bonds: + if atoms[a.index].residue.chain == atoms[b.index].residue.chain: + topology.add_bond(atoms[a.index], atoms[b.index]) + split = md.Trajectory(original.xyz.copy(), topology) + expected, _ = bb.backbone_angles_traj(original) + actual, names = bb.backbone_angles_traj(split) + assert names[3].endswith("_0") and names[4].endswith("_1") + expected[:, 3, 4:6] = np.nan + expected[:, 4, 0] = np.nan + np.testing.assert_array_equal(actual, expected) + expected, _ = bb.jcouplings_traj(original, couplings=EPSILON_COUPLINGS) + actual, _ = bb.jcouplings_traj(split, couplings=EPSILON_COUPLINGS) + expected[:, 3, :] = np.nan + np.testing.assert_array_equal(actual, expected)