feat(metadata): sync up core scripts

cw/align.py

@@ -16,12 +16,32 @@
 
 from cw.io import load_protein, write_transformed_cif
 
+GAP_PENALTY = (-10, -1)
+
+
+def paired_alignment_trace(
+    seq_a: bseq.ProteinSequence, seq_b: bseq.ProteinSequence
+) -> np.ndarray:
+    """Global BLOSUM62 alignment of two protein sequences.
+
+    Returns the (n, 2) array of trace rows where neither sequence has a gap;
+    column 0 indexes seq_a, column 1 indexes seq_b. Terminal gaps are free, so
+    ragged termini (missing density) are not penalized.
+    """
+    matrix = bseq_align.SubstitutionMatrix.std_protein_matrix()
+    alignment = bseq_align.align_optimal(  # type: ignore
+        seq_a, seq_b, matrix, gap_penalty=GAP_PENALTY, terminal_penalty=False
+    )[0]
+    trace = alignment.trace
+    return trace[(trace[:, 0] != -1) & (trace[:, 1] != -1)]
+
 
 def get_ca_coords_and_sequence(protein: struc.AtomArray) -> tuple[np.ndarray, str]:
     """Extract Cα positions and one-letter sequence in matching order.
 
     Ported verbatim from align_pdbs.py::get_ca_coords_and_sequence.
     Non-standard residues without a known one-letter code are skipped.
+    All protein chains are concatenated (deliberate: superposition uses every Cα).
     """
     ca = protein[protein.atom_name == "CA"]
     coords, seq_chars = [], []
@@ -46,15 +66,9 @@ def get_paired_ca_positions(
     mob_coords, mob_seq_str = get_ca_coords_and_sequence(mobile)
     ref_coords, ref_seq_str = get_ca_coords_and_sequence(ref)
 
-    mob_seq = bseq.ProteinSequence(mob_seq_str)
-    ref_seq = bseq.ProteinSequence(ref_seq_str)
-
-    matrix = bseq_align.SubstitutionMatrix.std_protein_matrix()
-    alignments = bseq_align.align_optimal(  # type: ignore
-        mob_seq, ref_seq, matrix, gap_penalty=(-10, -1), terminal_penalty=False
+    paired = paired_alignment_trace(
+        bseq.ProteinSequence(mob_seq_str), bseq.ProteinSequence(ref_seq_str)
     )
-    trace = alignments[0].trace
-    paired = trace[(trace[:, 0] != -1) & (trace[:, 1] != -1)]
     return mob_coords[paired[:, 0]], ref_coords[paired[:, 1]], len(paired)
 
 

cw/metadata.py

@@ -167,6 +167,65 @@ def _first(key: str) -> float | str:
     return _first("ls_R_factor_R_work"), _first("ls_R_factor_R_free")
 
 
+def resolve_diffrn_temp(entry: dict | None) -> float | str:
+    """Lowest ``_diffrn.ambient_temp`` (K) on an RCSB entry, or ``'<missing>'``.
+
+    The re-refined local CIFs mostly drop this category, so it comes from the same
+    entry JSON already fetched for experiment_condition. Lowest wins when an entry
+    lists several datasets: it is the temperature that determines whether the
+    solvent is frozen, which is what a water analysis cares about.
+    """
+    if entry is None:
+        return "<missing>"
+    temps = []
+    for block in entry.get("diffrn") or []:
+        value = block.get("ambient_temp")
+        if value is None:
+            continue
+        try:
+            temps.append(float(value))
+        except (TypeError, ValueError):
+            continue
+    return min(temps) if temps else "<missing>"
+
+
+def resolve_crystal_grow(entry: dict | None) -> tuple[float | str, float | str]:
+    """``(ph, crystal_grow_temp)`` from an RCSB entry's ``exptl_crystal_grow`` array.
+
+    These are crystallisation conditions — the same category ``experiment_condition``
+    is parsed from — and are distinct from ``diffrn_temp``, which is the temperature
+    the diffraction data were collected at. Returns the first parseable value across
+    blocks for each field, else ``'<missing>'``.
+
+    ``ph`` prefers the scalar ``pH``; when an entry gives only ``pdbx_pH_range`` the
+    raw range string is passed through verbatim ('4.4-4.8'). Both fields are optional
+    and roughly half of a typical cohort supplies neither, so the column holds floats
+    and range strings side by side — read it as object dtype, not float.
+    """
+    if entry is None:
+        return "<missing>", "<missing>"
+    blocks = entry.get("exptl_crystal_grow") or []
+
+    def _first(key: str) -> float | str:
+        for block in blocks:
+            value = block.get(key)
+            if value is not None:
+                try:
+                    return float(value)
+                except (TypeError, ValueError):
+                    pass
+        return "<missing>"
+
+    ph = _first("pH")
+    if ph == "<missing>":
+        for block in blocks:
+            ph_range = str(block.get("pdbx_pH_range") or "").strip()
+            if ph_range:
+                ph = ph_range
+                break
+    return ph, _first("temp")
+
+
 def metadata_row(cif_path: Path) -> dict:
     """Extract one metadata CSV row from a local mmCIF file + RCSB API.
 
@@ -175,7 +234,8 @@ def metadata_row(cif_path: Path) -> dict:
       - gemmi.cif.read       → r_work, r_free (local re-refined CIF has these)
                                ligand_names (from _pdbx_entity_nonpoly.comp_id loop)
       - RCSB Data API        → experiment_condition, starting_model,
-                               deposited_r_work, deposited_r_free
+                               deposited_r_work, deposited_r_free, ph,
+                               crystal_grow_temp, diffrn_temp
                                (absent from re-refined local CIFs)
     """
     pdb_id = cif_path.stem.removesuffix("_final")
@@ -215,6 +275,7 @@ def metadata_row(cif_path: Path) -> dict:
     entry = _fetch_rcsb_entry(pdb_id)
     _codes, _status, starting_model = resolve_starting_model(entry, pdb_id)
     deposited_r_work, deposited_r_free = resolve_deposited_r_factors(entry)
+    ph, crystal_grow_temp = resolve_crystal_grow(entry)
     if entry is not None:
         grow_blocks = entry.get("exptl_crystal_grow") or []
         grow_details = [
@@ -242,5 +303,8 @@ def metadata_row(cif_path: Path) -> dict:
         "num_water": num_water,
         "ligand_names": ligand_names,
         "experiment_condition": experiment_condition,
+        "ph": ph,
+        "crystal_grow_temp": crystal_grow_temp,
+        "diffrn_temp": resolve_diffrn_temp(entry),
         "starting_model": starting_model,
     }

scripts/build_metadata.py

@@ -3,8 +3,9 @@
 Usage:
     uv run scripts/build_metadata.py <cohort.txt> [-o metadata.csv]
 
-Fetches experiment_condition and starting_model from the RCSB Data API;
-all other fields come from local mmCIF files.
+Fetches experiment_condition, starting_model, the deposited R-factors, and the
+crystallisation/data-collection conditions (ph, crystal_grow_temp, diffrn_temp)
+from the RCSB Data API; all other fields come from local mmCIF files.
 Writes to data/<cohort_id>/metadata.csv by default.
 """