diff --git a/docs/changelog.rst b/docs/changelog.rst index 615ec9aa..32780de4 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -2,6 +2,14 @@ Changelog ========= +Unreleased +---------- + +Fixes that change results +~~~~~~~~~~~~~~~~~~~~~~~~~ +- Disulfide detection no longer bonds conformer ``A`` of one cysteine to conformer ``B`` of its partner: an ``SG`` with an altloc pairs only with the same altloc or a blank one + + Version 0.7.0 ------------- diff --git a/tests/unit/topology/test_links.py b/tests/unit/topology/test_links.py index 5c136b0d..fda075be 100644 --- a/tests/unit/topology/test_links.py +++ b/tests/unit/topology/test_links.py @@ -264,3 +264,59 @@ def test_split_sg_disulfide_restrains_each_conformer(pdb_dir): cb53b = _row(model, "A", 53, "CB", "B") vdw = restraints.restraints["vdw"]["indices"].tolist() assert not any({a, b} == {cb53b, sg21} for a, b in vdw) + + +def _split_cys21(src, dst): + """Copy of 3A5V with CYS21's CB and SG split into altlocs A and B. + + B is A shifted by 0.3 Å along x, so every SG of CYS21 lies within bonding distance + of every SG of CYS53. + """ + lines = [] + for line in src.read_text().splitlines(): + split = ( + line.startswith("ATOM") + and line[17:20] == "CYS" + and line[22:26] == " 21" + and line[12:16].strip() in ("CB", "SG") + ) + if not split: + lines.append(line) + continue + shifted = f"{float(line[30:38]) + 0.3:8.3f}" + conformer_b = line[:16] + "B" + line[17:30] + shifted + line[38:54] + lines.append(line[:16] + "A" + line[17:54] + " 0.50" + line[60:]) + lines.append(conformer_b + " 0.50" + line[60:]) + dst.write_text("\n".join(lines) + "\n") + return dst + + +@pytest.mark.unit +def test_disulfide_bonds_stay_within_one_conformer(pdb_dir, tmp_path): + """Two split cysteines bond A to A and B to B, never A to B. + + Each bond carries its own conformer's angles and torsion, as for one split + cysteine against an unsplit partner. + """ + model = _load(_split_cys21(pdb_dir / "3A5V.pdb", tmp_path / "3A5V_split21.pdb")) + restraints = model.restraints.restraints + sg = {alt: _row(model, "A", 21, "SG", alt) for alt in ("A", "B")} + sg53 = {alt: _row(model, "A", 53, "SG", alt) for alt in ("A", "B")} + sg_rows = set(sg.values()) | set(sg53.values()) + + bonds = restraints["bond"]["disulfide"]["indices"].tolist() + found = {frozenset(row) for row in bonds if set(row) <= sg_rows} + assert found == { + frozenset((sg["A"], sg53["A"])), + frozenset((sg["B"], sg53["B"])), + } + + for alt in ("A", "B"): + cb21, cb53 = _row(model, "A", 21, "CB", alt), _row(model, "A", 53, "CB", alt) + for edge_type, expected in ( + ("angle", {(cb21, sg[alt], sg53[alt]), (sg[alt], sg53[alt], cb53)}), + ("torsion", {(cb21, sg[alt], sg53[alt], cb53)}), + ): + rows = restraints[edge_type]["disulfide"]["indices"].tolist() + found = [tuple(row) for row in rows if {sg[alt], sg53[alt]} <= set(row)] + assert sorted(found) == sorted(expected), (edge_type, alt) diff --git a/torchref/topology/build.py b/torchref/topology/build.py index c0cd86cd..cbf576a1 100644 --- a/torchref/topology/build.py +++ b/torchref/topology/build.py @@ -1002,7 +1002,7 @@ def build_topology_with_values( residue_of_row[row] = r sg_rows = [row for row in range(len(cols["name"])) if cols["name"][row] == "SG"] - disulfide_pairs = find_disulfide_links(sg_rows, residue_of_row, xyz) + disulfide_pairs = find_disulfide_links(sg_rows, residue_of_row, xyz, cols["altloc"]) disulfide, disulfide_values = _disulfide_edges( topology, residue_of_row, disulfide_pairs, link_dict, verbose ) diff --git a/torchref/topology/residue_graph.py b/torchref/topology/residue_graph.py index 81be493b..7ead1a6f 100644 --- a/torchref/topology/residue_graph.py +++ b/torchref/topology/residue_graph.py @@ -312,6 +312,7 @@ def find_disulfide_links( sg_rows: Sequence[int], residue_of_row: Dict[int, int], xyz: torch.Tensor, + altloc: Optional[Sequence[str]] = None, ) -> List[Tuple[int, int]]: """``SG`` atom pairs within :data:`DISULFIDE_MAX_DISTANCE` in different residues. @@ -321,6 +322,11 @@ def find_disulfide_links( residue -- its own alternative conformers -- are within bonding distance of each other and are excluded by the differing-residue test. + A bond joins one state of the model, so an ``SG`` in conformer ``A`` pairs only + with an ``SG`` in conformer ``A`` or one without an altloc. Conformer ``A`` of one + cysteine and conformer ``B`` of its partner are often within bonding distance of + each other too, and restraining that pair would pull the two conformers together. + Parameters ---------- sg_rows : sequence of int @@ -329,6 +335,9 @@ def find_disulfide_links( ``{atom row: residue index}`` for those rows. xyz : torch.Tensor Cartesian coordinates, shape ``(N, 3)``. + altloc : sequence of str, optional + Per-atom altloc labels indexed by row, ``' '`` (or empty) for an atom in no + alternative conformation. Omitted, every pair is taken to share a state. Returns ------- @@ -347,8 +356,13 @@ def find_disulfide_links( pairs = [] for i, j in zip(a[hit].cpu().tolist(), b[hit].cpu().tolist()): row_i, row_j = rows[i], rows[j] - if residue_of_row[row_i] != residue_of_row[row_j]: - pairs.append((row_i, row_j)) + if residue_of_row[row_i] == residue_of_row[row_j]: + continue + if altloc is not None: + alt_i, alt_j = str(altloc[row_i]).strip(), str(altloc[row_j]).strip() + if alt_i and alt_j and alt_i != alt_j: + continue + pairs.append((row_i, row_j)) return pairs