Dbsnp Database

pip install requests pandas matplotlib
# xml.etree.ElementTree is part of Python stdlib — no additional install needed

import requests
import json

EMAIL = "your@email.com"      # required by NCBI policy
BASE_EUTILS = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"
BASE_VARIATION = "https://api.ncbi.nlm.nih.gov/variation/v0"

def fetch_snp_by_rsid(rsid: str) -> dict:
    """Fetch a dbSNP record by rsID using the NCBI Variation Services API (structured JSON)."""
    rs_num = str(rsid).lstrip("rs")
    r = requests.get(f"{BASE_VARIATION}/refsnp/{rs_num}", timeout=15)
    r.raise_for_status()
    return r.json()

record = fetch_snp_by_rsid("rs1800497")  # DRD2 Taq1A
print(f"rsID: rs{record['refsnp_id']}")
print(f"Variant type: {record['primary_snapshot_data'].get('variant_type')}")
# Top-level keys: citations, create_date, dbsnp1_merges, last_update_build_id,
# last_update_date, lost_obs_movements, mane_select_ids, present_obs_movements,
# primary_snapshot_data, refsnp_id. (No top-level `organism` field.)
# rsID: rs1800497
# Variant type: snv

import requests
import xml.etree.ElementTree as ET

EMAIL = "your@email.com"
BASE = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def efetch_snp_xml(rsid: str) -> ET.Element:
    """Fetch dbSNP XML record for a single rsID via the docsum rettype.
    Note: rettype="xml" returns a namespaced ExchangeSet; rettype="docsum"
    returns the simpler eSummaryResult/DocumentSummary tree without namespaces."""
    rs_num = str(rsid).lstrip("rs")
    r = requests.get(f"{BASE}/efetch.fcgi",
                     params={"db": "snp", "id": rs_num,
                             "rettype": "docsum", "retmode": "xml",
                             "email": EMAIL},
                     timeout=20)
    r.raise_for_status()
    return ET.fromstring(r.text)

root = efetch_snp_xml("rs80357906")

# Parse the DocumentSummary record (MAF/MAFALLELE were removed in 2024;
# GLOBAL_MAFS is a sub-tree — use ESummary JSON below for easier access)
for docsum in root.iter("DocumentSummary"):
    rs_id = docsum.get("uid")
    snp_class = docsum.findtext("SNP_CLASS", "Unknown")
    chr_pos = docsum.findtext("CHRPOS", "N/A")
    clin_sig = docsum.findtext("CLINICAL_SIGNIFICANCE", "N/A")
    print(f"rs{rs_id} | Class: {snp_class} | Position: {chr_pos}")
    print(f"  ClinSig: {clin_sig}")
# rs80357906 | Class: delins | Position: 17:43057062
#   ClinSig: pathogenic,risk-factor,uncertain-significance

# Fetch using ESummary for structured JSON (preferred for batch)
def esummary_snp(rsid: str) -> dict:
    rs_num = str(rsid).lstrip("rs")
    r = requests.get(f"{BASE}/esummary.fcgi",
                     params={"db": "snp", "id": rs_num,
                             "retmode": "json", "email": EMAIL},
                     timeout=15)
    r.raise_for_status()
    result = r.json()["result"]
    return result.get(rs_num, {})

rec = esummary_snp("rs80357906")
print(f"rs{rec.get('snp_id')}:")
print(f"  Class       : {rec.get('snp_class')}")            # e.g., 'delins'
# `maf`/`mafallele` were removed from ESummary in 2024 — use `global_mafs`
# (list of {study, freq}) and pick a study (e.g., 'GnomAD_genomes') or the
# global aggregate ('TOPMED'/'1000Genomes').
for m in rec.get('global_mafs', [])[:4]:
    print(f"  MAF[{m['study']:18s}]: {m['freq']}")
print(f"  ChrPos      : {rec.get('chrpos')}")               # 17:43057062
print(f"  ClinSig     : {rec.get('clinical_significance')}")
print(f"  FxnClass    : {rec.get('fxn_class')}")

import requests

