← All skills

gwas-catalog

data

Retrieve variant-trait associations from the NHGRI-EBI GWAS Catalog REST API — by trait, study accession, variant or gene — together with the discovery and replication cohort ancestry, sample sizes and effect estimates that decide whether an association transfers to anyone else. Handles the HAL response shape, the 20-row default page, and p-values that underflow a float to 0.0.

GWAS Catalog — variant-trait associations with study context

The NHGRI-EBI GWAS Catalog is the curated record of what published genome-wide association studies actually reported: on 2026-08-27, 229,952 studies and 1,191,948 variant-trait associations, each carrying the paper it came from, the cohort it was measured in, and the ontology term a curator assigned to the trait.

This skill is about getting associations out of the REST API and onto disk with the study record attached, because an association without its cohort is not a result. It is a number estimated in some particular set of people, at a threshold the Catalog sets lower than genome-wide significance, and both of those facts live in the study record rather than in the association.

Nothing here needs an account, a key or a click-through. The REST API and the FTP release files are open over HTTPS. EMBL-EBI's Terms of Use state that it "places no additional restrictions on the use or redistribution of the data available via its Data Resources and Tools other than those provided by the original data owners", and the Catalog's own summary-statistics page says "The majority of data are made available either through 'CC0' or 'EMBL-EBI's standard terms of use' with a small number of exceptions" — those exceptions are per-study, flagged as a Usage License on the download tables, and matter when you redistribute a study's summary statistics rather than when you query the catalogue. Cite the Catalog papers under Sources.

Four mechanical traps cause most wrong GWAS Catalog results, and all four return a well-formed HTTP 200 rather than an error. They are the next four sections; two more — what counts as one association, and what a "mapped gene" actually is — follow after them.

Trap 1 — this is HAL, and there are three different response shapes

Every response is application/hal+json. Records are never at the top level; they are nested under _embedded.<collectionName>, and the collection name is the plural resource name, not results or data. Code written against a plain JSON list finds nothing, raises nothing, and reports no hits.

The part that is easy to miss is that the same API answers in three shapes, and a single unwrap helper written against the first one breaks on the other two:

Shape Looks like Returned by
Paged collection _embedded + _links + page /studies, /associations, /efoTraits, and the searches whose URI template lists {?…,page,size,…}
Unpaged collection _embedded + _links, and no page key at all /studies/{acc}/associations, /efoTraits/{id}/studies, and the searches whose template lists only {?…,projection}
Single resource a bare object with _links, and no _embedded /studies/{acc}, /singleNucleotidePolymorphisms/{rsId}, and findByAccessionId
curl -s "https://www.ebi.ac.uk/gwas/rest/api/studies?size=1" \
| python3 -c "import json,sys
d = json.load(sys.stdin)
print('top-level keys :', sorted(d))
print('collection name:', list(d['_embedded']))
print('page           :', d['page'])
print('link rels      :', sorted(d['_links']))"
top-level keys : ['_embedded', '_links', 'page']
collection name: ['studies']
page           : {'size': 1, 'totalElements': 229952, 'totalPages': 229952, 'number': 0}
link rels      : ['first', 'last', 'next', 'profile', 'search', 'self']

profile and search are always there; first, last, next and prev appear only as the position in the result set allows, which is what makes next a usable end-of-pages signal.

Write one reader that copes with all three, and make it say which shape it got rather than guessing:

import json, urllib.error, urllib.parse, urllib.request

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


def fetch(path, **params):
    """One GET against the API. Returns the parsed HAL document."""
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    with urllib.request.urlopen(url, timeout=300) as fh:
        return json.loads(fh.read())


def shape(doc):
    """Which of the three shapes came back."""
    if "_embedded" not in doc:
        return "single"
    return "paged" if "page" in doc else "unpaged"


def records(doc, name):
    """Rows out of a collection document, or an exception naming what arrived.

    An empty result still carries `_embedded`, so a missing key means the shape
    is wrong -- usually a single resource where a collection was expected."""
    if "_embedded" not in doc:
        raise TypeError(f"single resource, not a {name} collection: {sorted(doc)}")
    if name not in doc["_embedded"]:
        raise KeyError(f"no '{name}' in _embedded: {list(doc['_embedded'])}")
    return doc["_embedded"][name]


# 1. Paged collection.
d = fetch("studies", size=1)
print(f"{'/studies':52} {shape(d):8} {len(records(d, 'studies'))} row(s), "
      f"{d['page']['totalElements']} total")

# 2. Unpaged collection -- no `page`, so `d['page']` is a KeyError, and any
#    loop keyed on totalElements never runs.
d = fetch("studies/GCST012219/associations")
print(f"{'/studies/GCST012219/associations':52} {shape(d):8} "
      f"{len(records(d, 'associations'))} row(s), page key present: {'page' in d}")

# 3. Single resource -- no `_embedded` at all. This is the one that silently
#    returns zero rows if you reach for _embedded with .get().
d = fetch("studies/search/findByAccessionId", accessionId="GCST012219")
print(f"{'findByAccessionId':52} {shape(d):8} accessionId={d['accessionId']}")
try:
    records(d, "studies")
except TypeError as e:
    print("   reaching for a collection here:", e)
/studies                                             paged    1 row(s), 229952 total
/studies/GCST012219/associations                     unpaged  26 row(s), page key present: False
findByAccessionId                                    single   accessionId=GCST012219
   reaching for a collection here: single resource, not a studies collection: ['_links', 'accessionId', 'ancestries', 'diseaseTrait', 'fullPvalueSet', 'genotypingTechnologies', 'gxe', 'gxg', 'imputed', 'initialSampleSize', 'platforms', 'pooled', 'publicationInfo', 'qualifier', 'replicationSampleSize', 'snpCount', 'studyDesignComment', 'userRequested']

Read the URI template on /{resource}/search before writing a query. The root document and each /search document advertise exactly which parameters an endpoint takes, and whether page,size is among them. That is the difference between shape one and shape two, and it is the only reliable way to know:

