CYP2C19 and SLCO1B1 From a Raw DNA File in Python

The Short Answer

You can go from a consumer raw DNA file to a CPIC phenotype for CYP2C19 and SLCO1B1 in about sixty lines of Python, with no dependencies beyond requests. The file gives you genotypes at a handful of positions; a small table turns them into star alleles; a rule turns star alleles into a phenotype; and a cited source tells you what that phenotype means for a drug.

TL;DR: Read the file into a {rsid: genotype} dict, handling both the 23andMe and AncestryDNA layouts and their no-call markers. Count the variant alleles at rs4244285 (*2), rs4986893 (*3) and rs12248560 (*17) for CYP2C19, and at rs4149056 for SLCO1B1. Map the counts to CPIC phenotypes. Fetch the guidance from a source that cites CPIC instead of hard-coding it — the code below uses DeepDNA's free knowledge endpoints, which need no key. Then write down, in the output, everything an array cannot tell you: phase, rare alleles and whether a position was genotyped at all.

These two genes are a good first project because they are the friendliest pharmacogenes to read from a SNP array. Neither has the copy-number problems that make CYP2D6 unreliable on a chip, their most important alleles are tagged by single SNPs, and both have peer-reviewed CPIC guidelines behind them: clopidogrel for CYP2C19, statins for SLCO1B1.

Everything here is for building and learning. A raw file is not a clinical test, and the output of this script is not a prescribing decision; the last sections are about exactly why.

What Is Actually in a Raw DNA File

A raw DNA file from a consumer genotyping service is a plain-text table with one row per position on the chip — roughly 600,000 to 700,000 rows. Two layouts cover most of what you will meet:

23andMe — tab-separated, four columns, comment lines starting with #:

# rsid	chromosome	position	genotype
rs4244285	10	96541616	GA
rs12248560	10	96521657	CC
rs4149056	12	21331549	TT
i6010053	1	1234567	--

AncestryDNA — tab-separated, five columns, a header row, and the two alleles in separate columns:

rsid	chromosome	position	allele1	allele2
rs4244285	10	96541616	G	A
rs4149056	12	21331549	T	C

Four details break naive parsers, and every one of them fails silently:

• No-calls are not empty. 23andMe writes --; AncestryDNA writes 0 in each allele column. A parser that treats -- as a genotype will happily count zero variant alleles and call the person a normal metabolizer.
• Not every row is an rsID. 23andMe mixes in internal identifiers starting with i. Looking things up by rsID simply skips them, which is what you want here.
• Coordinates are GRCh37. Both layouts report positions on the older genome build. If you join on position against a GRCh38 source, you will match the wrong base. Joining on rsID sidesteps that — which is why the code below never touches the position column.
• Genotypes are reported on the forward (plus) strand of the reference. For the four SNPs in this article that is convenient: CYP2C19 and SLCO1B1 are both annotated on the plus strand of chromosomes 10 and 12 (check it in Ensembl or NCBI Gene), so the alleles in the file are the same letters as the gene-level notation (G>A for *2, T>C for rs4149056). For genes on the minus strand you would have to complement the bases before comparing. Check the strand of every gene you add.

MyHeritage exports a comma-separated file with quoted fields; the same logic applies once you switch the delimiter and strip the quotes.

Parse It in Fifteen Lines

NO_CALLS = {"--", "00", "0", ""}

def read_raw(path):
    """Return {rsid: genotype} from a 23andMe or AncestryDNA raw file."""
    genotypes = {}
    with open(path) as f:
        for line in f:
            if line.startswith("#") or line.lower().startswith("rsid"):
                continue
            cols = line.rstrip("\r\n").split("\t")
            if len(cols) == 4:                       # 23andMe
                rsid, genotype = cols[0], cols[3]
            elif len(cols) == 5:                     # AncestryDNA
                rsid, genotype = cols[0], cols[3] + cols[4]
            else:
                continue
            if genotype not in NO_CALLS and "0" not in genotype:
                genotypes[rsid] = genotype.upper()
    return genotypes

Two design choices matter more than they look. No-calls are dropped, not stored, so "the position was not genotyped" and "the position was not on the chip" both come out as a missing key — and the caller has to deal with a missing key explicitly. And there is no pandas: for four SNPs out of 700,000 rows, a dict built in one pass is simpler and fast enough.

From Genotypes to Star Alleles

Pharmacogenes are described in star alleles: named haplotypes, where *1 is the reference and each other number is a defined combination of variants. For CYP2C19 on an array, three SNPs do almost all the work:

*2 is the most common loss-of-function allele worldwide — around 15% of European chromosomes and roughly 30% of East Asian ones. *3 is rare in Europeans but common in East Asians. *17 is carried on about a fifth of European chromosomes.

