Skip to content

Commit 15db852

Browse files
committed
chore: drop dead SSBOND cross-check + trim CysteineHit fields
`cross_check_ssbonds` was called for its return value but the caller discarded it and the function had no side effects (no stderr, no raise). That made the call a pure no-op, and the supporting code (`Ssbond` dataclass, `parsed.ssbonds` list, SSBOND parsing block, helper key functions) was a 180-line dead chain. The audit also showed nine `CysteineHit` fields that were only read inside `cross_check_ssbonds`: - resName, chainRole, bondingState - sasa, sidechainSasa, rsasa - partnerChainId, partnerResSeq, partnerIcode plus chainId / resSeq / iCode which had a single read inside the same dead function. With the per-cys PFrame out of scope per the refreshed spec, all downstream consumers (`scoring._cys_class_bump`, `scoring._has_integrity_issue`, the four count loops in main.py) only ever read `cysClass` and `sidechainRsasa`. Trimming the dataclass to those two fields drops construction boilerplate from 14 args to 2 and pulls `cysteines.py` from 212 to 154 lines. Net change: -181 lines across cysteines / structure / main / test_scoring, no behavior delta (65 pytest tests still pass).
1 parent 62b339c commit 15db852

5 files changed

Lines changed: 30 additions & 211 deletions

File tree

README.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -101,7 +101,7 @@ domain on every output.
101101
clonotype filter (R1 subset).
102102
3. **Python** (`software/liabilities-script/`):
103103
- Parses PDB (`structure.py`): v3.30 fixed-column parser; pulls
104-
ATOM / HETATM / SSBOND / `REMARK 99 PLATFORMA CDR*` records.
104+
ATOM / HETATM / `REMARK 99 PLATFORMA CDR*` records.
105105
- Runs FreeSASA against in-block Ala-X-Ala heavy-atom references
106106
(`data/heavy_atom_max_sasa.tsv`, spec R11).
107107
- Detects motifs (`motifs.py`): 11 patterns from

software/liabilities-script/cysteines.py

Lines changed: 23 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,15 @@
1111
# - disulfide : both canonical Cys present + bonded
1212
# - disulfide_broken : both canonical Cys present + unbonded
1313
# - disulfide_missing : a canonical position has no Cys (a phantom hit
14-
# is emitted at the expected position with cysClass
15-
# set to disulfide_missing)
14+
# is emitted at that role with cysClass set to
15+
# disulfide_missing)
1616
# - cys_extra : a Cys at any non-canonical position, regardless
1717
# of bonding state
18+
#
19+
# Per-cys drill-down is out of scope per the refreshed spec, so each hit
20+
# only carries the two fields actually consumed downstream: the class label
21+
# (for counts + the developability bump) and the side-chain rSASA (for the
22+
# "exposed extra" count + the buried-vs-exposed gate in the scoring bump).
1823

1924
import math
2025
from dataclasses import dataclass
@@ -29,20 +34,8 @@
2934

3035
@dataclass
3136
class CysteineHit:
32-
chainId: str
33-
resSeq: int
34-
iCode: str
35-
resName: str # always "CYS" for real hits; "-" for phantom missing entries
36-
chainRole: str # "H" / "L" / "-" when chain isn't mapped to either
3737
cysClass: str # "disulfide" | "disulfide_broken" | "disulfide_missing" | "cys_extra" | "bonded" | "unbonded"
38-
bondingState: str # "bonded" | "unbonded" | "missing"
39-
sasa: Optional[float]
40-
rsasa: Optional[float]
41-
sidechainSasa: Optional[float]
4238
sidechainRsasa: Optional[float]
43-
partnerChainId: str
44-
partnerResSeq: Optional[int]
45-
partnerIcode: str
4639

4740

