PGX Pipeline — an end-to-end pharmacogenomics platform for WES cohorts

FASTQ → aligned reads → SNV/indel + CNV calls → VEP annotation → pharmacogenomic database enrichment → cohort analytics → an interactive variant browser a researcher can actually use.

Built solo during my M1 GENIOMHE (AI & Bioinformatics) internship at Université Paris-Saclay / Université Évry, hosted in a clinical pharmacogenomics lab (Grupo Duponte).

This is a showcase repository. It contains documentation, results, figures, two live interactive demos, and selected source excerpts. The full pipeline (~1.5 MB of Python across 40+ modules) is private. Everything shown here runs on synthetic data — no patient data appears anywhere in this repo.

These are real, unmodified outputs of the pipeline, generated from a synthetic 128-sample cohort. Both are self-contained single-page apps: search, filter, click points, open variant cards.

Pharmacogenomic cohort work has a tooling gap.

Getting from raw sequencing reads to something a researcher can actually interpret means chaining together alignment, several variant callers, annotation, and evidence drawn from a dozen separate knowledge bases — and then repeating the database lookups variant by variant, across an entire cohort, largely by hand. Copy-number variants are frequently left out altogether, because calling them reliably from exome data is hard enough that many workflows simply skip it. In pharmacogenomics that omission matters: CNVs in genes like CYP2D6 change metaboliser status outright.

During a research internship in clinical pharmacogenomics I needed that entire path to exist, to be reproducible, and to run locally. So I built it — a pipeline taking a cohort from FASTQ to interpretation-ready evidence, with CNVs treated as first-class variants and patient data never leaving the machine it was sequenced on.

flowchart TD
    A[FASTQ] -->|Bowtie2 / BWA-MEM2| B[BAM]
    A2[BAM/CRAM] --> B
    B --> C{Variant calling}

    C -->|Exome-wide backbone| D[GATK HaplotypeCaller<br/>+ Strelka2]
    C -->|PGx genes| E[DeepVariant<br/>+ Octopus]
    C -->|Read-depth CNV| F[5-caller CNV ensemble<br/>CNVkit · ExomeDepth · GATK-gCNV<br/>CODEX2 · panelcn.MOPS]

    D --> G[PASS / REVIEW / FAIL tiering]
    E --> G
    G --> H[ML triage model<br/>calibrated gradient-boosted trees]
    H --> I[Joint cohort VCF]
    F --> I

    V1[Raw VCF] --> I
    V2[Filtered VCF] --> I

    I --> J[VEP annotation<br/>offline cache + hit cache]
    J --> K[Database enrichment<br/>ClinPGX · CPIC · PharmVar · DPWG<br/>ClinVar · gnomAD · GeneBe]
    K --> L[MegaVCF<br/>whole cohort, one file]
    K --> M[Cohort analytics<br/>PCA · IBS · ROH · HWE · LD<br/>PLINK2 Firth association]

    L --> N[MegaVCF Explorer]
    M --> O[Analytics Dashboard]
    L --> P[Per-patient reports · BED/TSV exports]

Four entry points — FASTQ, BAM/CRAM, raw VCF, or filtered VCF — so a lab that already has aligned data or existing calls can jump straight to annotation and analytics.

Full detail in docs/ARCHITECTURE.md and docs/METHODS.md.

Two independent validations: a synthetic cohort with a known truth set, and NIST/GIAB HG002. Full breakdown in docs/RESULTS.md.

F1 on GIAB HG002, whole-exome:

The 2 missed SNVs sat inside a homozygous deletion — correct non-calls. The missed CNV was a heterozygous duplication, the hardest class to detect from WES read depth.

Zero false positives in both validations is the design goal, not a coincidence: clinical research conclusions must never rest on variants that don't exist.

Ensemble calling alone gives 0.907. The triage model takes it to 0.929.

AUROC 0.80 with AUPRC 0.97 is the interesting shape: the model's ranking is decent, but the variants it actually promotes are almost never artifacts. That's exactly the behaviour you want when the cost of a false positive is a wrong biological conclusion.

WES CNV calling is genuinely hard: you only see captured exons, so a deletion starting in an intron has invisible edges. Single callers disagree constantly.

The ensemble (design notes):

1. Merge within-caller first. Fragmented adjacent bins from one caller collapse into one event before voting — otherwise a high-resolution caller casts five votes for what is biologically one CNV and dominates the result.
2. Cluster across callers by sample, chromosome, direction, interval overlap, and shared gene/target evidence.
3. One vote per caller per event. 5/5 → consensus PASS. 3–4/5 → REVIEW. <3 → dropped.
4. Report breakpoint uncertainty instead of faking precision — reported interval, target-confirmed core, outer supported span, affected exons, and an explicit precision flag are all separate fields.

Deliberately mixing bin-based callers (ExomeDepth, CODEX2, panelcn.MOPS at 250 bp) with target-based ones (CNVkit, GATK-gCNV, sometimes whole-gene) means the two failure modes don't correlate.

Population frequency (gnomAD-CNV, DGV, DECIPHER) is attached as annotation only and never filters a CNV out — a common CNV can still differ between responders and non-responders.

→ code-excerpts/review_triage_model.py — complete, runnable, ~300 lines.

