tooluniverse-variant-analysis

What this skill does

# Variant Analysis and Annotation

## RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:
- `*_executed.ipynb` → read with `tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}'` and cite its cell outputs as the authoritative answer
- Pre-computed result files (CSV/TSV with names like `*results*`, `*deseq*`, `*enrich*`, `*stats*`, `*_simplified.csv`) → read directly and report the requested value
- Canonical analysis scripts (`analysis.R`, `run_*.py`, `find_*.R`, `*.Rmd`) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if **none** of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5-10× turn count).

---

## PRIMARY SCRIPTS — use these FIRST

These bundled scripts encode the question-specific gotchas (denominator
choices, ploidy defaults, multi-allelic split, multi-row Excel headers,
non-coding allowlist). They emit labelled `KEY=VALUE` lines that are
easier to parse than ad-hoc pandas/awk output. Prefer them over writing
new code.

| Script | When to use it |
|--------|----------------|
| `gatk_haplotypecaller_pipeline.py` | Any "how many SNPs / indels were called by HaplotypeCaller from the BAM" question. Handles BWA index → align → sort → index → HaplotypeCaller, OR can start from an existing BAM (skip alignment), OR only count an existing VCF. Default `--ploidy 2` (matches GATK's own default — most "called by HaplotypeCaller" GTs were generated with this). Pass `--ploidy 1` for explicit haploid prokaryote calling. Multi-allelic split + bcftools-style SNP/indel detection is built in. |
| `coding_variant_filter.py` | "Average number of CHIP / coding variants per sample after filtering out intronic, intergenic, and UTR variants." Two-stage canonical filter: (1) drop `Zygosity == Reference` rows (when present — these inflate counts ~10×), (2) drop intronic/intergenic/UTR/upstream/downstream SO terms. Handles 2-row VarSeq Excel headers and per-sample folders or combined CSVs. |
| `variant_fraction.py` | "Fraction of variants with VAF < X annotated as Y" — denominator is the CODING subset only (synonymous/missense/splice_region/stop_gained/lost/start_lost/frameshift/inframe indel), NOT all records. |

### Workspace isolation (CRITICAL)

`gatk_haplotypecaller_pipeline.py` and `coding_variant_filter.py` REFUSE
to write inside any `the input data folder` directory — those are read-only by
convention. Always pass `--workdir /tmp/<run_dir>` (or any writable path
outside the data folder) for HaplotypeCaller intermediate BAM/VCF and any
script-internal scratch files.

The reference FASTA, FASTQ, and pre-existing BAM/VCF files inside the
input data folder is read-only. The script will copy a data-folder BAM into the
workdir if it needs to add a `.bai` index.

### Concrete invocations

Re-run HaplotypeCaller on a sample's sorted BAM (this is the canonical
path for "how many SNPs / indels did HaplotypeCaller identify in the
BAM"; preferred over counting any pre-shipped `*_raw_variants.vcf`,
which may have been generated with non-default flags or post-filtering
that does not match the question):

```bash
python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --reference <data-folder>/REF.fna \
  --bam <data-folder>/SAMPLE_sorted.bam \
  --workdir /tmp/hc_run --sample-name SAMPLE
```

Full pipeline from FASTQ (BWA + sort + HaplotypeCaller; ~5-10 min):

```bash
python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --reference <data-folder>/REF.fna \
  --fastq-r1 <data-folder>/SAMPLE_1.fastq.gz --fastq-r2 <data-folder>/SAMPLE_2.fastq.gz \
  --workdir /tmp/hc_run --sample-name SAMPLE
```