CYP2C19_MARKERS = {
    "rs4244285":  ("A", "*2",  "no_function"),
    "rs4986893":  ("A", "*3",  "no_function"),
    "rs12248560": ("T", "*17", "increased"),
}

def call_cyp2c19(genotypes):
    """Return (diplotype, phenotype, warnings)."""
    alleles, warnings = [], []
    for rsid, (alt, star, function) in CYP2C19_MARKERS.items():
        genotype = genotypes.get(rsid)
        if genotype is None:
            warnings.append(f"{rsid} ({star}) not genotyped: cannot rule out {star}")
            continue
        alleles += [(star, function)] * genotype.count(alt)

    if len(alleles) > 2:
        return None, "Indeterminate", warnings + ["more than two variant alleles: needs phasing"]

    no_function = sum(1 for _, f in alleles if f == "no_function")
    increased = sum(1 for _, f in alleles if f == "increased")
    stars = [s for s, _ in alleles] + ["*1"] * (2 - len(alleles))

    if no_function == 2:
        phenotype = "Poor metabolizer"
    elif no_function == 1:
        phenotype = "Intermediate metabolizer"      # *1/*2, and also *2/*17
    elif increased == 2:
        phenotype = "Ultrarapid metabolizer"
    elif increased == 1:
        phenotype = "Rapid metabolizer"
    else:
        phenotype = "Normal metabolizer"
    return "/".join(sorted(stars, key=lambda s: int(s[1:]))), phenotype, warnings

The phenotype rules follow CPIC's 2022 CYP2C19 terms: two no-function alleles is a poor metabolizer; one no-function allele is an intermediate metabolizer **whether its partner is 1 or 17; two *17 alleles is ultrarapid; one *17 with a *1 is rapid. The *2/*17 case is the one people most often get wrong — the increased-function allele does not cancel the no-function one.

