78% of 2,503 people carry a gene variant that would change their prescription.

def resolve_coords(rsids):
    """GRCh38 coordinates from Ensembl -- never hardcoded from memory."""
    req = urllib.request.Request(
        "https://rest.ensembl.org/variation/homo_sapiens?content-type=application/json",
        data=json.dumps({"ids": rsids}).encode(),
        headers={"Content-Type": "application/json", "Accept": "application/json"})
    data = json.load(urllib.request.urlopen(req, timeout=120))
    out = {}
    for rs in rsids:
        v = data.get(rs)
        # drop patch/alt scaffolds; keep the primary GRCh38 assembly only
        m = [x for x in v["mappings"]
             if x["assembly_name"] == "GRCh38"
             and "_" not in str(x["seq_region_name"])][0]
        out[rs] = {"chrom": str(m["seq_region_name"]), "pos": int(m["start"])}
    return out

VCF = ("https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/"
       "1000_genomes_project/release/20190312_biallelic_SNV_and_INDEL/"
       "ALL.chr{c}.shapeit2_integrated_snvindels_v2a_27022019.GRCh38.phased.vcf.gz")

vf = pysam.VariantFile(VCF.format(c=chrom))   # reads the remote .tbi index only
for r in vf.fetch(chrom, pos - 1, pos):
    if r.id == rs or r.pos == pos:
        rec = r
        break

# keep the two haplotypes separately -- see step 04 for why
h1, h2 = [], []
for s in samples:
    gt = rec.samples[s]["GT"]
    h1.append(gt[0]); h2.append(gt[1])

def call_haplotype(gene, hap_alleles):
    """Name one haplotype from its marker alleles.
    TPMT*3A = *3B and *3C in cis -- which only phased data can
    distinguish from *3B/*3C in trans."""
    var = {MARKERS[gene][rs] for rs, a in hap_alleles.items() if a == 1}
    if not var:
        return "*1"
    if gene == "TPMT":
        if {"*3B", "*3C"} <= var:
            return "*3A"
        return sorted(var)[0]
    if gene == "CYP2C19":
        # loss-of-function dominates a haplotype that also carries *17
        for a in ("*3", "*2", "*17"):
            if a in var:
                return a
    return sorted(var)[0]

ACTIVITY = {
    "CYP2C9": {"*1": 1.0, "*2": 0.5, "*3": 0.0},
    "DPYD":   {"*1": 1.0, "c.1905+1G>A": 0.0, "c.2846A>T": 0.5},
}

def phenotype(gene, d1, d2):
    if gene == "CYP2C19":                  # diplotype lookup, not a score
        lof = [a for a in (d1, d2) if a in CYP2C19_LOF]
        if len(lof) == 2: return "Poor metaboliser", None
        if len(lof) == 1: return "Intermediate metaboliser", None
        return CYP2C19_PHENO.get(frozenset([d1, d2]), "Normal metaboliser"), None
    if gene in ACTIVITY:                   # activity score
        score = ACTIVITY[gene].get(d1, 1.0) + ACTIVITY[gene].get(d2, 1.0)
        ...