diff --git a/cw/align.py b/cw/align.py index 2fe77e1..706b9cb 100644 --- a/cw/align.py +++ b/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) diff --git a/cw/metadata.py b/cw/metadata.py index 17b3fa1..a259f8e 100644 --- a/cw/metadata.py +++ b/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 ``''``. + + 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 "" + 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 "" + + +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 ``''``. + + ``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 "", "" + 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 "" + + ph = _first("pH") + if ph == "": + 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, } diff --git a/scripts/build_metadata.py b/scripts/build_metadata.py index 34df79e..b5b0c5f 100644 --- a/scripts/build_metadata.py +++ b/scripts/build_metadata.py @@ -3,8 +3,9 @@ Usage: uv run scripts/build_metadata.py [-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//metadata.csv by default. """