EMAIL = "your@email.com"
BASE = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def esearch_snp(query: str, retmax: int = 100) -> tuple[list, int]:
    """Search dbSNP using a query string. Returns (id_list, total_count)."""
    r = requests.get(f"{BASE}/esearch.fcgi",
                     params={"db": "snp", "term": query,
                             "retmax": retmax, "retmode": "json",
                             "email": EMAIL},
                     timeout=15)
    r.raise_for_status()
    result = r.json()["esearchresult"]
    return result["idlist"], int(result["count"])

# All variants in BRCA1
ids, total = esearch_snp("BRCA1[gene] AND human[orgn]", retmax=20)
print(f"BRCA1 variants in dbSNP: {total:,} total")
print(f"First 5 rsIDs: {['rs' + i for i in ids[:5]]}")

# Only clinical variants (linked to ClinVar)
ids_clin, total_clin = esearch_snp(
    "BRCA1[gene] AND human[orgn] AND clinsig[filter]", retmax=50)
print(f"BRCA1 variants with clinical significance: {total_clin:,}")

import requests

EMAIL = "your@email.com"
BASE = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def search_region(chrom: str, start: int, stop: int,
                  assembly: str = "GRCh38", retmax: int = 200) -> tuple[list, int]:
    """Find all dbSNP variants in a chromosomal region."""
    query = f"{chrom}[CHR] AND {start}:{stop}[CHRPOS37]" if assembly == "GRCh37" else \
            f"{chrom}[CHR] AND {start}:{stop}[CHRPOS]"
    r = requests.get(f"{BASE}/esearch.fcgi",
                     params={"db": "snp", "term": query,
                             "retmax": retmax, "retmode": "json",
                             "email": EMAIL},
                     timeout=20)
    r.raise_for_status()
    result = r.json()["esearchresult"]
    return result["idlist"], int(result["count"])

# PCSK9 exon 4 region (GRCh38)
ids, total = search_region("1", 55039700, 55040200)
print(f"Variants in chr1:55039700-55040200: {total:,} total")
print(f"Retrieved {len(ids)} rsIDs: {['rs' + i for i in ids[:5]]}")

import requests, json

EMAIL = "your@email.com"
BASE = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def fetch_snp_summaries(rsids: list) -> dict:
    """Fetch ESummary records for a list of rsIDs. Returns dict keyed by rs number."""
    ids_str = ",".join(str(r).lstrip("rs") for r in rsids)
    r = requests.post(f"{BASE}/esummary.fcgi",
                      data={"db": "snp", "id": ids_str,
                            "retmode": "json", "email": EMAIL},
                      timeout=20)
    r.raise_for_status()
    return r.json()["result"]

rsids = ["rs80357906", "rs80357220", "rs28897672", "rs1801133"]
result = fetch_snp_summaries(rsids)

for rs_num in rsids:
    uid = str(rs_num).lstrip("rs")
    rec = result.get(uid, {})
    # `global_mafs` is a list of {study, freq}; pick the first or filter by study
    gmafs = rec.get("global_mafs", [])
    maf_str = gmafs[0]["freq"] if gmafs else "N/A"
    maf_study = gmafs[0]["study"] if gmafs else ""
    print(f"\n{rs_num}:")
    print(f"  Class        : {rec.get('snp_class', 'N/A')}")
    print(f"  MAF          : {maf_str} (from {maf_study})")
    print(f"  Location     : {rec.get('chrpos', 'N/A')}")
    print(f"  ClinSig      : {rec.get('clinical_significance', 'N/A')}")
    print(f"  Function     : {rec.get('fxn_class', 'N/A')}")

import requests, time, pandas as pd

EMAIL = "your@email.com"
BASE = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def epost_ids(id_list: list) -> tuple[str, str]:
    """Upload rsIDs to NCBI history server. Returns (WebEnv, query_key)."""
    ids_str = ",".join(str(i).lstrip("rs") for i in id_list)
    r = requests.post(f"{BASE}/epost.fcgi",
                      data={"db": "snp", "id": ids_str, "email": EMAIL},
                      timeout=30)
    r.raise_for_status()
    import xml.etree.ElementTree as ET
    root = ET.fromstring(r.text)
    webenv = root.findtext("WebEnv")
    query_key = root.findtext("QueryKey")
    return webenv, query_key

