dbsnp-database

---
name: dbsnp-database
description: >
license: Apache-2.0
---

# dbSNP Database Integration

## Prerequisites

1.  **`uv`**: Read the `uv` skill and follow its Setup instructions to ensure
    `uv` is installed and on PATH.
2.  **User Notification**: If .licenses/dbsnp_database_LICENSE.txt does not
    already exist in the workspace root directory then (1) prominently notify
    the user to check the terms at https://www.ncbi.nlm.nih.gov/snp/, then (2)
    create the file recording the notification text and timestamp.
3.  **`.env` file**: Make sure the `.env` file exists in your home directory.
    Create one if it does not exist.
4.  **`NCBI_API_KEY`** (optional): Raises the NCBI rate limit from 3 to 10
    requests/second. The skill works without it, but a key is recommended if the
    user plans many queries or encounters a 429 error. You can register for a
    key for free at https://www.ncbi.nlm.nih.gov/account/settings/. You **MUST**
    use the safe credentials protocol in the `credentials` skill to check for
    and request this key if this skill looks relevant to the user's request.

## Core Rules

-   **Use the Wrapper**: ALWAYS execute the provided wrapper script
    `scripts/dbsnp_cli.py` to query the database rather than constructing custom
    HTTP or curl requests. The script automatically handles rate limiting,
    retries, and JSON parsing.
-   **Command Choice**: Do NOT use `search-region` to find the rsID of a
    specific variant; use `resolve-variant` instead.
-   **Output Size**: Avoid using `--full` on `get-variant` unless specifically
    needed, as raw payloads can exceed 1 MB.
-   **Shell Safety**: Always wrap HGVS strings in single quotes to prevent shell
    expansion errors.
-   **Notification**: If this skill is used, ensure this is mentioned in the
    output.

## When to Use

**Use this skill when you need to:**

-   Map a genomic variant to its canonical rsID (from VCF coordinates or HGVS
    notation).
-   Retrieve summary data for an rsID: variant type, gene associations, clinical
    significance, and population allele frequencies.
-   Convert an rsID back to genomic coordinates on a specific assembly.
-   Find all known variants within a chromosomal region.

**Do NOT use when you need to:**

-   Obtain clinical pathogenicity classifications with submitter rationales (use
    **clinvar-database**).
-   Get precise population-level allele frequencies stratified by ancestry (use
    **gnomad-database**).
-   Predict the functional effect of a novel mutation (use
    **alphagenome-single-variant-analysis**).
-   View 3D protein structures affected by a variant (use
    **alphafold-database-fetch-and-analyze / pdb-database**).

## Command Selection Guide

**Pick the right command on the first try.** Match the user's input to the
correct subcommand below — one command call is almost always sufficient.

-   User gives you…: Run this command
-   An rsID (e.g. `rs7412`, `rs268`): `get-variant`
-   Genomic coordinates: chrom pos ref alt (e.g. `8 19962213 C T`):
    `resolve-variant`
-   An HGVS string (e.g. `NC_000008.11:g.19962213del`): `resolve-hgvs`
-   An rsID and they want coordinates back: `resolve-rsid`
-   A chromosomal region (chrom start end): `search-region`

> [!CAUTION] **Do NOT use `search-region` to find the rsID of a specific
> variant.** If the user provides a chromosome, position, reference allele, and
> alternate allele (four values), use `resolve-variant` — it is a direct,
> single-API-call lookup. `search-region` is only for surveying all variants
> within a positional range and returns hundreds/thousands of results.

## Quick Start

```bash
# Look up variant rs7412: type, gene, clinical significance, MAF
uv run scripts/dbsnp_cli.py get-variant rs7412 --output /tmp/rs7412.json

# Find the rsID for a variant at chr8:19962213 C>T
uv run scripts/dbsnp_cli.py resolve-variant 8 19962213 C T \
  --output /tmp/resolve.json
```

All subcommands write JSON to disk. Always save output in the `/tmp/` directory.
The `--output` flag is required.

## Commands

### 1. `get-variant` — Fetch Variant Record

Retrieve the RefSNP record for one rsID. By default the output is abbreviated to
the most useful fields. Both `rs268` and `268` are accepted.

```bash
uv run scripts/dbsnp_cli.py get-variant rs268 --output /tmp/rs268.json
uv run scripts/dbsnp_cli.py get-variant 268 --assembly GCF_000001405.40 \
  --output /tmp/rs268.json
```

*Arguments:*

-   `rsid` (positional, required): The RefSNP identifier.
-   `--assembly`: RefSeq assembly accession (default: `GCF_000001405.40` =
    GRCh38).
