From 8f8d728eddf6bb6e79cb1af0d8af63a167d39a8e Mon Sep 17 00:00:00 2001 From: Matiq Date: Mon, 7 Sep 2026 01:00:58 +0300 Subject: [PATCH] pick best transcript per variant (1:1) default, --all-transcripts for expansion - offline GTF now picks MANE>appris1>canonical>basic>protein_coding (318 vs 3018 rows at VAF>10%) - online Ensembl REST likewise picks best (mane_select/canonical) - --all-transcripts restores previous all-transcripts behavior --- annotate_clean.py | 121 +++++++++++++++++++++++++++++++--------------- 1 file changed, 81 insertions(+), 40 deletions(-) diff --git a/annotate_clean.py b/annotate_clean.py index c7071a9..cce6343 100644 --- a/annotate_clean.py +++ b/annotate_clean.py @@ -54,13 +54,12 @@ def load_token(path="~/.config/oncokb/token"): return p.read_text().strip() return os.environ.get("ONCOKB_TOKEN", "").strip() -def vep_annotate(df, reference, gtf_path=None, offline=False): +def vep_annotate(df, reference, gtf_path=None, offline=False, all_transcripts=False): """Try VEP if installed and cache exists, else GTF offline, else Ensembl REST, else fallback.""" import shutil vep_bin = shutil.which("vep") if vep_bin and gtf_path and Path(gtf_path).expanduser().is_file(): pass - # Load GTF for offline all-transcript expansion gtf_transcripts = None if offline and gtf_path: gtf_file = Path(gtf_path).expanduser() @@ -89,22 +88,37 @@ def vep_annotate(df, reference, gtf_path=None, offline=False): }) continue if offline and gtf_transcripts is not None: - # GTF offline: all transcripts overlapping pos hits = [t for t in gtf_transcripts if t["chrom"] == chrom and t["start"] <= pos1 <= t["end"]] if hits: - for t in hits: - # Simple HGVS: use transcript ID + positional offset (placeholder, VEP would give exact c./p.) - hgvs_c = f"{t['tx']}:c.{pos1}{ref}>{alt}" - hgvs_p = "p.(?)" - effect = t["biotype"] or "transcript_variant" - rows.append({ - "_orig_idx": r.name, - "Ген": t["gene"], - "HGVS_c": hgvs_c, - "HGVS_p": hgvs_p, - "HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", - "Тип варианта и эффект": effect, - }) + if all_transcripts: + for t in hits: + hgvs_c = f"{t['tx']}:c.{pos1}{ref}>{alt}" + hgvs_p = "p.(?)" + effect = t["biotype"] or "transcript_variant" + rows.append({ + "_orig_idx": r.name, + "Ген": t["gene"], + "HGVS_c": hgvs_c, + "HGVS_p": hgvs_p, + "HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", + "Тип варианта и эффект": effect, + }) + continue + # pick best transcript (MANE > appris1 > canonical > basic > protein_coding) + def _score(t): + return (t["is_mane"], t["appris"] == 1, t["appris"] == 2, t["is_canonical"], t["is_basic"], t["is_ccds"], t["biotype"] == "protein_coding") + best = max(hits, key=_score) + hgvs_c = f"{best['tx']}:c.{pos1}{ref}>{alt}" + hgvs_p = "p.(?)" + effect = best["biotype"] or "transcript_variant" + rows.append({ + "_orig_idx": r.name, + "Ген": best["gene"], + "HGVS_c": hgvs_c, + "HGVS_p": hgvs_p, + "HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", + "Тип варианта и эффект": effect, + }) continue rows.append({ "_orig_idx": r.name, @@ -125,7 +139,7 @@ def vep_annotate(df, reference, gtf_path=None, offline=False): "Тип варианта и эффект": "SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel", }) continue - # Online: Ensembl REST, expand all transcripts + # Online: Ensembl REST try: hgvs_ens = f"{chrom.replace('chr','')}:g.{pos1}{ref}>{alt}" url = f"https://rest.ensembl.org/vep/homo_sapiens/hgvs/{hgvs_ens}?content-type=application/json" @@ -134,24 +148,46 @@ def vep_annotate(df, reference, gtf_path=None, offline=False): j = resp.json() tcs = j[0].get("transcript_consequences", []) if j and isinstance(j, list) and j[0] else [] if tcs: - for tc in tcs: - gene = tc.get("gene_symbol") or "intergenic" - hgvs_c = tc.get("hgvsc") or f"c.{pos1}{ref}>{alt}" - hgvs_p = tc.get("hgvsp") or "p.(?)" - cons = tc.get("consequence_terms", []) - effect = ", ".join(cons) if cons else ("SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel") - if ":" in hgvs_c: - hgvs_c = hgvs_c.split(":")[-1] - if ":" in hgvs_p: - hgvs_p = hgvs_p.split(":")[-1] - rows.append({ - "_orig_idx": r.name, - "Ген": gene, - "HGVS_c": hgvs_c, - "HGVS_p": hgvs_p, - "HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", - "Тип варианта и эффект": effect, - }) + if all_transcripts: + for tc in tcs: + gene = tc.get("gene_symbol") or "intergenic" + hgvs_c = tc.get("hgvsc") or f"c.{pos1}{ref}>{alt}" + hgvs_p = tc.get("hgvsp") or "p.(?)" + cons = tc.get("consequence_terms", []) + effect = ", ".join(cons) if cons else ("SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel") + if ":" in hgvs_c: + hgvs_c = hgvs_c.split(":")[-1] + if ":" in hgvs_p: + hgvs_p = hgvs_p.split(":")[-1] + rows.append({ + "_orig_idx": r.name, + "Ген": gene, + "HGVS_c": hgvs_c, + "HGVS_p": hgvs_p, + "HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", + "Тип варианта и эффект": effect, + }) + continue + def _score_tc(tc): + return (tc.get("mane_select") is not None, tc.get("canonical") == 1, tc.get("biotype") == "protein_coding", tc.get("impact") == "HIGH") + tc = max(tcs, key=_score_tc) + gene = tc.get("gene_symbol") or "intergenic" + hgvs_c = tc.get("hgvsc") or f"c.{pos1}{ref}>{alt}" + hgvs_p = tc.get("hgvsp") or "p.(?)" + cons = tc.get("consequence_terms", []) + effect = ", ".join(cons) if cons else ("SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel") + if ":" in hgvs_c: + hgvs_c = hgvs_c.split(":")[-1] + if ":" in hgvs_p: + hgvs_p = hgvs_p.split(":")[-1] + rows.append({ + "_orig_idx": r.name, + "Ген": gene, + "HGVS_c": hgvs_c, + "HGVS_p": hgvs_p, + "HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", + "Тип варианта и эффект": effect, + }) continue except Exception: pass @@ -295,6 +331,7 @@ def build_excel(df_clean, df_annot, out_path): def _load_gtf(gtf_path): import gzip + import re transcripts = [] try: opener = gzip.open if str(gtf_path).endswith(".gz") else open @@ -308,13 +345,11 @@ def _load_gtf(gtf_path): chrom = parts[0] start = int(parts[3]); end = int(parts[4]) attr = parts[8] - # parse gene_name, transcript_id - import re m_gene = re.search(r'gene_name "([^"]+)"', attr) m_tx = re.search(r'transcript_id "([^"]+)"', attr) m_biotype = re.search(r'transcript_type "([^"]+)"', attr) if m_gene and m_tx: - # Normalize chr prefix: GTF uses chr1, our clean uses chr1 -> keep as is + tags = set(re.findall(r'tag "([^"]+)"', attr)) transcripts.append({ "chrom": chrom, "start": start, @@ -322,6 +357,11 @@ def _load_gtf(gtf_path): "gene": m_gene.group(1), "tx": m_tx.group(1), "biotype": m_biotype.group(1) if m_biotype else "", + "is_mane": "MANE_Select" in tags, + "is_canonical": "Ensembl_canonical" in tags, + "is_basic": "basic" in tags, + "is_ccds": "CCDS" in tags, + "appris": next((int(t.split("_")[-1]) for t in tags if t.startswith("appris_principal_")), 99), }) except Exception as e: print(f"[gtf] failed to load {gtf_path}: {e}", file=sys.stderr) @@ -337,10 +377,11 @@ def main(): ap.add_argument("--tumor-type", default="All Solid Tumors", help="OncoKB tumor type (pan-cancer default)") ap.add_argument("--token", default="~/.config/oncokb/token", help="OncoKB token file or env") ap.add_argument("--offline", action="store_true", help="skip Ensembl/OncoKB network calls (fast, offline)") + ap.add_argument("--all-transcripts", action="store_true", help="expand all transcripts (default: pick best MANE/canonical per variant)") ap.add_argument("--min-vaf", type=float, default=None, help="filter VAF > threshold (e.g. 0.10 for 10%%)") ap.add_argument("--min-depth", type=int, default=None, help="filter depth >= threshold") ap.add_argument("--gtf", default="~/Projects/ffpe_damage/references/gencode.v44.annotation.gtf.gz", - help="GENCODE GTF for offline all-transcript expansion") + help="GENCODE GTF for offline expansion (pick best when not --all-transcripts)") args = ap.parse_args() clean_path = Path(args.clean) @@ -365,7 +406,7 @@ def main(): print(f"[warn] {len(df_clean)} variants - network annotation will be slow. Use --offline for fast placeholder.", flush=True) # HGVS / gene / effect - df_annot = vep_annotate(df_clean, args.reference, gtf_path=args.gtf, offline=args.offline) + df_annot = vep_annotate(df_clean, args.reference, gtf_path=args.gtf, offline=args.offline, all_transcripts=args.all_transcripts) # OncoKB token = load_token(args.token)