From 9f4312ba3a3b275e5ea301502bf8121e336793b6 Mon Sep 17 00:00:00 2001 From: Matiq Date: Mon, 7 Sep 2026 00:47:51 +0300 Subject: [PATCH] offline all-transcripts via GENCODE GTF + VAF>10% filter - GENCODE v44 50M (252k transcripts) for offline gene expansion (1 variant -> N rows) - --min-vaf/--min-depth filter before annotation (3501 -> 318 at VAF>10%) - --gtf support, --offline now uses GTF (no network, no VEP cache needed) - pan-cancer OncoKB retained for online mode --- annotate_clean.py | 93 ++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 88 insertions(+), 5 deletions(-) diff --git a/annotate_clean.py b/annotate_clean.py index b2e1ed7..c7071a9 100644 --- a/annotate_clean.py +++ b/annotate_clean.py @@ -55,11 +55,19 @@ def load_token(path="~/.config/oncokb/token"): return os.environ.get("ONCOKB_TOKEN", "").strip() def vep_annotate(df, reference, gtf_path=None, offline=False): - """Try VEP if installed and cache exists, else Ensembl REST, else fallback.""" + """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).is_file(): + 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() + if gtf_file.is_file(): + print(f" loading GTF {gtf_file} ...", flush=True) + gtf_transcripts = _load_gtf(gtf_file) + print(f" GTF loaded: {len(gtf_transcripts)} transcripts", flush=True) rows = [] total = len(df) for idx, (_, r) in enumerate(df.iterrows()): @@ -80,6 +88,33 @@ 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 + 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, + }) + continue + rows.append({ + "_orig_idx": r.name, + "Ген": "intergenic", + "HGVS_c": f"c.{pos1}{ref}>{alt}", + "HGVS_p": "p.(?)", + "HGVS": f"c.{pos1}{ref}>{alt} p.(?) ({hgvs_g})", + "Тип варианта и эффект": "SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel", + }) + continue if offline: rows.append({ "_orig_idx": r.name, @@ -90,7 +125,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 - # Real hg38 - try Ensembl REST, expand all transcripts + # Online: Ensembl REST, expand all transcripts 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" @@ -120,7 +155,6 @@ def vep_annotate(df, reference, gtf_path=None, offline=False): continue except Exception: pass - # Fallback single row rows.append({ "_orig_idx": r.name, "Ген": "intergenic", @@ -259,6 +293,42 @@ def build_excel(df_clean, df_annot, out_path): wb.save(out_path) print(f"written {out_path} ({ws.max_row-1} rows)") +def _load_gtf(gtf_path): + import gzip + transcripts = [] + try: + opener = gzip.open if str(gtf_path).endswith(".gz") else open + with opener(gtf_path, "rt") as fh: + for line in fh: + if line.startswith("#"): + continue + parts = line.rstrip("\n").split("\t") + if len(parts) < 9 or parts[2] != "transcript": + continue + 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 + transcripts.append({ + "chrom": chrom, + "start": start, + "end": end, + "gene": m_gene.group(1), + "tx": m_tx.group(1), + "biotype": m_biotype.group(1) if m_biotype else "", + }) + except Exception as e: + print(f"[gtf] failed to load {gtf_path}: {e}", file=sys.stderr) + return [] + return transcripts + + def main(): ap = argparse.ArgumentParser(description="Annotate clean_variants.csv -> Excel") ap.add_argument("--clean", required=True, help="clean_variants.csv from ffpe_damage_v2.py") @@ -267,10 +337,23 @@ 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("--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") args = ap.parse_args() clean_path = Path(args.clean) df_clean = pd.read_csv(clean_path) + # VAF/depth filter before annotation (VAF is alt_count/depth fraction) + n_before = len(df_clean) + if args.min_vaf is not None: + df_clean = df_clean[df_clean["VAF"] > args.min_vaf] + if args.min_depth is not None: + df_clean = df_clean[df_clean["depth"] >= args.min_depth] + if n_before != len(df_clean): + print(f"filter VAF>{args.min_vaf} depth>={args.min_depth}: {n_before} -> {len(df_clean)} variants", flush=True) + df_clean = df_clean.reset_index(drop=True) if df_clean.empty: print("clean_variants.csv is empty (or only header) - nothing to annotate") out = args.out or str(clean_path).replace(".csv", ".annotated.xlsx") @@ -282,7 +365,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, offline=args.offline) + df_annot = vep_annotate(df_clean, args.reference, gtf_path=args.gtf, offline=args.offline) # OncoKB token = load_token(args.token)