Fix disulfide bond detection for split cysteines with altlocs

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
 -------------
 

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)

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
     )

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