-   `--full`: Return the complete raw JSON payload — see warning below.
-   `--output`: Output file path (default: `/tmp/dbsnp_output.json`).

*Abbreviated output fields:*

-   `refsnp_id`: Numeric rsID
-   `variant_type`: e.g. `snv`, `ins`, `del`, `delins`
-   `genes`: Sorted list of gene symbols (locus names)
-   `clinical_significances`: List of clinical significance labels
-   `minor_allele_frequencies`: Study name, allele count, total count
-   `placements`: Genomic placements for the requested assembly

> [!WARNING] **About `--full`:** The raw RefSNP payload is typically 50–500 KB
> and can exceed 1 MB for clinically significant variants with many submissions.
> Only use `--full` when you specifically need data absent from the abbreviated
> output — for example:
>
> -   The complete HGVS nomenclature across every transcript and protein
>     isoform.
> -   Full submission history with individual submitter details and timestamps.
> -   Population-level allele frequency breakdowns by sub-population within a
>     study (e.g. per-population gnomAD counts).
> -   The full set of genomic placements across multiple assemblies (GRCh37 and
>     GRCh38 simultaneously).
> -   Merge history showing which older rsIDs were merged into this one.

### 2. `resolve-variant` — Genomic Coordinates → rsID

Determine the rsID(s) for a variant given its genomic coordinates (chromosome,
position, reference allele, alternate allele). **This is the command to use when
the user provides a variant as space-separated coordinates** like `8 19962213 C
T`.

```bash
uv run scripts/dbsnp_cli.py resolve-variant 8 19962213 C T \
  --output /tmp/resolve.json
```

*Arguments:*

-   `chrom` (positional): Chromosome number (e.g. `8`) or RefSeq sequence
    accession (e.g. `NC_000008.11`). **Chromosomes X and Y must be passed as
    their numeric equivalents: `23` for X and `24` for Y.**
-   `pos` (positional): 1-based genomic position.
-   `ref` (positional): Reference allele (e.g. `C`).
-   `alts` (positional): Alternate allele(s), comma-separated (e.g. `T`).
-   `--assembly`: RefSeq assembly accession (default: `GCF_000001405.40`).
-   `--output`: Output file path (default: `/tmp/dbsnp_output.json`).

*Output:* `{"rsids": ["12345", "67890"]}`

### 3. `resolve-rsid` — rsID → Genomic Coordinates

Get the genomic placement (sequence ID and allele details) for a known rsID on a
specific assembly.

```bash
uv run scripts/dbsnp_cli.py resolve-rsid rs7412 --output /tmp/coords.json
```

*Arguments:*

-   `rsid` (positional): The RefSNP identifier.
-   `--assembly`: RefSeq assembly accession (default: `GCF_000001405.40`).
-   `--output`: Output file path (default: `/tmp/dbsnp_output.json`).

*Output:* `{"rsid": "7412", "assembly": "...", "placements": [...]}`

### 4. `resolve-hgvs` — HGVS → rsID

Find the rsID(s) corresponding to an HGVS expression.

```bash
uv run scripts/dbsnp_cli.py resolve-hgvs 'NC_000008.11:g.19962213del' \
  --output /tmp/hgvs.json
```

*Arguments:*

-   `hgvs` (positional): The HGVS string.
-   `--assembly`: RefSeq assembly accession (default: `GCF_000001405.40`).
-   `--output`: Output file path (default: `/tmp/dbsnp_output.json`).

*Output:* `{"rsids": ["12345"]}`

> [!TIP] HGVS strings often contain characters that shells interpret (colons,
> greater-than signs). Always wrap them in **single quotes** to prevent shell
> expansion.

### 5. `search-region` — Regional Variant Search

Find all rsIDs within a bounded chromosomal region.

```bash
uv run scripts/dbsnp_cli.py search-region 7 117100000 117300000 \
  --output /tmp/region.json
```

*Arguments:*

-   `chrom` (positional): Chromosome (e.g. `7`). **Use `23` for chromosome X and
    `24` for chromosome Y.**
-   `start` (positional): Start position.
-   `end` (positional): End position.
-   `--retmax`: Maximum rsIDs to return (default: 500, ceiling: 5 000).
-   `--output`: Output file path (default: `/tmp/dbsnp_output.json`).

*Output:*

```json
{
  "rsids": ["12345", "67890", "..."],
  "returned": 500,
  "total_available": 1423,
  "truncated": true,
  "note": "Only 500 of 1423 variants returned.  Increase --retmax ..."
}
```

When `total_available` exceeds the returned count, the output includes a
`truncated` flag and a `note`. Increase `--retmax` to retrieve more (up to 5
000).

