How to Look Up Your Own Variants in NCBI dbSNP

By the end of this you will have your own variant calls annotated with rsIDs and population allele frequencies from a local copy of the dbSNP release VCF, a scripted way to resolve any single rsID through NCBI’s Variation Services API, and a clear sense of which rsIDs in your file are trustworthy identifiers versus artifacts of a 25-year-old submission archive. You need your variants as a VCF (from whole-genome or exome sequencing) or a genotyping array export. For software you need bcftools 1.17 or newer, htslib/tabix, curl, jq, and Python 3. You need the GRCh38 reference FASTA you called against. You also need about 60 GB of free disk if you want the full dbSNP VCF locally, or none if you only need single-variant lookups. Everything below is measurement and interpretation. It does not tell you what any variant means for your health, and for anything that looks clinically relevant you want a clinician or a certified genetic counselor reading it with you.

1. Understand what an rsID identifies

An rsID is an accession assigned to a cluster of submitted variant records (ss numbers) that dbSNP believes describe the same change at the same place. An rsID is not a variant. Submitters send ss records. Then dbSNP maps them to the reference assembly, clusters co-located submissions, and mints an rs accession for the cluster. The Japanese Millennium Genome Project is an example. It submitted 190,562 gene-based variations discovered by resequencing, which then became ss records that clustered into rs accessions.¹ The early public map of 1.42 million SNPs came from a similar pooling of independently discovered sites.² That history is why dbSNP behaves the way it does.

Three consequences you will hit in practice:

• An rsID can cover more than one allele. rs334 is the classic sickle-cell position, but plenty of rs accessions are multi-allelic clusters. An rsID plus a genotype string is ambiguous unless you also know the reference and the alternate allele on a named assembly.
• rsIDs get merged and withdrawn. When dbSNP decides two clusters are the same, the higher-numbered rs is merged into the lower-numbered one. Old array manifests and old papers are full of rsIDs that now redirect.
• dbSNP is an archive. No one curates it as a truth set. A variant having an rsID tells you someone once submitted it. It does not mean the variant is common, real, or validated. Build 156 contains over a billion rs accessions, the vast majority of them rare or singleton observations from large sequencing cohorts.

If you want a stable, unambiguous identifier for a variant, use a normalized CHROM POS REF ALT tuple on a named assembly. The alternative is NCBI’s SPDI notation, which stands for sequence, position, deletion and insertion. As for rsIDs, they are a convenient join key and a way into the literature. They are a poor primary key.

2. Normalize your variants first

Annotation is a coordinate join. If your representation of an indel differs from dbSNP’s by one base of left-alignment, you get no match and no warning. Normalize before you annotate, always.

bcftools norm \
  -f GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  -m -both \
  -c s \
  --threads 4 \
  -Oz -o mine.norm.vcf.gz \
  mine.vcf.gz

bcftools index -t mine.norm.vcf.gz

-m -both splits multi-allelic records into one ALT per line, which is what you want for a per-allele join. -c s sets a mismatching REF base to the reference and warns, rather than silently proceeding. Check the stderr summary: if bcftools norm reports a large number of REF mismatches, you are on the wrong assembly or the wrong reference build (the no-ALT analysis set versus the full GRCh38 with patches matters for a handful of regions, and hs37d5 versus GRCh37 matters everywhere).

Confirm the assembly before going further:

bcftools view -h mine.norm.vcf.gz | grep -E '^##(reference|contig=<ID=(chr)?1,)'

If the length of chromosome 1 is 249,250,621 you are on GRCh37 and need to lift over. If it is 248,956,422 you are on GRCh38.

For GRCh37 input, CrossMap is the tool we reach for, because it handles VCF properly. It will flip REF/ALT when the strand changes, and drop records it cannot map:

CrossMap vcf hg19ToHg38.over.chain.gz mine.norm.vcf.gz \
  GCA_000001405.15_GRCh38_no_alt_analysis_set.fna \
  mine.hg38.vcf

Expect to lose roughly 0.1–0.5% of sites as unmapped, concentrated in segmental duplications and the pseudoautosomal boundaries. Re-run bcftools norm after lifting.

3. Download the dbSNP release VCF and fix the chromosome names

NCBI publishes the per-build VCF keyed by RefSeq assembly accession:

