Gwas Database

pip install requests matplotlib numpy

import requests
import time

BASE = "https://www.ebi.ac.uk/gwas/rest/api"

def gwas_get(endpoint, params=None):
    """GWAS Catalog REST API helper with rate limiting and pagination support."""
    url = f"{BASE}/{endpoint}"
    resp = requests.get(url, params=params or {})
    resp.raise_for_status()
    time.sleep(0.2)
    return resp.json()

# Find studies for a trait keyword. Study records have no top-level `title`
# — the publication title lives at `publicationInfo.title`; the trait label
# lives at `diseaseTrait.trait`.
data = gwas_get("studies/search/findByDiseaseTrait", {"diseaseTrait": "diabetes"})
studies = data["_embedded"]["studies"]
print(f"Found {len(studies)} studies for 'diabetes'")
for s in studies[:3]:
    title = (s.get("publicationInfo") or {}).get("title", "N/A")
    trait = (s.get("diseaseTrait") or {}).get("trait", "N/A")
    print(f"  {s['accessionId']} | {trait[:40]:<40} | {title[:60]}")

# Search studies by disease trait
data = gwas_get("studies/search/findByDiseaseTrait", {"diseaseTrait": "breast cancer"})
studies = data["_embedded"]["studies"]
for s in studies[:5]:
    pi = s.get("publicationInfo") or {}
    print(f"  {s['accessionId']} | PMID:{pi.get('pubmedId','N/A')} | {pi.get('title','')[:60]}")

time.sleep(0.2)

# Search by PubMed ID. NOTE: the older `findByPubmedId` 404s on /studies/;
# the working endpoint is `findByPublicationIdPubmedId`.
data = gwas_get("studies/search/findByPublicationIdPubmedId", {"pubmedId": "25673413"})
studies = data["_embedded"]["studies"]
print(f"Studies from PMID 25673413: {len(studies)}")
for s in studies:
    trait = (s.get("diseaseTrait") or {}).get("trait", "N/A")
    print(f"  {s['accessionId']}: {trait}")

# Associations by EFO trait. The old path `efoTraits/{shortForm}/associations`
# also works *if* you have the current shortForm — but trait shortForms have
# been re-mapped to MONDO (e.g. EFO_0000249 → MONDO_0004975). The most reliable
# path is `associations/search/findByEfoTrait?efoTrait=<canonical trait name>`.
data = gwas_get("associations/search/findByEfoTrait",
                {"efoTrait": "type 2 diabetes mellitus", "size": 50})
assocs = data["_embedded"]["associations"]
print(f"Associations for 'type 2 diabetes mellitus': {len(assocs)}")

for a in assocs[:5]:
    pval = a.get("pvalue", None)
    genes = []
    for locus in a.get("loci", []) or []:
        for gene in locus.get("authorReportedGenes", []) or []:
            genes.append(gene.get("geneName", ""))
    loci = a.get("loci") or [{}]
    snps = [r.get("snps", [{}])[0].get("rsId", "N/A")
            for r in (loci[0].get("strongestRiskAlleles") or [])]
    print(f"  rs={snps} | p={pval} | genes={genes}")

# Associations for a specific variant. NOTE: association records do not embed
# `efoTraits` inline — they expose them via the `_links.efoTraits.href`
# HAL link. Follow the link (cached if needed) to resolve trait names.
data = gwas_get("singleNucleotidePolymorphisms/rs7903146/associations", {"size": 5})
assocs = data["_embedded"]["associations"]
print(f"Associations for rs7903146 (first page): {len(assocs)}")

def association_traits(assoc):
    """Resolve efoTraits via the HAL link on an association record."""
    href = (assoc.get("_links") or {}).get("efoTraits", {}).get("href")
    if not href:
        return []
    r = requests.get(href, timeout=15)
    if not r.ok:
        return []
    return [t.get("trait") for t in r.json().get("_embedded", {}).get("efoTraits", [])]

for a in assocs[:5]:
    traits = association_traits(a)
    print(f"  p={a.get('pvalue')} | OR={a.get('orPerCopyNum', 'N/A')} | traits={traits}")
    time.sleep(0.1)

# Lookup single variant
data = gwas_get("singleNucleotidePolymorphisms/rs7903146")
loc = data.get("locations", [{}])[0]
print(f"rs7903146: chr{loc.get('chromosomeName', '?')}:{loc.get('chromosomePosition', '?')}")
print(f"  Functional class: {data.get('functionalClass', 'N/A')}")
print(f"  Merged into: {data.get('merged', 'N/A')}")

time.sleep(0.2)

