""" Pharmacogenomic star-allele calling from 1000 Genomes (GRCh38, phased). Fetches only the pharmacogene variant positions over HTTPS (no full download), calls diplotypes, assigns CPIC phenotypes, and writes tidy CSVs for the notebook. """ import json, urllib.request, sys from pathlib import Path import pysam import pandas as pd OUT = Path(__file__).parent / "pgx_out" OUT.mkdir(exist_ok=True) 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") PANEL = ("https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/" "integrated_call_samples_v3.20130502.ALL.panel") # Star alleles defined by single SNVs. Each entry: the rsID and which allele # marks the variant haplotype. Deliberately restricted to alleles that CAN be # called from biallelic short-read SNVs -- see NOT_CALLABLE below. MARKERS = { "CYP2C19": {"rs4244285": "*2", "rs4986893": "*3", "rs12248560": "*17"}, "CYP2C9": {"rs1799853": "*2", "rs1057910": "*3"}, "TPMT": {"rs1800462": "*2", "rs1800460": "*3B", "rs1142345": "*3C"}, "DPYD": {"rs3918290": "c.1905+1G>A", "rs67376798": "c.2846A>T"}, "SLCO1B1": {"rs4149056": "*5"}, "VKORC1": {"rs9923231": "-1639G>A"}, } ALL_RS = [r for g in MARKERS.values() for r in g] def resolve_coords(rsids): """GRCh38 coordinates from Ensembl -- never hardcoded from memory.""" cache = OUT / "coords.json" if cache.exists(): return json.loads(cache.read_text()) 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) if not v: raise SystemExit(f"Ensembl did not return {rs}") 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"]), "alleles": m["allele_string"], "consequence": v.get("most_severe_consequence", "")} cache.write_text(json.dumps(out, indent=2)) return out def load_panel(): cache = OUT / "panel.csv" if cache.exists(): return pd.read_csv(cache) # The published panel has trailing empty fields in its header row, so let # pandas infer rather than zipping header names to columns by hand. df = pd.read_csv(PANEL, sep="\t") df = df[[c for c in df.columns if not str(c).startswith("Unnamed")]] df = df[["sample", "pop", "super_pop"]] df.to_csv(cache, index=False) return df def fetch_genotypes(coords): """Pull the haplotype pair at each marker for every sample.""" cache = OUT / "haplotypes.csv" if cache.exists(): return pd.read_csv(cache) by_chrom = {} for rs, c in coords.items(): by_chrom.setdefault(c["chrom"], []).append(rs) frames, meta = [], [] for chrom, rsids in sorted(by_chrom.items()): print(f" chr{chrom}: {len(rsids)} marker(s)", flush=True) vf = pysam.VariantFile(VCF.format(c=chrom)) samples = list(vf.header.samples) for rs in rsids: pos = coords[rs]["pos"] rec = None for r in vf.fetch(chrom, pos - 1, pos): if r.id == rs or r.pos == pos: rec = r break if rec is None: print(f" !! {rs} not found at {chrom}:{pos}", flush=True) continue meta.append({"rsid": rs, "chrom": chrom, "pos": rec.pos, "vcf_id": rec.id, "ref": rec.ref, "alt": rec.alts[0]}) # haplotype 1 / haplotype 2, as 0 (ref) or 1 (alt) h1, h2, phased = [], [], [] for s in samples: gt = rec.samples[s]["GT"] h1.append(gt[0]); h2.append(gt[1]) phased.append(bool(rec.samples[s].phased)) frames.append(pd.DataFrame({"sample": samples, f"{rs}_h1": h1, f"{rs}_h2": h2})) print(f" {rs} -> {chrom}:{rec.pos} {rec.ref}>{rec.alts[0]} " f"(phased: {sum(phased)}/{len(phased)})", flush=True) df = frames[0] for f in frames[1:]: df = df.merge(f, on="sample") df.to_csv(cache, index=False) pd.DataFrame(meta).to_csv(OUT / "markers.csv", index=False) return df # ---------------------------------------------------------------- calling --- def call_haplotype(gene, hap_alleles): """Name one haplotype from its marker alleles. hap_alleles: {rsid: 0|1}. 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": # *2/*3 are loss-of-function and dominate a haplotype carrying *17 for a in ("*3", "*2", "*17"): if a in var: return a return sorted(var)[0] CYP2C19_PHENO = { frozenset(["*17", "*17"]): "Ultrarapid metaboliser", frozenset(["*1", "*17"]): "Rapid metaboliser", frozenset(["*1", "*1"]): "Normal metaboliser", } CYP2C19_LOF = {"*2", "*3"} 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": 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: score = ACTIVITY[gene].get(d1, 1.0) + ACTIVITY[gene].get(d2, 1.0) if gene == "CYP2C9": ph = ("Normal metaboliser" if score == 2.0 else "Intermediate metaboliser" if score >= 1.0 else "Poor metaboliser") else: # DPYD ph = ("Normal metaboliser" if score == 2.0 else "Intermediate metaboliser" if score >= 1.0 else "Poor metaboliser") return ph, score if gene == "TPMT": nf = sum(1 for a in (d1, d2) if a != "*1") return ("Normal metaboliser" if nf == 0 else "Intermediate metaboliser" if nf == 1 else "Poor metaboliser"), None if gene == "SLCO1B1": n = sum(1 for a in (d1, d2) if a == "*5") return ("Normal function" if n == 0 else "Decreased function" if n == 1 else "Poor function"), None if gene == "VKORC1": n = sum(1 for a in (d1, d2) if a != "*1") return ("Low warfarin sensitivity" if n == 0 else "Intermediate warfarin sensitivity" if n == 1 else "High warfarin sensitivity"), None return "Indeterminate", None def main(): print("1. Resolving coordinates from Ensembl ...", flush=True) coords = resolve_coords(ALL_RS) print("2. Loading 1000G population panel ...", flush=True) panel = load_panel() print("3. Fetching genotypes from 1000 Genomes over HTTPS ...", flush=True) hap = fetch_genotypes(coords) print("4. Calling diplotypes ...", flush=True) rows = [] for _, r in hap.iterrows(): rec = {"sample": r["sample"]} for gene, markers in MARKERS.items(): present = [rs for rs in markers if f"{rs}_h1" in hap.columns] if not present: continue a1 = call_haplotype(gene, {rs: r[f"{rs}_h1"] for rs in present}) a2 = call_haplotype(gene, {rs: r[f"{rs}_h2"] for rs in present}) dip = "/".join(sorted([a1, a2])) ph, score = phenotype(gene, a1, a2) rec[f"{gene}_diplotype"] = dip rec[f"{gene}_phenotype"] = ph if score is not None: rec[f"{gene}_activity_score"] = score rows.append(rec) calls = pd.DataFrame(rows).merge(panel, on="sample", how="left") calls = calls[calls["super_pop"].notna()] calls.to_csv(OUT / "calls.csv", index=False) print(f" {len(calls)} samples with population labels", flush=True) for gene in MARKERS: col = f"{gene}_phenotype" if col not in calls.columns: continue t = (calls.groupby(["super_pop", col]).size() .unstack(fill_value=0)) t = (t.T / t.T.sum()).T * 100 t.round(2).to_csv(OUT / f"freq_{gene}.csv") print(f"\n{gene} phenotype % by superpopulation") print(t.round(1).to_string()) print("\nDone ->", OUT, flush=True) if __name__ == "__main__": main()