BASE=https://ftp.ncbi.nlm.nih.gov/snp/latest_release/VCF
curl -O $BASE/GCF_000001405.40.gz        # GRCh38
curl -O $BASE/GCF_000001405.40.gz.tbi
curl -O $BASE/GCF_000001405.40.gz.md5
md5sum -c GCF_000001405.40.gz.md5

GCF_000001405.40 is GRCh38.p14. GCF_000001405.25 is GRCh37.p13 if you must stay on the older assembly. The compressed file is in the tens of gigabytes. Download it once, keep it, and do not decompress it.

The catch: this VCF uses RefSeq sequence accessions in CHROM (NC_000001.11, NC_000023.11, NT_187361.1 for unplaced contigs), while your VCF almost certainly uses chr1 or 1. bcftools annotate will happily return zero annotations and exit 0. Build a rename map from the assembly report:

curl -sO https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.40_GRCh38.p14/GCF_000001405.40_GRCh38.p14_assembly_report.txt

awk -F'\t' '!/^#/ && $10 != "na" { print $7"\t"$10 }' \
  GCF_000001405.40_GRCh38.p14_assembly_report.txt > refseq2ucsc.txt
head -3 refseq2ucsc.txt
# NC_000001.11    chr1
# NC_000002.12    chr2

Column 7 is the RefSeq accession, column 10 is the UCSC-style name. If your VCF uses bare 1, 2 and X, use column 1 (the sequence name from the assembly) instead of column 10.

Rather than rewriting the 40 GB annotation file, extract only the regions you need. For a whole genome you do need the whole thing, so rename once and cache it:

bcftools annotate --rename-chrs refseq2ucsc.txt \
  --threads 8 -Oz -o dbsnp156.hg38.vcf.gz GCF_000001405.40.gz
bcftools index -t --threads 8 dbsnp156.hg38.vcf.gz

This takes a few hours on a laptop and produces a file the same size again. If disk is tight, rename your VCF into RefSeq space instead (it is thousands of times smaller) and rename back afterward.

4. Annotate with rsIDs and frequencies

Look at the dbSNP header before you pick fields, because INFO tags change between builds:

bcftools view -h dbsnp156.hg38.vcf.gz | grep '^##INFO'

You will typically find these fields.

• RS (the numeric accession)
• VC (variation class: SNV, DEL, INS or MNV)
• GENEINFO (gene symbol:GeneID pairs)
• FREQ (per-study allele frequencies, pipe-delimited by study name)
• COMMON
• ClinVar-derived fields when the record has a ClinVar submission

Then:

bcftools annotate \
  -a dbsnp156.hg38.vcf.gz \
  -c ID,INFO/RS,INFO/VC,INFO/GENEINFO,INFO/FREQ \
  --threads 4 \
  -Oz -o mine.dbsnp.vcf.gz \
  mine.norm.vcf.gz

bcftools index -t mine.dbsnp.vcf.gz

bcftools annotate matches on CHROM, POS, REF, and ALT by default when transferring ID, which is the behavior you want. If you pass -c ID and see rsIDs appearing on variants whose ALT allele is not in dbSNP, you have a normalization mismatch. The cause is not a bug.

Measure your match rate. For a 30x short-read genome, a plausible result is 97–99% of SNVs carrying an rsID and 85–95% of indels, with the indel gap driven by representation differences and by genuinely novel calls:

bcftools query() { bcftools view -H "$@" | wc -l; }
total=$(bcftools view -H -v snps mine.dbsnp.vcf.gz | wc -l)
known=$(bcftools view -H -v snps -i 'ID!="."' mine.dbsnp.vcf.gz | wc -l)
python3 -c "print(f'SNV rsID rate: {$known/$total:.4f}')"

A rate below about 95% for SNVs in a genome from a well-sampled population usually means a coordinate or normalization problem, or a high false-positive rate in the calls themselves. The original demonstration of accurate whole-genome short-read sequencing used concordance with known polymorphic sites as one of its primary quality checks, and the logic still holds: novel-variant fraction is a sensitive detector of calling artifacts.³

Pull out the FREQ field for the sites you care about. It is awkward to parse because it packs multiple studies into one string:

bcftools query -f '%CHROM\t%POS\t%ID\t%REF\t%ALT\t%INFO/FREQ\n' \
  -r chr11:5225000-5230000 mine.dbsnp.vcf.gz

We treat dbSNP FREQ as a rough orientation check rather than the frequency source of record. It aggregates submitted studies with uneven ancestry coverage. The studies are ALFA, TOPMED, 1000Genomes and GnomAD projections. For frequency work, annotate separately against a gnomAD release VCF and use its per-population AF fields with their own filter flags.

