Matiq/ffpe_damage
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

pick best transcript per variant (1:1) default, --all-transcripts for expansion

Browse Source 
- 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
This commit is contained in:
Matiq
2026-09-07 01:00:58 +03:00
parent 9f4312ba3a
commit 8f8d728edd
1 changed files with 81 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
annotate_clean.py
+81 -40
View File
@@ -54,13 +54,12 @@ def load_token(path="~/.config/oncokb/token"):
return p.read_text().strip() return p.read_text().strip()
return os.environ.get("ONCOKB_TOKEN", "").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.""" """Try VEP if installed and cache exists, else GTF offline, else Ensembl REST, else fallback."""
import shutil import shutil
vep_bin = shutil.which("vep") vep_bin = shutil.which("vep")
if vep_bin and gtf_path and Path(gtf_path).expanduser().is_file(): if vep_bin and gtf_path and Path(gtf_path).expanduser().is_file():
pass pass
# Load GTF for offline all-transcript expansion
gtf_transcripts = None gtf_transcripts = None
if offline and gtf_path: if offline and gtf_path:
gtf_file = Path(gtf_path).expanduser() gtf_file = Path(gtf_path).expanduser()
@@ -89,22 +88,37 @@ def vep_annotate(df, reference, gtf_path=None, offline=False):
}) })
continue continue
if offline and gtf_transcripts is not None: 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"]] hits = [t for t in gtf_transcripts if t["chrom"] == chrom and t["start"] <= pos1 <= t["end"]]
if hits: if hits:
for t in hits: if all_transcripts:
# Simple HGVS: use transcript ID + positional offset (placeholder, VEP would give exact c./p.) for t in hits:
hgvs_c = f"{t['tx']}:c.{pos1}{ref}>{alt}" hgvs_c = f"{t['tx']}:c.{pos1}{ref}>{alt}"
hgvs_p = "p.(?)" hgvs_p = "p.(?)"
effect = t["biotype"] or "transcript_variant" effect = t["biotype"] or "transcript_variant"
rows.append({ rows.append({
"_orig_idx": r.name, "_orig_idx": r.name,
"Ген": t["gene"], "Ген": t["gene"],
"HGVS_c": hgvs_c, "HGVS_c": hgvs_c,
"HGVS_p": hgvs_p, "HGVS_p": hgvs_p,
"HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", "HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})",
"Тип варианта и эффект": effect, "Тип варианта и эффект": 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 continue
rows.append({ rows.append({
"_orig_idx": r.name, "_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", "Тип варианта и эффект": "SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel",
}) })
continue continue
# Online: Ensembl REST, expand all transcripts # Online: Ensembl REST
try: try:
hgvs_ens = f"{chrom.replace('chr','')}:g.{pos1}{ref}>{alt}" 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" 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() j = resp.json()
tcs = j[0].get("transcript_consequences", []) if j and isinstance(j, list) and j[0] else [] tcs = j[0].get("transcript_consequences", []) if j and isinstance(j, list) and j[0] else []
if tcs: if tcs:
for tc in tcs: if all_transcripts:
gene = tc.get("gene_symbol") or "intergenic" for tc in tcs:
hgvs_c = tc.get("hgvsc") or f"c.{pos1}{ref}>{alt}" gene = tc.get("gene_symbol") or "intergenic"
hgvs_p = tc.get("hgvsp") or "p.(?)" hgvs_c = tc.get("hgvsc") or f"c.{pos1}{ref}>{alt}"
cons = tc.get("consequence_terms", []) hgvs_p = tc.get("hgvsp") or "p.(?)"
effect = ", ".join(cons) if cons else ("SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel") cons = tc.get("consequence_terms", [])
if ":" in hgvs_c: effect = ", ".join(cons) if cons else ("SNV, missense_variant (predicted)" if len(ref)==1 and len(alt)==1 else "indel")
hgvs_c = hgvs_c.split(":")[-1] if ":" in hgvs_c:
if ":" in hgvs_p: hgvs_c = hgvs_c.split(":")[-1]
hgvs_p = hgvs_p.split(":")[-1] if ":" in hgvs_p:
rows.append({ hgvs_p = hgvs_p.split(":")[-1]
"_orig_idx": r.name, rows.append({
"Ген": gene, "_orig_idx": r.name,
"HGVS_c": hgvs_c, "Ген": gene,
"HGVS_p": hgvs_p, "HGVS_c": hgvs_c,
"HGVS": f"{hgvs_c} {hgvs_p} ({hgvs_g})", "HGVS_p": hgvs_p,
"Тип варианта и эффект": effect, "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 continue
except Exception: except Exception:
pass pass
@@ -295,6 +331,7 @@ def build_excel(df_clean, df_annot, out_path):
def _load_gtf(gtf_path): def _load_gtf(gtf_path):
import gzip import gzip
import re
transcripts = [] transcripts = []
try: try:
opener = gzip.open if str(gtf_path).endswith(".gz") else open opener = gzip.open if str(gtf_path).endswith(".gz") else open
@@ -308,13 +345,11 @@ def _load_gtf(gtf_path):
chrom = parts[0] chrom = parts[0]
start = int(parts[3]); end = int(parts[4]) start = int(parts[3]); end = int(parts[4])
attr = parts[8] attr = parts[8]
# parse gene_name, transcript_id
import re
m_gene = re.search(r'gene_name "([^"]+)"', attr) m_gene = re.search(r'gene_name "([^"]+)"', attr)
m_tx = re.search(r'transcript_id "([^"]+)"', attr) m_tx = re.search(r'transcript_id "([^"]+)"', attr)
m_biotype = re.search(r'transcript_type "([^"]+)"', attr) m_biotype = re.search(r'transcript_type "([^"]+)"', attr)
if m_gene and m_tx: 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({ transcripts.append({
"chrom": chrom, "chrom": chrom,
"start": start, "start": start,
@@ -322,6 +357,11 @@ def _load_gtf(gtf_path):
"gene": m_gene.group(1), "gene": m_gene.group(1),
"tx": m_tx.group(1), "tx": m_tx.group(1),
"biotype": m_biotype.group(1) if m_biotype else "", "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: except Exception as e:
print(f"[gtf] failed to load {gtf_path}: {e}", file=sys.stderr) 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("--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("--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("--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-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("--min-depth", type=int, default=None, help="filter depth >= threshold")
ap.add_argument("--gtf", default="~/Projects/ffpe_damage/references/gencode.v44.annotation.gtf.gz", 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() args = ap.parse_args()
clean_path = Path(args.clean) 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) print(f"[warn] {len(df_clean)} variants - network annotation will be slow. Use --offline for fast placeholder.", flush=True)
# HGVS / gene / effect # 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 # OncoKB
token = load_token(args.token) token = load_token(args.token)