Spaces:
Running on Zero
Running on Zero
| """ | |
| 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) | |
| 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 | |
| def aa_len(self) -> int: | |
| return len(self.protein) | |
| def nt_len(self) -> int: | |
| return len(self.nt_seq) | |
| 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) | |
| 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)) | |