5. Query single rsIDs through the Variation Services API

For one-off lookups, skip the local file. The Variation Services API returns the full JSON refsnp record:

curl -s 'https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/334' | jq '
  {
    rs: .refsnp_id,
    status: .merge_status,
    placements: [
      .primary_snapshot_data.placements_with_allele[]
      | select(.seq_id | startswith("NC_"))
      | {seq: .seq_id,
         hgvs: [.alleles[].hgvs]}
    ]
  }'

Note the endpoint takes the bare number, not rs334. Useful fields in the JSON:

• primary_snapshot_data.placements_with_allele[] — one entry per sequence the variant is placed on, including RefSeqGene and transcript sequences. Filter to NC_ accessions for genomic coordinates, and match the assembly by accession version.
• .alleles[].allele.spdi holds SPDI tuples. Each tuple has seq_id, position as 0-based interbase, deleted_sequence and inserted_sequence. This is the representation to store if you want something that survives a dbSNP build change.
• primary_snapshot_data.allele_annotations[].frequency[] holds per-study counts, with study_name, allele_count and total_count. The counts let you judge how much to trust a number. Frequencies alone do not.
• primary_snapshot_data.allele_annotations[].clinical[] — ClinVar links with clinical_significances and review_status.

For search rather than direct retrieval, E-utilities still works against db=snp:

curl -s 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=HBB[GENE]+AND+missense_variant[FXN]&retmax=200&retmode=json' \
  | jq -r '.esearchresult.idlist[]'

Rate limits: 3 requests/second without an API key, 10/second with one. Pass &api_key=$NCBI_API_KEY and put a time.sleep(0.11) between calls. If you need more than a few thousand lookups, use the local VCF or the JSON refsnp dumps under /snp/latest_release/JSON/ instead of hammering the API.

6. Resolve merged and withdrawn rsIDs

Old array manifests, old GWAS catalogs, and old papers cite rsIDs that no longer exist as primary accessions. The API tells you where they went:

rs_resolve () {
  curl -s "https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/$1" \
    | jq -r 'if .merged_snapshot_data then
               "rs\(.refsnp_id) MERGED_INTO rs\(.merged_snapshot_data.merged_into[0])"
             elif .present_obs_movements then
               "rs\(.refsnp_id) MOVED"
             else
               "rs\(.refsnp_id) LIVE \(.merge_status // "n/a")"
             end'
}
rs_resolve 1799853

A record can also come back with "status": "withdrawn" or a snapshot_release_datetime and no primary_snapshot_data, which means dbSNP retired it (usually because the underlying submission turned out to be an assembly artifact or a false mapping). Treat withdrawn rsIDs the way you would treat a retracted paper: drop them, and be suspicious of any downstream result that depended on one.

Batch resolution against a list of legacy IDs, with rate limiting:

import json, time, urllib.request

API = "https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/{}"

def resolve(rsnum, key=None):
    url = API.format(rsnum) + (f"?api_key={key}" if key else "")
    with urllib.request.urlopen(url, timeout=30) as fh:
        rec = json.load(fh)
    if "merged_snapshot_data" in rec:
        return ("merged", rec["merged_snapshot_data"]["merged_into"])
    if "primary_snapshot_data" not in rec:
        return ("withdrawn_or_unsupported", [])
    return ("live", [rec["refsnp_id"]])

for rs in open("legacy_rsids.txt"):
    rs = rs.strip().lstrip("rs")
    print(rs, resolve(rs))
    time.sleep(0.12)

7. Join array genotypes to dbSNP without flipping strands

A consumer array export is tab-separated, with columns rsid, chromosome, position, genotype, and is on GRCh37 for most vendors. If that is your input, you have two extra problems.

First, the file’s rsIDs come from a manifest that may be several dbSNP builds old. Resolve merges before joining, or you will silently lose sites.

Second, strand. Array genotypes are reported on a strand defined by the probe design, and while most vendors report on the reference plus strand, the exported allele letters for A/T and C/G sites are ambiguous: you cannot tell A/T on the plus strand from A/T on the minus strand by inspection. The unambiguous sites are A/G, A/C, T/G and T/C. For those, a mismatch between your genotype letters and the dbSNP REF/ALT pair means a flip, and you can correct by complementing. For the ambiguous A/T and C/G sites, the only reliable resolution is allele frequency comparison against a matched reference panel, and it fails when the minor allele frequency is near 0.5. We drop ambiguous-strand sites with MAF above roughly 0.4 rather than guess. Roughly 15–20% of array markers are A/T or C/G, so this is not a rounding error.

