How to Check a Single SNP in Your Own Genome Data

Oak
snp wgs vcf bcftools annotation dbsnp
A tripod-mounted lens clamps and brightly illuminates one bead in a long braid of faintly glowing seed-pods on a forest trunk at night.

By the end of this you will have, for any rsID you care about, a genotype call with the evidence behind it: read depth, allele depth, genotype quality, the reference allele on the plus strand of a named genome build, and the predicted consequence on a MANE Select transcript. You will also know when the answer is “no call” rather than “reference,” which is the single most common way people get this wrong. You need: your own variant calls (a .vcf.gz or .g.vcf.gz with an index), ideally the aligned reads (.cram or .bam plus index), the matching reference FASTA, and bcftools/htslib 1.17+, samtools, and optionally Ensembl VEP with an offline cache. Everything below assumes GRCh38. If you are on hg19, read step 1 before doing anything else.

1. Establish which build and which file type you have

A position means nothing without a build. chr19:44908684 is APOE on GRCh38 and points somewhere unrelated on GRCh37. Check the header first:

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

If the ##reference line is missing or unhelpful, use contig lengths. chr1 is 248,956,422 bp on GRCh38 and 249,250,621 bp on GRCh37. That one number settles it.

Then determine whether you have a variant-only VCF or a gVCF:

bcftools view -H sample.vcf.gz | head -3

A gVCF has <NON_REF> in the ALT column and END= in INFO, and it contains records for every base of the genome. A plain VCF contains only sites where the caller emitted a variant. The difference matters enormously: in a variant-only VCF, absence of a record at your SNP position means either you are homozygous reference or the caller never had enough evidence to say anything. Those two states are not the same, and the file will not distinguish them. A gVCF will, via the reference-block GQ.

Consumer array exports (the 23andMe-style .txt with rsid/chromosome/position/genotype) are a third category and are handled in step 5.

If your data is on GRCh37 and your annotation resources are GRCh38, lift over the small set of positions you care about rather than the whole callset:

CrossMap vcf hg19ToHg38.over.chain.gz query.b37.vcf GRCh38.fa query.b38.vcf

Lift the query, not the genome. Liftover of an entire WGS callset silently drops and inverts records in segmental duplications, and you will not notice.

2. Resolve the rsID to coordinates and reference alleles

Do not trust a coordinate you copied from a blog post. Resolve it against dbSNP. The full dbSNP VCF for GRCh38 is GCF_000001405.40.gz from the NCBI FTP site. It uses RefSeq accessions as chromosome names (NC_000019.10, not chr19), which trips up every first attempt:

tabix GCF_000001405.40.gz NC_000019.10:44908684-44908684

Going the other direction, from rsID to position, is the annoying case, because the VCF is indexed by coordinate and not by ID. For a handful of IDs, use the RefSNP API:

curl -s https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/429358 \
  | jq '.primary_snapshot_data.placements_with_allele[]
        | select(.seq_id | startswith("NC_000019"))
        | {seq_id, alleles: [.alleles[].allele.spdi]}'

For thousands of IDs, build a local lookup once:

bcftools query -f '%ID\t%CHROM\t%POS\t%REF\t%ALT\n' GCF_000001405.40.gz \
  | sort -k1,1 > rsid_to_pos.tsv

That is roughly an hour of I/O and about 30 GB of output, and it is worth it if you do this regularly.

Two properties of rsIDs to internalize. First, they are identifiers for a position and allele set, not for a specific allele, so an rsID alone never tells you what you carry 1. Second, they get merged and retired: when dbSNP determines two submitted variants are the same, the higher-numbered rsID is merged into the lower and stops appearing in current VCFs. The RefSNP API returns the merge target, so an ID that “does not exist” in a current release is usually alive under a different number.

3. Query your own VCF

With the coordinate in hand:

bcftools view -r chr19:44908684 -Ov sample.vcf.gz \
  | bcftools norm -m -any -f GRCh38.fa -Ov \
  | bcftools query -f '%CHROM\t%POS\t%ID\t%REF\t%ALT[\t%GT\t%AD\t%DP\t%GQ]\n'

The norm step is not optional. Callers emit multiallelic records (ALT=A,T) and left-align indels inconsistently, so a naive string comparison of your ALT against dbSNP’s ALT fails on perfectly good data. -m -any splits multiallelics into one record per ALT and -f re-normalizes against the reference.

Note the region syntax: -r chr19:44908684 is 1-based and inclusive. If you pass a BED file with -R, that file is 0-based half-open. Mixing them shifts every query by one base, which produces plausible-looking wrong answers.

