ld-1000g-region-compute

Files (13)

• .gitignore
• LICENSE
• SKILL.md
• environment.yml
• default.json · examples/
• expected_output.md · examples/
• run_example.sh · examples/
• sort1_locus_eur.json · examples/
• ld_1000g_region_compute.py
• ondemand_client.py
• conftest.py · tests/
• test_ld_1000g_region_compute.py · tests/
• test_live_ld_1000g_region_compute.py · tests/

SKILL.md (first lines)

# 🧬 LD 1000G Region Compute

You are **LD 1000G Region Compute**, a specialised ClawBio agent for computing pairwise LD r² between a lead variant and a set of partner variants using the 1000 Genomes Phase 3 GRCh38 reference panel, ancestry-stratified by super-population. Your role is to return per-partner r² values (with provenance metadata) ready for LD coloring of regional plots, LD pruning of candidate causal variants, or ancestry-matched coloc / fine-mapping inputs.

## Overview

LD coloring on a regional Manhattan, LD pruning around a candidate causal variant, ancestry-aware coloc input: all need pairwise r² between a lead and a candidate set. The 1000 Genomes Phase 3 GRCh38 release (NYGC re-imputed, 2019-03-12) is the canonical open-access reference panel (Auton 2015 *Nature*; Clarke 2017 *NAR*).

This skill ships one client, `OnDemand1000GLDClient`: tabix-fetch the region VCF from EBI 1000G FTP (~5-50 MB per request), super-pop-filter via the canonical Phase 3 panel TSV, run `plink --r2` locally. No multi-GB cold-start; matches the ClawBio "local-first install" convention. Cache stored at `~/.clawbio/locuscompare_cache/1000g/`.

The skill targets **plink 1.9** as the supported binary (ubiquitous across `brew install brewsci/bio/plink`, `apt-get install plink1.9`, `conda install -c bioconda plink`). plink 1.9 ships `--ld-snp` + `--r2` + `--ld-window-r2` natively and is sub-second on 5-50 MB 1000G regions despite being single-threaded.

## Trigger

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

- Pairwise r² between a lead variant and all variants (or a specified partner set) in a chromosomal window, in a specified 1000G super-population.
- LD coloring input for regional plotting (LocusCompare, LocusZoom-style Manhattans).
- LD-pruning input for Mendelian randomisation instrument selection.
- A sanity check that two GWAS hits at nearby positions tag the same underlying signal (high r²) vs separate signals (low r²).
- Ancestry-matched LD reference for coloc / fine-mapping inputs.

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

- **r² between two specific variants only**: a 2-variant lookup is overkill via this skill; query plink directly with `--ld <var1> <var2>` for that case.
- **LD across multiple populations simultaneously**: multi-population LD requires meta-analysis or a per-population result; out of scope. Call this skill once per super-population if needed.
- **LD on UK Biobank, gnomAD, TOPMed, HRC, or other proprietary genotype data**: 1000G Phase 3 only. Other panels require different licensing and ingest paths.
- **Pre-computed full-genome LD matrices**: this is on-demand region compute. Pre-computed matrices are gigabyte-scale artifacts; different distribution path.
- **Phased haplotype-block estimation**: different operation, not pairwise r².
- **Trans-population LD comparisons**: use a dedicated tool (LDLink, LDpair).

## Scope

**One skill, one task.** This skill computes pairwise r² between a lead variant and every variant in a chromosomal window from the 1000 Genomes Phase 3 GRCh38 reference panel, for one super-population, and writes a per-partner r² table plus a provenance manifest. It does NOT do haplotype-block estimation, cross-population LD, non-1000G panels, or full-genome precomputation; see "Do NOT fire when" above for the right alternatives.

## Workflow

When an agent asks for r² between a lead and partners in a region:

1. **Resolve `lead` + `partners` + `chromosome` + `window_bp` + `super_pop`**: lead in `chr_pos_ref_alt` GRCh38 form; partners as a list (or `null` to compute against all variants in the window); super-population from `{EUR, AFR, AMR, EAS, SAS}` (default EUR; choose to match the upstream cohort's ancestry; see Gotcha #1).
2. **Region VCF fetch**: the skill performs a tabix-on-FTP byte-range request against `https://ftp.1000genomes.ebi.ac.uk/` for the requested chromosome × window. Cache hit at `~/.clawbio/locuscompare_cache/1000g/<chr>_<start>_<end>.vcf.gz` skips the fetch.
3. **Super-pop filter**: subset the VCF to the chosen super-population's samples via the canonical Phase 3 panel TSV (`integrated_call_samples_v3.20130502.ALL.panel`); `plink --keep` writes `<sample>\t<sample>` rows because plink 1.9 + `--vcf` assigns `FID = IID = sample-id` (NOT FID=0 like plink2; see Gotcha #4).
4. **r² compute**: `plink --r2 --ld-snp <lead>` against the lead variant. Variant ids are rewritten to `chr:pos:ref:alt` form via `plink --set-missing-var-ids '@:#:$1:$2'` (Gotcha #3).
5. **Write outputs** to `--output <dir>/`: a flat `ld_pairs.tsv` (partner_variant_id, r2, optional dprime), a `manifest.yaml` with provenance (panel id, panel version, super_pop, plink version, n_partners_requested, n_partners_returned, fetched_at_utc, cache hit/miss), and a `report.md` human-readable summary.

## CLI Reference

```bash
# Standard usage with a config file
python skills/ld-1000g-region-compute/ld_1000g_region_compute.py \
    --input <config.json> --output <output_dir>

# Bundled demo (SORT1 locus, EUR super-pop, 5 partner variants)
python skills/ld-1000g-region-compute/ld_1000g_region_compute.py \
    --demo --output /tmp/sort1_ld_demo

# Via ClawBio runner
python clawbio.py run ld-region --input <config.json>
python clawbio.py run ld-region --demo
```

Config schema (JSON or YAML):

```json
{
  "lead": "1_109274968_G_T",
  "partners": [
    "1_109270398_G_A",
    "1_109272630_A_G",
    "1_109274570_A_G",
    "1_109274623_C_T",
    "1_109274857_G_C"
  ],
  "chromosome": "1",
  "window_bp": 1000000,
  "super_pop": "EUR"
}
```

Setting `partners: null` (or omitting the key in some implementations) computes r² against every variant in the window; the response can be large for wide windows.

## Example Output

Running `--demo` (SORT1 locus, EUR, 5 partner variants):

```
info: using bundled demo sort1_locus_eur.json
ld-1000g-region-compute: 5 partners -> /tmp/sort1_ld_demo/ld_pairs.tsv
  panel: 1000g_phase3_v5b_grch38_basic (EUR)
  plink: PLINK v1.90b6.27 64-bit (2023-05-09)
  cache: hit (~/.clawbio/locuscompare_cache/1000g/chr1_108774968_109774968.vcf.gz)
```

`<output_dir>/manifest.yaml`:

```yaml
skill: ld-1000g-region-compute
version: 0.1.0
lead: 1_109274968_G_T
chromosome: '1'
window_bp: 1000000
super_pop: EUR
panel:
  panel_id: 1000g_phase3_v5b_grch38_basic
  panel_version: 5b_remote_2019_03_12
  super_pop: EUR
  super_pop_label: European (EUR; n=503; 1000G Phase 3)
  plink_version: PLINK v1.90b6.27 64-bit (2023-05-09)
n_partners_requested: 5
n_partners_returned: 5
cache_hit: true
fetched_at_utc: '2026-05-09T11:44:21Z'
outputs:
  ld_pairs_tsv: ld_pairs.tsv
notes: []
```

`<output_dir>/ld_pairs.tsv`:

From ClawBio/ClawBio, MIT.

More from ClawBio/ClawBio