4841
def _dist(a, b) -> float:
@@ -93,41 +86,20 @@ def _scan_disulfides(cys_records) -> dict[int, int]:
9386
def _classify_cys(
9487
res_seq: int,
9588
role: Optional[str],
96-
bonding_state: str,
89+
bonded: bool,
9790
canonical_keys: set[tuple[str, int]],
9891
canonical_positions: dict,
9992
) -> str:
10093
"""R23 four-state classification. Falls back to raw bonding state
10194
when numbering isn't wired (no role / no canonical positions for the
102-
scheme), so the table is still useful in the auto-detect-only case."""
95+
scheme), so the hit is still classified usefully in the auto-detect-only
96+
case (though the resulting "bonded"/"unbonded" values are inert
97+
downstream)."""
10398
if role is not None and (role, res_seq) in canonical_keys:
104-
return "disulfide" if bonding_state == "bonded" else "disulfide_broken"
99+
return "disulfide" if bonded else "disulfide_broken"
105100
if role is not None and canonical_positions:
106101
return "cys_extra"
107-
return bonding_state
108-
109-
110-
def _phantom_missing_cys(role: str, pos: int, chain_id_for_role: str) -> CysteineHit:
111-
"""Emit a placeholder row for a canonical Cys position the structure
112-
doesn't actually contain (resName "-" sentinel). Spec R23: missing
113-
canonical Cys is a developability red flag and needs to appear in
114-
the cys table even without a real residue to attach to."""
115-
return CysteineHit(
116-
chainId=chain_id_for_role,
117-
resSeq=pos,
118-
iCode="-",
119-
resName="-",
120-
chainRole=role,
121-
cysClass="disulfide_missing",
122-
bondingState="missing",
123-
sasa=None,
124-
rsasa=None,
125-
sidechainSasa=None,
126-
sidechainRsasa=None,
127-
partnerChainId="-",
128-
partnerResSeq=None,
129-
partnerIcode="-",
130-
)
102+
return "bonded" if bonded else "unbonded"
131103

132104

133105
def detect_cysteines(
@@ -144,69 +116,42 @@ def detect_cysteines(
144116
cys_records = _collect_cys_records(parsed)
145117
partner_of = _scan_disulfides(cys_records)
146118

147-
# Per-scheme canonical positions + the set we test each real Cys against.
148119
canonical_positions: dict[str, tuple[int, int]] = (
149120
CANONICAL_CYS_POSITIONS.get(numbering_scheme, {}) if numbering_scheme else {}
150121
)
151122
canonical_keys: set[tuple[str, int]] = {
152123
(role, pos) for role, positions in canonical_positions.items() for pos in positions
153124
}
154-
# Used below to skip phantom rows for canonical positions that ARE
155-
# filled , keyed by (role, res_seq).
156-
cys_by_role_pos: dict[tuple[str, int], int] = {}
157-
for idx, (chain_id, r, _ca, _sg) in enumerate(cys_records):
125+
# Tracks which canonical (role, res_seq) slots are actually filled, so
126+
# missing-Cys phantom rows below don't duplicate a real hit.
127+
filled_canonical: set[tuple[str, int]] = set()
128+
for chain_id, r, _ca, _sg in cys_records:
158129
role = role_of_chain(chain_id, heavy_chain_id, light_chain_id)
159130
if role is not None:
160-
cys_by_role_pos[(role, r.res_seq)] = idx
131+
filled_canonical.add((role, r.res_seq))
161132

162133
hits: list[CysteineHit] = []
163134
for idx, (chain_id, r, _ca, _sg) in enumerate(cys_records):
164135
key = (chain_id, f"{r.res_seq}{r.i_code}".strip())
165136
sasa_info = sasa_lookup.get(key, {})
166137
role = role_of_chain(chain_id, heavy_chain_id, light_chain_id)
167-
168-
partner_idx = partner_of.get(idx)
169-
if partner_idx is not None:
170-
p_chain, p_res, _, _ = cys_records[partner_idx]
171-
bonding_state = "bonded"
172-
partner_chain = p_chain
173-
partner_res_seq: Optional[int] = p_res.res_seq
174-
partner_icode = p_res.i_code or "-"
175-
else:
176-
bonding_state = "unbonded"
177-
partner_chain = "-"
178-
partner_res_seq = None
179-
partner_icode = "-"
180-
138+
bonded = partner_of.get(idx) is not None
181139
hits.append(
182140
CysteineHit(
183-
chainId=chain_id,
184-
resSeq=r.res_seq,
185-
iCode=r.i_code or "-",
186-
resName="CYS",
187-
chainRole=role or "-",
188141
cysClass=_classify_cys(
189-
r.res_seq, role, bonding_state, canonical_keys, canonical_positions
142+
r.res_seq, role, bonded, canonical_keys, canonical_positions
190143
),
191-
bondingState=bonding_state,
192-
sasa=sasa_info.get("sasa"),
193-
rsasa=sasa_info.get("rsasa"),
194-
sidechainSasa=sasa_info.get("sideChainSasa"),
195144
sidechainRsasa=sasa_info.get("sideChainRsasa"),
196-
partnerChainId=partner_chain,
197-
partnerResSeq=partner_res_seq,
198-
partnerIcode=partner_icode,
199145
)
200146
)
201147

202148
# Phantom rows for missing canonical Cys (R23).
203149
for role, (p1, p2) in canonical_positions.items():
204-
chain_id_for_role = heavy_chain_id if role == "H" else light_chain_id
205-
if not chain_id_for_role:
150+
if not (heavy_chain_id if role == "H" else light_chain_id):
206151
continue
207152
for pos in (p1, p2):
208-
if (role, pos) in cys_by_role_pos:
153+
if (role, pos) in filled_canonical:
209154
continue
210-
hits.append(_phantom_missing_cys(role, pos, chain_id_for_role))
155+
hits.append(CysteineHit(cysClass="disulfide_missing", sidechainRsasa=None))
211156

212157
return hits

software/liabilities-script/main.py

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
from metrics import compute_metrics
2424
from motifs import detect_motifs
2525
from scoring import compute_developability
26-
from structure import check_hallmark_tetrad, cross_check_ssbonds, parse_pdb
26+
from structure import check_hallmark_tetrad, parse_pdb
2727

2828

2929
# Heavy-atom Ala-X-Ala SASA references (R11), loaded from
@@ -253,9 +253,7 @@ def analyze_pdb(
253253
)
254254
mode = "TAP" if n_chains == 2 else "TNP"
255255

256-
# R21 SSBOND cross-check + R33 hallmark tetrad. Hallmark mismatch
257-
# surfaces in the TSV; SSBOND mismatches stay stderr-only.
258-
cross_check_ssbonds(parsed.ssbonds, cys_hits)
256+
# R33 hallmark tetrad re-check; mismatch surfaces as a TSV warning column.
259257
hallmark = check_hallmark_tetrad(
260258
parsed, numbering_scheme, heavy_chain_id, chain_count_mode=mode
261259
)

software/liabilities-script/structure.py

Lines changed: 3 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,8 @@
44
PLATFORMA CDR record extraction, FreeSASA wrapper output normalization
55
at the caller, B-factor read-through, the numbering-scheme constants
66
(R14: canonical Cys positions, CDR ranges, hallmark tetrad, CDRH3
7-
compactness anchors), `region_for` resolution, plus the spec defensive
8-
checks (R21 SSBOND cross-check, R33 hallmark tetrad re-check).
7+
compactness anchors), `region_for` resolution, plus the R33 hallmark
8+
tetrad defensive check.
99
"""
1010