curl -s "https://www.ebi.ac.uk/gwas/rest/api/studies/search" \
| python3 -c "import json,sys
for rel, link in sorted(json.load(sys.stdin)['_links'].items()):
    print(f\"{rel:30} {link['href'].split('/')[-1]}\")"
findByAccessionId              findByAccessionId{?accessionId,projection}
findByDiseaseTrait             findByDiseaseTrait{?diseaseTrait,page,size,sort,projection}
findByEfoTrait                 findByEfoTrait{?efoTrait,page,size,sort,projection}
findByEfoUri                   findByEfoUri{?uri,page,size,sort,projection}
findByFullPvalueSet            findByFullPvalueSet{?fullPvalueSet,page,size,sort,projection}
findByPublicationIdPubmedId    findByPublicationIdPubmedId{?pubmedId,page,size,sort,projection}
findByUserRequested            findByUserRequested{?userRequested,page,size,sort,projection}
self                           search

Trap 2 — the default page is 20 rows, and elsewhere size does nothing at all

Two opposite failures, and which one you get depends on the endpoint you happened to pick.

On a paged endpoint the default is 20 and there is no warning. Ask for the studies of a common trait and you get twenty of them, in a complete-looking document, with the evidence of truncation only in page.totalElements and in whether _links carries a next.

On an unpaged endpoint size is accepted and ignored — the whole result set comes back regardless. Probing such an endpoint with size=1 to "look at one record" fetches everything; associations/search/findByRsId?rsId=rs7412&size=1 returned all 1,182 associations for APOE, and took minutes.

And size on a paged endpoint is capped at 1000: ?size=5000 answers HTTP 200 with page.size quietly rewritten to 1000, and page.totalPages recomputed to match. Code that asks for 5000 and then loops total // 5000 times gets 1000 rows a page and stops after a fifth of them — with no error, and with every row it did read being correct.

import json, urllib.parse, urllib.request

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


def fetch_url(url):
    with urllib.request.urlopen(url, timeout=300) as fh:
        return json.loads(fh.read())


def fetch(path, **params):
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    return fetch_url(url)


def walk(path, name, per_page=100, cap=100_000, **params):
    """Follow _links.next to exhaustion, then prove nothing was left behind.

    Paging is driven by the `next` link, not by incrementing an offset: the
    server owns the page size and rewrites it if you ask for too much. On the
    last page `next` is simply absent, which is what ends the loop."""
    doc = fetch(path, size=per_page, **params)
    if "page" not in doc:
        return doc["_embedded"][name]          # unpaged endpoint: this is all of it
    total = doc["page"]["totalElements"]
    if total > cap:
        raise RuntimeError(f"{path}: {total} rows exceeds cap {cap}")
    rows = list(doc["_embedded"][name])
    while "next" in doc["_links"]:
        doc = fetch_url(doc["_links"]["next"]["href"])
        rows += doc["_embedded"][name]
    assert len(rows) == total, f"{path}: collected {len(rows)} of {total}"
    return rows


# a. The 20-row default, on a trait with far more than 20 studies.
d = fetch("studies/search/findByEfoTrait", efoTrait="type 2 diabetes mellitus")
print(f"default page   : {len(d['_embedded']['studies'])} studies returned, "
      f"{d['page']['totalElements']} exist, next link: {'next' in d['_links']}")

# b. `size` is honoured up to 1000 and silently rewritten above it.
for asked in (100, 1000, 5000):
    p = fetch("efoTraits", size=asked)["page"]
    print(f"asked size={asked:<5} -> page.size={p['size']:<5} "
          f"totalPages={p['totalPages']}")

# c. Following `next` collects the lot. Small trait, so this is four quick pages.
rows = walk("studies/search/findByEfoTrait", "studies", per_page=5,
            efoTrait="eye color")
print(f"walked         : {len(rows)} eye-colour studies, "
      f"{len({r['accessionId'] for r in rows})} distinct accessions")

# d. An unpaged endpoint ignores `size` completely -- both calls cost the same.
one = fetch("studies/GCST012219/associations", size=1)["_embedded"]["associations"]
all_ = fetch("studies/GCST012219/associations")["_embedded"]["associations"]
print(f"size=1 ignored : {len(one)} rows with size=1, {len(all_)} rows without")
default page   : 20 studies returned, 343 exist, next link: True
asked size=100   -> page.size=100   totalPages=247
asked size=1000  -> page.size=1000  totalPages=25
asked size=5000  -> page.size=1000  totalPages=25
walked         : 16 eye-colour studies, 16 distinct accessions
size=1 ignored : 26 rows with size=1, 26 rows without

Rules that fall out of this:

  • Never read _embedded without reading page in the same breath. len(rows) == page.totalElements is the assertion; a next link is the alarm.
  • Page by following _links.next, not by incrementing your own offset. The server decides the page size, and it will rewrite yours without saying so.
  • Do not probe an unpaged endpoint with a small size. It costs the full result set. Check the URI template first, or expect the whole thing.
  • The API is not fast: the global /studies collection returns roughly 100 study records per 40 seconds, so a 343-study trait is about three minutes and the full 229,952 studies is not a REST job at all. See Get the files for the bulk route.

Trap 3 — p-values are a mantissa and an exponent, and the float underflows to zero

Every association carries three p-value fields, and they do not agree at the extremes:

Field Type At p = 2e-28539
pvalueMantissa integer, 1-9 2
pvalueExponent negative integer -28539
pvalue IEEE-754 double 0.0

(The same association: rs1129038-T in GCST012219.)

An IEEE-754 double cannot represent anything below about 5e-324, and the Catalog's strongest associations are far below that — meta-analyses of hundreds of thousands of people routinely report exponents in the hundreds or thousands. So pvalue is a convenience field that silently becomes 0.0 exactly for the associations you most wanted to find.

The damage is rarely a crash. 0.0 < 5e-8 is True, so a significance filter still passes. What breaks is everything that treats zero as a value: -log10(pvalue) raises, ranking ties every strong hit at the same place, and a "drop rows with no p-value" step keyed on falsiness throws away the top of the table.

Reconstruct from the mantissa and exponent, and keep them as the sort key.

import json, math, urllib.parse, urllib.request

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


def fetch(path, **params):
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    with urllib.request.urlopen(url, timeout=300) as fh:
        return json.loads(fh.read())


def pkey(assoc):
    """Sortable, printable p-value that survives underflow.

    Returns (exponent, mantissa) -- ascending puts the strongest first, and it
    is exact, because both fields are integers. Never round-trip through the
    `pvalue` double to compare two associations."""
    return (assoc["pvalueExponent"], assoc["pvalueMantissa"])


def neglog10p(assoc):
    """-log10(p) from the integer fields. math.log10(assoc['pvalue']) raises
    ValueError on exactly the rows this exists for."""
    return -(math.log10(assoc["pvalueMantissa"]) + assoc["pvalueExponent"])


rows = fetch("studies/GCST012219/associations")["_embedded"]["associations"]
rows.sort(key=pkey)

underflowed = [a for a in rows if a["pvalue"] == 0.0]
print(f"GCST012219 (eye colour, Simcoe 2021): {len(rows)} associations, "
      f"{len(underflowed)} with pvalue == 0.0")
print(f"{'risk allele':22} {'mantissa':>9} {'exponent':>9} "
      f"{'pvalue (double)':>16} {'-log10(p)':>10}")
for a in rows[:3] + rows[-2:]:
    allele = a["loci"][0]["strongestRiskAlleles"][0]["riskAlleleName"]
    print(f"{allele:22} {a['pvalueMantissa']:>9} {a['pvalueExponent']:>9} "
          f"{a['pvalue']:>16} {neglog10p(a):>10.1f}")

try:
    math.log10(rows[0]["pvalue"])
except ValueError as e:
    print("\nmath.log10(pvalue) on the strongest association:", e)
GCST012219 (eye colour, Simcoe 2021): 26 associations, 2 with pvalue == 0.0
risk allele             mantissa  exponent  pvalue (double)  -log10(p)
rs1129038-T                    2    -28539              0.0    28538.7
rs4778218-G                    1      -582              0.0      582.0
rs12203592-T                   2      -321       1.996e-321      320.7
rs2351061-A                    4        -8            4e-08        7.4
rs80308281-T                   5        -8            5e-08        7.3

math.log10(pvalue) on the strongest association: math domain error

Note the third row. 2e-321 is representable, but only as a subnormal double, so pvalue comes back as 1.996e-321 — the nearest double, three significant figures short. Between roughly 1e-308 and 5e-324 the float is present and lossy rather than absent, which is worse than zero because nothing about it looks wrong.

Sorting /associations server-side is fragile — sort locally instead. Two separate failures, both reproduced on 2026-08-27:

  • sort=…,desc on a p-value field returns HTTP 500. /associations?sort=pvalue,desc, sort=pvalueExponent,desc and sort=pvalueMantissa,desc all fail with a server error. The same fields ascending are fine, as is desc on other fields (snpType, standardError) and on /studies and /efoTraits.
  • A sorted page can arrive truncated, at HTTP 200. The last page of /associations under sort=pvalueExponent,asc came back as a 200 with the body cut off mid-record and unparseable — three times out of three at size=100&page=11919, and twice out of twice at size=1000&page=1191. The same page without sort parsed cleanly, 48 rows. So the status code is not the check: json.loads raising is the only signal, and a client that swallows it loses the tail of the collection silently.

Sort the rows you have in hand instead — pkey above is exact and costs nothing.

Effect sizes come in two mutually exclusive forms, and both can be null at once. orPerCopyNum is an odds ratio, for case-control traits; betaNum with betaUnit and betaDirection is a linear effect, for quantitative traits, and the direction is a word ("increase" / "decrease") rather than a sign, so a beta is unsigned until you read two fields. range is the reported confidence interval as a string, standardError the SE. Never coerce a beta into an OR — and do not assume every row has either: GCST012219 carries a beta on 24 of its 26 associations and neither an OR nor a beta on the other two, because the paper did not report one.

Trap 4 — traits are ontology terms, matched exactly, and they are not all EFO

The trait endpoints are named efoTraits, and the name is now the smaller half of the story. Counting the distinct terms in gwas-efo-trait-mappings.tsv from the 2026-08-24 release — 20,047 of them, over 110,338 curated trait strings — gives:

Prefix Distinct terms Ontology
OBA 10,027 Ontology of Biological Attributes — measurements
EFO 9,032 Experimental Factor Ontology
MONDO 491 Mondo Disease Ontology
HP 397 Human Phenotype Ontology
GO, Orphanet, PATO, NCIT, UBERON, GSSO, MP, HANCESTRO 100 between them

Most terms are no longer EFO, and the uri switches namespace with the prefix: http://www.ebi.ac.uk/efo/EFO_0003949, http://purl.obolibrary.org/obo/MONDO_0005148, http://www.orpha.net/ORDO/Orphanet_654. Three namespaces, twelve prefixes. Any code that recovers an id by stripping http://www.ebi.ac.uk/efo/ hands back a full URI for roughly a third of the mappings and every one of the 10,027 OBA terms.

The label search is an exact match on the whole string — case is ignored, and nothing else is. Not a prefix, not a substring, no synonyms, no spelling variants, and a trailing space is fatal. And the two "not found" answers differ: findByEfoTrait returns HTTP 200 with an empty list, findByShortForm returns HTTP 404.

import json, urllib.error, urllib.parse, urllib.request

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


def fetch(path, **params):
    """Returns (status, document). 404 is a real answer here, not a failure."""
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    try:
        with urllib.request.urlopen(url, timeout=300) as fh:
            return fh.status, json.loads(fh.read())
    except urllib.error.HTTPError as e:
        return e.code, None


def resolve_trait(label):
    """Exact trait label -> (shortForm, uri). Raises rather than guessing."""
    status, doc = fetch("efoTraits/search/findByEfoTrait", trait=label)
    hits = doc["_embedded"]["efoTraits"] if doc else []
    if not hits:
        raise LookupError(f"{label!r}: no trait with that exact label "
                          f"(HTTP {status} with an empty list, not an error)")
    if len(hits) > 1:
        raise LookupError(f"{label!r}: {[h['shortForm'] for h in hits]}")
    return hits[0]["shortForm"], hits[0]["uri"]


for label in ("type 2 diabetes",              # a substring is not a match
              "type II diabetes mellitus",    # the wording of the retired EFO term
              "eye colour",                   # British spelling
              "type 2 diabetes mellitus ",    # one trailing space
              "TYPE 2 DIABETES MELLITUS",     # case, and only case, is ignored
              "type 2 diabetes mellitus",
              "body mass index",
              "eye color"):
    try:
        sf, uri = resolve_trait(label)
        print(f"{label!r:29} -> {sf:16} {uri}")
    except LookupError:
        print(f"{label!r:29} -> no match (HTTP 200, empty list)")

# The retired EFO id for type 2 diabetes 404s; the live term is MONDO.
for sf in ("EFO_0001360", "MONDO_0005148", "EFO_0003949"):
    status, doc = fetch("efoTraits/search/findByShortForm", shortForm=sf)
    label = doc["trait"] if doc else "-"
    print(f"shortForm {sf:16} HTTP {status}  {label}")
'type 2 diabetes'             -> no match (HTTP 200, empty list)
'type II diabetes mellitus'   -> no match (HTTP 200, empty list)
'eye colour'                  -> no match (HTTP 200, empty list)
'type 2 diabetes mellitus '   -> no match (HTTP 200, empty list)
'TYPE 2 DIABETES MELLITUS'    -> MONDO_0005148    http://purl.obolibrary.org/obo/MONDO_0005148
'type 2 diabetes mellitus'    -> MONDO_0005148    http://purl.obolibrary.org/obo/MONDO_0005148
'body mass index'             -> EFO_0004340      http://www.ebi.ac.uk/efo/EFO_0004340
'eye color'                   -> EFO_0003949      http://www.ebi.ac.uk/efo/EFO_0003949
shortForm EFO_0001360      HTTP 404  -
shortForm MONDO_0005148    HTTP 200  type 2 diabetes mellitus
shortForm EFO_0003949      HTTP 200  eye color

EFO_0001360 — the id every pre-2023 pipeline used for type 2 diabetes — is gone. So:

  • Resolve the label to a short form once, then key on the short form, and record which one you resolved alongside the results. A hardcoded ontology id has a shelf life.
  • There is no fuzzy trait search in this API. When you do not already know the exact label, get the vocabulary rather than guessing at it — page /efoTraits (24,618 terms, ~1.3 s per 100-row page) and match locally, or read gwas-efo-trait-mappings.tsv from the release, below.
  • Do not construct short forms by string-surgery on uri. Read the shortForm field, which the API supplies for both namespaces.
  • The Catalog also keeps a curator-written free-text trait on the study (diseaseTrait.trait, e.g. "Eye color", capitalised) distinct from the ontology term on the association ("eye color"). They are different fields with different vocabularies; the file gwas-efo-trait-mappings.tsv in each release maps one to the other offline, for all 110,338 curated trait strings.
  • findByDiseaseTrait searches that free-text study trait, findByEfoTrait the ontology label. Picking the wrong one is a silent empty result.

An association is not a result

Everything above is plumbing. This is the part that decides whether the number you return means anything, and it is why the study record has to travel with the association rather than being fetched only when someone asks.

Three facts change the interpretation, and all three are in the study record:

  1. The Catalog's curation threshold is 1×10⁻⁵, not genome-wide significance. Associations weaker than 5×10⁻⁸ are in the table by design — the weakest association in the whole catalogue is exactly 1×10⁻⁵, and nothing beyond it is curated. Filtering at 5×10⁻⁸ is the reader's job, and the Catalog does not do it.
  2. Ancestry of the discovery cohort. ancestries is a list with type of "initial" (discovery) or "replication", each with numberOfIndividuals, ancestralGroups and countryOfRecruitment. An effect estimated in Europeans has no guaranteed transfer — allele frequency, LD structure and effect size all differ — and the great majority of the Catalog is European.
  3. Whether anything replicated. A study with no "replication" ancestry entry and replicationSampleSize of NA reported a discovery result and stopped there.

projection=associationByStudy and projection=associationByEfoTrait inline the whole study record — ancestries included — into each association, which turns N+1 requests into one:

import json, urllib.parse, urllib.request
from collections import Counter

BASE = "https://www.ebi.ac.uk/gwas/rest/api"
GWAS_SIG = 5e-8


def fetch(path, **params):
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    with urllib.request.urlopen(url, timeout=300) as fh:
        return json.loads(fh.read())


def p_of(a):
    """Exact p-value as a float where representable, else the true magnitude
    as a (mantissa, exponent) pair. Only ever used for a threshold test."""
    return a["pvalueMantissa"] * 10.0 ** a["pvalueExponent"]


def cohort(study, stage):
    """(individuals, ancestral groups) for one stage of a study."""
    entries = [a for a in study["ancestries"] if a["type"] == stage]
    n = sum(a["numberOfIndividuals"] or 0 for a in entries)
    groups = {g["ancestralGroup"] for a in entries for g in a["ancestralGroups"]}
    return n, groups


# One request: every association for the trait, each with its study inlined.
rows = fetch("efoTraits/EFO_0003949/associations",
             projection="associationByEfoTrait")["_embedded"]["associations"]
studies = {a["study"]["accessionId"]: a["study"] for a in rows}
print(f"eye color (EFO_0003949): {len(rows)} associations across "
      f"{len(studies)} studies, in ONE request")

weak = [a for a in rows if p_of(a) > GWAS_SIG]
print(f"weaker than 5e-8       : {len(weak)} of {len(rows)} "
      f"({100 * len(weak) / len(rows):.0f}%) -- curated at p < 1e-5, not 5e-8")

disc = Counter()
for s in studies.values():
    n, groups = cohort(s, "initial")
    for g in groups:
        disc[g] += n
total = sum(disc.values())
print("\ndiscovery individuals by ancestral group")
for g, n in disc.most_common():
    print(f"  {g:24} {n:>9,}  {100 * n / total:5.1f}%")

no_rep = [acc for acc, s in sorted(studies.items()) if not cohort(s, "replication")[0]]
print(f"\nstudies with no replication stage: {len(no_rep)} of {len(studies)}")
print("  " + ", ".join(no_rep))

# The sub-phenotype hides in a free-text field on the ASSOCIATION, not the trait.
descs = Counter(repr(a["pvalueDescription"]) for a in rows)
print("\npvalueDescription, top 6 of "
      f"{len(descs)} distinct values on one ontology trait")
for d, n in descs.most_common(6):
    print(f"  {d:36} {n:>4}")
eye color (EFO_0003949): 201 associations across 16 studies, in ONE request
weaker than 5e-8       : 41 of 201 (20%) -- curated at p < 1e-5, not 5e-8

discovery individuals by ancestral group
  European                   199,091   98.7%
  East Asian                   2,616    1.3%

studies with no replication stage: 8 of 16
  GCST000710, GCST005093, GCST005094, GCST005095, GCST005096, GCST007489, GCST90026489, GCST90255684

pvalueDescription, top 6 of 13 distinct values on one ontology trait
  None                                  138
  '(EA)'                                 22
  '(Eye color)'                           8
  ' '                                     6
  '(eye color)'                           5
  '(Eye color phototype score)'           5

What that output is actually telling you, and what to carry into a write-up:

  • A fifth of the rows are below genome-wide significance. Report the threshold you applied, every time. p_of() above is the only safe way to test it, because a["pvalue"] is 0.0 for the strongest rows.
  • 98.7% of the discovery sample is European, and the entire non-European remainder is one East Asian cohort. Any statement of the form "variant X raises risk of Y" inherits that, and a polygenic score built on it is known to lose accuracy out of ancestry. Where an ancestralGroup is "NR" it means not reported — treat it as missing, not as a category, and do not let it dilute a percentage silently.
  • Eight of sixteen studies have no replication stage. That is not a defect; it is what a discovery paper looks like. It is a defect to report it as though it had one.
  • One ontology trait is several phenotypes. pvalueDescription is free text on the association, and thirteen distinct values sit under the single term "eye color" — '(EA)', '(Eye color phototype score)', and elsewhere in the same 201 rows '(Blue eye color)', '(Brown eye color)' and '(Intermediate/Green eye color)', which take opposite effect directions at the same locus. It is inconsistently cased ('(Eye color)' and '(eye color)' are both present), None on 138 rows, and on six rows it is a single space rather than empty. Strip it, read it, and never pool effect sizes across distinct values.

Two further study-record fields worth carrying: fullPvalueSet (whether full summary statistics were deposited, which decides if fine-mapping or Mendelian randomisation is even possible) and snpCount with imputed (how many variants were tested, which sets the real multiple-testing burden).

Variants, mapped genes and genomic context

/singleNucleotidePolymorphisms/{rsId} gives the variant's mapped position, cytogenetic region, functional class and the genes around it. Two things about genomicContexts are easy to get wrong:

  • It is mapped twice, once from Ensembl and once from NCBI, so every gene appears more than once and the two sources disagree on distances. Counting rows counts sources.
  • isClosestGene is not unique, and it is not the gene the variant is in. It is flagged per source and per direction — nearest upstream and nearest downstream both get it — so a single variant can have four rows marked closest. The gene the variant actually sits inside is the one with isIntergenic: false, and it is flagged isClosestGene: false.
import json, urllib.parse, urllib.request

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


def fetch(path, **params):
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    with urllib.request.urlopen(url, timeout=300) as fh:
        return json.loads(fh.read())


snp = fetch("singleNucleotidePolymorphisms/rs1129038")
loc = snp["locations"][0]
print(f"{snp['rsId']}  chr{loc['chromosomeName']}:{loc['chromosomePosition']}  "
      f"{loc['region']['name']}  {snp['functionalClass']}  merged={snp['merged']}")

ctx = snp["genomicContexts"]
print(f"\n{len(ctx)} genomicContexts rows, "
      f"{len({c['gene']['geneName'] for c in ctx})} distinct genes, "
      f"sources {sorted({c['source'] for c in ctx})}")
print(f"{'gene':10} {'source':9} {'distance':>9} {'intergenic':>11} "
      f"{'upstream':>9} {'closest':>8}")
for c in sorted(ctx, key=lambda c: (c["gene"]["geneName"], c["source"])):
    print(f"{c['gene']['geneName']:10} {c['source']:9} {c['distance']:>9} "
          f"{str(c['isIntergenic']):>11} {str(c['isUpstream']):>9} "
          f"{str(c['isClosestGene']):>8}")

# Genes the variant is actually inside, deduplicated across the two sources.
inside = sorted({c["gene"]["geneName"] for c in ctx if not c["isIntergenic"]})
closest = sorted({c["gene"]["geneName"] for c in ctx if c["isClosestGene"]})
print(f"\nvariant lies within : {inside}")
print(f"flagged closest     : {closest}")
rs1129038  chr15:28111713  15q13.1  3_prime_UTR_variant  merged=0

6 genomicContexts rows, 3 distinct genes, sources ['Ensembl', 'NCBI']
gene       source     distance  intergenic  upstream  closest
HERC2      Ensembl           0       False     False    False
HERC2      NCBI              0       False     False    False
OCA2       Ensembl       12049        True      True     True
OCA2       NCBI          12398        True      True     True
RPL41P2    Ensembl      203938        True     False     True
RPL41P2    NCBI         203940        True     False     True
variant lies within : ['HERC2']
flagged closest     : ['OCA2', 'RPL41P2']

That is the textbook eye-colour locus, and every part of the trap is visible in six rows. The variant is in the 3′ UTR of HERC2; the gene it regulates is OCA2, about 12 kb upstream; and RPL41P2 — a ribosomal protein pseudogene 204 kb downstream — is flagged isClosestGene just as OCA2 is, because each direction gets its own nearest gene. Ensembl and NCBI put OCA2 at 12,049 and 12,398 bp respectively, so even the distance depends on which row you read. Taking genomicContexts[0], or filtering on isClosestGene alone, returns a pseudogene as the answer to "which gene is this?"

Lookups by variant, gene and region all hang off the same resource, and they split across both collection shapes — which is why the template, not the name, tells you what you are getting. On /singleNucleotidePolymorphisms/search, findByRsId, findByBpLocation, findByEfoTrait and findByDiseaseTrait are unpaged and return the whole result set; findByGene, findByChromBpLocationRange, findIdsByLocations… and findByPubmedId are paged and default to 20.

And the two "not found" answers differ again: an rsID the Catalog does not hold gives HTTP 200 with an empty list from associations/search/findByRsId?rsId=…, but 404 from /singleNucleotidePolymorphisms/{rsId}. Same fact, two signals, and only one of them raises.

Not every association is one variant

loci[0]["strongestRiskAlleles"][0] is the obvious way to read a risk allele, and it is right for almost every row and silently wrong for two whole classes of them. They are rare — of 1,000 consecutive /associations rows sampled on 2026-08-27, 999 were single variants and one was a 3-SNP haplotype — which is exactly why nobody writes the code that handles them, and why the error survives review. The Catalog records three forms, distinguished by two booleans on the association:

Form multiSnpHaplotype snpInteraction Shape
Single variant false false one locus, one risk allele
Haplotype true false one locus, N risk alleles, haplotypeSnpCount = N
SNP × SNP interaction false true N loci, one risk allele each

So the two failure modes point in opposite directions: the haplotype hides its extra variants inside loci[0], and the interaction hides them in the loci you never looked at. Both come back as a perfectly ordinary association row.

import json, textwrap, urllib.parse, urllib.request
from collections import Counter

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


def fetch(path, **params):
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    with urllib.request.urlopen(url, timeout=300) as fh:
        return json.loads(fh.read())


def risk_alleles(assoc):
    """Every risk allele in an association, flattened across loci, in order."""
    return [a["riskAlleleName"] for locus in assoc["loci"]
            for a in locus["strongestRiskAlleles"]]


def form(assoc):
    if assoc["snpInteraction"]:
        return "interaction"
    return "haplotype" if assoc["multiSnpHaplotype"] else "single"


for acc in ("GCST012219", "GCST000071", "GCST006494"):
    rows = fetch(f"studies/{acc}/associations")["_embedded"]["associations"]
    counts = Counter(form(r) for r in rows)
    # Ties broken on the allele list itself, so this picks the same row every
    # run -- the API does not promise an order and max() keeps the first max.
    widest = max(rows, key=lambda r: (len(risk_alleles(r)), risk_alleles(r)))
    alleles = risk_alleles(widest)
    locus = widest["loci"][0]
    print(f"{acc}: {len(rows)} associations {dict(counts)}")
    print(f"  widest row      : {len(widest['loci'])} loci, {len(alleles)} risk "
          f"alleles, description {locus['description']!r}, "
          f"haplotypeSnpCount {locus['haplotypeSnpCount']}")
    print(f"  the naive read  : {alleles[0]}")
    print("  what is there   :")
    print(textwrap.fill(", ".join(alleles), width=70,
                        initial_indent="    ", subsequent_indent="    "))
GCST012219: 26 associations {'single': 26}
  widest row      : 1 loci, 1 risk alleles, description 'Single variant', haplotypeSnpCount None
  the naive read  : rs9971729-C
  what is there   :
    rs9971729-C
GCST000071: 3 associations {'haplotype': 2, 'single': 1}
  widest row      : 1 loci, 16 risk alleles, description '17 marker haplotype-1', haplotypeSnpCount 16
  the naive read  : rs11209003-?
  what is there   :
    rs11209003-?, rs11209002-?, rs2064689-?, rs1004819-?, rs2902440-?,
    rs11465802-?, rs2201841-?, rs11465804-?, rs11209026-?,
    rs1343151-?, rs10889676-?, rs10889677-?, rs9988642-?,
    rs12567232-?, rs6669582-?, rs10789230-?
GCST006494: 7 associations {'interaction': 7}
  widest row      : 2 loci, 2 risk alleles, description 'SNP x SNP interaction', haplotypeSnpCount None
  the naive read  : rs10830963-G
  what is there   :
    rs10830963-G, rs73659517-G

Reporting rs11209003-? as the variant for that row turns a 16-marker haplotype into a single-SNP claim, and reporting rs10830963-G alone turns an interaction into a main effect. Three further details visible in that output:

  • -? means the risk allele was not reported. The suffix on riskAlleleName is the allele, and ? is a real value in the data — do not parse it as a base.
  • description and haplotypeSnpCount can disagree. That row is labelled '17 marker haplotype-1' and carries 16 alleles with haplotypeSnpCount: 16. Trust the array you can count, not the curator's free-text label.
  • authorReportedGenes lives on the locus, not the association, so an interaction has one gene list per locus. Flatten across loci, as risk_alleles does above.

Get the files

Everything above prints. This writes, because a target dossier, a figure or a meta-analysis needs files — and because the catalogue grows every fortnight, so a directory of TSVs with no release stamp cannot be compared against a later pull.

Route A — the REST API, for one trait. One request for the associations with the study record inlined, then flatten to two TSVs plus a manifest recording what was resolved and when.

import csv, datetime, json, os, urllib.parse, urllib.request

BASE = "https://www.ebi.ac.uk/gwas/rest/api"
TRAIT_LABEL = "eye color"
OUT = "Data/gwas-catalog"
GWAS_SIG = 5e-8


def fetch(path, **params):
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    with urllib.request.urlopen(url, timeout=300) as fh:
        return url, json.loads(fh.read())


def write_tsv(path, rows, cols):
    with open(path, "w", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=cols, delimiter="\t",
                           extrasaction="ignore")
        w.writeheader()
        w.writerows(rows)
    return {"path": path, "rows": len(rows), "bytes": os.path.getsize(path)}


os.makedirs(OUT, exist_ok=True)
sources, files = {}, []

# 1. Exact label -> ontology short form. Never hardcode the id.
url, doc = fetch("efoTraits/search/findByEfoTrait", trait=TRAIT_LABEL)
hits = doc["_embedded"]["efoTraits"]
if len(hits) != 1:
    raise LookupError(f"{TRAIT_LABEL!r}: {[h['shortForm'] for h in hits]}")
trait = hits[0]
sources["trait"] = url
print("resolved:", trait["shortForm"], "|", trait["uri"])

# 2. One request for every association, study record inlined.
url, doc = fetch(f"efoTraits/{trait['shortForm']}/associations",
                 projection="associationByEfoTrait")
assocs = doc["_embedded"]["associations"]
sources["associations"] = url

# 3. Associations, flat. Mantissa and exponent are the p-value of record; the
#    double is carried too, so a reader can see where it underflowed.
arows = []
for a in assocs:
    # Flatten across every locus: a haplotype is one locus with N alleles, an
    # interaction is N loci with one each, and taking loci[0][0] loses both.
    alleles = [al["riskAlleleName"] for L in a["loci"]
               for al in L["strongestRiskAlleles"]]
    genes = sorted({g["geneName"] for L in a["loci"]
                    for g in L["authorReportedGenes"]})
    p = a["pvalueMantissa"] * 10.0 ** a["pvalueExponent"]
    arows.append({
        "studyAccession": a["study"]["accessionId"],
        "associationForm": ("interaction" if a["snpInteraction"] else
                            "haplotype" if a["multiSnpHaplotype"] else "single"),
        "lociCount": len(a["loci"]),
        "riskAlleleCount": len(alleles),
        "riskAlleles": ";".join(alleles),
        "riskFrequency": a["riskFrequency"],
        "authorReportedGenes": ",".join(genes),
        "pvalueMantissa": a["pvalueMantissa"],
        "pvalueExponent": a["pvalueExponent"],
        "pvalueDouble": a["pvalue"],
        "genomeWideSignificant": p <= GWAS_SIG,
        "pvalueDescription": (a["pvalueDescription"] or "").strip(),
        "orPerCopyNum": a["orPerCopyNum"],
        "betaNum": a["betaNum"],
        "betaUnit": a["betaUnit"],
        "betaDirection": a["betaDirection"],
        "standardError": a["standardError"],
        "range": a["range"],
        "snpType": a["snpType"],
    })
arows.sort(key=lambda r: (r["pvalueExponent"], r["pvalueMantissa"]))
files.append(write_tsv(f"{OUT}/{trait['shortForm']}_associations.tsv", arows,
                       list(arows[0])))

# 4. Studies, one row each, with the cohort that decides interpretation.
srows = []
for acc, s in sorted({a["study"]["accessionId"]: a["study"]
                      for a in assocs}.items()):
    stage = {}
    for st in ("initial", "replication"):
        e = [x for x in s["ancestries"] if x["type"] == st]
        stage[st] = (sum(x["numberOfIndividuals"] or 0 for x in e),
                     sorted({g["ancestralGroup"] for x in e
                             for g in x["ancestralGroups"]}))
    srows.append({
        "studyAccession": acc,
        "pubmedId": s["publicationInfo"]["pubmedId"],
        "firstAuthor": s["publicationInfo"]["author"]["fullname"],
        "publicationDate": s["publicationInfo"]["publicationDate"],
        "journal": s["publicationInfo"]["publication"],
        "diseaseTrait": s["diseaseTrait"]["trait"],
        "discoveryN": stage["initial"][0],
        "discoveryAncestry": "|".join(stage["initial"][1]),
        "replicationN": stage["replication"][0],
        "replicationAncestry": "|".join(stage["replication"][1]),
        "initialSampleSize": s["initialSampleSize"],
        "replicationSampleSize": s["replicationSampleSize"],
        "snpCount": s["snpCount"],
        "imputed": s["imputed"],
        "fullPvalueSet": s["fullPvalueSet"],
        "genotypingTechnology": "|".join(
            t["genotypingTechnology"] for t in s["genotypingTechnologies"]),
    })
files.append(write_tsv(f"{OUT}/{trait['shortForm']}_studies.tsv", srows,
                       list(srows[0])))

# 5. The manifest. Without it these numbers cannot be compared to a later pull.
_, root = fetch("studies", size=1)
manifest = {
    "retrieved": datetime.date.today().isoformat(),
    "traitLabel": TRAIT_LABEL,
    "traitShortForm": trait["shortForm"],
    "traitUri": trait["uri"],
    "catalogueStudiesTotal": root["page"]["totalElements"],
    "significanceThreshold": GWAS_SIG,
    "associationsRetrieved": len(arows),
    "associationsAtThreshold": sum(r["genomeWideSignificant"] for r in arows),
    "studiesRetrieved": len(srows),
    "sources": sources,
    "files": files,
}
with open(f"{OUT}/manifest.json", "w") as fh:
    json.dump(manifest, fh, indent=2)

for f in files:
    print(f"{f['path']:52} {f['rows']:>4} rows  {f['bytes']:>7,} bytes")
print(f"{OUT + '/manifest.json':52} "
      f"{manifest['associationsAtThreshold']} of "
      f"{manifest['associationsRetrieved']} at p <= 5e-8")
resolved: EFO_0003949 | http://www.ebi.ac.uk/efo/EFO_0003949
Data/gwas-catalog/EFO_0003949_associations.tsv        201 rows   24,094 bytes
Data/gwas-catalog/EFO_0003949_studies.tsv              16 rows    3,492 bytes
Data/gwas-catalog/manifest.json                      160 of 201 at p <= 5e-8

Route B — the release files, for anything catalogue-wide. At roughly half a second per study record the REST API cannot answer questions about all 229,952 studies. Each fortnightly release publishes the whole catalogue as flat files over anonymous HTTPS at https://ftp.ebi.ac.uk/pub/databases/gwas/releases/latest/:

File Size What it carries
gwas-catalog-download-studies-v1.0.3.1.txt 123 MB every study, with STUDY ACCESSION, MAPPED_TRAIT_URI, COHORT, FULL SUMMARY STATISTICS, SUMMARY STATS LOCATION
gwas-catalog-ancestry.tsv 49 MB one row per study per ancestry stage — the cohort table
gwas-catalog-associations_ontology-annotated-full.zip 70 MB every association, ontology-annotated
gwas-efo-trait-mappings.tsv 19 MB free-text study trait → ontology term and parent

Stream it rather than loading it, and join on the accessions Route A already collected. This block downloads about 49 MB and takes a few minutes on a slow link:

import csv, io, json, os, urllib.request
from collections import Counter

FTP = "https://ftp.ebi.ac.uk/pub/databases/gwas/releases/latest"
OUT = "Data/gwas-catalog"

wanted = {r["studyAccession"] for r in
          csv.DictReader(open(f"{OUT}/EFO_0003949_studies.tsv"), delimiter="\t")}
print(f"{len(wanted)} accessions from Route A")

# The header is misspelled upstream -- NUMBER OF INDIVDUALS, ADDITONAL ANCESTRY
# DESCRIPTION. Read the header row; never hand-type these names.
url = f"{FTP}/gwas-catalog-ancestry.tsv"
rows, seen, stages = [], 0, Counter()
with urllib.request.urlopen(url, timeout=600) as fh:
    reader = csv.DictReader(io.TextIOWrapper(fh, encoding="utf-8"), delimiter="\t")
    header = reader.fieldnames
    for rec in reader:
        seen += 1
        stages[rec["STAGE"]] += 1
        if rec["STUDY ACCESSION"] in wanted:
            rows.append(rec)

print(f"streamed {seen:,} ancestry rows for the whole catalogue")
print(f"stages: {dict(stages.most_common())}")
print(f"kept {len(rows)} rows for {len({r['STUDY ACCESSION'] for r in rows})} studies")

with open(f"{OUT}/EFO_0003949_ancestry.tsv", "w", newline="") as fh:
    w = csv.DictWriter(fh, fieldnames=header, delimiter="\t")
    w.writeheader()
    w.writerows(sorted(rows, key=lambda r: (r["STUDY ACCESSION"], r["STAGE"])))

with open(f"{OUT}/manifest.json") as fh:
    manifest = json.load(fh)
manifest["bulkRelease"] = {"url": url, "ancestryRowsInRelease": seen,
                           "ancestryRowsKept": len(rows), "columns": header}
with open(f"{OUT}/manifest.json", "w") as fh:
    json.dump(manifest, fh, indent=2)

path = f"{OUT}/EFO_0003949_ancestry.tsv"
print(f"{path:52} {len(rows):>4} rows  {os.path.getsize(path):>7,} bytes")
16 accessions from Route A
streamed 298,269 ancestry rows for the whole catalogue
stages: {'initial': 260878, 'replication': 37391}
kept 26 rows for 16 studies
Data/gwas-catalog/EFO_0003949_ancestry.tsv             26 rows    4,452 bytes

That output settles something the per-trait view only hints at. Across all 229,952 studies the release carries 260,878 initial rows and 37,391 replication rows — and since a study with a replication stage contributes at least one row, at most one study in six replicated anything at all. The eight-of-sixteen seen above for eye colour is the ordinary case, not an unlucky draw, and "was it replicated?" is a question worth asking of every association you plan to report.

Full summary statistics, when a study deposited them. Studies with fullPvalueSet: true have a directory under https://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/, bucketed a thousand accessions at a time — GCST012001-GCST013000/GCST012219/. Do not build the filename from the accession: submissions since the GWAS-SSF standard are {ACCESSION}.tsv.gz with a -meta.yaml beside them, but legacy studies are not (GCST000028/Saxena-17463246.txt, and GCST000510 holds fifteen separate .assoc.gz files, one per blood trait). List the directory, read the -meta.yaml for the genome assembly and imputation panel, and check md5sum.txt. A study with fullPvalueSet: false has no directory at all, and the URL 404s.

Try it

A self-contained check that this skill still holds. Public data, no account, no key, Python standard library only — no pip install step, and nothing written to disk.

Data — the GWAS Catalog REST API, reached through nine endpoints:

https://www.ebi.ac.uk/gwas/rest/api/studies?size=1
https://www.ebi.ac.uk/gwas/rest/api/studies/search/findByAccessionId?accessionId=GCST012219
https://www.ebi.ac.uk/gwas/rest/api/studies/GCST012219/associations
https://www.ebi.ac.uk/gwas/rest/api/studies/GCST000071/associations
https://www.ebi.ac.uk/gwas/rest/api/studies/GCST006494/associations
https://www.ebi.ac.uk/gwas/rest/api/studies/search/findByEfoTrait?efoTrait=eye%20color&size=5
https://www.ebi.ac.uk/gwas/rest/api/efoTraits/search/findByEfoTrait?trait=type%202%20diabetes%20mellitus
https://www.ebi.ac.uk/gwas/rest/api/efoTraits/search/findByShortForm?shortForm=EFO_0001360
https://www.ebi.ac.uk/gwas/rest/api/efoTraits/EFO_0003949/associations?projection=associationByEfoTrait

The Catalog is open EMBL-EBI infrastructure; no account or licence acceptance is needed to read it. Each probe is chosen for a trap rather than for being typical. GCST012219 (Simcoe et al. 2021, PMID 33692100, eye colour in 192,986 people) exercises three at once: it is reached through an unpaged sub-resource, its strongest association is 2×10⁻²⁸⁵³⁹ and therefore underflows the pvalue double to 0.0, and it has no replication stage. EFO_0003949 and MONDO_0005148 probe the trait vocabulary, EFO_0001360 is the retired id that now 404s, and GCST000071 (a 16-marker haplotype) and GCST006494 (SNP × SNP interactions) are the two association shapes that are not one variant. Last confirmed reachable 2026-08-27.

Budget about two minutes. Most calls answer in a second or two; the last one — every association for a trait, with the study record inlined — took roughly 70 seconds, which is the price of the one request that replaces seventeen.

Run — this block asserts the shape of the API, not only its values, so a renamed field or restructured pagination fails it loudly:

import json, math, urllib.error, urllib.parse, urllib.request

BASE = "https://www.ebi.ac.uk/gwas/rest/api"
GWAS_SIG = 5e-8
DOUBLE_FLOOR = -324           # below this an IEEE-754 double is exactly 0.0


def fetch(url):
    """Returns (status, document). 404 is an answer here, not a failure."""
    try:
        with urllib.request.urlopen(url, timeout=300) as fh:
            return fh.status, json.loads(fh.read())
    except urllib.error.HTTPError as e:
        return e.code, None


def get(path, **params):
    url = f"{BASE}/{path}"
    if params:
        url += "?" + urllib.parse.urlencode(params, doseq=True)
    return fetch(url)


# 1. SHAPE: a paged collection is _embedded + _links + page, and nothing else.
status, doc = get("studies", size=1)
assert status == 200, status
assert sorted(doc) == ["_embedded", "_links", "page"], sorted(doc)
assert list(doc["_embedded"]) == ["studies"], list(doc["_embedded"])
assert sorted(doc["page"]) == ["number", "size", "totalElements", "totalPages"], \
    sorted(doc["page"])
assert doc["page"]["size"] == 1 and doc["page"]["number"] == 0, doc["page"]
assert "next" in doc["_links"], sorted(doc["_links"])
n_studies = doc["page"]["totalElements"]
print(f"1 paged shape      : _embedded.studies + page, {n_studies:,} studies")

# 2. SHAPE: findByAccessionId is a SINGLE resource -- no _embedded at all.
status, doc = get("studies/search/findByAccessionId", accessionId="GCST012219")
assert status == 200 and "_embedded" not in doc, sorted(doc)
assert doc["accessionId"] == "GCST012219", doc["accessionId"]
assert doc["fullPvalueSet"] is False and doc["imputed"] is True
print(f"2 single resource  : GCST012219, {doc['snpCount']:,} variants tested, "
      f"fullPvalueSet={doc['fullPvalueSet']}")

# 3. SHAPE: a sub-resource is UNPAGED -- no `page` key, and `size` is ignored.
status, one = get("studies/GCST012219/associations", size=1)
status2, all_ = get("studies/GCST012219/associations")
assert status == status2 == 200
assert "page" not in one and "page" not in all_, "sub-resource grew a page key"
rows = all_["_embedded"]["associations"]
assert len(one["_embedded"]["associations"]) == len(rows), "size is now honoured"
print(f"3 unpaged shape    : size=1 returned all {len(rows)} associations")

# 4. INVARIANT: mantissa x 10^exponent is the p-value of record, and the
#    `pvalue` double is exactly 0.0 for everything below the double floor.
for a in rows:
    assert 1 <= a["pvalueMantissa"] <= 9, a["pvalueMantissa"]
    assert a["pvalueExponent"] < 0, a["pvalueExponent"]
    if a["pvalueExponent"] <= DOUBLE_FLOOR:
        assert a["pvalue"] == 0.0, (a["pvalueMantissa"], a["pvalueExponent"])
    else:
        assert math.isclose(a["pvalue"],
                            a["pvalueMantissa"] * 10.0 ** a["pvalueExponent"],
                            rel_tol=1e-9), a["pvalue"]
best = min(rows, key=lambda a: (a["pvalueExponent"], a["pvalueMantissa"]))
under = [a for a in rows if a["pvalue"] == 0.0]
assert best["pvalue"] == 0.0 and best in under
print(f"4 p-value underflow: strongest {best['pvalueMantissa']}e"
      f"{best['pvalueExponent']} -> pvalue {best['pvalue']}, "
      f"{len(under)} of {len(rows)} underflowed")

# 5. INVARIANT: paging is driven by _links.next, and the walk is exhaustive.
status, doc = get("studies/search/findByEfoTrait", efoTrait="eye color", size=5)
total, seen, pages = doc["page"]["totalElements"], [], 0
while True:
    pages += 1
    seen += doc["_embedded"]["studies"]
    if "next" not in doc["_links"]:
        break
    status, doc = fetch(doc["_links"]["next"]["href"])
    assert status == 200, status
assert len(seen) == total, f"collected {len(seen)} of {total}"
assert pages == math.ceil(total / 5), pages
print(f"5 next-link paging : {pages} pages of 5 -> {len(seen)} of {total} studies")

# 6. INVARIANT: the trait label search is exact, and the two misses differ.
status, doc = get("efoTraits/search/findByEfoTrait", trait="type 2 diabetes")
assert status == 200 and doc["_embedded"]["efoTraits"] == [], "prefix match now?"
status, doc = get("efoTraits/search/findByEfoTrait",
                  trait="type 2 diabetes mellitus")
hit = doc["_embedded"]["efoTraits"][0]
assert hit["shortForm"] == "MONDO_0005148", hit["shortForm"]
assert hit["uri"].endswith("MONDO_0005148") and "purl.obolibrary.org" in hit["uri"]
retired, _ = get("efoTraits/search/findByShortForm", shortForm="EFO_0001360")
live, _ = get("efoTraits/search/findByShortForm", shortForm="EFO_0003949")
assert (retired, live) == (404, 200), (retired, live)
print(f"6 trait vocabulary : 'type 2 diabetes' -> 0 hits; exact label -> "
      f"{hit['shortForm']}; EFO_0001360 -> HTTP {retired}")

# 7. INVARIANT: the projection inlines the study, so ancestry travels with the
#    association -- and the Catalog curates below genome-wide significance.
status, doc = get("efoTraits/EFO_0003949/associations",
                  projection="associationByEfoTrait")
assoc = doc["_embedded"]["associations"]
assert all("study" in a and "ancestries" in a["study"] for a in assoc)
studies = {a["study"]["accessionId"]: a["study"] for a in assoc}
weak = [a for a in assoc
        if a["pvalueMantissa"] * 10.0 ** a["pvalueExponent"] > GWAS_SIG]
assert weak, "no sub-threshold associations -- curation threshold changed?"
euro = sum(x["numberOfIndividuals"] or 0
           for s in studies.values() for x in s["ancestries"]
           if x["type"] == "initial"
           and any(g["ancestralGroup"] == "European" for g in x["ancestralGroups"]))
disc = sum(x["numberOfIndividuals"] or 0
           for s in studies.values() for x in s["ancestries"]
           if x["type"] == "initial")
print(f"7 study context    : {len(assoc)} associations / {len(studies)} studies, "
      f"{len(weak)} weaker than 5e-8, {100 * euro / disc:.0f}% European discovery")

# 8. INVARIANT: an association is not always one variant, and the two
#    multi-variant forms hide their extra alleles in different places.
status, doc = get("studies/GCST000071/associations")
hap = [a for a in doc["_embedded"]["associations"] if a["multiSnpHaplotype"]]
assert hap, "GCST000071 no longer carries a haplotype association"
wide = max(hap, key=lambda a: len(a["loci"][0]["strongestRiskAlleles"]))
n_alleles = len(wide["loci"][0]["strongestRiskAlleles"])
assert len(wide["loci"]) == 1 and n_alleles > 1, (len(wide["loci"]), n_alleles)
assert wide["loci"][0]["haplotypeSnpCount"] == n_alleles, wide["loci"][0]
status, doc = get("studies/GCST006494/associations")
inter = [a for a in doc["_embedded"]["associations"] if a["snpInteraction"]]
assert inter, "GCST006494 no longer carries an interaction association"
assert all(len(a["loci"]) > 1 for a in inter), "interaction is no longer multi-locus"
assert all(len(L["strongestRiskAlleles"]) == 1
           for a in inter for L in a["loci"]), "interaction locus grew alleles"
print(f"8 multi-variant    : haplotype = 1 locus x {n_alleles} alleles; "
      f"interaction = {len(inter[0]['loci'])} loci x 1 allele")

Expect — two kinds of line, and they fail for different reasons.

Invariants. These hold whatever the release, so a failure means the skill is wrong, not stale: three response shapes distinguished by the presence of page and _embedded; page carrying exactly number, size, totalElements, totalPages; size ignored on a sub-resource; pvalueMantissa in 1–9 with a negative exponent; pvalue exactly 0.0 at or below exponent −324 and matching the reconstruction above it; a _links.next walk collecting exactly totalElements; exact-match trait labels; HTTP 404 from findByShortForm versus HTTP 200 with an empty list from findByEfoTrait; EFO_0001360 retired and type 2 diabetes living at MONDO_0005148; every associationByEfoTrait row carrying an inlined study with ancestries; at least one association weaker than 5×10⁻⁸; and a haplotype association being one locus with many risk alleles while an interaction is many loci with one each.

Observed values, 2026-08-27. These move as the Catalog grows — a mismatch is drift to investigate, not a bug. Study and association counts rise every fortnight; the eye colour trait had 201 associations across 16 studies; GCST012219 had 26 associations, 2 of them underflowed; the widest haplotype in GCST000071 carried 16 risk alleles.

1 paged shape      : _embedded.studies + page, 229,952 studies
2 single resource  : GCST012219, 11,532,091 variants tested, fullPvalueSet=False
3 unpaged shape    : size=1 returned all 26 associations
4 p-value underflow: strongest 2e-28539 -> pvalue 0.0, 2 of 26 underflowed
5 next-link paging : 4 pages of 5 -> 16 of 16 studies
6 trait vocabulary : 'type 2 diabetes' -> 0 hits; exact label -> MONDO_0005148; EFO_0001360 -> HTTP 404
7 study context    : 201 associations / 16 studies, 41 weaker than 5e-8, 99% European discovery
8 multi-variant    : haplotype = 1 locus x 16 alleles; interaction = 2 loci x 1 allele

Limits worth stating in a write-up

  • The Catalog records what was published, not what is true. It is a curated literature index. An association appears because a paper reported it at 1×10⁻⁵ or stronger, not because anyone re-analysed the data.
  • It is overwhelmingly European. Say so with a number from the study records you actually retrieved, not as a generic caveat.
  • Curated associations are the reported lead variant, not a fine-mapped causal one. The risk allele is in linkage disequilibrium with whatever is causal, and authorReportedGenes and genomicContexts are proximity, not mechanism.
  • Effect sizes are not comparable across studies without care — different covariates, different phenotype definitions, different reference alleles, and pvalueDescription sub-phenotypes pooled under one ontology term.
  • Sample sizes overlap. Many studies draw on the same biobanks, so summing numberOfIndividuals across studies of a trait counts people more than once. The 99% figure above is a composition estimate, not a headcount.
  • State the release. A new release lands roughly every fortnight and the totals move with it — 229,952 studies and 1,191,948 associations on 2026-08-27, against a release dated 2026-08-24. Any count you report is a dated observation, not a constant, which is what the manifest in Get the files exists to record.

Sources

The GWAS Catalog is a collaboration between the National Human Genome Research Institute and EMBL-EBI. Cite Cerezo et al. (2025) for the resource, and cite the original study behind any association you report — the PubMed id travels with every study record for exactly that reason.