python code-excerpts/review_triage_model.py --demo

Hard quality thresholds throw away real variants. GATK's own filters cap out at F1 0.902 on GIAB — many of the calls they discard are genuine, just borderline. That's especially damaging in pharmacogenes, which sit next to pseudogenes and near-identical paralogs and therefore always look borderline.

So instead of one hard cut, calls are tiered PASS / REVIEW / FAIL, and a model re-scores the REVIEW tier into P(real):

• Features: QUAL, QD, FS, SOR, MQ, MQRankSum, ReadPosRankSum, BaseQRankSum, DP, GQ, allele balance from AD, is_indel, indel length.
• Model: HistGradientBoostingClassifier(max_depth=3, learning_rate=0.05, max_iter=300, l2_regularization=1.0) wrapped in CalibratedClassifierCV(method="isotonic") — calibration matters because the thresholds are probabilities, not scores.
• Labels come free: every call in a GIAB run is TP/FP-labelable against the NIST truth set.
• Decision: P ≥ 0.70 → rescue to PASS. P < 0.30 → drop to FAIL. Otherwise stay REVIEW.
• Validation splits by whole chromosome, not randomly, so nearby variants can't leak between train and test. The rescue threshold is re-tuned on the training fold each time and applied blind to the held-out fold — no threshold-peeking.

It's measured against two deliberately dumb baselines — keep-all-REVIEW and drop-all-REVIEW — because "the model helps" only means something relative to doing nothing.

The point of all this is that a researcher opens one file and has everything.

MegaVCF joins every variant across the whole cohort — SNVs, indels, and symbolic CNVs (<DEL>/<DUP>) — into a single VCF, then renders a self-contained HTML explorer over it. No server, no install; it opens in a browser and works offline.

Per variant you get consequence and impact, ACMG classification with an editor that recomputes the class live when you add evidence, CPIC/DPWG/PharmGKB drug guidance, PharmVar star alleles, gnomAD frequencies by population, ClinVar significance, which callers supported the call, and which samples carry it. Reviewer comments and filtered exports are built in.

The analytics dashboard covers population structure (PCA, IBS, ROH), QC (heterozygosity, call rate, depth and breadth at 1×/10×/20×, Ti/Tv, HWE, allele-frequency spectrum), association (PLINK2 logistic regression with Firth penalisation → Manhattan, QQ, forest, gene-signal), and LD plus a PGx-readiness heatmap showing which pharmacogenes are actually covered well enough to trust.

Interactive plots are Plotly; publication-ready static figures are matplotlib; statistics are SciPy and PLINK2. More figures →

Language Python 3.11 · pandas · NumPy · scikit-learn · SciPy · Plotly · matplotlib · PyTorch (optional)

Genomics Bowtie2 · BWA-MEM2 · samtools/bcftools · GATK4 · Strelka2 · DeepVariant · Octopus · CNVkit · ExomeDepth · GATK-gCNV · CODEX2 · panelcn.MOPS · Ensembl VEP · PLINK2 · PyPGx

Knowledge bases ClinPGX · CPIC · PharmGKB · PharmVar · DPWG · ClinVar · gnomAD · GeneBe · DGV · DECIPHER · OMIM · Reactome · DrugBank

Delivery local web UI (server + single-page app) · conda environment · Docker · per-run manifests with tool versions and content hashes for reproducibility · stage checkpointing with resume/retry · preflight validation that fails in seconds instead of four hours in

README.md                    you are here
docs/
  ARCHITECTURE.md            pipeline stages, CNV ensemble design, module layout
  METHODS.md                 calling, tiering, annotation, enrichment, analytics
  RESULTS.md                 full validation numbers and how they were produced
  PUBLISHING.md              how to enable Pages and refresh the demos
  demo/                      the two live interactive demos + figure dashboard
figures/                     13 annotated screenshots (synthetic cohort)
code-excerpts/
  review_triage_model.py     the ML triage model, verbatim and runnable
  README.md                  what's excerpted and what isn't
scripts/refresh_demo.sh      regenerate docs/demo from a completed pipeline run

This is research-use-only software, and it isn't finished.

• No clinical validation has been performed. Nothing here is a diagnostic device, and every output needs expert review before it informs anything.
• WES CNV breakpoints are read-depth estimates. Single-exon events can be real but need caution.
• Small cohorts give weak CNV reference structure.
• Association results are research signals — ancestry structure, capture bias, missingness and multiple testing all have to be weighed.
• Star-allele calling is limited by WES capture, paralogy and phasing.

In progress: star-allele/haplotype calling validated against CDC GeT-RM consensus; a deep learning model to predict function and dosing for novel pharmacogene variants (known ones just look up CPIC); and extending toward multi-omics, starting with transcriptomics.

The showcase content in this repository is © 2026 Bruno Young de Castro — see LICENSE. The full pipeline source is private and not licensed for redistribution.

Happy to walk through the architecture, the validation design, or the parts I'd rebuild differently.

Bruno Young de Castro — LinkedIn · brunoyc@icloud.com

Developed during a research internship in clinical pharmacogenomics, as part of the M1 GENIOMHE programme at Université Paris-Saclay / Université Évry.