1111
import re
@@ -46,21 +46,10 @@ def atom(self, name: str) -> Optional[Atom]:
4646
return None
4747

4848

49-
@dataclass
50-
class Ssbond:
51-
chain1: str
52-
res1: int
53-
i_code1: str
54-
chain2: str
55-
res2: int
56-
i_code2: str
57-
58-
5949
@dataclass
6050
class Parsed:
6151
chain_order: List[str] = field(default_factory=list)
6252
residues_by_chain: Dict[str, List[Residue]] = field(default_factory=dict)
63-
ssbonds: List[Ssbond] = field(default_factory=list)
6453
# Spec R10 , CDR ranges from `REMARK 99 PLATFORMA CDR*` records emitted
6554
# by the Structure Prediction block. Shape: {"H": {"CDR1": (start, end),
6655
# "CDR2": (...), "CDR3": (...)}, "L": {...}}. Empty when not present;
@@ -125,29 +114,6 @@ def parse_pdb(text: str) -> Parsed:
125114
out.chain_role_to_pdb_chain[role] = chain_start
126115
elif existing.upper() != chain_start.upper():
127116
out.chain_role_to_pdb_chain.pop(role, None)
128-
elif tag == "SSBOND":
129-
# SSBOND record fixed offsets (PDB v3.30):
130-
# col 15 = chainID1, 17-20 = resSeq1, 21 = iCode1
131-
# col 29 = chainID2, 31-34 = resSeq2, 35 = iCode2
132-
try:
133-
chain1 = line[15:16]
134-
res1 = int(line[17:21].strip())
135-
i_code1 = line[21:22].strip()
136-
chain2 = line[29:30]
137-
res2 = int(line[31:35].strip())
138-
i_code2 = line[35:36].strip()
139-
except ValueError:
140-
continue
141-
out.ssbonds.append(
142-
Ssbond(
143-
chain1=chain1,
144-
res1=res1,
145-
i_code1=i_code1,
146-
chain2=chain2,
147-
res2=res2,
148-
i_code2=i_code2,
149-
)
150-
)
151117
elif tag in ("ATOM", "HETATM"):
152118
if not in_first_model:
153119
continue
@@ -338,7 +304,7 @@ def role_of_chain(
338304

339305

340306
# ---------------------------------------------------------------------------
341-
# Section 3: defensive checks (R21 SSBOND cross-check, R33 hallmark tetrad)
307+
# Section 3: R33 hallmark tetrad defensive check
342308
# ---------------------------------------------------------------------------
343309

344310

@@ -375,81 +341,6 @@ def _match_score(observed_one_letter: list[Optional[str]], canonical: tuple) ->
375341
return score
376342

377343

378-
def _norm_icode(s: Optional[str]) -> str:
379-
"""Normalize iCode for comparison. Both the empty-string convention (PDB
380-
SSBOND records) and the `"-"` sentinel (cysteines.py emits this so PColumn
381-
String axes accept the value) collapse to empty-string here so the two
382-
sides of the cross-check actually meet."""
383-
if not s:
384-
return ""
385-
t = s.strip()
386-
if t in ("", "-"):
387-
return ""
388-
return t
389-
390-
391-
def _pair_key(chain_a, res_a, icode_a, chain_b, res_b, icode_b):
392-
"""Order-independent key for a Cys pair (sorted by chain+resSeq)."""
393-
a = (chain_a, res_a, _norm_icode(icode_a))
394-
b = (chain_b, res_b, _norm_icode(icode_b))
395-
return tuple(sorted([a, b]))
396-
397-
398-
def cross_check_ssbonds(ssbonds, cys_hits) -> dict:
399-
"""Compare PDB SSBOND records with geometry-detected disulfide bonds.
400-
401-
Returns a dict suitable for JSON emission:
402-
{
403-
"headerBondCount": int,
404-
"geometryBondCount": int,
405-
"matched": int,
406-
"headerOnly": [...], # SSBOND declared, geometry rejected
407-
"geometryOnly": [...], # geometry found, no SSBOND
408-
}
409-
"""
410-
header_pairs: set = set()
411-
for s in ssbonds:
412-
header_pairs.add(
413-
_pair_key(s.chain1, s.res1, s.i_code1, s.chain2, s.res2, s.i_code2)
414-
)
415-
416-
# Each disulfide is represented twice in cys_hits (once per partner). De-dup.
417-
geom_pairs: set = set()
418-
for h in cys_hits:
419-
if h.cysClass not in ("disulfide", "disulfide_broken"):
420-
continue
421-
if h.partnerChainId is None or h.partnerResSeq is None:
422-
continue
423-
if h.cysClass == "disulfide_broken":
424-
# geometry rejected this; it stays under headerOnly, not geometryOnly.
425-
continue
426-
geom_pairs.add(
427-
_pair_key(
428-
h.chainId, h.resSeq, h.iCode,
429-
h.partnerChainId, h.partnerResSeq, h.partnerIcode,
430-
)
431-
)
432-
433-
matched = header_pairs & geom_pairs
434-
header_only = header_pairs - geom_pairs
435-
geometry_only = geom_pairs - header_pairs
436-
437-
def _fmt(pair):
438-
a, b = pair
439-
return {
440-
"chain1": a[0], "resSeq1": a[1], "iCode1": a[2] or "",
441-
"chain2": b[0], "resSeq2": b[1], "iCode2": b[2] or "",
442-
}
443-
444-
return {
445-
"headerBondCount": len(header_pairs),
446-
"geometryBondCount": len(geom_pairs),
447-
"matched": len(matched),
448-
"headerOnly": [_fmt(p) for p in sorted(header_only)],
449-
"geometryOnly": [_fmt(p) for p in sorted(geometry_only)],
450-
}
451-
452-
453344
def check_hallmark_tetrad(
454345
parsed,
455346
numbering_scheme: Optional[str],

software/tests/test_scoring.py

Lines changed: 1 addition & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -36,22 +36,7 @@ def _motif(
3636

3737

3838
def _cys(cys_class="cys_extra", sidechain_rsasa=0.5) -> CysteineHit:
39-
return CysteineHit(
40-
chainId="H",
41-
resSeq=100,
42-
iCode="",
43-
resName="CYS",
44-
chainRole="H",
45-
cysClass=cys_class,
46-
bondingState="unbonded",
47-
sasa=80.0,
48-
rsasa=0.5,
49-
sidechainSasa=40.0,
50-
sidechainRsasa=sidechain_rsasa,
51-
partnerChainId="",
52-
partnerResSeq=None,
53-
partnerIcode="",
54-
)
39+
return CysteineHit(cysClass=cys_class, sidechainRsasa=sidechain_rsasa)
5540

5641

5742
class TestComputeFlagsFv:

0 commit comments

Comments
 (0)