def efetch_history(webenv: str, query_key: str,
                   retstart: int = 0, retmax: int = 100) -> dict:
    """Retrieve records from NCBI history server using ESummary."""
    r = requests.get(f"{BASE}/esummary.fcgi",
                     params={"db": "snp", "WebEnv": webenv,
                             "query_key": query_key, "retstart": retstart,
                             "retmax": retmax, "retmode": "json",
                             "email": EMAIL},
                     timeout=30)
    r.raise_for_status()
    return r.json()["result"]

rsid_batch = ["rs80357906", "rs80357220", "rs28897672", "rs1801133",
              "rs429358", "rs7412", "rs1800497"]

webenv, query_key = epost_ids(rsid_batch)
print(f"Posted {len(rsid_batch)} IDs | WebEnv: {webenv[:40]}...")

records = []
for start in range(0, len(rsid_batch), 100):
    result = efetch_history(webenv, query_key, retstart=start, retmax=100)
    for uid in result.get("uids", []):
        rec = result[uid]
        # global_mafs replaced maf/mafallele in 2024 — flatten the first entry
        gmafs = rec.get("global_mafs", [])
        records.append({
            "rsid": f"rs{uid}",
            "snp_class": rec.get("snp_class"),
            "maf": gmafs[0]["freq"] if gmafs else None,
            "maf_study": gmafs[0]["study"] if gmafs else None,
            "chrpos": rec.get("chrpos"),
            "clinical_sig": rec.get("clinical_significance"),
        })
    time.sleep(0.5)

df = pd.DataFrame(records)
print(f"\nRetrieved {len(df)} records:")
print(df.to_string(index=False))

import requests
import pandas as pd

BASE_VARIATION = "https://api.ncbi.nlm.nih.gov/variation/v0"

def fetch_refsnp(rsid: str) -> dict:
    """Fetch structured JSON from NCBI Variation Services API."""
    rs_num = str(rsid).lstrip("rs")
    r = requests.get(f"{BASE_VARIATION}/refsnp/{rs_num}", timeout=15)
    r.raise_for_status()
    return r.json()

def parse_allele_frequencies(record: dict) -> list:
    """Extract allele frequencies from a Variation Services record.
    freq[i] keys: study_name, study_version, local_row_id, observation, allele_count, total_count.
    The asserted allele identifier lives in freq.observation.inserted_sequence (NOT at ann.allele)."""
    freqs = []
    snapshot = record.get("primary_snapshot_data", {})
    for ann in snapshot.get("allele_annotations", []):
        for f in ann.get("frequency", []):
            obs = f.get("observation", {})
            freqs.append({
                "allele": obs.get("inserted_sequence"),
                "study": f.get("study_name"),
                "allele_count": f.get("allele_count"),
                "total_count": f.get("total_count"),
                "freq": (f["allele_count"] / f["total_count"]
                         if f.get("total_count") else None),
            })
    return freqs

record = fetch_refsnp("rs1800497")   # DRD2 Taq1A variant
print(f"rsID: rs{record['refsnp_id']}")
print(f"Variant type: {record['primary_snapshot_data'].get('variant_type')}")  # 'snv'

placements = record["primary_snapshot_data"].get("placements_with_allele", [])
for placement in placements[:2]:
    seq_id = placement.get("seq_id")
    for allele in placement.get("alleles", []):
        spdi = allele.get("allele", {}).get("spdi", {})
        print(f"  Placement: {seq_id} | SPDI: {spdi.get('inserted_sequence')}")

freqs = parse_allele_frequencies(record)
if freqs:
    df_freq = pd.DataFrame(freqs[:5])
    print(f"\nAllele frequencies ({len(freqs)} entries):")
    print(df_freq.to_string(index=False))

import requests, time, pandas as pd

EMAIL = "your@email.com"
BASE = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def epost_snp(rsids):
    ids = ",".join(str(r).lstrip("rs") for r in rsids)
    r = requests.post(f"{BASE}/epost.fcgi",
                      data={"db": "snp", "id": ids, "email": EMAIL}, timeout=30)
    r.raise_for_status()
    import xml.etree.ElementTree as ET
    root = ET.fromstring(r.text)
    return root.findtext("WebEnv"), root.findtext("QueryKey")

