Fix overlapping assembly copies and invalid protein interfaces

docs/dataset.md

@@ -175,6 +175,11 @@ for each side, the number of contacting residue pairs, and PRODIGY-cryst
 features. Full deposited sequences and source taxonomies are in the chain and
 entry tables.
 
+Interfaces require at least seven contacting residues per side by default and
+at least one known residue in each chain. Near-coincident assembly copies are
+removed; pairs with at least ten backbone-clashing residues per side within
+1.5 Å are excluded while their chains remain available for other interfaces.
+
 `interface_representatives` stores the full-interface representatives,
 while `interface_half_representatives` stores individual interface
 sides. The latter is useful when one wants diverse protein surfaces without

src/plinder/data/annotations/cif_utils.py

@@ -776,7 +776,9 @@ def drop_self_clashing_symmetry_copies(
     C-alpha atoms: a chain is removed if at least five and more than 10%, or
     at least twenty regardless of chain length, lie within 2 Å of previously
     retained protein chains. Such C-alpha distances are nonphysical for
-    separate chains. Requires the ``sym_id`` and ``label_asym_id`` annotations.
+    separate chains. Near-coincident protein copies are also removed when
+    more than half of their C-alpha atoms lie within 2.5 Å of a single kept
+    chain. Requires the ``sym_id`` and ``label_asym_id`` annotations.
     """
     categories = set(assembly.get_annotation_categories())
     if not {"sym_id", "label_asym_id"}.issubset(categories):
@@ -835,11 +837,11 @@ def drop_self_clashing_symmetry_copies(
         )
         unique_chains, chain_index = np.unique(chain_names, return_inverse=True)
         chain_sizes = np.bincount(chain_index)
-        eligible = [index for index, size in enumerate(chain_sizes) if size >= 10]
+        eligible = [index for index, size in enumerate(chain_sizes) if size >= 5]
         if len(eligible) > 1:
             cell_list = struc.CellList(ca, cell_size=3.0)
             retained: list[int] = []
-            dropped_chains: list[tuple[str, int, float]] = []
+            dropped_chains: list[tuple[str, int, float, float]] = []
             for index in sorted(
                 eligible,
                 key=lambda value: (
@@ -849,33 +851,62 @@ def drop_self_clashing_symmetry_copies(
             ):
                 query_mask = chain_index == index
                 if retained:
-                    neighbors = cell_list.get_atoms(ca.coord[query_mask], radius=2.0)
-                    clashing = (
-                        (neighbors >= 0)
-                        & np.isin(chain_index[np.clip(neighbors, 0, None)], retained)
-                    ).any(axis=1)
+                    query_coords = ca.coord[query_mask]
+                    neighbors = cell_list.get_atoms(query_coords, radius=2.5)
+                    retained_neighbors = (neighbors >= 0) & np.isin(
+                        chain_index[np.clip(neighbors, 0, None)], retained
+                    )
+                    query_rows, neighbor_columns = np.nonzero(retained_neighbors)
+                    neighbor_indices = neighbors[query_rows, neighbor_columns]
+                    clashing = np.zeros(len(query_coords), dtype=bool)
+                    distances = np.linalg.norm(
+                        query_coords[query_rows] - ca.coord[neighbor_indices], axis=1
+                    )
+                    clashing[query_rows[distances <= 2.0]] = True
+                    # Count residues per partner, not multiple nearby atoms or
+                    # the combined contacts to several distinct chains.
+                    residue_partners = np.unique(
+                        np.column_stack((query_rows, chain_index[neighbor_indices])),
+                        axis=0,
+                    )
+                    _, partner_counts = np.unique(
+                        residue_partners[:, 1], return_counts=True
+                    )
+                    overlap_fraction = float(
+                        partner_counts.max(initial=0) / len(query_coords)
+                    )
+                    near_copy = overlap_fraction > 0.5
                     fraction = float(clashing.mean())
                     clash_count = int(clashing.sum())
-                    if clash_count >= 20 or (clash_count >= 5 and fraction > 0.1):
+                    if (
+                        near_copy
+                        or clash_count >= 20
+                        or (clash_count >= 5 and fraction > 0.1)
+                    ):
                         sym_id, asym_id = unique_chains[index].split(".", 1)
                         keep &= ~(
                             (assembly.sym_id == int(sym_id))
                             & (assembly.label_asym_id == asym_id)
                         )
                         dropped_chains.append(
-                            (str(unique_chains[index]), clash_count, fraction)
+                            (
+                                str(unique_chains[index]),
+                                clash_count,
+                                fraction,
+                                overlap_fraction,
+                            )
                         )
                         continue
                 retained.append(index)
             if dropped_chains:
                 LOG.info(
                     "%s: dropped %d clashing protein chains (first 10 chain, "
-                    "C-alpha clash count, fraction): %s",
+                    "2 Å C-alpha clash count, fraction, 2.5 Å partner overlap): %s",
                     context,
                     len(dropped_chains),
                     [
-                        (chain, count, round(fraction, 2))
-                        for chain, count, fraction in dropped_chains[:10]
+                        (chain, count, round(fraction, 2), round(overlap, 2))
+                        for chain, count, fraction, overlap in dropped_chains[:10]
                     ],
                 )
     if keep.all():

src/plinder/data/annotations/interface_utils.py

@@ -19,6 +19,8 @@
 
 PROTEIN_BACKBONE_ATOMS = frozenset({"N", "CA", "C", "O"})
 DEFAULT_MIN_INTERFACE_RESIDUES = 7
+BACKBONE_CLASH_DISTANCE = 1.5
+MIN_BACKBONE_CLASH_RESIDUES = 10
 MIN_INTERFACE_RESIDUES_METADATA_KEY = b"plinder.interface.min_interface_residues"
 PRODIGY_CONTACT_RADIUS = 5.0
 ProteinChainContacts = dict[tuple[str, str], set[tuple[int, int]]]
@@ -400,6 +402,18 @@ def find_protein_chain_contacts(
     }
     if len(protein_chains) < 2:
         return {}
+    return _backbone_residue_contacts(
+        atoms, protein_chains, contact_radius, spatial_index
+    )
+
+
+def _backbone_residue_contacts(
+    atoms: struc.AtomArray,
+    protein_chains: set[str],
+    radius: float,
+    spatial_index: BiounitSpatialIndex,
+) -> ProteinChainContacts:
+    """Collect backbone residue pairs using the assembly's spatial index."""
 
     protein_backbone = np.fromiter(
         (
@@ -413,9 +427,7 @@ def find_protein_chain_contacts(
     contacts: ProteinChainContacts = {}
     for atom_index in np.flatnonzero(protein_backbone):
         neighbors = np.asarray(
-            spatial_index.cell_list.get_atoms(
-                atoms.coord[atom_index], radius=contact_radius
-            ),
+            spatial_index.cell_list.get_atoms(atoms.coord[atom_index], radius=radius),
             dtype=int,
         ).reshape(-1)
         neighbors = neighbors[(neighbors > atom_index) & (neighbors < len(atoms))]
@@ -458,7 +470,10 @@ def detect_protein_interfaces(
     ``N``, ``CA``, ``C`` or ``O`` atoms lies within ``contact_radius`` of a
     backbone atom in another eligible assembly-chain instance.  Interfaces
     are retained only when both sides contain at least
-    ``min_interface_residues`` resolved residues.
+    ``min_interface_residues`` resolved residues and neither chain is modeled
+    entirely as ``UNK``. Pairs with at least ten distinct residues on each
+    side involved in backbone clashes within 1.5 Å are excluded; their chains
+    remain available for other interfaces and for protein-contact counts.
 
     ``chain_pair_contact_areas`` are the assembly-wide Voronota-LT areas keyed
     by sorted chain pair; when given, each interface records its pair's area
@@ -478,12 +493,19 @@ def detect_protein_interfaces(
             "interface contact radius exceeds the biological-assembly index radius"
         )
 
+    known_backbone_chains = set(
+        atoms.chain_id[
+            np.isin(atoms.atom_name, list(PROTEIN_BACKBONE_ATOMS))
+            & (atoms.res_name != "UNK")
+        ]
+    )
     eligible_chains = {
         instance_chain
         for instance_chain in spatial_index.chain_ids
         if (chain := chains.get(_asym_id(instance_chain))) is not None
         and "polypeptide" in chain.chain_type_str.lower()
         and chain.length >= min_chain_length
+        and instance_chain in known_backbone_chains
     }
     if len(eligible_chains) < 2:
         return []
@@ -498,11 +520,25 @@ def detect_protein_interfaces(
         if chain_contacts is None
         else chain_contacts
     )
+    clashes = _backbone_residue_contacts(
+        atoms,
+        eligible_chains,
+        min(BACKBONE_CLASH_DISTANCE, contact_radius),
+        spatial_index,
+    )
+    clashing_pairs = {
+        pair
+        for pair, residues in clashes.items()
+        if len({left for left, _ in residues}) >= MIN_BACKBONE_CLASH_RESIDUES
+        and len({right for _, right in residues}) >= MIN_BACKBONE_CLASH_RESIDUES
+    }
 
     interfaces: list[ProteinInterface] = []
     for (chain_1, chain_2), residue_pairs in sorted(contacts.items()):
         if chain_1 not in eligible_chains or chain_2 not in eligible_chains:
             continue
+        if (chain_1, chain_2) in clashing_pairs:
+            continue
         chain_1_numbers, chain_1_indices = _residue_mappings(
             instance_chain=chain_1,
             residue_numbers=(pair[0] for pair in residue_pairs),

tests/test_annotations.py

@@ -171,6 +171,90 @@ def test_protein_contacts_include_short_polypeptide_partners() -> None:
     )
 
 
+@pytest.mark.parametrize("clashing_residues", [9, 10])
+@pytest.mark.parametrize("cached_contacts", [False, True])
+def test_protein_interface_clashes_exclude_pairs_not_chains(
+    clashing_residues, cached_contacts
+):
+    atoms = struc.AtomArray(36)
+    atoms.chain_id = np.repeat(["1.A", "1.B", "1.C"], 12)
+    atoms.res_id = np.tile(np.arange(1, 13), 3)
+    atoms.atom_name[:] = "CA"
+    atoms.element[:] = "C"
+    atoms.res_name[:] = "ALA"
+    base = np.column_stack((np.arange(12) * 3.8, np.zeros(12), np.zeros(12)))
+    close = base + [0.0, 1.0, 0.0]
+    close[clashing_residues:, 1] = 3.0
+    atoms.coord = np.concatenate((base, close, base + [0.0, 6.0, 0.0]))
+    chains = {name: _interface_test_chain(list(range(1, 13))) for name in "ABC"}
+    spatial_index = BiounitSpatialIndex.from_atoms(atoms, 10.0)
+    contacts = find_protein_chain_contacts(
+        atoms, chains=chains, contact_radius=10.0, spatial_index=spatial_index
+    )
+
+    interfaces = detect_protein_interfaces(
+        atoms,
+        pdb_id="1abc",
+        biounit_id="1",
+        chains=chains,
+        spatial_index=spatial_index,
+        chain_contacts=contacts if cached_contacts else None,
+    )
+
+    pairs = {(row.chain_1, row.chain_2) for row in interfaces}
+    expected = {("1.A", "1.C"), ("1.B", "1.C")}
+    if clashing_residues < 10:
+        expected.add(("1.A", "1.B"))
+    assert pairs == expected
+    assert ("1.A", "1.B") in contacts
+
+
+def test_protein_interface_clash_threshold_counts_residues_not_atoms():
+    atoms = struc.AtomArray(72)
+    atoms.chain_id = np.repeat(["1.A", "1.B"], 36)
+    atoms.res_id = np.tile(np.repeat(np.arange(1, 10), 4), 2)
+    atoms.atom_name = np.tile(["N", "CA", "C", "O"], 18)
+    atoms.element = np.tile(["N", "C", "C", "O"], 18)
+    atoms.res_name[:] = "ALA"
+    base = np.column_stack(
+        (np.repeat(np.arange(9) * 3.8, 4), np.zeros(36), np.zeros(36))
+    )
+    atoms.coord = np.concatenate((base, base + [0.0, 1.0, 0.0]))
+
+    interfaces = detect_protein_interfaces(
+        atoms,
+        pdb_id="1abc",
+        biounit_id="1",
+        chains={name: _interface_test_chain(list(range(1, 10))) for name in "AB"},
+    )
+
+    assert len(interfaces) == 1
+
+
+@pytest.mark.parametrize("known_residues", [0, 1, 12])
+@pytest.mark.parametrize("atom_name", ["CA", "N"])
+def test_protein_interfaces_require_some_known_residues(known_residues, atom_name):
+    atoms = struc.AtomArray(24)
+    atoms.chain_id = np.repeat(["1.A", "1.B"], 12)
+    atoms.res_id = np.tile(np.arange(1, 13), 2)
+    atoms.atom_name[:] = atom_name
+    atoms.element[:] = "C" if atom_name == "CA" else "N"
+    atoms.res_name[:] = "ALA"
+    atoms.res_name[12:] = "UNK"
+    atoms.res_name[12 : 12 + known_residues] = "ALA"
+    base = np.column_stack((np.arange(12) * 3.8, np.zeros(12), np.zeros(12)))
+    atoms.coord = np.concatenate((base, base + [0.0, 4.5, 0.0]))
+
+    interfaces = detect_protein_interfaces(
+        atoms,
+        pdb_id="1abc",
+        biounit_id="1",
+        chains={name: _interface_test_chain(list(range(1, 13))) for name in "AB"},
+    )
+
+    assert len(interfaces) == int(known_residues > 0)
+
+
 def test_interface_system_id_is_unordered_and_rejects_self_interfaces():
     assert interface_system_id("1ABC", "2", "2.B", "1.A") == interface_system_id(
         "1abc", "2", "1.A", "2.B"

tests/test_cif_utils.py

@@ -1443,3 +1443,52 @@ def test_drop_long_protein_chain_with_many_absolute_ca_clashes():
     kept = drop_self_clashing_symmetry_copies(atoms)
 
     assert set(kept.label_asym_id) == {"A", "C"}
+
+
+@pytest.mark.parametrize("length", [5, 9, 20])
+@pytest.mark.parametrize("offset,removed", [(2.0, True), (2.5, True), (2.6, False)])
+def test_drop_near_coincident_protein_copies(length, offset, removed):
+    import biotite.structure as struc
+
+    from plinder.data.annotations.cif_utils import drop_self_clashing_symmetry_copies
+
+    base = np.column_stack(
+        (np.arange(length) * 3.8, np.zeros(length), np.zeros(length))
+    )
+    atoms = struc.AtomArray(2 * length)
+    atoms.coord = np.concatenate((base, base + [0.0, offset, 0.0]))
+    atoms.set_annotation("label_asym_id", np.repeat(["A", "B"], length))
+    atoms.set_annotation("sym_id", np.zeros(len(atoms), dtype=int))
+    atoms.chain_id = atoms.label_asym_id.copy()
+    atoms.atom_name[:] = "CA"
+    atoms.element[:] = "C"
+    atoms.res_name[:] = "ALA"
+    atoms.res_id = np.tile(np.arange(1, length + 1), 2)
+
+    kept = drop_self_clashing_symmetry_copies(atoms)
+
+    assert set(kept.label_asym_id) == ({"A"} if removed else {"A", "B"})
+
+
+def test_near_copy_requires_majority_overlap_with_one_partner():
+    import biotite.structure as struc
+
+    from plinder.data.annotations.cif_utils import drop_self_clashing_symmetry_copies
+
+    base = np.column_stack((np.arange(20) * 3.8, np.zeros(20), np.zeros(20)))
+    other = base + [0.0, 20.0, 0.0]
+    mixed = base + [0.0, 2.5, 0.0]
+    mixed[10:] = other[10:] + [0.0, 2.5, 0.0]
+    atoms = struc.AtomArray(60)
+    atoms.coord = np.concatenate((base, other, mixed))
+    atoms.set_annotation("label_asym_id", np.repeat(["A", "B", "C"], 20))
+    atoms.set_annotation("sym_id", np.zeros(len(atoms), dtype=int))
+    atoms.chain_id = atoms.label_asym_id.copy()
+    atoms.atom_name[:] = "CA"
+    atoms.element[:] = "C"
+    atoms.res_name[:] = "ALA"
+    atoms.res_id = np.tile(np.arange(1, 21), 3)
+
+    kept = drop_self_clashing_symmetry_copies(atoms)
+
+    assert set(kept.label_asym_id) == {"A", "B", "C"}