On a gVCF, the record covering your position may start thousands of bases earlier with END= extending over it. htslib’s index accounts for INFO/END, so -r normally returns the enclosing block, but verify this once on a position you know sits mid-block before you rely on it across a batch.

4. Go back to the reads when the VCF says nothing

If step 3 returns no record, you have not yet learned anything. Pile up the reads:

samtools mpileup -f GRCh38.fa -r chr19:44908684-44908684 \
  -q 20 -Q 20 sample.cram

Output columns are chromosome, position, reference base, depth, the read bases, and base qualities. Dots and commas are reference-matching reads on the forward and reverse strands. Letters are mismatches. If you see 32 reads and all of them are dots, you are homozygous reference at adequate depth. If you see 3 reads, you have no information and the caller was right to stay quiet.

To get a formal call including a homozygous-reference genotype:

bcftools mpileup -f GRCh38.fa -r chr19:44908684 -a AD,DP,SP \
    -q 20 -Q 20 sample.cram \
  | bcftools call -m -Ov \
  | bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%GT\t%AD\t%DP]\n'

Two things to watch. samtools mpileup caps depth (-d, default 8000) and applies base alignment quality recalibration by default, which can suppress mismatches near indels. Disable it with -B if a site looks suspiciously clean next to a known indel. And a 30x genome has a Poisson-distributed depth: roughly 1-2% of positions land below 10x even in a good library, and GC-extreme or repetitive regions do much worse. Low depth at one specific SNP is normal, not evidence of a problem with the sequencing.

We would always look at the reads for any single SNP we plan to act on. It takes thirty seconds and it catches the strand and no-call failure modes that a table of genotypes hides.

5. Get the alleles and the strand right

dbSNP reports alleles on the plus strand of the reference assembly. Your VCF is also plus-strand. So far so good. The trouble comes from three directions.

Genotyping arrays historically used TOP/BOT or A/B allele encodings, and consumer exports sometimes report on the gene’s strand rather than the assembly’s. For a SNP on a minus-strand gene, an array export saying CT and a VCF saying G/A can describe the same genotype. For unambiguous SNPs (A/G, C/T) you can resolve this by complementing and checking which orientation contains the reference base. For strand-ambiguous SNPs (A/T and C/G) you cannot resolve it from the genotype alone, and allele-frequency-based inference only works when the minor allele frequency is far from 0.5. Treat A/T and C/G calls from array data as unresolved unless you have the manifest.

Second, the dbSNP reference allele is a property of the current assembly and occasionally flips between builds when the reference itself is corrected. If your annotation source and your VCF disagree about REF at a position, stop and check the build before assuming a swap.

Third, some rsIDs cover more than two alleles. Human SNPs are overwhelmingly biallelic, but triallelic sites exist and a REF=C ALT=T record tells you nothing about a third allele the caller did not observe 2.

6. Annotate the consequence

Now attach biology to the call. With an offline VEP cache:

vep -i query.vcf --cache --offline --assembly GRCh38 \
    --dir_cache ~/.vep --fasta GRCh38.fa \
    --tab -o out.tsv --fork 4 \
    --mane --canonical --symbol --hgvsc --hgvsp \
    --af_gnomadg --pick_allele_gene

--pick_allele_gene gives you one consequence per allele per gene, which is what you want in a table. Drop it and you get one row per transcript, which is what you want when the transcript matters. --mane flags the MANE Select transcript, which is the one to quote when you write down an HGVS notation for a clinician.

This is also where the “types of SNP” question gets answered, because there are three different classification axes and people conflate them.

By nucleotide change, a SNP is a transition (purine to purine, A↔G, or pyrimidine to pyrimidine, C↔T) or a transversion (purine to pyrimidine or vice versa). Transitions outnumber transversions roughly two to one genome-wide, largely because methylated CpG cytosines deaminate to thymine. The transition/transversion ratio of your callset is a quality metric: expect around 2.0-2.1 genome-wide and around 3.0 in exons, and suspect contamination or a filtering problem if you see 1.5.

By genomic location, SNPs are coding, intronic, untranslated-region, regulatory, or intergenic. The great majority are non-coding, since coding sequence is about 1-2% of the genome. Across the genome, common variation sits on the order of one SNP per kilobase between two chromosomes 1, and the first genome-scale maps were built by shotgun sequencing a reduced representation of pooled samples to find exactly these positions 3.

By effect on protein, coding SNPs are synonymous (no amino acid change), missense (one amino acid substituted), or nonsense (a stop codon introduced). Splice-site SNPs sit slightly outside this taxonomy and are frequently the more consequential ones. These categories describe what a variant does to a sequence, not what it does to a person.

