How to Look Up an rsID in Your Own DNA Data

Oak
rsid dbsnp vcf annotation clinvar genomics
A giant glowing fern frond in a dark forest where one leaflet shines amber among thousands of blue ones.

By the end of this guide you will have a normalized, indexed VCF of your own genotypes on a known reference build, a command that answers “what is my genotype at rs671?” in under a second, and an annotated copy of that file carrying dbSNP identifiers, predicted consequence, gnomAD allele frequencies, and ClinVar review status on every variant. You need your raw data (either a whole-genome VCF or an array export from a consumer testing service), roughly 60 GB of disk for the reference annotation sets, and bcftools, tabix, and samtools on your path. If you plan to annotate the whole file rather than a handful of positions, install Ensembl VEP with a local cache; the online endpoints are fine for one variant and hopeless for five million.

1. What an rsID is, precisely

An rsID (Reference SNP cluster ID, written rs followed by digits) is an accession number assigned by NCBI’s dbSNP to a cluster of submitted variant observations that map to the same position and represent the same change on the same reference assembly. The important consequence of that definition is that an rsID names a location and a set of alleles, not your genotype. rs671 refers to the G>A change at chr12:111,803,962 on GRCh38 in ALDH2; it says nothing about whether you carry zero, one, or two copies of the A allele. Half the confusion in rsID lookups comes from forgetting this.

Two further properties matter for anyone writing code against these identifiers. First, rsIDs are not stable in the way a UUID is: when dbSNP determines that two clusters describe the same variant it merges them, retiring one accession and redirecting it to the survivor, and occasionally an rsID is withdrawn entirely when the underlying submission turns out to be an artifact. Any lookup table you build from a single dbSNP build will silently fail on retired identifiers. Second, an rsID can carry more than two alleles, so a single rsID may correspond to several rows in a properly decomposed VCF. Despite these wrinkles the accession remains the common key across the literature, clinical trial registries, and variant databases, which is why tools that mine those sources index on rsIDs directly 1.

2. Get your genotypes into a normalized VCF on a known build

Everything downstream depends on knowing which reference assembly your coordinates refer to. Consumer array exports from 23andMe and AncestryDNA have historically been on GRCh37 (build 37, also called hg19), while most sequencing pipelines delivered in the last few years use GRCh38. The two differ by millions of bases of offset in places, so a GRCh37 position looked up in a GRCh38 database returns a confidently wrong answer rather than an error.

If you have an array text export, convert it and let bcftools fill in the reference allele from the FASTA:

bcftools convert --tsv2vcf genome.txt \
  -f GRCh37.primary_assembly.fa \
  -s ME --haploid2diploid \
  -Oz -o me.grch37.vcf.gz
tabix -p vcf me.grch37.vcf.gz

Then lift to GRCh38 if you want to use current annotation sets:

CrossMap vcf hg19ToHg38.over.chain.gz me.grch37.vcf.gz \
  GRCh38.primary_assembly.fa me.grch38.vcf

Expect CrossMap to drop a fraction of a percent of sites where the chain file has no unique mapping. Keep the unmapped file it writes; those positions are disproportionately in segmental duplications, which is exactly where array probes misbehave anyway.

If you have a whole-genome VCF, normalize it before anything else. Split multiallelic records into one alternate allele per row and left-align indels against the reference:

bcftools norm -m -any -f GRCh38.primary_assembly.fa \
  -Oz -o me.norm.vcf.gz me.vcf.gz
bcftools index -t me.norm.vcf.gz

Skipping norm is the single most common reason an annotation join comes back empty. A deletion represented as CTT>C in your file and TT>. in the database will never match on a naive key.

3. Look up a single rsID

For one identifier, the fastest correct route is the Ensembl REST API, because it follows dbSNP merge history and returns coordinates on both assemblies:

curl -s 'https://rest.ensembl.org/variation/human/rs671?content-type=application/json' \
  | jq '{name, synonyms, mappings: [.mappings[] | {assembly_name, seq_region_name, start, allele_string, strand}]}'

If the accession you asked about was merged, the response comes back under the current name with your query listed in synonyms. That is the behavior you want and the reason to prefer this over grepping a static file.

Once you have the coordinate, query your own VCF by region rather than by ID. Region queries use the tabix index and return immediately, while filtering on ID forces a full scan of the file:

bcftools query -r chr12:111803962 \
  -f '%CHROM\t%POS\t%ID\t%REF\t%ALT[\t%GT\t%DP\t%GQ]\n' \
  me.norm.vcf.gz

A result of chr12 111803962 rs671 G A 0/1 34 99 means one reference copy and one alternate copy, called from 34 reads with a genotype quality of 99. If the position returns nothing at all, read step 5 before concluding you are homozygous reference.

For a handful of variants at once, write the positions to a BED file and pass -R regions.bed. For thousands, move to step 4.