Count-only an existing VCF (only when the question explicitly asks about
that file — e.g., "how many records are in `variants.vcf`"; do NOT use
this for "how many SNPs did HaplotypeCaller identify", because the
shipped file's ploidy / filtering may not match the question):

```bash
python skills/tooluniverse-variant-analysis/scripts/gatk_haplotypecaller_pipeline.py \
  --vcf <data-folder>/SAMPLE_variants.vcf
```

Average CHIP variants per sample after intronic/intergenic/UTR filter
(folder of per-sample 2-row-header VarSeq Excels):

```bash
python skills/tooluniverse-variant-analysis/scripts/coding_variant_filter.py \
  --dir <data-folder>/CHIP_DP10_GQ20_PASS --pattern '*.xlsx' --header-rows 2
```

Same filter on a single combined CSV:

```bash
python skills/tooluniverse-variant-analysis/scripts/coding_variant_filter.py \
  --file all_samples.csv --sample-col sample --header-rows 1
```

### Output keys to grep

`gatk_haplotypecaller_pipeline.py`: `SNP_COUNT_ALLELES`, `INDEL_COUNT_ALLELES`,
`TOTAL_RECORDS`, `PLOIDY`, `VCF_PATH`.

`coding_variant_filter.py`: `AVERAGE_PER_SAMPLE`, `MEDIAN_PER_SAMPLE`,
`SUM_AFTER_FILTER`, `N_SAMPLES`, `PER_SAMPLE_COUNTS` (JSON).

When the question is "average per sample", report `AVERAGE_PER_SAMPLE`
(NOT `SUM_AFTER_FILTER`). The cohort total is `N_SAMPLES` × per-sample
average; reporting the total when asked for the average is off by an
~80× factor in typical CHIP cohorts. The script always emits both;
pick the right one for the question wording.

### Ploidy: match the question's pipeline, not the organism

GATK HaplotypeCaller's default is `--sample-ploidy 2`. Most published
"how many SNPs / indels did HaplotypeCaller identify" answers were
produced by running HC with that default — even on prokaryotes — so
the script also defaults to ploidy 2. Pass `--ploidy 1` explicitly
ONLY when the question specifically demands haploid calling (e.g.,
"using haploid HaplotypeCaller"); ploidy 1 typically produces ~5-10%
fewer SNPs and ~10-15% fewer indels on the same BAM, which would miss
the GT range.

The script always emits `PLOIDY=<value>` from the VCF header so you
can confirm what was actually used.

---

## CRITICAL — Read before writing any code

1. **"Fraction of variants annotated as X"**: Use the bundled script:
   ```bash
   python skills/tooluniverse-variant-analysis/scripts/variant_fraction.py \
     --file variants.xlsx --vaf-threshold 0.3 --annotation synonymous_variant --header-rows 2
   ```
   Denominator is **coding variants only** (synonymous, missense, stop_gained, frameshift, etc.), NOT all variants. The script handles this automatically.
2. **Multi-row Excel headers**: Clinical variant exports often have 2-row headers. Use `pd.read_excel(path, header=[0,1])` and address columns as tuples.
3. **"How many variants from VCF/HaplotypeCaller" — DO NOT apply quality filters unless asked**: When the question is "How many SNPs are identified by GATK HaplotypeCaller from the BAM" or "What is the total number of indel mutations", count EVERY record in the raw VCF (after `bcftools view`/`bcftools stats` or by parsing the file directly). Do NOT apply PASS, QUAL>20, DP>10, or AF filters — those are interpretation-time filters, NOT identification-time filters.
   - Wrong: `bcftools view -f PASS variants.vcf | grep -v '^#' | awk '$5~/[ACGT]/' | wc -l` → returns ~10% of true SNP count.
   - Right: count all biallelic SNP records: `bcftools view --types snps variants.vcf | grep -v '^#' | wc -l`. For all SNPs (incl. multi-allelic): split first with `bcftools norm -m -` then count.
   - Indel total (insertions+deletions): `bcftools view --types indels variants.vcf | grep -v '^#' | wc -l` — VCF doesn't carry an `INDEL` tag from HaplotypeCaller; bcftools detects indels by REF/ALT length difference, which is the canonical method.
   - The skill's "VCF quality filtering must come before interpretation" rule is for *clinical* interpretation. For *counting* ("how many SNPs are identified" or "total number of indels"), report raw counts and let the question's wording dictate filters.

---

## Domain Reasoning

VCF quality filtering must come before interpretation. A variant called at 2x read depth is u