""" ORF extraction and translation. The glue between TIS prediction and MHC binding prediction: given a transcript sequence and a predicted translation initiation site, read off the open reading frame, translate it, and chop the result into the peptide lengths that MHC class I molecules actually present. Nothing here is novel or uncertain -- it is the standard genetic code and some bookkeeping. It is a separate module precisely because it is the one part of the pipeline that can be tested exactly. """ import re from dataclasses import dataclass, field from typing import Iterable # Standard genetic code (NCBI translation table 1). CODON_TABLE = { "TTT": "F", "TTC": "F", "TTA": "L", "TTG": "L", "CTT": "L", "CTC": "L", "CTA": "L", "CTG": "L", "ATT": "I", "ATC": "I", "ATA": "I", "ATG": "M", "GTT": "V", "GTC": "V", "GTA": "V", "GTG": "V", "TCT": "S", "TCC": "S", "TCA": "S", "TCG": "S", "CCT": "P", "CCC": "P", "CCA": "P", "CCG": "P", "ACT": "T", "ACC": "T", "ACA": "T", "ACG": "T", "GCT": "A", "GCC": "A", "GCA": "A", "GCG": "A", "TAT": "Y", "TAC": "Y", "TAA": "*", "TAG": "*", "CAT": "H", "CAC": "H", "CAA": "Q", "CAG": "Q", "AAT": "N", "AAC": "N", "AAA": "K", "AAG": "K", "GAT": "D", "GAC": "D", "GAA": "E", "GAG": "E", "TGT": "C", "TGC": "C", "TGA": "*", "TGG": "W", "CGT": "R", "CGC": "R", "CGA": "R", "CGG": "R", "AGT": "S", "AGC": "S", "AGA": "R", "AGG": "R", "GGT": "G", "GGC": "G", "GGA": "G", "GGG": "G", } STOP_CODONS = {"TAA", "TAG", "TGA"} # MHC class I peptides are overwhelmingly 8-11mers, with 9 the mode. MHC1_PEPTIDE_LENGTHS = (8, 9, 10, 11) def translate(nt: str, stop_at_stop: bool = True) -> str: """Translate a nucleotide sequence in frame 0. Codons containing N (or anything not in the table) become 'X'. A stop codon is rendered '*' and, by default, ends translation. """ nt = nt.upper().replace("U", "T") out = [] for i in range(0, len(nt) - len(nt) % 3, 3): codon = nt[i : i + 3] aa = CODON_TABLE.get(codon, "X") if aa == "*": if stop_at_stop: break out.append(aa) continue out.append(aa) return "".join(out) @dataclass class ORF: """One open reading frame read off from a predicted initiation site.""" transcript: str tis_pos: int # 0-based position of the first base of the start codon start_codon: str stop_pos: int | None # 0-based position of first base of the stop codon nt_seq: str # start codon through stop codon (stop included if found) protein: str # translated, stop excluded has_stop: bool # False => ORF runs off the end of the transcript orf_type: str = "unknown" tis_score: float | None = None @property def aa_len(self) -> int: return len(self.protein) @property def nt_len(self) -> int: return len(self.nt_seq) @property def is_microprotein(self) -> bool: """The field's usual, arbitrary cutoff: under 100 aa.""" return self.aa_len < 100 def extract_orf( seq: str, tis_pos: int, transcript: str = "query", tis_score: float | None = None, ) -> ORF: """Read the ORF beginning at `tis_pos` until an in-frame stop or the end. `tis_pos` is 0-based and points at the first base of the start codon. """ seq = seq.upper().replace("U", "T") if not 0 <= tis_pos < len(seq): raise ValueError(f"TIS position {tis_pos} outside transcript of {len(seq)} nt") start_codon = seq[tis_pos : tis_pos + 3] stop_pos = None for i in range(tis_pos, len(seq) - 2, 3): if seq[i : i + 3] in STOP_CODONS: stop_pos = i break if stop_pos is None: # Runs off the 3' end: keep whole codons only. usable = (len(seq) - tis_pos) // 3 * 3 nt_seq = seq[tis_pos : tis_pos + usable] has_stop = False else: nt_seq = seq[tis_pos : stop_pos + 3] has_stop = True return ORF( transcript=transcript, tis_pos=tis_pos, start_codon=start_codon, stop_pos=stop_pos, nt_seq=nt_seq, protein=translate(nt_seq), has_stop=has_stop, tis_score=tis_score, ) def classify_orf( orf: ORF, cds_start: int | None = None, cds_end: int | None = None, ) -> str: """Name the ORF relative to the transcript's annotated CDS, if known. `cds_start` / `cds_end` are 0-based; `cds_end` is exclusive and includes the stop codon. With no annotation we can only say whether it looks canonical. Types follow the usual translatome vocabulary: CDS -- matches the annotated start uORF -- entirely within the 5'UTR uoORF -- starts in the 5'UTR, overlaps the CDS out of frame N-ext -- starts upstream, in frame with the CDS (N-terminal extension) intORF -- starts inside the CDS, out of frame N-trunc -- starts inside the CDS, in frame (N-terminal truncation) dORF -- entirely within the 3'UTR ncORF -- no annotated CDS (e.g. lncRNA, processed pseudogene) """ if cds_start is None or cds_end is None: return "ncORF" tis = orf.tis_pos orf_end = (orf.stop_pos + 3) if orf.stop_pos is not None else len(orf.nt_seq) + tis in_frame = (tis - cds_start) % 3 == 0 if tis == cds_start: return "CDS" if tis < cds_start: if orf_end <= cds_start: return "uORF" return "N-ext" if in_frame else "uoORF" if tis < cds_end: return "N-trunc" if in_frame else "intORF" return "dORF" def tile_peptides( protein: str, lengths: Iterable[int] = MHC1_PEPTIDE_LENGTHS, drop_initiator_met: bool = True, ) -> list[str]: """Every subsequence of the given lengths -- the candidate MHC-I ligands. Methionine aminopeptidase usually removes the initiator Met, so by default we also tile the Met-less form. Peptides containing X are dropped because binding predictors cannot score them. """ forms = {protein} if drop_initiator_met and protein.startswith("M") and len(protein) > 1: forms.add(protein[1:]) peptides: set[str] = set() for form in forms: for n in lengths: if len(form) < n: continue for i in range(len(form) - n + 1): pep = form[i : i + n] if "X" not in pep and "*" not in pep: peptides.add(pep) return sorted(peptides) @dataclass class ScoredORF: """An ORF plus whatever the downstream predictors said about it.""" orf: ORF peptides: list[str] = field(default_factory=list) best_peptide: str | None = None best_allele: str | None = None best_affinity_nm: float | None = None best_percentile: float | None = None # --- Interop with TIS Transformer ------------------------------------------- # Its FASTA-path CSV reports TIS_pos 1-based (the source does `idx + 1`), while # everything in this module is 0-based. Converting in one named place rather # than inline is deliberate: an off-by-one here would shift every reading frame # and the output would still look plausible. def tis_pos_to_zero_based(tis_pos: int) -> int: """Convert TIS Transformer's 1-based TIS_pos to a 0-based offset.""" return int(tis_pos) - 1 # --- Annotation from FASTA headers ------------------------------------------ # Headers written by the demo fetcher carry the CDS layout, e.g. # >INS_ENST00000381330.5 ... len=465nt 5UTR=59nt CDS=333nt 3UTR=73nt # Parsing it lets us name ORF types (uORF / N-ext / intORF ...) without a GTF. # Users can annotate their own sequences the same way; absent the tags we fall # back to "ncORF", which is honest rather than wrong. _UTR_RE = re.compile(r"5UTR=(\d+)\s*nt", re.I) _CDS_RE = re.compile(r"\bCDS=(\d+)\s*nt", re.I) def parse_cds_annotation(header: str) -> tuple[int | None, int | None]: """Pull (cds_start, cds_end) 0-based/exclusive out of a FASTA header.""" utr = _UTR_RE.search(header) cds = _CDS_RE.search(header) if not utr or not cds: return None, None start = int(utr.group(1)) return start, start + int(cds.group(1))