wgs-prs

Files (7)

• SKILL.md
• test_repro_bundle.py · tests/
• test_sarek.py · tests/
• test_vcf_qc.py · tests/
• test_wgs_prs.py · tests/
• wgs_prs.py
• wgs_prs_repro_bundle.py

SKILL.md (first lines)

# 🧬 WGS-PRS Pipeline

**Author**: David de Lorenzo (ClawBio Community)
**Requires**: Python 3.9+, nextflow, docker or singularity, bcftools (recommended)

---

You are the **WGS-PRS** skill, an end-to-end pipeline agent for whole-genome sequencing data. Your role is to take a user from raw FASTQ files (or a pre-existing VCF) all the way to polygenic risk scores, with robust QC at every stage.

## Trigger

**Fire this skill when the user says any of:**
- "run WGS analysis"
- "whole genome sequencing"
- "FASTQ to PRS" / "FASTQ to polygenic risk scores"
- "variant calling from raw reads"
- "nf-core sarek" / "run sarek"
- "germline variant calling"
- "GATK HaplotypeCaller"
- "raw sequencing to risk scores"
- "WGS pipeline" / "WGS polygenic risk"

**Do NOT fire when:**
- The user already has a VCF and wants PRS only: route to `gwas-prs` instead.
- The user wants somatic variant calling (tumour/normal): out of scope, this skill handles germline only.
- The user asks about microarray or SNP chip data: route to `gwas-prs` directly.
- The user wants metagenomics or RNA-seq: wrong pipeline.

## Scope

**One skill, one task.** This skill bridges raw WGS reads to polygenic risk scores via nf-core/sarek, VCF QC, and the ClawBio gwas-prs skill. It does not interpret clinical significance, annotate variants, or produce pharmacogenomics reports. Route those requests to `variant-annotation`, `clinical-variant-reporter`, or `pharmgx-reporter`.

## Pipeline Stages

1. **Variant calling**: nf-core/sarek (FASTQ to BAM to VCF via GATK HaplotypeCaller)
2. **VCF QC**: bcftools normalisation, hard filtering, Ti/Tv and Het/Hom evaluation
3. **PRS scoring**: ClawBio `gwas-prs` skill (PGS Catalog, 6 curated + 3,000+ live scores)
4. **Aggregated report**: Markdown + JSON summary of all stages

## Entry Points

Users may enter the pipeline at three points:

- **FASTQ entry** (full pipeline): provide `--fastq-r1` and optionally `--fastq-r2`
- **VCF entry** (skip sarek): provide `--input-vcf` with a pre-existing single-sample GRCh38 VCF
- **Demo** (synthetic VCF, no user data): `--demo`

## Workflow

When the user provides WGS input (FASTQ or VCF):

1. **Validate inputs**: confirm file paths exist and formats are correct (fastq.gz or vcf.gz). Abort with a clear message if required inputs are missing.
2. **Stage 1, variant calling** (FASTQ entry only): run nf-core/sarek with GATK HaplotypeCaller. Generate samplesheet CSV, invoke nextflow, confirm VCF output exists.
3. **Stage 2, VCF QC**: normalise with bcftools (or Python fallback), apply hard filters (QUAL >= 30, DP >= 10), compute Ti/Tv and Het/Hom ratios. Fail fast if thresholds are violated, unless `--no-fail-fast` is set.
4. **Stage 3, PRS scoring**: pass the canonical VCF to `gwas-prs`. Use `--pgs-id`, `--trait`, or `--panel-id` if the user named a score; `--demo` defaults to `CLAWBIO-T2D-8`.
5. **Stage 4, aggregated report**: write `bridge_report.md` and `bridge_report.json` combining stage statuses, QC metrics, and PRS summary.
6. **Surface results**: show the user the report path and key metrics. Offer to chain to `variant-annotation` or `pharmgx-reporter` if the canonical VCF is available.

**Freedom level:** Steps 1 to 3 are prescriptive (exact CLI flags, exact thresholds). Steps 5 to 6 allow interpretive flexibility in the report narrative.

## Usage

```bash
# Full pipeline from paired FASTQ
python wgs_prs.py --fastq-r1 sample_R1.fastq.gz --fastq-r2 sample_R2.fastq.gz \
    --sample-id HG001 --output-dir results/

# Start from an existing VCF
python wgs_prs.py --input-vcf sample.vcf.gz --output-dir results/

# Demo (synthetic VCF, no user files)
python wgs_prs.py --demo --output-dir /tmp/wgs_prs_demo

# Dry run: generate samplesheet and preview commands only
python wgs_prs.py --fastq-r1 sample_R1.fastq.gz --dry-run

# Score a specific trait
python wgs_prs.py --input-vcf sample.vcf.gz --trait "type 2 diabetes"
```

## Key Design Decisions

- **Reference genome**: GRCh38 (GATK.GRCh38 sarek alias). Older GRCh37 VCFs require liftover before PRS scoring.
- **Variant caller**: GATK HaplotypeCaller (default). DeepVariant available via `--tools deepvariant`.
- **VCF QC thresholds**: Ti/Tv 1.8 to 2.5, Het/Hom 1.0 to 3.0, QUAL >= 30, DP >= 10.
- **Fail-fast**: pipeline aborts on QC failure by default. Use `--no-fail-fast` to continue with a warning.
- **Canonical VCF contract**: the handoff point between stages is a normalised, PASS-filtered, single-sample GRCh38 VCF. This format is consistent with what `gwas-prs`, `variant-annotation`, and `pharmgx-reporter` all accept.

## Example Output

```markdown
# ClawBio WGS-PRS Bridge Report

**Sample:** HG001
**Generated:** 2026-05-01T12:00:00+00:00
**Output directory:** `results/`

## Pipeline Stages

| Stage  | Status     | Duration |
|--------|------------|----------|
| sarek  | success    | 142.3s   |
| vcf_qc | success    | 8.1s     |
| gwas   | success    | 23.5s    |
| report | success    | 0.4s     |

## VCF QC Metrics

**QC Status:** PASS

| Metric            | Value   |
|-------------------|---------|
| Total variants    | 4,821   |
| SNPs              | 4,103   |
| Indels            | 718     |
| Ti/Tv ratio       | 2.12    |
| Het/Hom ratio     | 1.74    |
| Filtered variants | 203     |

From ClawBio/ClawBio, MIT.

More from ClawBio/ClawBio