# Search variants by chromosomal region
data = gwas_get("singleNucleotidePolymorphisms/search/findByChromBpLocationRange",
                {"chrom": "10", "bpStart": "114750000", "bpEnd": "114800000", "size": 50})
snps = data["_embedded"]["singleNucleotidePolymorphisms"]
print(f"Variants in chr10:114750000-114800000: {len(snps)}")
for v in snps[:5]:
    print(f"  {v['rsId']}: {v.get('functionalClass', 'N/A')}")

# Search variants by gene name (the cytogenetic-band endpoint
# `findByCytogeneticBand` was removed — use gene or chromosome-range instead).
data = gwas_get("singleNucleotidePolymorphisms/search/findByGene",
                {"geneName": "TCF7L2", "size": 5})
snps = data["_embedded"]["singleNucleotidePolymorphisms"]
print(f"Variants in TCF7L2: {len(snps)}")

# Search traits by exact name (the older `findByDescription` endpoint was
# removed — search/efoTrait now expects the canonical trait label).
data = gwas_get("efoTraits/search/findByEfoTrait", {"trait": "Alzheimer disease"})
traits = data["_embedded"]["efoTraits"]
print(f"Traits matching 'Alzheimer disease': {len(traits)}")
for t in traits[:5]:
    print(f"  {t['shortForm']}: {t['trait']}  (uri={t['uri']})")

time.sleep(0.2)

# Get specific trait by shortForm. NOTE: many legacy EFO IDs have been
# re-mapped to MONDO (e.g. old `EFO_0000249` for Alzheimer is now
# `MONDO_0004975` — `efoTraits/EFO_0000249` returns 404). Resolve via search
# above first, then use the current shortForm:
short_form = traits[0]["shortForm"]   # e.g. 'MONDO_0004975'
data = gwas_get(f"efoTraits/{short_form}")
print(f"Trait: {data['trait']}")
print(f"  URI       : {data['uri']}")
print(f"  shortForm : {data['shortForm']}")

# List studies with available summary statistics
data = gwas_get("studies/search/findByFullPvalueSet", {"fullPvalueSet": True, "size": 10})
studies = data["_embedded"]["studies"]
print(f"Studies with summary stats (first page): {len(studies)}")
for s in studies[:5]:
    trait = (s.get("diseaseTrait") or {}).get("trait", "N/A")
    print(f"  {s['accessionId']}: {trait[:50]}")

time.sleep(0.2)

# Summary statistics metadata is NOT exposed via the REST API
# (`studies/{acc}/summaryStatistics` returns 404). Use the GWAS Catalog FTP
# directly — paths are predictable by study accession:
ftp_base = "http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics"
acc = studies[0]["accessionId"]
ftp_url = f"{ftp_base}/{acc[:-3]}001-{acc[:-3]}999/{acc}/"
print(f"Summary stats FTP directory: {ftp_url}")

# Download summary statistics file (FTP)
# Summary statistics are hosted on the GWAS Catalog FTP, not the REST API
import urllib.request

study_id = "GCST006867"  # Example study
ftp_base = "http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics"
# Actual paths vary by study; check the study page for the download link
# Example: ftp_base/{study_id}/{study_id}.tsv.gz
url = f"{ftp_base}/{study_id}"
print(f"Summary stats FTP directory: {url}")
# Use requests.get() or urllib to download the .tsv.gz file

# Search associations by gene
data = gwas_get("singleNucleotidePolymorphisms/search/findByGene",
                {"geneName": "BRCA1", "size": 50})
snps = data["_embedded"]["singleNucleotidePolymorphisms"]
print(f"Variants near BRCA1: {len(snps)}")
for v in snps[:5]:
    locs = v.get("locations", [{}])
    pos = locs[0].get("chromosomePosition", "?") if locs else "?"
    print(f"  {v['rsId']}: chr{locs[0].get('chromosomeName', '?')}:{pos}")

time.sleep(0.2)

# Get study publication details
data = gwas_get("studies/GCST000392")
pub = data.get("publicationInfo", {})
print(f"Study: {data['accessionId']}")
print(f"  Author: {pub.get('author', {}).get('fullname', 'N/A')}")
print(f"  Journal: {pub.get('publication', 'N/A')}")
print(f"  PMID: {pub.get('pubmedId', 'N/A')}")
print(f"  Date: {pub.get('publicationDate', 'N/A')}")

import requests, time

BASE = "https://www.ebi.ac.uk/gwas/rest/api"

def gwas_get(endpoint, params=None):
    url = f"{BASE}/{endpoint}"
    resp = requests.get(url, params=params or {})
    resp.raise_for_status()
    time.sleep(0.2)
    return resp.json()

