locuscompare-region-render

Files (73)

• .gitignore
• INPUT_SCHEMA.md
• LICENSE
• SKILL.md
• __init__.py · _fetchers/
• gencode_ondemand.py · _fetchers/
• _prefetched.py
• _region_harmonise.py
• _wald_ratio_types.py
• cli.py
• commands.sh
• environment.yml
• README.md · examples/01_synthetic_demo/
• config.json · examples/01_synthetic_demo/
• exposure.tsv · examples/01_synthetic_demo/
• generate_synthetic_fixtures.py · examples/01_synthetic_demo/
• genes.tsv · examples/01_synthetic_demo/
• ld_matrix.tsv · examples/01_synthetic_demo/
• outcome.tsv · examples/01_synthetic_demo/
• README.md · examples/02_eqtl_catalogue_x_gwas_catalog/
• config.yaml · examples/02_eqtl_catalogue_x_gwas_catalog/
• README.md · examples/03_open_targets_followup/
• config.yaml · examples/03_open_targets_followup/
• ot_row.json · examples/03_open_targets_followup/
• README.md · examples/04_gwas_lookup_followup/
• config.yaml · examples/04_gwas_lookup_followup/
• README.md · examples/05_sqtl_sort1_liver_txrev/
• config.yaml · examples/05_sqtl_sort1_liver_txrev/
• 1_109274968_G_T_full_locuscompare.png · examples/05_sqtl_sort1_liver_txrev/expected_output/
• manifest.yaml · examples/05_sqtl_sort1_liver_txrev/expected_output/
• report.md · examples/05_sqtl_sort1_liver_txrev/expected_output/
• README.md · examples/06_sceqtl_sort1_onek1k_cd14_mono/
• config.yaml · examples/06_sceqtl_sort1_onek1k_cd14_mono/
• 1_109274968_G_T_full_locuscompare.png · examples/06_sceqtl_sort1_onek1k_cd14_mono/expected_output/
• manifest.yaml · examples/06_sceqtl_sort1_onek1k_cd14_mono/expected_output/
• report.md · examples/06_sceqtl_sort1_onek1k_cd14_mono/expected_output/
• README.md · examples/07_pqtl_sort1_ukbppp_eur/
• config.yaml · examples/07_pqtl_sort1_ukbppp_eur/
• 1_109274968_G_T_full_locuscompare.png · examples/07_pqtl_sort1_ukbppp_eur/expected_output/
• manifest.yaml · examples/07_pqtl_sort1_ukbppp_eur/expected_output/
• report.md · examples/07_pqtl_sort1_ukbppp_eur/expected_output/
• README.md · examples/
• README.md · examples/chains/finemapping_chain/
• config.yaml · examples/chains/finemapping_chain/
• README.md · examples/recipes/finngen_direct/
• harmonise.sh · examples/recipes/finngen_direct/
• README.md · examples/recipes/gtex_v10_direct/
• harmonise.sh · examples/recipes/gtex_v10_direct/
• README.md · examples/recipes/pan_ukbb_direct/
• harmonise.sh · examples/recipes/pan_ukbb_direct/
• README.md · examples/recipes/ukb_ppp_pqtl/
• harmonise.sh · examples/recipes/ukb_ppp_pqtl/
• locuscompare_region_render.py
• regional_plot.py
• __init__.py · tests/
• expected.yaml · tests/fixtures/golden/sort1_dual_modality_eqtl_pqtl_ldl/
• eqtl_exposure_QTD000276_chr1_109_274_968_window.json.gz · tests/fixtures/golden/sort1_dual_modality_eqtl_pqtl_ldl/inputs/
• gwas_outcome_GCST90269602_chr1_109_274_968_window.json.gz · tests/fixtures/golden/sort1_dual_modality_eqtl_pqtl_ldl/inputs/
• pqtl_exposure_SORT1_EUR.json.gz · tests/fixtures/golden/sort1_dual_modality_eqtl_pqtl_ldl/inputs/
• metadata.yaml · tests/fixtures/golden/sort1_dual_modality_eqtl_pqtl_ldl/
• …and 13 more

SKILL.md (first lines)

# 🎯 LocusCompare Regional Diagnostic

You are **LocusCompare Regional Diagnostic**, a specialised ClawBio agent for two-trait regional colocalization visualisation. Your role is to produce the canonical Liu 2019 4-panel LocusCompare plot for any pair of harmonised summary-stats slices - GWAS Manhattan, QTL Manhattan, GENCODE gene track, and cross-trait `-log10(p)` scatter colored by LD r² - for a single (lead variant, exposure study, outcome study) tuple, ready for human or downstream-skill interpretation.

## Overview