4. Put rsIDs on every variant in your file

Rather than looking up identifiers one at a time forever, annotate once. Download the dbSNP VCF for your assembly from NCBI. The file uses RefSeq accessions (NC_000001.11) as contig names, which will not match a VCF using chr1, so rename first:

# refseq2chr.txt: two columns, e.g. "NC_000001.11  chr1"
bcftools annotate --rename-chrs refseq2chr.txt \
  -Oz -o dbsnp.b156.chr.vcf.gz GCF_000001405.40.gz
bcftools index -t dbsnp.b156.chr.vcf.gz

bcftools annotate -a dbsnp.b156.chr.vcf.gz -c ID \
  -Oz -o me.rsid.vcf.gz me.norm.vcf.gz
bcftools index -t me.rsid.vcf.gz

bcftools annotate matches on CHROM, POS, REF, and ALT, which is why the normalization in step 2 was not optional. Expect this to take twenty to forty minutes on a whole-genome file and to leave a meaningful minority of your variants without an rsID: rare private variants and most novel indels simply are not in dbSNP.

You now have a file where ID is populated, so you can build a small local index of only the identifiers you personally carry, which is a few million rows rather than the billion-plus in dbSNP:

bcftools query -f '%ID\t%CHROM\t%POS\t%REF\t%ALT\n' me.rsid.vcf.gz \
  | grep -v '^\.' | sqlite3 -cmd '.mode tabs' me.db '.import /dev/stdin rsid'
sqlite3 me.db 'CREATE INDEX i_rsid ON rsid(ID);'

5. Read the genotype correctly

A returned genotype is not automatically a true genotype, and three failure modes account for most misreadings.

The first is strand. Array platforms report alleles relative to a design strand that may be the reverse complement of the reference. For a C/T variant this is detectable, because a reported G/A at a C/T site can only be the complement. For A/T and C/G variants it is ambiguous by inspection, and the only reliable fixes are to trust the vendor’s stated orientation or to resolve the strand using allele frequency, which fails for variants near 50 percent frequency. We treat A/T and C/G array calls as unresolved unless there is a specific reason to believe otherwise. Sequencing data does not have this problem, since reads are aligned to the reference and alleles are reported against it.

The second is absence. A position missing from your VCF can mean homozygous reference, or it can mean no coverage. If your file is a variants-only VCF, you cannot distinguish the two. Check coverage at the position directly:

samtools depth -r chr12:111803962-111803962 -a me.cram \
  --reference GRCh38.primary_assembly.fa

Fewer than about ten reads at a site in a 30x genome usually means a mapping problem rather than a real deletion. For array data, missing means the probe failed or the variant was never on the chip, and consumer chips genotype well under a million sites out of the roughly five million variants a person carries.

The third is call quality. Filter before you interpret. A 0/1 call with GQ below 20, or with an allele balance far from 0.5, is a candidate for a sequencing artifact rather than a heterozygous site. bcftools query -f '[%AD]' gives the per-allele read depths you need to check balance.

6. Add consequence and population frequency

Knowing your genotype at an rsID tells you little without knowing what the variant does to the protein and how common it is. Run VEP offline over the annotated file:

vep -i me.rsid.vcf.gz --cache --offline --assembly GRCh38 \
  --fork 8 --vcf --compress_output bgzip \
  --everything --check_existing --af_gnomade --af_gnomadg \
  --pick_allele_gene \
  -o me.vep.vcf.gz

--everything turns on SIFT, PolyPhen, canonical-transcript flags, and HGVS notation. --pick_allele_gene collapses the many transcript-level consequences to one per gene per allele, which makes the output tractable at the cost of hiding transcript-specific effects. If you are examining a specific gene, drop --pick_allele_gene and read all transcripts.

Allele frequency is the strongest single filter for triaging a variant list, and it is also where ancestry matters most. gnomAD’s frequency estimates are far better powered for European and African ancestries than for many East and South Asian populations, so a variant that looks alarmingly rare globally may be common in your own background. Population-specific resources exist to fill this in: the Human Genetic Variation Database catalogs variation from Japanese exomes and genomes specifically for this purpose 2. A careful per-gene example of how frequency, predicted effect, and functional annotation combine is the systematic annotation of roughly 1,350 common variants across the 19-member ALDH gene family from gnomAD, which sorts them by predicted functional consequence rather than treating all coding changes alike 3.

7. Check clinical interpretation, and what “VUS” means

ClinVar aggregates submitted clinical assertions about variants. Join it in the same way, and keep the review status, not only the significance:

bcftools annotate -a clinvar_GRCh38.vcf.gz \
  -c INFO/CLNSIG,INFO/CLNREVSTAT,INFO/CLNDN,INFO/CLNVI \
  -Oz -o me.clinvar.vcf.gz me.vep.vcf.gz