# Step 1: Resolve trait → current shortForm via findByEfoTrait
# (the older `findByDescription` endpoint was removed).
traits = gwas_get("efoTraits/search/findByEfoTrait", {"trait": "schizophrenia"})
efo_id = traits["_embedded"]["efoTraits"][0]["shortForm"]
print(f"Using EFO: {efo_id}")

# Step 2: Get all associations for this trait — nested path
# `efoTraits/{shortForm}/associations` works once you have the canonical shortForm.
all_assocs = []
page = 0
while True:
    data = gwas_get(f"efoTraits/{efo_id}/associations",
                    {"size": 500, "page": page})
    assocs = data["_embedded"]["associations"]
    all_assocs.extend(assocs)
    if page >= data["page"]["totalPages"] - 1:
        break
    page += 1

# Step 3: Filter genome-wide significant
significant = [a for a in all_assocs if a.get("pvalue") and a["pvalue"] < 5e-8]
print(f"Total associations: {len(all_assocs)}, genome-wide significant: {len(significant)}")

# Step 4: Extract variant and effect details
for a in significant[:10]:
    risk_alleles = a.get("loci", [{}])[0].get("strongestRiskAlleles", [])
    snp = risk_alleles[0].get("snps", [{}])[0].get("rsId", "N/A") if risk_alleles else "N/A"
    or_val = a.get("orPerCopyNum", "N/A")
    beta = a.get("betaNum", "N/A")
    print(f"  {snp} | p={a['pvalue']:.2e} | OR={or_val} | beta={beta}")

import requests, time

BASE = "https://www.ebi.ac.uk/gwas/rest/api"

def gwas_get(endpoint, params=None):
    url = f"{BASE}/{endpoint}"
    resp = requests.get(url, params=params or {})
    resp.raise_for_status()
    time.sleep(0.2)
    return resp.json()

rs_id = "rs7903146"  # Well-known pleiotropic variant

# Get all associations for this variant
data = gwas_get(f"singleNucleotidePolymorphisms/{rs_id}/associations", {"size": 500})
assocs = data["_embedded"]["associations"]

# Association records don't embed efoTraits inline — follow the HAL
# `_links.efoTraits.href` link per association. Cache per-href to avoid
# duplicate fetches.
trait_set = {}
href_cache = {}
def fetch_traits(href):
    if href in href_cache:
        return href_cache[href]
    r = requests.get(href, timeout=15)
    href_cache[href] = (r.json().get("_embedded", {}).get("efoTraits", [])) if r.ok else []
    return href_cache[href]

for a in assocs:
    href = (a.get("_links") or {}).get("efoTraits", {}).get("href")
    if not href:
        continue
    for t in fetch_traits(href):
        tid = t.get("shortForm", "unknown")
        if tid not in trait_set:
            trait_set[tid] = {"trait": t.get("trait", "N/A"),
                              "best_pval": a.get("pvalue", 1), "count": 0}
        trait_set[tid]["count"] += 1
        if a.get("pvalue") and a["pvalue"] < trait_set[tid]["best_pval"]:
            trait_set[tid]["best_pval"] = a["pvalue"]
    time.sleep(0.1)

print(f"{rs_id} is associated with {len(trait_set)} distinct traits:")
for tid, info in sorted(trait_set.items(), key=lambda x: x[1]["best_pval"]):
    print(f"  {tid}: {info['trait']} (best p={info['best_pval']:.2e}, n={info['count']})")

import numpy as np
import matplotlib.pyplot as plt
import csv, gzip, io, requests

# Simulated summary statistics for demonstration
# In practice: download from GWAS Catalog FTP
# url = "http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCSTXXXXXX/..."
# resp = requests.get(url); data = gzip.decompress(resp.content)

np.random.seed(42)
n_snps = 5000
chroms = np.random.choice(range(1, 23), size=n_snps)
positions = np.random.randint(1, 250_000_000, size=n_snps)
pvalues = np.random.uniform(0, 1, size=n_snps)
# Add some "hits"
pvalues[:20] = 10 ** np.random.uniform(-15, -8, size=20)

# Manhattan plot
log_p = -np.log10(pvalues)
chrom_offsets = {}
running_offset = 0
for c in range(1, 23):
    chrom_offsets[c] = running_offset
    running_offset += 250_000_000

x_positions = [positions[i] + chrom_offsets[chroms[i]] for i in range(n_snps)]
colors = ["#1f77b4" if c % 2 == 0 else "#ff7f0e" for c in chroms]

fig, ax = plt.subplots(figsize=(14, 5))
ax.scatter(x_positions, log_p, c=colors, s=4, alpha=0.6)
ax.axhline(y=-np.log10(5e-8), color="red", linestyle="--", linewidth=0.8,
           label="Genome-wide significance (p=5e-8)")