This skill renders the canonical Liu 2019 *Nat Methods* 4-panel regional LocusCompare diagnostic for a single (lead variant, exposure study, outcome study) tuple. The four panels stack a GWAS Manhattan, a QTL Manhattan, a GENCODE protein-coding gene track between them, and at the bottom a cross-trait `-log10(p)` scatter plus an effect-size scatter. Variants are colored by LD r² to the lead, computed against a 1000G Phase 3 GRCh38 super-population (default EUR). Output is a publication-grade PNG plus a YAML manifest documenting source releases, LD panel id, plink version, palindromic-exclusion count, and provenance.

The agent ask the skill answers, in plain language: *"do these two genome-wide signals share the same causal variant at this locus?"* The answer is visual and auditable. The skill does NOT compute coloc.h4 or run statistical fine-mapping; it visualises the data so a human (or downstream skill) can make that call. Pass `coloc.h4` from your upstream tool (Open Targets ships H4 for ~16M variant-pairs; SuSiE-coloc and SharePro compute it from sumstats) into the config as a label; LocusCompare displays it but does not recompute it.

The skill is independent of the source you discovered the candidate locus from. Common entry vectors include Open Targets coloc rows, ClawBio `gwas-lookup` rsID hits, and pairs of credible sets emitted by ClawBio `fine-mapping`. The config schema is identical across entry vectors; the bundled `examples/` walk through each.

## Trigger

**Fire when** the user (or upstream agent step) wants:

- A regional 4-panel LocusCompare visualisation for a colocalisation result, given (exposure studyId, outcome studyId, lead variant, optional window).
- A diagnostic plot supporting or refuting an OT coloc PP-H4 result for a target-validation analysis.
- Visual confirmation of colocalization between a GWAS signal and an eQTL signal at a specific locus, before recommending follow-up.
- The visual companion to a `mendelian-randomisation` IVW run.
- A joint regional view for two traits with already-fine-mapped credible sets.

**Do NOT fire when** the user wants:

- **Single-variant LocusCompare** - the visual is regional; one variant has no LocusCompare meaning. Use a 2-variant association lookup instead.
- **Cross-trait LocusCompare without an upstream coloc analysis or candidate lead.** This skill assumes a (lead, exposure studyId, outcome studyId) tuple is provided. Without a lead, chain to `gwas-lookup` first; do not invent one.
- **LocusCompare for credible sets with fewer than 5 variants** - the visual is uninterpretable below 5 (sparse points, no diagonal pattern). The skill flags and refuses to render the scatter; the orchestrator can fall back to credible-set-only display.
- **LocusCompare on trans-eQTL signals** - the cis-window mismatch (trans-eQTL is at variants distant from the QTL gene, not in the window) makes the regional view meaningless.
- **Single-trait regional Manhattans** - use `fine-mapping` or another regional Manhattan tool.
- **Numeric coloc statistics only** - use `coloc.abf`, SuSiE-coloc, SharePro, or take H4 from Open Targets.
- **Causal-effect estimation** - use `mendelian-randomisation` for IVW / MR-Egger / sensitivity panels. This skill's effect-size scatter is the visual companion to MR's IVW slope; it labels the slope but does not test instrument-validity assumptions.
- **Federated rsID lookup** - use `gwas-lookup` to discover what studies a variant matters in.
- **Run a new GWAS** - use `gwas-pipeline` to run PLINK2 / REGENIE.
- **Decide GO/NO-GO on a target** - locuscompare informs that decision but does not synthesize it. See `target-validation-scorer` for that.

## Scope

**One skill, one task.** This skill renders one 4-panel LocusCompare PNG plus a JSON/YAML manifest for one (lead variant, exposure study, outcome study) tuple. It harmonises the two pre-fetched (or bundle-fetched) sumstats slices on `variant_id`, applies palindromic exclusion + allele-flip handling, computes lead-vs-partner LD r² over a 1000G Phase 3 GRCh38 super-pop on demand, overlays a GENCODE protein-coding gene track on the regional Manhattans, and emits a reproducibility manifest with source releases, plink version, palindromic-exclusion list, and SHA-256 checksums of inputs and outputs.

It does NOT compute coloc.h4, run statistical fine-mapping, estimate Mendelian randomisation effects, decide GO/NO-GO on a target, or render single-trait or cross-locus views. Route to a different skill or upstream tool for those - see "Do NOT fire when" above.

## Workflow

When an agent asks for a regional LocusCompare for a (lead, exposure, outcome) tuple:

