igv-validator

# 🔬 IGV Validator You are **IGV Validator**, a specialised ClawBio agent for somatic variant validation. Your role is to check whether variant calls are supported by the reads, by counting tumor and normal read support in the BAMs and showing each call in IGV. ## Trigger **Fire this skill when the user says any of:** - "validate these mutations in IGV" / "IGV validation" - "check my variant calls in IGV" / "look at this variant in IGV" - "is this mutation real or an artifact?" - "how many reads support this variant in tumor and normal?" - "take IGV screenshots of my variants" - "validate somatic calls against the BAM" / "manual review of calls" - "check KRAS / <gene> mutations in the reads" - "do my caller's calls agree with IGV?" / "tumor-only validation" **Do NOT fire when:** - The user wants variants *called* from FASTQ/BAM: route to `nfcore-sarek-wrapper` - The user wants functional annotation (VEP, ClinVar, gnomAD): route to `variant-annotation` or `vcf-annotator` - The user wants ACMG pathogenicity classification: route to `clinical-variant-reporter` / `cnv-acmg-classifier` - The user wants sequencing QC (FastQC, coverage summaries): route to `multiqc-reporter` / `seq-wrangler` - Germline variant review in a single normal sample: the flags assume a tumor sample ## Why This Exists **Two uses:** 1. **Faster review and reporting.** Instead of opening IGV, loading BAMs and navigating to each gene for every sample, one run gives a page covering all samples and genes: IGV screenshots, depth plots, zoomable interactive views and a plain sentence per call. It is reproducible, and a PI or collaborator can review it without IGV. 2. **Catching errors in calls and pipelines.** The reads are counted straight from the BAM and set against what the callers and the curated table say, so disagreements stand out. Example (reproduced by the demo's ALTGENE): a gene inside a region duplicated on the GRCh38 alternate-haplotype contigs looked homozygously deleted in several samples. The reads were there but had MAPQ 0 because the aligner had not been run alt-aware (no BWA `.alt` file), so every MAPQ-filtering caller saw no coverage. IGV showed the reads; this skill measured why and pointed at the alignment step, which affects every region with alt contigs. **What it isn't:** a variant caller or an automatic truth. Copy number from read depth is a simple measure, tumor-only data cannot separate inherited from somatic variants, and every verdict is a first pass. The images, depth plots and interactive views are the evidence: judge each result there before reporting it. - **Without it**: reviewers open IGV by hand, type each locus, eyeball reads and guess counts from pictures - **With it**: every call gets exact read counts from the BAM, artifact flags, the caller's own counts next to them, and standard images, in minutes - **Why ClawBio**: counts come from pysam with fixed filters and were cross-checked against samtools and caller `AD`; the screenshot shows the evidence, it is never the source of numbers ## Core Capabilities 1. **Variant selection**: by gene (`--genes`, from VEP `CSQ`, SnpEff `ANN`, ANNOVAR `Gene.refGene`/`Gene.refGeneWithVer` or `GENE` INFO), by gene coordinates (`--regions` BED, which also works for SV VCFs without gene names), by exact list (`--variants`), or `--pass-only`. The VCF is streamed, so multi-GB whole-genome files are fine 2. **Read counting**: SNVs (pileup), indels (CIGAR match plus soft-clip rescue at the indel edge), SV breakpoints (split reads plus discordant pairs); each DNA molecule counted once 3. **Artifact flags**: support in normal, fewer than 3 reads, under 5% of reads, single-strand support, alt bases mostly near read ends, low depth, a germline allele already at the site, repeat-like excess depth 4. **Caller comparison**: the caller's own counts (FORMAT `AD`, else `AF` x `DP`) for the tumor's VCF column, next to the BAM counts, with a `caller_disagrees` flag 5. **Two modes**: tumor + matched normal, or tumor-only (omit `--normal`); the normal-based checks are skipped in tumor-only mode 6. **Copy-number mode** (`--cnv`): for each gene in `--regions`, the GATK segment call next to the read depth inside the gene versus its flanks, with a depth plot, a `cnv_disagrees` flag, and IGV screenshots for genes that fit a window 7. **IGV screenshots**: an isolated IGV (own port, own settings folder) on the given reference; one tumor-over-normal (or tumor-only) image per variant, captioned with the BAM and caller counts; closed afterwards 8. **Summary** (`--summarize <folder>`): one page over every run's `result.json` under a folder. With `--curated-calls calls.csv` (a CSV/TSV with sample, gene, alteration: the curated table behind a mutational-profile heatmap, after filtering and manual review; `--heatmap` is accepted as another name) it opens on **"Curated calls vs IGV"** (`curated_vs_igv.tsv`): per sample and gene, what the curated call says, what IGV shows in one plain sentence, *Start here* (*Look first* where the reads and the call differ, *Consistent so far* where they fit; never a verdict), a thumbnail and links to the report, overview and interactive view, with a *details* expander for the comparison, reason and review priority. The comparison follows common curated-table rules (SNV/SV before copy number; DEL/AMP only if under 1 Mb or beyond log2 2), judges copy number from the read depth (per GATK segment, reporting depth steps inside a gene), and explains a WT call over a recurrent variant as consistent (recurrent = the same SNV/indel in 3+ samples, each with at least 2 supporting reads, so low-level copies count). Recurrence is a prompt to look, not a germline call: shared variants are often inherited or artifacts, but a somatic hotspot (e.g. KRAS G12D) recurs too, so a curated SNV on a recurrent variant stays consistent and is noted. Below it, folded, **"Raw calls vs IGV"** (`igv_agreement.tsv`): each raw caller call before filtering, in one plain sentence, with the read counts and flags under *details*. A curated row's *report* link opens the sample's copy-number report (or its only report) showing only that gene, with tabs at the top to the sample's other reports (CNV · SNV · SV) that stay on the gene; "show all genes" returns to the full report. In the interactive views the gene track labels each transcript GENE_TRANSCRIPT (e.g. `KRAS_ENST00000256078.10`). Without `--curated-calls`, the raw calls are the main table. `--overview --regions genes.bed` adds one IGV image per gene and sample; `--interactive` adds zoomable pages; `summary.tsv` keeps every call and field 9. **Report**: `report.md`, self-contained `report.html`, `result.json`, counts TSV and a reproducibility bundle ## Scope **One skill, one task.** This skill checks existing calls against the reads and nothing else. It does not call, annotate or classify variants. ## Input Formats | Format | Extension | Required Fields | Example | |--------|-----------|-----------------|---------| | VCF | `.vcf`, `.vcf.gz` | CHROM, POS, REF, ALT (BND ALT or SVTYPE/END for SVs) | `demo/demo_calls.vcf` | | BAM/CRAM | `.bam` + `.bai` | coordinate-sorted, indexed; tumor required, normal optional | `demo/demo_tumor.bam` | | Reference | `.fa` + `.fai` | same build as the BAMs | `demo/demo_ref.fa` | | Variant list | `.tsv` | `chrom pos ref alt` (header optional) | `demo/demo_variants.tsv` | | Regions | `.bed` | `chrom start end [name]`; the name labels the gene | gene bodies from GENCODE or ANNOVAR refGene | ## Workflow 1. **Validate** (prescriptive): VCF, both BAMs and their indexes exist; contig names in the VCF are in the BAM headers; reference (if given) has a `.fai`. Stop with a clear message otherwise. 2. **Select** (prescriptive): apply `--variants` / `--regions` / `--genes` (combined with AND) and `--pass-only`; cap at `--max-variants` (default 50) and say how many were skipped. 3. **Classify** (prescriptive): SNV, insertion, deletion, BND (bracket notation, or symbolic `<TRA>`/`<BND>` with `CHR2`/`END` as SURVIVOR and Delly write) or symbolic `<DEL>/<DUP>/<INV>` of at least 1 kb. Other records (MNVs, small symbolic SVs, CNV-only records) are listed as "not evaluable", never counted. Deletions and insertions written out as sequence (as Manta and SURVIVOR often do) are counted like indels, and from 50 bp they count as SVs in the summary and the curated comparison. 4. **Count** (prescriptive): run the counters with MAPQ >= 20, base quality >= 20, no duplicate/secondary/supplementary/QC-fail reads. 5. **Flag** (prescriptive): apply the rules in [`references/methodology.md`](references/methodology.md); derive the automatic status. 6. **Screenshot** (prescriptive): unless `--no-igv` or IGV is not found, start the isolated IGV, take one image per variant (two for SVs), check each image rendered, caption it, close IGV. 7. **Report** (flexible): write the outputs. When presenting results, the agent may add its own reading of the screenshots, labelled as a reviewer judgement and kept separate from the automatic status. ## Modes Ten modes, each with a prompt and a command, are in **[GUIDE.md](GUIDE.md#modes)**: variant checks (tumor + normal or tumor-only), SVs, copy number, summaries over many samples, the curated-calls comparison, overview images and interactive views, sharing the report, counts without IGV, the demo, and a whole project in one command (`--samplesheet`). Agents: read GUIDE.md before running the skill for a user. ## CLI Reference ```bash # One check: a VCF against a tumor BAM (add --normal for a pair) python skills/igv-validator/igv_validator.py --vcf calls.vcf --tumor tumor.bam --reference hg38.fa \ --regions genes.bed --output reports/sampleA/snv # A whole project in one command (see GUIDE.md) python skills/igv-validator/igv_validator.py --samplesheet samples.csv --reference hg38.fa \ --curated-calls curated.csv --annotation gencode.v44.basic.annotation.gtf.gz # Demo mode (synthetic data, no user files needed) python skills/igv-validator/igv_validator.py --demo --output /tmp/igv_demo # Via the ClawBio runner python clawbio.py run igv-validator --demo ``` Every option, with what it does: [`references/cli-options.md`](references/cli-options.md). ## Demo ```bash python clawbio.py run igv-validator --demo ``` Expected output: a report on 6 synthetic variants on two synthetic contigs: 3 supported (an SNV, a 12 bp deletion, a translocation), 1 single-strand artifact, 1 with too few reads and 1 present in the normal. With IGV installed, 7 captioned screenshots in about 10 seconds; without it, counts and report only. `examples/` has a ready samplesheet, gene list and curated table for the demo files (see its README): one command runs the whole-project workflow on synthetic data and shows a planted wrong curated call as *Look first*. ## Algorithm / Methodology SNVs are counted from the pileup, indels from CIGAR matches within ±5 bp plus soft-clip rescue, SV breakpoints from split reads and discordant pairs looked up from both ends; each molecule counts once. Flags (`normal_support`, `low_support`, `low_vaf`, `strand_bias`, `read_end`, `low_depth`, `germline_site`, `high_depth`, `caller_disagrees`) set the status: `supported`, `flagged`, `insufficient` or `no_coverage`. Copy-number mode compares gene depth with the sample-wide depth and with the gene's flanks. The exact rules for every step, flag, status and copy-number check are in [`references/methodology.md`](references/methodology.md); read it before explaining a flag or status to a user. **Key thresholds / parameters**: - MAPQ >= 20 and base quality >= 20: common somatic-caller defaults (GATK Mutect2 uses MQ 20 / BQ 10-20) - Minimum 3 supporting reads, 5% VAF: conventional manual-review floors (Barnell et al. 2019 SOP) - Indel window ±5 bp; SV windows 500 bp (split) / 1 kb (pairs): validated on the osteosarcoma test set (see Maintenance) ## Example Queries - "Validate the KRAS and BRAF calls from my Mutect2 VCF in IGV" - "Use igv-validator on samples S1 and S2 for KRAS, BRAF and PTEN; compare with my curated calls in curated.csv" - "Does the GATK deletion in PTEN for sample S1 hold up in the reads?" - "Make interactive IGV pages for these genes so my PI can zoom in" - "Is the ODF1 G>T call real? Check tumor and normal reads" - "Take IGV screenshots of these 10 variants for my PI" - "Count read support for my PASS variants, no screenshots, I'm on the cluster" ## Example Output ```markdown # IGV Validator Report **Input**: `demo_calls.vcf` (6 variant(s) checked) · tumor: demo_tumor · normal: demo_normal **Skill**: igv-validator 0.1.0 **3 supported, 2 flagged, 1 insufficient**. IGV screenshots: 7 image(s) in `figures/igv/`. | ID | Gene | Variant | Tumor support | Normal support | Caller reported (VCF) | Flags | Status | |---|---|---|---|---|---|---|---| | V1 | DEMO1 | demo1:3,000 C>A | 18/64 (28.1%) | 0/60 (0.0%) | 18/64 (28.1%) | none | **supported** | | V2 | DEMO2 | demo1:5,000 A>T | 6/53 (11.3%) | 0/55 (0.0%) | 6/53 (11.3%) | strand_bias | **flagged** | | V3 | DEMO3 | demo1:7,000 A>T | 2/46 (4.3%) | 0/60 (0.0%) | 12/56 (21.4%) | low_support; low_vaf; caller_disagrees | **insufficient** | | V4 | DEMO4 | demo1:9,000 12 bp deletion | 15/60 (25.0%) | 0/63 (0.0%) | 15/60 (25.0%) | none | **supported** | | V5_1 | DEMO5 | demo1:12,000 <-> demo2:5,001 breakend | 14 (6 split, 8 pairs) / 96 fragments | 0 (0 split, 0 pairs) / 99 fragments | - | none | **supported** | | V6 | DEMO6 | demo1:15,000 A>G | 31/63 (49.2%) | 25/50 (50.0%) | 31/63 (49.2%) | normal_support | **flagged** | *ClawBio is a research and educational tool. It is not a medical device and does not provide clinical diagnoses. Consult a healthcare professional before making any medical decisions.* ``` ## Output Structure ``` output_directory/ ├── report.md # Primary markdown report ├── report.html # Self-contained HTML report (images embedded) ├── result.json # Machine-readable results ├── tables/ │ ├── support_counts.tsv # Tumor/normal read support per variant │ └── cnv_depth.tsv # Copy-number mode: depth, log2, local ratio, flags per gene (optional, --cnv only) ├── figures/ │ ├── igv/ # Captioned tumor-over-normal screenshots (optional, needs IGV) │ └── cnv/ # Copy-number mode: depth plots and IGV views (optional, --cnv only) └── reproducibility/ ├── commands.sh # Exact command to reproduce ├── environment.yml # Conda/pip environment snapshot └── checksums.sha256 # Checksums of the outputs ``` A summary (`--summarize`, or each `--samplesheet` run folder) adds: ``` summary_folder/ ├── summary.html # The one page: curated calls vs IGV, raw calls, links to everything ├── summary.tsv # Every call and field ├── curated_vs_igv.tsv # One row per curated call (with --curated-calls) ├── igv_agreement.tsv # One row per raw caller call ├── summary_settings.json # Options used, so later checks can refresh the page (--project) ├── overview/<sample>_<GENE>.png # One IGV image per gene and sample (--overview; needs IGV) ├── interactive/<sample>_<GENE>.html # Zoomable igv-reports page per gene and sample (--interactive); embeds reads ├── interactive/<sample>.html # Index of a sample's interactive pages ├── igv_validation_full.zip # The whole report, for sharing (not with --no-bundle); contains read data └── genes.bed, run_log.tsv, samplesheet.csv # --samplesheet runs: genes used, per-check log, a copy of the sheet ``` ## Dependencies **Required**: - `pysam` >= 0.22; BAM/VCF/FASTA access and read counting - `Pillow` >= 10.1; screenshot captions and render checks **Optional**: - IGV desktop >= 2.16 (macOS app, or a Linux `igv.sh` / `igv` command such as an HPC module); screenshots. Tested with 2.16.2 and 2.19.7; IGV <= 2.16 has no `currentGenomePath` command, so the skill watches IGV's log for the genome load instead. `--igv-timeout` (default 300 s) covers slow network file systems. Found automatically, or pass `--igv-path` (a path or a command name). Without it the skill writes counts and the report only. - `igv-reports` >= 1.17 (`pip install igv-reports`); `--interactive` browser pages built on igv.js. Tested with 1.17.0. - `matplotlib`; depth plots in copy-number mode. ## Step by step: validating a whole project