"""Kyte-Doolittle hydropathy profiles and a proteome check of the 19-residue / 1.6 rule. Reproduces the numbers in content/guides/kyte-doolittle-hydropathy.mdx. python analyze.py [--cache DIR] Downloads (into --cache, default a temp folder, never into this directory): * UniProtKB entries P02724 (glycophorin A), P02945 (bacteriorhodopsin) and P02144 (myoglobin) as JSON; * reviewed human proteome UP000005640 as TSV with transmembrane, intramembrane, signal peptide and keyword columns. Writes next to this script: * profiles.csv window-19 and window-9 profiles of the three examples * window-sweep.csv TM-helix detection and false-positive rate per window * summary.json everything quoted in the guide, with the UniProt release Scale: Kyte J, Doolittle RF (1982) J Mol Biol 157:105-132, Table 2. Profile convention: mean of a full window, plotted at its centre residue (Kyte and Doolittle 1982; ExPASy ProtScale). Partial windows are not plotted. """ from __future__ import annotations import argparse import csv import gzip import json import re import tempfile import urllib.request from pathlib import Path KD = { 'I': 4.5, 'V': 4.2, 'L': 3.8, 'F': 2.8, 'C': 2.5, 'M': 1.9, 'A': 1.8, 'G': -0.4, 'T': -0.7, 'S': -0.8, 'W': -0.9, 'Y': -1.3, 'P': -1.6, 'H': -3.2, 'E': -3.5, 'Q': -3.5, 'D': -3.5, 'N': -3.5, 'K': -3.9, 'R': -4.5, } THRESHOLD = 1.6 WINDOWS = [7, 9, 11, 13, 15, 17, 19, 21] EXAMPLES = {'P02724': 'glycophorin A', 'P02945': 'bacteriorhodopsin', 'P02144': 'myoglobin'} HERE = Path(__file__).resolve().parent PROTEOME_URL = ( 'https://rest.uniprot.org/uniprotkb/stream?compressed=true&format=tsv' '&fields=accession,length,sequence,ft_transmem,ft_intramem,ft_signal,keywordid' '&query=%28proteome%3AUP000005640%29%20AND%20%28reviewed%3Atrue%29' ) # Evidence codes that rest on experiment or a deposited structure, not on # sequence-based prediction (ECO:0000255) or similarity (ECO:0000250). EXPERIMENTAL = {'ECO:0000269', 'ECO:0007744'} FEATURE = re.compile(r'(TRANSMEM|SIGNAL) [?]?(\d+)(?:; /note="([^"]*)")?(?:; /evidence="([^"]*)")?') def fetch(url: str, dest: Path) -> tuple[bytes, str]: if not dest.exists(): with urllib.request.urlopen(url) as response: dest.write_bytes(response.read()) dest.with_suffix('.release').write_text(response.headers.get('X-UniProt-Release', '')) release = dest.with_suffix('.release') return dest.read_bytes(), release.read_text() if release.exists() else '' def profile(seq: str, window: int) -> list[tuple[int, float]]: """(1-based centre position, window mean) for every full window.""" values = [KD[a] for a in seq] half = window // 2 total = sum(values[:window]) out = [(half + 1, total / window)] for i in range(window, len(values)): total += values[i] - values[i - window] out.append((i - half + 1, total / window)) return out def gravy(seq: str) -> float: return sum(KD[a] for a in seq) / len(seq) def runs_above(points: list[tuple[int, float]], threshold: float) -> list[list[int]]: """Contiguous runs of centre positions whose mean exceeds the threshold.""" runs: list[list[int]] = [] for pos, value in points: if value > threshold: if runs and runs[-1][1] == pos - 1: runs[-1][1] = pos else: runs.append([pos, pos]) return runs def parse_features(text: str, kind: str) -> list[tuple[int, int, str, set[str]]]: out = [] for m in FEATURE.finditer(text or ''): if m.group(1) != kind: continue out.append((int(m.group(2)), int(m.group(3)), (m.group(4) or '').split(';')[0], set(re.findall(r'ECO:\d+', m.group(5) or '')))) return out def examples(cache: Path) -> tuple[dict, list[dict]]: summary, rows = {}, [] for acc, name in EXAMPLES.items(): raw, release = fetch(f'https://rest.uniprot.org/uniprotkb/{acc}.json', cache / f'{acc}.json') entry = json.loads(raw) seq = entry['sequence']['value'] feats = { kind: [ (f['location']['start']['value'], f['location']['end']['value']) for f in entry.get('features', []) if f['type'] == kind ] for kind in ('Transmembrane', 'Signal', 'Chain', 'Initiator methionine') } p19, p9 = profile(seq, 19), profile(seq, 9) peak = max(p19, key=lambda t: t[1]) summary[acc] = { 'name': name, 'uniprotRelease': release, 'length': len(seq), 'gravy': round(gravy(seq), 3), 'features': feats, 'window19': { 'max': round(peak[1], 2), 'maxAt': peak[0], 'min': round(min(v for _, v in p19), 2), 'runsAbove1.6': runs_above(p19, THRESHOLD), }, 'window9': {'max': round(max(v for _, v in p9), 2), 'min': round(min(v for _, v in p9), 2)}, } for (pos, v19) in p19: rows.append({'protein': name, 'accession': acc, 'window': 19, 'position': pos, 'residue': seq[pos - 1], 'hydropathy': round(v19, 3)}) for (pos, v9) in p9: rows.append({'protein': name, 'accession': acc, 'window': 9, 'position': pos, 'residue': seq[pos - 1], 'hydropathy': round(v9, 3)}) return summary, rows def proteome(cache: Path) -> tuple[dict, list[dict]]: raw, release = fetch(PROTEOME_URL, cache / 'human-reviewed.tsv.gz') lines = gzip.decompress(raw).decode().splitlines()[1:] tm_segments: list[tuple[str, int, int]] = [] soluble: list[str] = [] secreted: list[tuple[str, int]] = [] skipped = 0 proteins = {} for line in lines: acc, _length, seq, transmem, intramem, signal, keywords = (line.split('\t') + [''] * 7)[:7] if set(seq) - KD.keys(): skipped += 1 continue proteins[acc] = seq tms = parse_features(transmem, 'TRANSMEM') sps = parse_features(signal, 'SIGNAL') kws = set(keywords.split(';')) if keywords else set() kws = {k.strip() for k in kws} membrane_kw = bool(kws & {'KW-0472', 'KW-0812'}) for start, end, note, ev in tms: if note == 'Helical' and ev & EXPERIMENTAL: tm_segments.append((acc, start, end)) if not tms and not intramem and not sps and not membrane_kw: soluble.append(acc) if sps and not tms and not intramem and not membrane_kw: secreted.append((acc, sps[0][1])) sweep = [] for window in WINDOWS: profiles: dict[str, dict[int, float]] = {} def prof(acc: str) -> dict[int, float]: if acc not in profiles: profiles[acc] = dict(profile(proteins[acc], window)) if len(proteins[acc]) >= window else {} return profiles[acc] hit = 0 for acc, start, end in tm_segments: p = prof(acc) if any(p.get(pos, -9) > THRESHOLD for pos in range(start, end + 1)): hit += 1 fp = sum(1 for acc in soluble if any(v > THRESHOLD for v in prof(acc).values())) sp = sum(1 for acc, sp_end in secreted if any(v > THRESHOLD for pos, v in prof(acc).items() if pos <= sp_end)) sweep.append({ 'window': window, 'tmHelices': len(tm_segments), 'tmHelicesDetected': hit, 'tmHelixDetectionPct': round(100 * hit / len(tm_segments), 1), 'solubleProteins': len(soluble), 'solubleWithWindowAbove1.6': fp, 'solubleFalsePositivePct': round(100 * fp / len(soluble), 1), 'secretedProteins': len(secreted), 'signalPeptideAbove1.6': sp, 'signalPeptideAbove1.6Pct': round(100 * sp / len(secreted), 1), }) def gravy_stats(accs: list[str]) -> dict: values = sorted(gravy(proteins[a]) for a in accs) pick = lambda q: round(values[min(len(values) - 1, int(q * len(values)))], 3) return { 'n': len(values), 'p10': pick(0.1), 'median': pick(0.5), 'p90': pick(0.9), 'belowZeroPct': round(100 * sum(v < 0 for v in values) / len(values), 1), } tm_proteins = sorted({a for a, _, _ in tm_segments}) summary = { 'uniprotRelease': release, 'gravy': { 'experimentalTmProteins': gravy_stats(tm_proteins), 'solubleProteins': gravy_stats(soluble), 'secretedSignalPeptideProteins': gravy_stats([a for a, _ in secreted]), }, 'reviewedHumanEntries': len(lines), 'skippedNonstandardResidues': skipped, 'experimentalHelicalTmSegments': len(tm_segments), 'experimentalTmProteins': len(tm_proteins), 'solubleProteins': len(soluble), 'secretedSignalPeptideProteins': len(secreted), 'thresholdRule': 'window mean > 1.6', 'definitions': { 'experimentalHelicalTmSegments': 'UniProt TRANSMEM features with note Helical and evidence ECO:0000269 or ECO:0007744', 'solubleProteins': 'no TRANSMEM, INTRAMEM or SIGNAL feature and no Membrane (KW-0472) or Transmembrane (KW-0812) keyword', 'secretedSignalPeptideProteins': 'SIGNAL feature, no TRANSMEM or INTRAMEM feature, no Membrane or Transmembrane keyword; hit = a window centred inside the signal peptide exceeds 1.6', 'tmHelixDetected': 'at least one window centred inside the annotated segment exceeds 1.6', }, } return summary, sweep def main() -> None: parser = argparse.ArgumentParser() parser.add_argument('--cache', type=Path, default=Path(tempfile.gettempdir()) / 'kd-hydropathy-cache') args = parser.parse_args() args.cache.mkdir(parents=True, exist_ok=True) example_summary, rows = examples(args.cache) proteome_summary, sweep = proteome(args.cache) with open(HERE / 'profiles.csv', 'w', newline='', encoding='utf-8') as fh: writer = csv.DictWriter(fh, fieldnames=list(rows[0])) writer.writeheader() writer.writerows(rows) with open(HERE / 'window-sweep.csv', 'w', newline='', encoding='utf-8') as fh: writer = csv.DictWriter(fh, fieldnames=list(sweep[0])) writer.writeheader() writer.writerows(sweep) (HERE / 'summary.json').write_text( json.dumps({'scale': KD, 'examples': example_summary, 'proteome': proteome_summary, 'windowSweep': sweep}, indent=2) + '\n', encoding='utf-8', ) print(json.dumps({'examples': example_summary, 'proteome': proteome_summary}, indent=1)) for row in sweep: print(row) if __name__ == '__main__': main()