From 39931f4bc0aa7686100cfc62e45d5e61f1b620a5 Mon Sep 17 00:00:00 2001 From: Ninjani Date: Fri, 9 Oct 2026 19:22:00 +0200 Subject: [PATCH] fix: filter overlapping assembly copies and invalid protein interfaces --- docs/dataset.md | 5 ++ src/plinder/data/annotations/cif_utils.py | 57 ++++++++++--- .../data/annotations/interface_utils.py | 44 +++++++++- tests/test_annotations.py | 84 +++++++++++++++++++ tests/test_cif_utils.py | 49 +++++++++++ 5 files changed, 222 insertions(+), 17 deletions(-) diff --git a/docs/dataset.md b/docs/dataset.md index d38441e4..bd22cd20 100644 --- a/docs/dataset.md +++ b/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 diff --git a/src/plinder/data/annotations/cif_utils.py b/src/plinder/data/annotations/cif_utils.py index 9730a899..66085cbc 100644 --- a/src/plinder/data/annotations/cif_utils.py +++ b/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(): diff --git a/src/plinder/data/annotations/interface_utils.py b/src/plinder/data/annotations/interface_utils.py index c558fdcd..8a096a38 100644 --- a/src/plinder/data/annotations/interface_utils.py +++ b/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), diff --git a/tests/test_annotations.py b/tests/test_annotations.py index 9cee0abe..6c020edb 100644 --- a/tests/test_annotations.py +++ b/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" diff --git a/tests/test_cif_utils.py b/tests/test_cif_utils.py index cbbae686..24df66f2 100644 --- a/tests/test_cif_utils.py +++ 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"}