## Typical Workflows

### Identify a known variant from coordinates

```bash
# Step 1: Map VCF coordinates to rsID
uv run scripts/dbsnp_cli.py resolve-variant 19 44908684 T C \
  --output /tmp/step1.json

# Step 2: Get the full details for the resolved rsID
uv run scripts/dbsnp_cli.py get-variant <rsid_from_step1> \
  --output /tmp/step2.json
```

### Survey variants in a gene region

```bash
# Step 1: Find all variants in a region spanning the CFTR gene
uv run scripts/dbsnp_cli.py search-region 7 117100000 117300000 \
  --retmax 1000 --output /tmp/region.json

# Step 2: Retrieve details on individual rsIDs of interest
uv run scripts/dbsnp_cli.py get-variant <rsid> --output /tmp/detail.json
```

### Translate HGVS notation to genomic coordinates

```bash
# Step 1: Get the rsID for an HGVS expression
uv run scripts/dbsnp_cli.py resolve-hgvs 'NC_000019.10:g.44908684T>C' \
  --output /tmp/hgvs.json

# Step 2: Resolve that rsID to VCF-style coordinates
uv run scripts/dbsnp_cli.py resolve-rsid <rsid> --output /tmp/coords.json
```

## Assembly Defaults and Automatic Fallback

The Variation Services endpoints (used by `get-variant`, `resolve-variant`,
`resolve-rsid`, `resolve-hgvs`) expect a **RefSeq assembly accession**. The
RefSeq accession for GRCh38 is `GCF_000001405.40`, and for GRCh37 it is
`GCF_000001405.25`.

The `search-region` subcommand always searches GRCh38 positions.

> [!IMPORTANT] **Automatic assembly fallback:** The `resolve-variant` and
> `resolve-hgvs` commands automatically try GRCh38 first. If no rsIDs are found,
> they retry with GRCh37 before reporting failure. When a fallback occurs the
> output JSON includes a `"note"` field explaining which assembly succeeded.
> **You do NOT need to manually retry with a different assembly** — the script
> handles this transparently.

You only need to override `--assembly` when you specifically want to
**restrict** the lookup to one assembly (e.g. because the user's coordinates are
known to be GRCh37).

## NCBI API Key and Rate Limiting

Without an API key the script is limited to **3 requests per second**. With a
key this increases to **10 requests per second**.

You can register for a key for free at
https://www.ncbi.nlm.nih.gov/account/settings/. You **MUST** use the safe
credentials protocol in the `credentials` skill to check for and request this
key if this skill looks relevant to the user's request.

## Troubleshooting HTTP 500 Errors

### Reference Allele Mismatch

If you receive an HTTP 500 error with a message detailing that the asserted
reference allele is not equal to the reference sequence:

**What it means:** The coordinate position is likely valid, but the reference
allele (`ref`) you provided does not match the base at that position in the
requested assembly.

**Action:** 1. **DO NOT RETRY** the exact same query mechanically. 2. **Check
the assembly**: Coordinates are assembly-specific. 3. **Switch assembly**: If
you were querying GRCh37, try GRCh38 (using `--assembly GCF_000001405.40`), or
if querying GRCh38, try GRCh37 (using `--assembly GCF_000001405.25`).

## Common Mistakes

-   **Mistake:** Forgetting to quote HGVS strings **Fix:** Wrap in single
    quotes: `'NC_000008.11:g.19962213del'`

-   **Mistake:** Passing a chromosome name to `resolve-variant` instead of a
    sequence accession **Fix:** Use the numeric chromosome ID (e.g. `8`) or a
    RefSeq accession like `NC_000008.11`

-   **Mistake:** Using `--full` on `get-variant` without needing it **Fix:** The
    abbreviated output covers most use cases; `--full` returns 50–500 KB+ of
    JSON

-   **Mistake:** Expecting `search-region` to return all results by default
    **Fix:** The default `--retmax` is 500; check `total_available` in the
    output to see if results were truncated

-   **Mistake:** Using GRCh37 coordinates with `search-region` **Fix:**
    `search-region` always uses GRCh38 positions; lift over coordinates first if
    starting from GRCh37

-   **Mistake:** Manually retrying `resolve-variant` or `resolve-hgvs` with a
    different `--assembly` when the first call fails **Fix:** The script
    automatically tries GRCh38 then GRCh37; a single call is sufficient

-   **Mistake:** Passing `X` or `Y` as the chromosome value **Fix:** Use the
    numeric equivalents: `23` for chromosome X and `24` for chromosome Y. The
    CLI treats chromosomes numerically by default.