Notice what the function does when a marker is missing: it still returns a call, but with a warning attached. That is deliberate. rs4244285 is on almost every consumer chip; rs4986893 is on some. Refusing to answer whenever *3 is missing would make the script useless for most European files, but calling someone *1/*1 without saying that *3 was never checked would be a silent error. The warning has to travel with the result.

SLCO1B1 and Statins

SLCO1B1 encodes OATP1B1, the transporter that moves statins into the liver. The key variant is rs4149056 (c.521T>C, the marker for the *5 and *15 haplotypes): the C allele reduces transport, so more of the statin stays in the blood and reaches muscle. It is the reason simvastatin myopathy has a genetic component — in the SEARCH trial, people with two C alleles on simvastatin 80 mg had roughly seventeen times the odds of myopathy.

def call_slco1b1(genotypes):
    genotype = genotypes.get("rs4149056")
    if genotype is None:
        return None, "Indeterminate", ["rs4149056 not genotyped"]
    phenotype = {0: "Normal function", 1: "Decreased function", 2: "Poor function"}
    return genotype, phenotype[genotype.count("C")], []

This is a simplification of CPIC's 2022 statin guideline, which defines SLCO1B1 function from full star-allele haplotypes; rs4149056 is the variant that carries the clinically important signal and the one a consumer chip reliably includes.

Fetch the Guidance Instead of Hard-Coding It

The tempting next step is a dictionary that maps "Poor metabolizer" to "avoid clopidogrel". Don't write it. Guidelines are revised — CPIC updated clopidogrel in 2022 and statins in 2022 — and a hard-coded recommendation has no source, no date and no way of telling you it went stale.

Fetch it from somewhere that cites its source. DeepDNA's knowledge endpoints return the CPIC guidance for a gene–drug pair, the phenotype definitions and the references, free and without a key:

import requests

API = "https://deepdna.ai/api/v1"

def guidance(gene, drug):
    r = requests.get(f"{API}/pgx/{gene}/{drug}", timeout=10)
    r.raise_for_status()
    body = r.json()
    return body["pair"], body["meta"]["disclaimer"]

pair, disclaimer = guidance("CYP2C19", "clopidogrel")
print(pair["guideline"], "—", pair["summary"])
for source in pair["sources"]:
    print("  source:", source["label"], source["url"] or "")

The response for clopidogrel carries the CPIC 2022 recommendation (an alternative antiplatelet for poor and intermediate metabolizers undergoing PCI for acute coronary syndromes), the FDA boxed warning, and links to CPIC, PharmGKB and PharmVar. The level field is null for now: the CPIC evidence level is not yet populated, so don't branch on it. For one variant rather than a drug pair, GET /variants/rs4149056 returns the alleles, the effect and the population frequencies, with sources. The full list of genes is at GET /genes, and every endpoint is described in the OpenAPI spec.

Put It Together

genotypes = read_raw("genome.txt")

diplotype, phenotype, warnings = call_cyp2c19(genotypes)
print(f"CYP2C19 {diplotype}: {phenotype}")
for w in warnings:
    print("  warning:", w)

genotype, function, warnings = call_slco1b1(genotypes)
print(f"SLCO1B1 rs4149056 {genotype}: {function}")

pair, disclaimer = guidance("SLCO1B1", "simvastatin")
print(pair["summary"])
print(disclaimer)

For a file carrying GA at rs4244285, CT at rs12248560 and TC at rs4149056, that prints a *2/*17 intermediate metabolizer, a decreased-function SLCO1B1 result, the simvastatin guidance, and a disclaimer that says, correctly, that none of this is medical advice.

What an Array Cannot Tell You

Everything above is correct as code and incomplete as genetics, and the output should say so. Four limits are worth printing next to every result:

1. Phase. A genotype is two letters with no record of which chromosome each came from. GA at rs4244285 plus CT at rs12248560 is almost always *2 on one copy and *17 on the other, which is how the code reads it. But the file cannot prove it, and a clinical call would.
2. Rare alleles. CYP2C19 has dozens of defined star alleles; a chip tags a few. A "normal metabolizer" call really means "none of the alleles we looked for" — which is why CPIC-aware laboratories genotype more positions.
3. Missing markers. If rs4986893 is not on the chip, *3 is not ruled out. In East Asian ancestry that matters.
4. Genes you should not try this with. CYP2D6 depends on whole-gene deletions, duplications and hybrids with a neighbouring pseudogene, none of which a SNP array can see. A CYP2D6 phenotype from a raw file is closer to a guess; it belongs in a different article.

There is also a limit that is not technical. A decision about clopidogrel after a stent, or about which statin to start, is made by a clinician with a validated test. What this script produces is a reason to have that conversation, not a replacement for it.

Wiring It Into an Agent

If you are building an agent rather than a script, the split of work is the same and it matters more. Let the agent keep the conversation and the context; let tools do the lookups. The knowledge endpoints above are already described in OpenAPI, so most agent frameworks can turn them into tools directly — a get_pgx_guidance(gene, drug) tool that the model calls instead of recalling CPIC from memory. There is no DeepDNA MCP server published yet; OpenAPI is the way in today.

The part that is harder to do safely inside an agent is the first step: reading the file itself. DeepDNA's POST /dna/parse and POST /dna/annotate are designed to do that — detect the layout, drop no-calls, state the genome build and return annotated variants with sources — and they are in private beta, not yet deployed. Their contract is in the same OpenAPI spec, so you can build against it now. If you are building something that reads genetic data, tell us what; a person reads every request.

The knowledge endpoints are live today. The API docs have the full list, and the gene pages for CYP2C19 and SLCO1B1 show the same records in a human-readable form. For the clinical background, our pharmacogenomics guide covers the guidelines in more depth.

Frequently Asked Questions

Can I get my CYP2C19 metabolizer status from a 23andMe raw file?

Usually, as an estimate. The *2 marker (rs4244285) is on nearly every 23andMe and AncestryDNA chip and the *17 marker (rs12248560) on most, so you can map the genotypes to a CPIC phenotype. It is not a clinical result: phase is unknown, rare alleles are not tested, and *3 may be missing. For a prescribing decision, a clinician should order a validated pharmacogenomic test.

Is *2/*17 an intermediate or a rapid metabolizer?

Under CPIC's 2022 terms, *2/*17 is an intermediate metabolizer. One no-function allele sets the phenotype to intermediate regardless of whether the other allele is normal (*1) or increased function (*17).

Why join on rsID instead of chromosome and position?

Because consumer files report positions on GRCh37 and many current annotation sources use GRCh38. Joining on position across builds matches the wrong base without any error. The rsID is stable across builds for the SNPs used here.

Can I do the same for CYP2D6?

Not reliably. CYP2D6 function depends on gene deletions, duplications and hybrid genes with a neighbouring pseudogene, which a SNP array cannot detect. Any CYP2D6 phenotype from a raw file should be treated as unknown.

Does the DeepDNA API need a key?

The knowledge endpoints (/genes, /variants, /pgx) are free and need no key. The endpoints that read whole DNA files are in private beta and not yet deployed.

This article is for educational purposes only and does not constitute medical advice. Pharmacogenomic results from consumer raw data are not clinically validated; decisions about medication should be made with a qualified clinician using a validated test.

This article was created with AI assistance and reviewed by the DeepDNA editorial team.