ax.axhline(y=-np.log10(1e-5), color="blue", linestyle=":", linewidth=0.5,
           label="Suggestive (p=1e-5)")
ax.set_xlabel("Chromosome")
ax.set_ylabel("-log10(p-value)")
ax.set_title("GWAS Manhattan Plot")
ax.legend(fontsize=8)

# Chromosome labels
tick_positions = [chrom_offsets[c] + 125_000_000 for c in range(1, 23)]
ax.set_xticks(tick_positions)
ax.set_xticklabels([str(c) for c in range(1, 23)], fontsize=7)
plt.tight_layout()
plt.savefig("manhattan_plot.png", dpi=300, bbox_inches="tight")
print("Saved manhattan_plot.png")

import requests, time

# Step 1: Get EFO/MONDO shortForm from GWAS Catalog
# (the old `findByDescription` endpoint was removed; use `findByEfoTrait`)
BASE = "https://www.ebi.ac.uk/gwas/rest/api"
traits = requests.get(f"{BASE}/efoTraits/search/findByEfoTrait",
                      params={"trait": "coronary artery disease"}).json()
efo_id = traits["_embedded"]["efoTraits"][0]["shortForm"]
time.sleep(0.2)

# Step 2: Query PGS Catalog for scores using same EFO
pgs_url = f"https://www.pgscatalog.org/rest/score/search?trait_id={efo_id}"
pgs_data = requests.get(pgs_url).json()
print(f"PGS scores for {efo_id}: {pgs_data.get('count', 0)}")
for score in pgs_data.get("results", [])[:5]:
    print(f"  {score['id']}: {score['name']} (variants: {score.get('variants_number', '?')})")

import requests, time

BASE = "https://www.ebi.ac.uk/gwas/rest/api"

# Query region: chr9:21,900,000-22,200,000 (CDKN2A/2B locus)
data = requests.get(f"{BASE}/singleNucleotidePolymorphisms/search/findByChromBpLocationRange",
                    params={"chrom": "9", "bpStart": "21900000",
                            "bpEnd": "22200000", "size": 200}).json()
time.sleep(0.2)
snps = data["_embedded"]["singleNucleotidePolymorphisms"]
print(f"Variants in CDKN2A/2B locus: {len(snps)}")

# Get associations for each variant
region_assocs = []
for v in snps[:10]:  # Limit for demo
    rs = v["rsId"]
    try:
        a_data = requests.get(f"{BASE}/singleNucleotidePolymorphisms/{rs}/associations",
                              params={"size": 50}).json()
        time.sleep(0.2)
        for a in a_data["_embedded"]["associations"]:
            for t in a.get("efoTraits", []):
                region_assocs.append({"rsId": rs, "trait": t["trait"],
                                      "pvalue": a.get("pvalue")})
    except Exception:
        continue

print(f"Total associations in region: {len(region_assocs)}")

import matplotlib.pyplot as plt
import numpy as np
import requests, time

BASE = "https://www.ebi.ac.uk/gwas/rest/api"
data = requests.get(f"{BASE}/singleNucleotidePolymorphisms/rs1801282/associations",
                    params={"size": 100}).json()
time.sleep(0.2)
assocs = data["_embedded"]["associations"]

# Extract OR values (filter to those with OR data)
entries = []
for a in assocs:
    or_val = a.get("orPerCopyNum")
    ci_text = a.get("range", "")
    trait = a.get("efoTraits", [{}])[0].get("trait", "N/A") if a.get("efoTraits") else "N/A"
    if or_val and or_val > 0:
        entries.append({"trait": trait[:30], "or": or_val, "ci": ci_text})

if entries:
    fig, ax = plt.subplots(figsize=(8, max(3, len(entries) * 0.4)))
    y_pos = range(len(entries))
    ors = [e["or"] for e in entries]
    labels = [f"{e['trait']} (OR={e['or']:.2f})" for e in entries]
    ax.barh(y_pos, [np.log(o) for o in ors], color=["#d62728" if o > 1 else "#2ca02c" for o in ors])
    ax.set_yticks(y_pos)
    ax.set_yticklabels(labels, fontsize=8)
    ax.axvline(x=0, color="black", linewidth=0.8)
    ax.set_xlabel("log(OR)")
    ax.set_title("rs1801282 (PPARG) Effect Sizes")
    plt.tight_layout()
    plt.savefig("forest_plot.png", dpi=300, bbox_inches="tight")
    print(f"Saved forest_plot.png with {len(entries)} entries")