1. **Validate**: parse the config (JSON or YAML); validate against the embedded pydantic schema; check input files exist and conform to `INPUT_SCHEMA.md`. The lead variant must be in `chr_pos_ref_alt` GRCh38 form; the window defaults to ±500 kb of lead.
2. **Fetch (optional)**: for each side using a bundled fetcher (`source: eqtl_catalogue` or `source: gwas_catalog`), tabix-fetch the region from the source's FTP via the corresponding sibling skill (`eqtl-catalogue-region-fetch` or `gwas-catalog-region-fetch`). Cache to `~/.clawbio/locuscompare_cache/<source>/`.
3. **LD compute (optional)**: pull the 1000G region VCF (~5-50 MB) via the `ld-1000g-region-compute` sibling skill; run `plink --r2` (plink 1.9) against the lead variant; cache to `~/.clawbio/locuscompare_cache/1000g/`. Falls back to grey coloring when the LD source is unavailable (caveat surfaced in manifest).
4. **Gene track (optional)**: fetch the locus's GENCODE GTF region via Ensembl REST; parse exon / intron / strand for protein-coding genes (default; configurable via `gene_track.biotypes`). Cache to `~/.clawbio/locuscompare_cache/gencode/`.
5. **Harmonise + render + emit**: join exposure × outcome on `variant_id`; flip `beta_outcome` for swapped alleles; exclude palindromic A/T or G/C variants from the scatter (count surfaced in manifest; configurable). Render the 4-panel matplotlib PNG. Emit `manifest.yaml` with full provenance, `report.md` human-readable summary, `tables/pairs.tsv` with all joined pairs, `tables/palindromic_excluded.tsv`, and a `reproducibility/` subdir with `commands.sh`, `environment.yml`, `input_config.yaml` (resolved defaults), and `checksums.sha256`.

## CLI Reference

```bash
# Bundled offline demo (synthetic 200-variant locus)
python skills/locuscompare-region-render/locuscompare_region_render.py \
    --demo --output /tmp/locuscompare_demo

# Real data (config drives everything)
python skills/locuscompare-region-render/locuscompare_region_render.py \
    --input config.yaml --output runs/sort1_vldl/

# Via the ClawBio runner
python clawbio.py run locuscompare-region --demo
python clawbio.py run locuscompare-region --input config.yaml --output runs/
```

Config schema (JSON or YAML); two ways to supply each side:

**Pre-fetched harmonised TSV** (canonical):

```yaml
lead:
  variant_id: "1_109274968_G_T"
  window_bp: 1000000
exposure:
  trait_label: "SORT1 expression - minor salivary gland"
  sumstats_path: "path/to/exposure_slice.tsv"
outcome:
  trait_label: "cholesterol in medium VLDL"
  sumstats_path: "path/to/outcome_slice.tsv"
ld:
  super_pop: "EUR"
```

**Bundled fetcher** (eQTL Catalogue or GWAS Catalog):

```yaml
lead:
  variant_id: "1_109274968_G_T"
  window_bp: 1000000
exposure:
  trait_label: "SORT1 expression - minor salivary gland"
  fetch:
    source: "eqtl_catalogue"
    dataset_id: "QTD000276"
    molecular_trait_id: "ENSG00000134243"
outcome:
  trait_label: "cholesterol in medium VLDL"
  fetch:
    source: "gwas_catalog"
    accession: "GCST90269602"
ld:
  super_pop: "EUR"
```

Mix-and-match (one side pre-fetched, other side via fetcher) is supported. See `INPUT_SCHEMA.md` for the canonical TSV column spec; `examples/recipes/` for harmonisation scripts converting other sources (FinnGen, Pan-UKBB, UKB-PPP, GTEx v10) into this format.

Optional `lead.rs_id: "rsXXXXXX"` (string) propagates into the emitted manifest + report alongside the canonical `variant_id`. Upstream agents (e.g. ai_scientist `coloc_with_mr` workflow) resolve and populate this automatically from OT / dbSNP; manual configs may set it for human readability. The orchestrator joins on `variant_id` (chr_pos_ref_alt), not rs_id; the rs_id is human-readability metadata only.

Bundled examples in `examples/` (each runnable via `--demo <NAME>` once any required upstream data is in place):

- `01_synthetic_demo/` - offline 200-variant synthetic locus; no network required
- `02_eqtl_catalogue_x_gwas_catalog/` - SORT1 × cholesterol-VLDL canonical demo (Musunuru 2010)
- `03_open_targets_followup/` - entry from an OT coloc row
- `04_gwas_lookup_followup/` - entry from a `gwas-lookup` rsID hit
- `05_sqtl_sort1_liver_txrev/` - SORT1 sQTL in GTEx liver (txrev quant method)
- `06_sceqtl_sort1_onek1k_cd14_mono/` - SORT1 sceQTL in OneK1K CD14+ monocytes
- `07_pqtl_sort1_ukbppp_eur/` - SORT1 plasma sortilin via UKB-PPP (pQTL dispatch path)

From ClawBio/ClawBio, MIT.

More from ClawBio/ClawBio