How to Call and Read Your Own Pharmacogenetic Results

Oak
pharmacogenomics pharmcat cyp2d6 whole-genome-sequencing cpic
Two parallel glowing vines on a dark tree trunk in a bioluminescent forest, one with a duplicated cluster of light-beads.

At the end of this you will have a directory containing a per-gene table of star-allele diplotypes (for example CYP2C19 *1/*2), a CPIC-standard phenotype and activity score for each, an HTML report listing the gene-drug pairs with published guidelines, and a separate high-confidence CYP2D6 call including copy number and hybrid alleles. You need aligned reads from whole-genome sequencing (a CRAM or BAM against GRCh38, 30x or better), the matching reference FASTA, Python 3.10+, Java 17, and about 20 minutes of compute. Exome data will work for most genes and will fail on the intronic and structural variation that defines several CYP2D6 alleles. A consumer SNP array will not work: array designs miss most of the low-frequency alleles and cannot distinguish a reference call from a no-call, which is the single most common way pharmacogenetic interpretation goes wrong.

Nothing below is medical advice. Diplotypes and phenotypes are measurements. Any decision about starting, stopping, or changing a medication belongs to a prescribing clinician or a pharmacist, and you should bring them the report rather than act on it.

1. Understand what you are calling before you run anything

Pharmacogenetic genes are not reported as individual variants. They are reported as haplotypes, called star alleles, defined by a specific combination of positions on the same chromosome copy. CYP2C19*2 is rs4244285 (c.681G>A, a splice defect). CYP2C19*17 is a promoter variant that increases expression. CYP2D6*4 is a splice variant, *5 is a whole-gene deletion, *10 and *41 are reduced-function haplotypes, and *1xN is a duplication of a functional copy.

Your two haplotypes together form a diplotype, and the diplotype maps to a phenotype through a lookup table maintained by CPIC. For CYP2D6 the mapping goes through an activity score: each allele is assigned a value (fully functional alleles like *1 and *2 score 1, no-function alleles like *4 and *5 score 0, reduced-function alleles score 0.25 or 0.5), the two are summed, and duplications multiply. Under the 2019 CPIC/DPWG consensus, a score of 0 is a poor metabolizer, 0.25 to 1.0 is an intermediate metabolizer, 1.25 to 2.25 is normal, and above 2.25 is ultrarapid. Roughly 5 to 10 percent of people of European ancestry are CYP2D6 poor metabolizers, and allele frequencies vary substantially by ancestry, which matters later.

“Poor metabolizer” describes enzyme activity, not a whole-body property. For a drug cleared by that enzyme, less clearance means higher exposure at the same dose. For a prodrug activated by that enzyme (codeine to morphine, tamoxifen to endoxifen, clopidogrel via CYP2C19), the same genotype means less active drug. The direction of the effect depends entirely on the drug, which is why there is no such thing as a general list of “drugs to avoid” that you can read off a genotype without the pharmacology attached.

2. Produce a VCF that includes reference calls

This is the step people skip, and it silently corrupts everything downstream. A standard variant-only VCF from a joint-calling pipeline contains rows only where you differ from the reference. A pharmacogenetic caller reading that file cannot tell whether a missing position is reference (you do not carry the variant) or simply uncovered (unknown). PharmCAT will treat absent positions as no-calls and refuse to assign a diplotype for the affected gene.

Two ways to fix it. If you have a gVCF from DeepVariant or GATK, use it directly. If you only have a variant-only VCF, regenerate calls across the PharmCAT position set from your CRAM:

# PharmCAT ships the exact positions it needs
wget https://github.com/PharmGKB/PharmCAT/releases/latest/download/pharmcat_positions.vcf.bgz
tabix -p vcf pharmcat_positions.vcf.bgz

bcftools mpileup \
  -f GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -R pharmcat_positions.vcf.bgz \
  -a FORMAT/AD,FORMAT/DP \
  -q 20 -Q 20 --max-depth 500 \
  -Ou sample.cram \
| bcftools call -m -Oz -o sample.pgx.vcf.gz
tabix -p vcf sample.pgx.vcf.gz

Then run the preprocessor, which normalizes representation, splits multiallelics, and checks your contig naming and build:

pharmcat_vcf_preprocessor \
  -vcf sample.pgx.vcf.gz \
  -refFna GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -refVcf pharmcat_positions.vcf.bgz \
  -o preprocessed/

Read the *.missing_pgx_var.vcf file it emits. It lists positions PharmCAT wanted and did not find. If a gene has a long list there, treat its call as provisional.

3. Run PharmCAT end to end

pharmcat_pipeline preprocessed/sample.preprocessed.vcf \
  -o results/ \
  -reporterJson \
  --research-mode CYP2D6

You get three artifacts per sample: sample.match.json (which haplotype combinations are consistent with your data, including ties), sample.phenotype.json (diplotype to activity score to phenotype), and sample.report.html (a readable summary organized by drug with CPIC and DPWG recommendation text).

Parse the JSON rather than reading the HTML if you want a table you can version:

import json, csv
d = json.load(open("results/sample.phenotype.json"))
rows = []
for gene, g in d["geneReports"]["CPIC"].items():
    for src in g.get("sourceDiplotypes", []):
        rows.append({
            "gene": gene,
            "diplotype": src.get("label"),
            "phenotype": "; ".join(src.get("phenotypes") or []),
            "activity_score": src.get("activityScore"),
            "called": g.get("called"),
        })
with open("pgx_calls.csv", "w", newline="") as fh:
    w = csv.DictWriter(fh, fieldnames=rows[0].keys()); w.writeheader(); w.writerows(rows)

PharmCAT covers roughly two dozen genes with CPIC guidelines, including CYP2C19, CYP2C9, CYP3A5, CYP2B6, DPYD, TPMT, NUDT15, SLCO1B1, VKORC1, UGT1A1, and G6PD. Its handling of CYP2D6 from a VCF alone is deliberately limited, which is the next step.

4. Call CYP2D6 separately, from reads

CYP2D6 sits next to CYP2D7, a near-identical pseudogene. Deletions, duplications, and CYP2D6/CYP2D7 hybrid alleles (*13, *36, *68) are common, and short-read aligners mismap across the two. Variant callers routinely get this wrong. Use a tool that models copy number from read depth against the CYP2D7 paralog.

Cyrius is the one we would reach for first on Illumina WGS:

echo "/data/sample.cram" > manifest.txt
python star_caller.py \
  --manifest manifest.txt \
  --genome 38 \
  --prefix sample \
  --outDir cyrius_out/ \
  --threads 8 \
  --reference GRCh38_full_analysis_set_plus_decoy_hla.fa

The output TSV gives a genotype string such as *1/*4 or *2x2/*4, plus a Filter column. Treat anything other than PASS as needing a second opinion.

Cross-check with a second caller. PyPGx handles 80+ genes and produces copy-number plots you can eyeball:

pypgx run-ngs-pipeline CYP2D6 cyp2d6-pipeline \
  --variants cohort.vcf.gz \
  --depth-of-coverage depth.zip \
  --control-statistics control-VDR.zip \
  --assembly GRCh38

Aldy is a third option and uses a different algorithm (integer linear programming over read evidence):

aldy genotype -p wgs -g CYP2D6 --genome hg38 -o aldy_cyp2d6.txt sample.bam

When two of three agree, you have a usable call. When they disagree, the disagreement is almost always about copy number or a hybrid allele, and the output is a range of possible activity scores rather than a single phenotype.

5. Map calls to guidelines, and know which pairs have evidence

A diplotype only matters where there is a gene-drug pair with published, graded evidence. Download the CPIC gene-drug pair list and the PharmGKB clinical annotation levels, and join them to your table:

curl -sL https://api.cpicpgx.org/v1/pair?select=genesymbol,drugname,cpiclevel \
  -H "Accept: application/json" -o cpic_pairs.json

Filter to CPIC level A and A/B. Those are the pairs where a guideline gives a specific prescribing action. Everything at level C or D is interesting biology and not a decision input. The high-value pairs in practice are DPYD with fluoropyrimidines, TPMT and NUDT15 with thiopurines, CYP2C19 with clopidogrel and with proton pump inhibitors, CYP2C9 plus VKORC1 with warfarin, SLCO1B1 with simvastatin, CYP2D6 with codeine and tramadol, and the HLA alleles (HLA-B*57:01, HLA-B*15:02, HLA-A*31:01) with specific hypersensitivity risks. HLA typing needs its own tool (HLA-LA or Optitype on your CRAM), since PharmCAT expects HLA alleles to be supplied as input rather than called from the VCF.

The important framing: a level A result tells a clinician that a standard dose is likely to produce a non-standard exposure. It does not tell you what to do instead. Do not change a dose or stop a medication based on a genotype without a clinician. Abrupt changes to psychiatric, anticoagulant, or antiepileptic medication can be dangerous independent of any genotype.

6. Read a commercial report against your own calls

If you have already paid for a panel (GeneSight, Genomind, OneOme, and similar), pull the diplotype table out of the PDF and diff it against your calls. The diplotypes should largely agree. Where they differ, the usual cause is CYP2D6 copy number or an allele the panel does not genotype. Sequencing detects pharmacogene variation that targeted genotyping panels miss by design, including rare and novel functional variants that never appear on a fixed array 1. Whole-genome data has been used to build full pharmacogenomic profiles across the drug-metabolizing gene set for exactly this reason 2.

What the commercial reports add on top of the diplotypes is a proprietary combinatorial algorithm that buckets drugs into colored tiers. That layer is where we would set the skepticism dial high. The underlying per-gene CPIC calls are reproducible and auditable. The tier assignment is a vendor-specific weighting you cannot inspect, and it rolls pharmacokinetic gene calls together with pharmacodynamic variants of much weaker evidence. A real-world analysis of psychiatric PGx test ordering found that a substantial share of tests were ordered for patients whose medications had no actionable gene-drug pair at all, which limits how much the report can change 3. Returning these results also carries documentation and duty-to-follow-up obligations for the ordering clinician that a direct-to-consumer framing tends to obscure 4.

On the ADHD question, which comes up constantly: no pharmacogenetic panel diagnoses ADHD. Diagnosis is clinical. The genetics of ADHD susceptibility is a separate research question from drug metabolism, involving rare and common variants across many genes, and studies of candidate genes such as BDNF report risk associations, not prescribing information 5. A PGx report can tell you that you carry two no-function CYP2D6 alleles, which is relevant to how certain medications are cleared. It cannot tell you whether you have the condition those medications treat.

7. Store the result so it survives allele-definition updates

Star-allele definitions change. New alleles get named, function assignments get revised, and PharmCAT ships new data releases several times a year. Keep the raw evidence, not just the conclusion.

pgx/
  input/sample.pgx.vcf.gz            # frozen, with reference calls
  runs/2026-09-14_pharmcat-3.0.0/
    sample.match.json
    sample.phenotype.json
    sample.report.html
    cyrius_out/sample.tsv
  pgx_calls.csv                      # joined table, with tool + version columns

Record the tool version and the CPIC data release in every row. Re-run the whole pipeline annually from the frozen VCF and CRAM. It costs minutes. The alternative is a PDF from 2023 with a phenotype that the current lookup table no longer assigns.

8. Hand it to a pharmacist, not to a prescription

The person best equipped to use this is a clinical pharmacist. Pharmacist-led PGx programs in ambulatory settings have been evaluated specifically for feasibility of turning test results into documented, actionable recommendations 6. Bring three things: the diplotype table with the calling tool named, your current and past medication list including over-the-counter drugs and supplements, and any history of a medication that did not work or caused a strong side effect at a standard dose.

Also mention inhibitors and inducers. A person with a normal-metabolizer CYP2D6 genotype who is taking a strong CYP2D6 inhibitor can behave like a poor metabolizer. This is phenoconversion, and it means the genotype sets a ceiling on activity rather than fixing it. Any interpretation that ignores concurrent medications is incomplete.

Common problems

Build mismatch. Coordinates on GRCh37 fed into a GRCh38 tool produce nonsense diplotypes rather than an error. Check ##contig lines and the reference MD5 in your CRAM header before anything else.

Variant-only VCF. Covered above and worth repeating, because it is the most frequent failure. If PharmCAT reports “not called” for most genes, you supplied a file without reference calls.

Chromosome naming. chr7 versus 7. bcftools annotate --rename-chrs fixes it.

Phasing ambiguity. With unphased data, two variants in the same gene can be on the same haplotype or on opposite ones, producing different diplotypes. PharmCAT lists all consistent combinations in match.json. Do not let a downstream script silently take the first one. If phase matters (it does for CYP2C9 and CYP2D6), long reads or trio data resolve it.

CYP2D6 hybrids. *13 and *68 are CYP2D7-derived and depth-based callers can report them as duplications or miss them entirely. If Cyrius filters the call, do not synthesize a phenotype from the raw variants yourself.

Ancestry bias in allele definitions. Star-allele catalogs were built largely from European cohorts. Alleles common in African ancestry populations (CYP2D6*17, *29, *45, *46) and in East Asian populations (*10, *36) are covered by current definitions but are systematically underrepresented on older genotyping panels. Sequencing-based calling narrows this gap because it does not depend on a pre-chosen probe set 1.

DPYD is not a star-allele gene in the usual sense. CPIC scores it by summing activity values for a small set of specific variants, and rare no-function variants matter a great deal clinically. Check for no-calls at the DPYD positions specifically.

Overreading pharmacodynamic genes. SLC6A4, HTR2A, COMT, and MTHFR appear on commercial reports. Their evidence for guiding drug selection is weak relative to the CYP genes, and CPIC does not issue prescribing guidance for most of them. We would report them as “measured, no actionable guideline” and leave it there.

Turnaround expectations. A commercial panel typically takes one to three weeks from sample to report. From a CRAM you already own, the full pipeline above runs in under 30 minutes on a laptop. Producing the numbers is the fast part. The interpretation conversation is the part that takes time and should.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Wenjian Yang, Gang Wu, Ulrich Broeckel, et al. Comparison of genome sequencing and clinical genotyping for pharmacogenes. Clinical Pharmacology & Therapeutics, 2016. https://doi.org/10.1002/cpt.411 ↩ ↩2

  2. Clint Mizzi, Brock A. Peters, Christina Mitropoulou, et al. Personalized Pharmacogenomics Profiling Using Whole-Genome Sequencing. Pharmacogenomics, 2014. https://doi.org/10.2217/pgs.14.102 ↩

  3. Lusi Zhang, Anthony Tholkes, Kaliya C. Jones, et al. Real‐World Characterization of Psychiatric Pharmacogenomic Test Ordering and Clinical Relevance in Adults and Children. Clinical and Translational Science, 2025. https://doi.org/10.1111/cts.70297 ↩

  4. Eric T. Ward, Kristin M. Kostick, Gabriel Lázaro‐Muñoz. Integrating Genomics into Psychiatric Practice: Ethical and Legal Challenges for Clinicians. Harvard Review of Psychiatry, 2019. https://doi.org/10.1097/hrp.0000000000000203 ↩

  5. Ziarih Hawi, Tarrant D.R. Cummins, Janette Tong, et al. Rare DNA variants in the brain-derived neurotrophic factor gene increase risk for attention-deficit hyperactivity disorder: a next-generation sequencing study. Molecular Psychiatry, 2016. https://doi.org/10.1038/mp.2016.117 ↩

  6. J. Oestreich, E. Palisoc, G. Bennett, et al. 122 Evaluating the feasibility and implementation of pharmacist-led pharmacogenomic testing in an ambulatory care setting. Journal of the American Pharmacists Association, 2026. https://doi.org/10.1016/j.japh.2026.103173 ↩