diff --git a/annotate_clean.py b/annotate_clean.py index cce6343..8b8c02b 100644 --- a/annotate_clean.py +++ b/annotate_clean.py @@ -239,14 +239,27 @@ def fetch_civic(): pass return None -def annotate_acmg(df): - """Placeholder ACMG: ClinVar lookup would go here. For MVP, all VUS unless known pathogenic.""" - # Real implementation would join ClinVar variant_summary + InterVar - # For synthetic data, mark as VUS +def annotate_acmg(df, clinvar_path=None): + """ACMG via ClinVar variant_summary if available, else VUS placeholder.""" + clinvar = {} + if clinvar_path: + cp = Path(clinvar_path).expanduser() + if cp.is_file(): + print(f" loading ClinVar {cp} ...", flush=True) + clinvar = _load_clinvar(cp) + print(f" ClinVar loaded: {len(clinvar)//2} variants", flush=True) out = [] for _, r in df.iterrows(): - # Heuristic: if VAF low and not in ClinVar, VUS - out.append("Uncertain significance") + chrom = str(r["chrom"]); pos = int(r["position"]) + 1 + ref = str(r["ref"]); alt = str(r["alt"]) + key = (chrom, pos, ref, alt) + key2 = (chrom.replace("chr",""), pos, ref, alt) + if key in clinvar: + out.append(clinvar[key]) + elif key2 in clinvar: + out.append(clinvar[key2]) + else: + out.append("Uncertain significance") return out def build_excel(df_clean, df_annot, out_path): @@ -369,6 +382,65 @@ def _load_gtf(gtf_path): return transcripts +def _load_clinvar(clinvar_path): + import gzip + m = {} + try: + opener = gzip.open if str(clinvar_path).endswith(".gz") else open + with opener(clinvar_path, "rt") as fh: + header = None + for line in fh: + if line.startswith("#"): + continue + if header is None: + header = line.rstrip("\n").split("\t") + # Find column indices + try: + ci_chr = header.index("Chromosome") + ci_start = header.index("Start") + ci_ref = header.index("ReferenceAllele") + ci_alt = header.index("AlternateAllele") + ci_sig = header.index("ClinicalSignificance") + except ValueError: + # Fallback for older format: try different names + continue + continue + parts = line.rstrip("\n").split("\t") + if len(parts) <= max(ci_chr, ci_start, ci_ref, ci_alt, ci_sig): + continue + chrom = parts[ci_chr] + # ClinVar Chromosome is 1,2.. not chr1 + chrom_norm = f"chr{chrom}" if not chrom.startswith("chr") else chrom + try: + pos = int(parts[ci_start]) + except ValueError: + continue + ref = parts[ci_ref]; alt = parts[ci_alt] + sig = parts[ci_sig] + # Normalize significance to ACMG + # ClinVar: Pathogenic, Likely pathogenic, Uncertain significance, Likely benign, Benign, etc. + # Map to our 5-tier + sig_lower = sig.lower() + if "pathogenic" in sig_lower and "likely" not in sig_lower: + acmg = "Pathogenic" + elif "likely pathogenic" in sig_lower: + acmg = "Likely pathogenic" + elif "benign" in sig_lower and "likely" not in sig_lower: + acmg = "Benign" + elif "likely benign" in sig_lower: + acmg = "Likely benign" + else: + acmg = "Uncertain significance" + key = (chrom_norm, pos, ref, alt) + # Also add without chr prefix for matching + key2 = (chrom, pos, ref, alt) + m[key] = acmg + m[key2] = acmg + except Exception as e: + print(f"[clinvar] failed to load {clinvar_path}: {e}", file=__import__("sys").stderr) + return {} + return m + 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") @@ -382,11 +454,17 @@ def main(): 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 expansion (pick best when not --all-transcripts)") + ap.add_argument("--clinvar", default="~/Projects/ffpe_damage/references/clinvar_variant_summary.txt.gz", + help="ClinVar variant_summary.txt.gz for ACMG P/L (auto if exists)") args = ap.parse_args() + # VAF>20% default for online (pan-cancer), leave --min-vaf configurable, offline keeps None + if args.min_vaf is None and not args.offline: + args.min_vaf = 0.20 + print(f"[info] online default --min-vaf 0.20 (use --min-vaf 0.05 to keep more)", flush=True) + 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] @@ -418,7 +496,7 @@ def main(): oncokb_map = fetch_oncokb(hgvs_g_for_oncokb, token, args.tumor_type, offline=args.offline) if token else {} # Attach ACMG/AMP/ONCO to annot rows - acmg_list = annotate_acmg(df_clean) + acmg_list = annotate_acmg(df_clean, clinvar_path=args.clinvar) # Map orig_idx -> acmg acmg_by_idx = {df_clean.index[i]: acmg_list[i] for i in range(len(acmg_list))} for idx, row in df_annot.iterrows():