7. Batch a list of rsIDs into a genotype table

# regions.txt: 1-based CHROM<TAB>POS, produced from the rsid lookup in step 2
bcftools view -R regions.txt -Ou sample.vcf.gz \
  | bcftools norm -m -any -f GRCh38.fa -Ou \
  | bcftools query -f '%CHROM\t%POS\t%ID\t%REF\t%ALT[\t%GT\t%AD\t%DP\t%GQ]\n' \
  > genotypes.tsv

Then do the important part: left-join this against your full rsID list and audit the misses. Every rsID with no row is a site you must resolve in step 4 before you call it reference. In our experience with a 30x WGS callset and a list of a few thousand rsIDs, 1-3% come back without a record, and a handful of those turn out to be genuine no-calls rather than reference.

A reasonable GQ floor for a single site you care about is 30 with DP of at least 15 and an allele balance between 0.3 and 0.7 for a heterozygous call. Outside those bounds, look at the pileup.

8. Put the call in context before you conclude anything

You now have a genotype. The gap between “I carry this allele” and “this means something for me” is where most SNP checking goes off the rails. Three things to do before drawing any conclusion.

Check the population frequency. gnomAD genome frequencies come through VEP with --af_gnomadg. An allele at 15% frequency in your ancestry group is not rare and is very unlikely to be individually informative.

Check the evidence for the association. Most reported genotype-phenotype links come from association studies measuring small average effects across large cohorts, and the effect sizes do not transfer cleanly to a single person. Whole-genome sequencing gives complete coverage of a person’s variation, but interpretation remains the bottleneck rather than data generation 4. A large fraction of variants found in any individual genome have no confident functional interpretation 5.

Check whether this is a question for a clinician. Variants in genes with established clinical significance, including APOE, BRCA1/2, LDLR, HFE, and the pharmacogenes, should be confirmed in a clinical laboratory and discussed with a genetic counselor or physician. Research-grade sequencing is not a clinical result, and the decision-making around a clinically significant finding depends on family history and personal context that no file contains 6. We do not interpret these calls for you and neither should a chatbot. Take the coordinate, the HGVS notation, and the read evidence to someone qualified to act on it.

Common problems

The rsID returns nothing from dbSNP. It was merged into a lower-numbered ID. Query the RefSNP API and follow the merge.

Chromosome names do not match. dbSNP VCFs use NC_000001.11, Ensembl uses 1, UCSC and most pipelines use chr1. Fix with bcftools annotate --rename-chrs chr_map.txt, where the map is two columns, old and new.

The position is right but the genotype is empty (./.). The caller emitted a record and then failed to genotype it. Go to the pileup.

Everything is shifted by one base. You passed a BED file to -r or a 1-based region to -R. They use different conventions.

A het call with 40 reads supporting ALT and 4 supporting REF. That is not a clean heterozygote. Look for a nearby indel, a paralogous region, or mapping quality zero reads. samtools view -c -q 1 versus -c at the locus tells you the fraction of multi-mapping reads.

An indel adjacent to your SNP changes the REF/ALT representation. Always bcftools norm -f before comparing to any external source.

The array export and the WGS call disagree. Check strand first, especially for A/T and C/G SNPs, then check whether the array probe sits over another variant you carry, which is a common cause of array no-calls and miscalls.

MT and chrY. Mitochondrial calls are haploid and often reported with heteroplasmy fractions rather than genotypes, and chrY has large regions where mapping is unreliable. Different rules apply to both. Do not reuse your autosomal thresholds there.

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. Anthony J. Brookes. The essence of SNPs. Gene, 1999. https://doi.org/10.1016/s0378-1119(99)00219-x ↩ ↩2

  2. Benjamin A. Salisbury, Manish Pungliya, Julie Y. Choi, et al. SNP and haplotype variation in the human genome. Mutation Research - Fundamental and Molecular Mechanisms of Mutagenesis, 2003. https://doi.org/10.1016/s0027-5107(03)00014-9 ↩

  3. David Altshuler, Victor J. Pollara, Chris R. Cowles, et al. An SNP map of the human genome generated by reduced representation shotgun sequencing. Nature, 2000. https://doi.org/10.1038/35035083 ↩

  4. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩

  5. Michael Snyder, Jiang Du, Mark Gerstein. Personal genome sequencing: current approaches and challenges. Genes & Development, 2010. https://doi.org/10.1101/gad.1864110 ↩

  6. A.J. Marian. Sequencing Your Genome: What Does it Mean?. Methodist DeBakey Cardiovascular Journal, 2014. https://doi.org/10.14797/mdcj-10-1-3 ↩