Converting the export to a VCF is best done with a tool that carries the manifest’s strand information, not by hand. plink2 --23file will read a 23andMe export and write a VCF against a reference, and bcftools +fixref can then check and correct orientation against your reference FASTA:

plink2 --23file genome.txt SAMPLE1 --recode vcf --out mine_array
bcftools +fixref mine_array.vcf -- -f GRCh37.fa -m top   # report, then fix

Run bcftools +fixref in report mode first (no -m). If it reports more than a percent or two of non-reference-allele mismatches, stop and fix the assembly before doing anything else.

8. Get ClinVar interpretation, and know where your job ends

dbSNP links to ClinVar but is not ClinVar. The significance label alone is not enough. For any variant where interpretation matters, query ClinVar directly and read the review status too:

curl -s 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=clinvar&term=rs334[VARIANT_ID]&retmode=json'

Or annotate in bulk from the ClinVar VCF (https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz), pulling CLNSIG, CLNREVSTAT, and CLNDN. A CLNREVSTAT of no_assertion_criteria_provided or criteria_provided,_single_submitter carries far less weight than reviewed_by_expert_panel or practice_guideline. Conflicting submissions are common.

This is the boundary. Annotating your genome against dbSNP and ClinVar is a data engineering exercise you can do well by yourself. Deciding what a variant means for you is clinical work. It depends on your phenotype, your family history, and the penetrance literature for that specific allele in your ancestry. It also depends on confirmation from a clinical-grade assay, because research-grade sequencing is not a diagnostic result. If an annotation in your own file worries you, take it to a clinical geneticist or a certified genetic counselor and ask for orthogonal confirmation before you conclude anything.

Common problems

Zero rsIDs annotated, exit code 0. Chromosome naming mismatch between your VCF and the dbSNP VCF. bcftools annotate does not warn. Check with bcftools index -s dbsnp156.hg38.vcf.gz | head and bcftools index -s mine.norm.vcf.gz | head and confirm the names match exactly.

Indels annotated at a much lower rate than SNVs. Either you skipped bcftools norm -f, or you normalized against a different reference FASTA than dbSNP’s. Left-alignment differences of a single base are enough to break the join.

REF mismatch warnings everywhere. Assembly or build mismatch. Compare contig lengths in your header against the reference. GRCh37 with hs37d5 decoys, GRCh38 with and without ALT contigs, and the T2T-CHM13 assembly all look superficially similar and will not join.

An rsID from a paper returns HTTP 404 or an empty record. Withdrawn accession, or you passed rs334 where the endpoint wants 334.

FREQ frequencies disagree with gnomAD. Expected. dbSNP aggregates submitted studies with different ancestry compositions and different filtering. Use one frequency source per analysis and state which.

Tabix errors on the dbSNP VCF. The published .tbi is built for the RefSeq-named file. If you renamed chromosomes, you must re-index. Also, the file is large enough that .tbi is near its coordinate limits on some contigs. If indexing fails, write .csi instead: bcftools index -c.

Memory blowup during annotation. bcftools annotate streams, so memory should stay under a gigabyte. If it does not, you are probably annotating an unsorted VCF. Run bcftools sort -m 4G first.

An rsID maps to multiple genomic positions. Real, and common in segmental duplications and the pseudoautosomal regions. Check placements_with_allele in the JSON for more than one NC_ placement on your assembly, and treat genotypes at those sites as unreliable from short-read data.

Footnotes

1. Hiroaki Haga, Ryo Yamada, Yoshifumi Ohnishi, et al. Gene-based SNP discovery as part of the Japanese Millennium Genome Project: identification of 190 562 genetic variations in the human genome. Journal of Human Genetics, 2002. https://doi.org/10.1007/s100380200092 ↩

2. The International SNP Map Working Group, Cold Spring Harbor Laboratories:, Ravi Sachidanandam, et al. A map of human genome sequence variation containing 1.42 million single nucleotide polymorphisms. Nature, 2001. https://doi.org/10.1038/35057149 ↩

3. David Bentley, Shankar Balasubramanian, Harold Swerdlow, et al. Accurate whole human genome sequencing using reversible terminator chemistry. Nature, 2008. https://doi.org/10.1038/nature07517 ↩