CLNREVSTAT is the field that separates a claim backed by an expert panel or practice guideline from one submitted by a single laboratory without assertion criteria. We weight anything below two-star review status as a hypothesis rather than a finding.

Most of what you find will be labeled “uncertain significance” (VUS): the variant has been observed, it plausibly alters a protein, and nobody has enough evidence to call it pathogenic or benign. This is the expected outcome, not an anomaly. A study of 1,000 healthy, ancestrally diverse participants sequenced and examined across 114 cancer-susceptibility genes found variants of uncertain significance in the large majority of individuals, while genuinely pathogenic findings were present in a small minority 4. VUS rates are also systematically higher in people of non-European ancestry, because the reference series used to establish benignity are thinner there.

The practical consequence is that a VUS should not change anything you do. Reclassification happens in both directions as evidence accumulates, and the population-scale frameworks developed for returning genomic results deliberately restrict reporting to variants meeting defined evidence thresholds in genes with actionable associations, precisely to avoid acting on uncertainty 5. If your file contains a pathogenic or likely pathogenic call in a gene with known clinical relevance, the correct next step is confirmation in a CLIA-certified laboratory using an orthogonal assay and a conversation with a genetic counselor or clinical geneticist. Research-grade sequencing and consumer array calls both produce false positives at rates that matter when a single variant is driving a decision, and interpretation in the context of your personal and family history is clinical work.

8. Build a stable key

If you are maintaining your own annotation store rather than re-deriving it each time, do not key on rsID. Key on the variant itself. VariantKey encodes chromosome, position, and alleles into a single reversible 64-bit integer, which gives you a sortable, joinable primary key that survives dbSNP releases and does not break when an accession is merged or withdrawn 6. Keep the rsID as an attribute for cross-referencing the literature, and refresh it whenever you update the dbSNP build.

Common problems

The annotation join returns nothing. Almost always a contig naming mismatch (1 versus chr1), a build mismatch, or unnormalized indels. Check with bcftools view -h file.vcf.gz | grep contig on both files before debugging anything else.

An rsID from a paper returns “not found”. The accession was merged or withdrawn. Query Ensembl or the dbSNP API rather than a static file, since both follow merge history.

One rsID maps to multiple rows. The variant is multiallelic, and bcftools norm -m -any split it correctly. Match on REF and ALT, not on ID alone.

The same rsID appears on two chromosomes. dbSNP occasionally clusters variants in paralogous or segmentally duplicated regions ambiguously. Treat genotypes in these regions as low confidence.

Your array says one thing and your sequence says another. Array probes fail in repetitive sequence and in the presence of a nearby second variant under the probe. When the two disagree, the sequencing call with adequate depth and balanced allele reads is the one to believe.

You are tempted to paste your file into a web lookup tool. A genotype profile is identifying. Genome-wide genotype data supports genealogical search that can locate distant relatives from a profile alone 7, which means an upload is a disclosure about your relatives as well as yourself. Run the annotation locally.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Shray Alag. Unique insights from ClinicalTrials.gov by mining protein mutations and RSids in addition to applying the Human Phenotype Ontology. PLOS ONE, 2020. https://doi.org/10.1371/journal.pone.0233438 ↩

  2. Koichiro Higasa, Noriko Miyake, Jun Yoshimura, et al. Human genetic variation database, a reference database of genetic variations in the Japanese population. Journal of Human Genetics, 2016. https://doi.org/10.1038/jhg.2016.12 ↩

  3. Che-Hong Chen, Benjamin R. Kraemer, Lucia Lee, et al. Annotation of 1350 Common Genetic Variants of the 19 ALDH Multigene Family from Global Human Genome Aggregation Database (gnomAD). Biomolecules, 2021. https://doi.org/10.3390/biom11101423 ↩

  4. Dale L. Bodian, Justine N. McCutcheon, Prachi Kothiyal, et al. Germline Variation in Cancer-Susceptibility Genes in a Healthy, Ancestrally Diverse Cohort: Implications for Individual Genome Sequencing. PLoS ONE, 2014. https://doi.org/10.1371/journal.pone.0094554 ↩

  5. Heather M McLaughlin, Ozge Ceyhan-Birsoy, Kurt D Christensen, et al. A systematic approach to the reporting of medically relevant findings from whole genome sequencing. BMC Medical Genetics, 2014. https://doi.org/10.1186/s12881-014-0134-1 ↩

  6. Nicola Asuni, Steven Wilder. VariantKey: A Reversible Numerical Representation of Human Genetic Variants. 2018. https://doi.org/10.1101/473744 ↩

  7. Yuan Wei, Ryan Lewis, Ardalan Naseri, et al. Genealogical search using whole-genome genotype profiles. Responsible Genomic Data Sharing, 2020. https://doi.org/10.1016/b978-0-12-816197-5.00004-8 ↩