def esummary_history(webenv, query_key, start, retmax=100):
    r = requests.get(f"{BASE}/esummary.fcgi",
                     params={"db": "snp", "WebEnv": webenv,
                             "query_key": query_key, "retstart": start,
                             "retmax": retmax, "retmode": "json", "email": EMAIL},
                     timeout=30)
    r.raise_for_status()
    return r.json()["result"]

# Example: VCF post-processing — annotate a list of called variants
variant_rsids = [
    "rs80357906", "rs80357220", "rs28897672", "rs1801133",
    "rs429358", "rs7412", "rs1800497", "rs2230199",
]

print(f"Posting {len(variant_rsids)} rsIDs to NCBI history server...")
webenv, query_key = epost_snp(variant_rsids)

records = []
for start in range(0, len(variant_rsids), 100):
    result = esummary_history(webenv, query_key, start=start, retmax=100)
    for uid in result.get("uids", []):
        rec = result[uid]
        gmafs = rec.get("global_mafs", [])  # 2024 schema (replaced maf/mafallele)
        records.append({
            "rsid": f"rs{uid}",
            "snp_class": rec.get("snp_class"),
            "maf": gmafs[0]["freq"] if gmafs else None,
            "maf_study": gmafs[0]["study"] if gmafs else None,
            "chrpos_grch38": rec.get("chrpos"),
            "gene": rec.get("genes", [{}])[0].get("name") if rec.get("genes") else None,
            "clinical_significance": rec.get("clinical_significance"),
            "fxn_class": rec.get("fxn_class"),
        })
    time.sleep(0.5)

df = pd.DataFrame(records)
df.to_csv("variant_annotations.csv", index=False)
print(f"Annotated {len(df)} variants → variant_annotations.csv")
print(df[["rsid", "snp_class", "maf", "chrpos_grch38", "clinical_significance"]].to_string(index=False))

import requests, time
import pandas as pd
import matplotlib.pyplot as plt

EMAIL = "your@email.com"
BASE = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils"

def search_gene_variants(gene: str, retmax: int = 500) -> list:
    r = requests.get(f"{BASE}/esearch.fcgi",
                     params={"db": "snp", "term": f"{gene}[gene] AND human[orgn]",
                             "retmax": retmax, "retmode": "json", "email": EMAIL},
                     timeout=20)
    r.raise_for_status()
    return r.json()["esearchresult"]["idlist"]

def fetch_summaries_batch(ids: list, batch_size: int = 100) -> list:
    records = []
    for i in range(0, len(ids), batch_size):
        batch = ids[i:i+batch_size]
        ids_str = ",".join(batch)
        r = requests.post(f"{BASE}/esummary.fcgi",
                          data={"db": "snp", "id": ids_str,
                                "retmode": "json", "email": EMAIL},
                          timeout=30)
        r.raise_for_status()
        result = r.json()["result"]
        for uid in result.get("uids", []):
            rec = result[uid]
            records.append({"rsid": f"rs{uid}", "snp_class": rec.get("snp_class", "unknown")})
        time.sleep(0.5)
    return records

gene = "CFTR"
print(f"Searching dbSNP for {gene} variants...")
ids = search_gene_variants(gene, retmax=300)
print(f"Found {len(ids)} IDs; fetching summaries...")

records = fetch_summaries_batch(ids)
df = pd.DataFrame(records)
class_counts = df["snp_class"].value_counts()

# Plot variant class distribution
fig, ax = plt.subplots(figsize=(8, 5))
colors = ["#4472C4", "#ED7D31", "#A9D18E", "#FF0000", "#FFC000", "#7030A0"]
bars = ax.bar(class_counts.index, class_counts.values,
              color=colors[:len(class_counts)], edgecolor="white")
ax.bar_label(bars, padding=3, fontsize=9)
ax.set_xlabel("Variant Class")
ax.set_ylabel("Count")
ax.set_title(f"dbSNP Variant Classes in {gene} (n={len(df)})")
plt.tight_layout()
plt.savefig(f"{gene}_variant_classes.png", dpi=150, bbox_inches="tight")
print(f"Saved {gene}_variant_classes.png")
print(class_counts.to_string())
# snv                 241
# indel                38
# mnv                  12
# del                   7
# ins                   2