#!/usr/bin/env python3 """Convert protein C-alpha coordinates from a PDB file to an Arduino header. The generated model is a connected 3-D backbone polyline. It is intended for small monochrome displays such as a 128x64 SSD1306 OLED. """ from __future__ import annotations import argparse import gzip from pathlib import Path import shlex AMINO_ACIDS = { "ALA", "ARG", "ASN", "ASP", "CYS", "GLN", "GLU", "GLY", "HIS", "ILE", "LEU", "LYS", "MET", "PHE", "PRO", "SER", "THR", "TRP", "TYR", "VAL", "MSE", } def read_pdb_ca_atoms(path: Path, chain: str | None) -> list[tuple[float, float, float]]: """Read one alpha-carbon point per residue from the selected chain. PDB is a fixed-width format, so coordinates are read by column position. Only ATOM records and common amino-acid residues are accepted. Alternate locations are reduced to the blank location or alternate location A. """ points: list[tuple[float, float, float]] = [] seen_residues: set[tuple[str, str, str, str]] = set() selected_chain = chain with open_text(path) as handle: for line in handle: if not line.startswith("ATOM ") or len(line) < 54: continue atom_name = line[12:16].strip() residue_name = line[17:20].strip().upper() alternate = line[16].strip() current_chain = line[21].strip() or "_" if atom_name != "CA" or residue_name not in AMINO_ACIDS: continue if alternate not in ("", "A"): continue if selected_chain is None: selected_chain = current_chain if current_chain != selected_chain: continue residue_number = line[22:26].strip() insertion_code = line[26].strip() residue_key = (current_chain, residue_number, insertion_code, residue_name) if residue_key in seen_residues: continue seen_residues.add(residue_key) try: x = float(line[30:38]) y = float(line[38:46]) z = float(line[46:54]) except ValueError: continue points.append((x, y, z)) return points def open_text(path: Path): """Open plain-text or gzip-compressed structure files.""" if path.suffix.lower() == ".gz": return gzip.open(path, "rt", encoding="utf-8", errors="replace") return path.open("r", encoding="utf-8", errors="replace") def read_mmcif_ca_atoms(path: Path, chain: str | None) -> list[tuple[float, float, float]]: """Read alpha-carbon points from a common PDBx/mmCIF atom_site loop. This intentionally handles the atom_site loop needed by this converter, rather than attempting to interpret every possible mmCIF data category. """ with open_text(path) as handle: lines = iter(handle) for line in lines: if line.strip().lower() != "loop_": continue headers: list[str] = [] pending = None for candidate in lines: stripped = candidate.strip() if stripped.startswith("_"): headers.append(stripped.split()[0]) continue pending = candidate break if not headers or not any(header.startswith("_atom_site.") for header in headers): continue rows: list[list[str]] = [] if pending is not None and pending.strip() and not pending.lstrip().startswith("#"): rows.append(shlex.split(pending, comments=False, posix=True)) for row_line in lines: stripped = row_line.strip() if not stripped: continue if stripped.startswith("#"): break if stripped.startswith("_") or stripped.lower() == "loop_": break rows.append(shlex.split(row_line, comments=False, posix=True)) required = { "group": "_atom_site.group_PDB", "atom": "_atom_site.auth_atom_id", "comp": "_atom_site.auth_comp_id", "chain": "_atom_site.auth_asym_id", "seq": "_atom_site.auth_seq_id", "x": "_atom_site.Cartn_x", "y": "_atom_site.Cartn_y", "z": "_atom_site.Cartn_z", } # Some files omit auth_* fields, so fall back to label_* fields. alternatives = { "atom": "_atom_site.label_atom_id", "comp": "_atom_site.label_comp_id", "chain": "_atom_site.label_asym_id", "seq": "_atom_site.label_seq_id", } indices: dict[str, int] = {} for key, header in required.items(): chosen = header if header in headers else alternatives.get(key) if chosen not in headers: indices = {} break indices[key] = headers.index(chosen) if not indices: continue alt_index = headers.index("_atom_site.label_alt_id") if "_atom_site.label_alt_id" in headers else None ins_index = headers.index("_atom_site.pdbx_PDB_ins_code") if "_atom_site.pdbx_PDB_ins_code" in headers else None points: list[tuple[float, float, float]] = [] seen_residues: set[tuple[str, str, str, str]] = set() selected_chain = chain for row in rows: if len(row) < len(headers): continue if row[indices["group"]].upper() != "ATOM": continue if row[indices["atom"]].strip().upper() != "CA": continue if row[indices["comp"]].strip().upper() not in AMINO_ACIDS: continue if alt_index is not None and row[alt_index] not in (".", "?", "A"): continue current_chain = row[indices["chain"]] if current_chain in (".", "?"): current_chain = "_" if selected_chain is None: selected_chain = current_chain if current_chain != selected_chain: continue insertion = row[ins_index] if ins_index is not None else "" if insertion in (".", "?"): insertion = "" residue_key = (current_chain, row[indices["seq"]], insertion, row[indices["comp"]]) if residue_key in seen_residues: continue try: point = ( float(row[indices["x"]]), float(row[indices["y"]]), float(row[indices["z"]]), ) except ValueError: continue seen_residues.add(residue_key) points.append(point) return points return [] def read_ca_atoms(path: Path, chain: str | None) -> list[tuple[float, float, float]]: name = path.name.lower() if name.endswith((".cif", ".mmcif", ".cif.gz", ".mmcif.gz")): return read_mmcif_ca_atoms(path, chain) return read_pdb_ca_atoms(path, chain) def downsample(points: list[tuple[float, float, float]], maximum: int) -> list[tuple[float, float, float]]: if maximum < 2: raise ValueError("maximum must be at least 2") if len(points) <= maximum: return points # Keep the first and last points and choose evenly spaced intermediate ones. result = [] for index in range(maximum): source_index = round(index * (len(points) - 1) / (maximum - 1)) result.append(points[source_index]) return result def normalize(points: list[tuple[float, float, float]], size: float) -> list[tuple[float, float, float]]: center = tuple(sum(point[i] for point in points) / len(points) for i in range(3)) centered = [tuple(point[i] - center[i] for i in range(3)) for point in points] largest = max(max(abs(value) for value in point) for point in centered) if largest == 0: return centered scale = size / largest return [tuple(value * scale for value in point) for point in centered] def write_header(path: Path, points: list[tuple[float, float, float]], source: Path, chain: str | None) -> None: with path.open("w", encoding="utf-8") as handle: handle.write("#ifndef PROTEIN_MODEL_H\n#define PROTEIN_MODEL_H\n\n") handle.write("// Generated by pdb_to_oled.py\n") handle.write(f"// Source: {source.name}; chain: {chain or 'first protein chain'}\n\n") handle.write("struct Point3D { float x; float y; float z; };\n\n") handle.write(f"const int NUM_POINTS = {len(points)};\n") handle.write("const Point3D protein[NUM_POINTS] = {\n") for x, y, z in points: handle.write(f" {{{x:8.3f}f, {y:8.3f}f, {z:8.3f}f}},\n") handle.write("};\n\n#endif\n") def main() -> None: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("pdb", type=Path, help="Input PDB file") parser.add_argument("-o", "--output", type=Path, default=Path("protein_model.h")) parser.add_argument("--chain", help="Chain ID; default is the first protein chain") parser.add_argument("--max-points", type=int, default=180, help="Maximum C-alpha points to export (default: 180)") parser.add_argument("--size", type=float, default=27.0, help="Largest centered coordinate after scaling (default: 27)") args = parser.parse_args() points = read_ca_atoms(args.pdb, args.chain) if len(points) < 2: raise SystemExit("No usable protein C-alpha atoms found. Check the PDB and chain ID.") points = downsample(points, args.max_points) points = normalize(points, args.size) write_header(args.output, points, args.pdb, args.chain) print(f"Wrote {len(points)} points to {args.output}") if __name__ == "__main__": main()