dbSNP — reference variant identifiers (rsID) and the COMMON flag

dbSNP157/GRCh38 is NCBI's catalog of short genetic variation — every submitted, deduplicated SNP/indel gets a stable rs# accession. iris reads it as a single bgzipped, tabix-indexed VCF, GCF_000001405.40.gz, at /mnt/cephfs/hot_nvme/dbsnp/dbSNP157_hg38/ (29.55 GB / 27.52 GiB, .tbi index ~3.0 MB). Unlike every other tabix source iris reads, this file is one file for the whole genome, not split per chromosome, and its #CHROM values are NCBI RefSeq accessions (NC_000001.11, …) rather than chr1/1.

This page documents the file itself — contigs, INFO schema, a couple of real parsing gotchas discovered while building the adapter — and shows how iris reads it via DbSnpAnnotator (src/iris/adapters/driven/annotators/dbsnp.py). Every number and example below ran against the real file.

##fileformat=VCFv4.2
##fileDate=20241205
##source=dbSNP
##dbSNP_BUILD_ID=157
##reference=GRCh38.p14

Scale

698 contigs total, but only 25 are usable by iris (see Contigs below). Full-genome record counts weren't exhaustively tallied (that's a multi-minute scan even with an index), but real per-contig counts, fetched end-to-end:

Contig	Records	Fetch time
NC_000001.11 (chr1, 248.9 Mb)	93,124,941	59.8s
NC_000021.9 (chr21, 46.7 Mb)	15,323,112	10.1s
NC_000024.10 (chrY, 57.2 Mb)	3,047,929	1.8s
NC_012920.1 (chrMT, 16.6 kb)	9,229	<0.1s

Density varies a lot by contig — chr1 and chr21 are both ~0.33–0.37 records/bp, chrY drops to ~0.053 records/bp (poorly callable, repetitive), and MT is ~0.56 records/bp (mitochondrial hypervariability). Don't extrapolate a genome-wide total from one contig's density.

Contigs: only 25/698 are reachable from iris

NC_ (chromosomes)        25   NC_000001.11 … NC_000022.11, NC_000023.11 (X),
                               NC_000024.10 (Y), NC_012920.1 (MT)
NW_ (unplaced scaffolds) 323
NT_ (alt loci / patches) 350

DbSnpAnnotator resolves contigs via resolve_human_chromosome(), which only knows the 25 primary NC_* accessions (see HUMAN_CHROMOSOMES). The 673 NW_/NT_ scaffold contigs are silently skipped — any dbSNP record that exists only on an alt haplotype or unplaced scaffold is unreachable through this annotator. This mirrors iris' general policy of working against the primary assembly; it's not a bug, but worth knowing before assuming "no match" means "not in dbSNP" for a coordinate on a patch contig.

INFO schema

31 INFO fields, grouped by what DbSnpAnnotator actually consumes vs. not:

Consumed by DbSnpAnnotator

Field	Meaning
ID (column, not INFO)	The rs# accession itself — this is dbsnp_id
COMMON	Flag: ≥1% minor allele frequency in ≥1 1000Genomes population, with ≥2 founders contributing to that frequency (dbSNP's own definition, from the VCF header)

Present but not consumed

Field	Meaning
RS	The rs# as an integer (redundant with ID)
GENEINFO	symbol:gene_id pairs, \|-delimited for multiple genes, e.g. HES4:57801
PSEUDOGENEINFO	Same format, for overlapping pseudogenes
dbSNPBuildID	First dbSNP build this rs# appeared in — a rough recency/provenance signal
SAO	Allele origin: 0 unspecified, 1 germline, 2 somatic, 3 both
SSR	Suspect reason bitmask: paralog, byEST, oldAlign, Para_EST, 1kg_failed, other
VC	Variant class — observed values: SNV, INDEL, INS, DEL, MNV
PM, NSF, NSM, NSN, SYN, U3, U5, ASS, DSS, INT, R3, R5, GNO, PUB	Boolean flags for publication status and coarse functional class (non-synonymous frameshift/missense/nonsense, synonymous, UTR, splice site, intron, gene-flanking, genotypes-available) — all superseded by GencodeVariantAnnotator/VepFlightAnnotator's actual transcript-consequence calls
FREQ	Per-study allele frequencies — see below, not currently extracted into AlleleFrequency
CLNHGVS, CLNVI, CLNORIGIN, CLNSIG, CLNDISDB, CLNDN, CLNREVSTAT, CLNACC	Legacy embedded ClinVar mirror — see below

The FREQ field: a pysam parsing gotcha

FREQ packs a two-level structure into one string: studies separated by |, each as STUDY:freq_ref,freq_alt1,freq_alt2,…:

FREQ=KOREAN:0.9891,0.0109,.|SGDP_PRJ:0,1,.|dbGaP_PopFreq:1,.,0

Its header declares Number=. (variable-length), which tells pysam "split on commas" — pysam has no idea the string also uses | and : internally, so it mangles the pipe-delimited studies into meaningless comma fragments:

>>> record.info["FREQ"]
('KOREAN:0.9891', '0.0109', '.|SGDP_PRJ:0', '1', '.|dbGaP_PopFreq:1', '.', '0')

The correct reconstruction rejoins with , before splitting on |:

>>> ",".join(record.info["FREQ"]).split("|")
['KOREAN:0.9891,0.0109,.', 'SGDP_PRJ:0,1,.', 'dbGaP_PopFreq:1,.,0']

Sampling ~140K records across 4 chromosomes found 32 distinct study cohorts (not an exhaustive count — dbSNP157 draws from more across the whole genome):

1000Genomes, 1000Genomes_30X, ALSPAC, Chileans, Daghestan, Estonian, ExAC,
FINRISK, GENOGRAPHIC, GENOME_DK, GnomAD_exomes, GnomAD_genomes, GoESP, GoNL,
HGDP_Stanford, HapMap, KOREAN, Korea1K, Korea4K, MGP, NorthernSweden,
PAGE_STUDY, PRJEB36033, PRJEB37584, Qatari, SGDP_PRJ, Siberian, TOMMO,
TOPMED, TWINSUK, Vietnamese, dbGaP_PopFreq

DbSnpAnnotator doesn't extract any of these into AlleleFrequency — they're per-study cohorts of wildly varying size and ancestry composition, not the standardized population buckets AlleleFrequency.population expects elsewhere (compare MCPS, which has exactly that: consistent labeled populations). If a future need calls for gnomAD- or TOPMed-specific frequencies from dbSNP's FREQ, use the reconstruction above — GnomAD_exomes, GnomAD_genomes, and TOPMED are all present as study labels.

The legacy CLN* block: stale, not empty

DbSnpAnnotator's docstring used to claim these fields were "no longer populated by recent builds" — that turned out to be wrong. In a 286K-record sample, CLNSIG was populated on 4,337 records (~1.5%):

NC_000001.11  1013466  rs3841266     CLNSIG=(.,2|2)  CLNDN=(.,Mendelian_susceptibility_to_mycobacterial_diseases_due_to_complete_ISG15_deficiency|not_specified)
NC_000001.11  1013549  rs553176213   CLNSIG=(.,3,.,.,.)  CLNDN=(.,ISG15-related_disorder,.,.,.)

The catch: CLNSIG here uses dbSNP's own numeric code table (from this file's own header — 2 = Benign, 3 = Likely benign, 5 = Pathogenic, …), not ClinVar's current string vocabulary ("Pathogenic", "Likely_benign", …) that ClinvarAnnotator parses from the authoritative ClinVar TSV. It's a snapshot embedded at some earlier ClinVar release and not kept in sync with dbSNP's own build cadence — a stale, differently-encoded duplicate, which is why DbSnpAnnotator deliberately doesn't parse it. Also positionally aligned to ALT (. placeholder per allele with no data, |-joined when an allele has multiple ClinVar submissions) rather than a single scalar — another reason to prefer ClinvarAnnotator's purpose-built parsing over reimplementing this.

Access patterns from iris

Registering the source

# ~/.iris/config.toml
[sources.dbsnp]
host = "local"
path = "/mnt/cephfs/hot_nvme/dbsnp/dbSNP157_hg38/GCF_000001405.40.gz"

from iris.adapters.driven.annotators.dbsnp import DbSnpAnnotator

annotator = DbSnpAnnotator.from_registry()          # resolves [sources.dbsnp]
annotator = DbSnpAnnotator.from_registry(workers=8)  # sets the default workers

annotate() — point-lookup enrichment

Matches by Variant.as_key(), the same (chrom, pos, ref, alt) contract every VariantAnnotator uses. Live run against two real, COMMON-flagged records on different chromosomes:

from iris.domain.genomics.entities import Variant, VariantAnnotation
from iris.domain.genomics.objects import Allele, GenomicPosition, VariantNomenclature, VariantType

def va(chrom, pos, ref, alt):
    key = (chrom, pos, ref, alt)
    return VariantAnnotation(variant=Variant(
        variant_type=VariantType.SNV,
        nomenclature=VariantNomenclature.from_key(key),
        alleles=Allele(reference=ref, alternate=(alt,)),
        position=GenomicPosition(chromosome=chrom, start=pos, end=pos),
    ))

items = [va("1", 10_001, "T", "A"), va("2", 10_026, "A", "G")]
result = annotator.annotate(items)

1:10001:T:A → dbsnp_id=rs1570391677  dbsnp_common score=1.0 (version=157)
2:10026:A:G → dbsnp_id=rs1366167113  dbsnp_common score=1.0 (version=157)

Both are real NC_000001.11/NC_000002.12 records, resolved transparently from the canonical "1"/"2" chromosome names the caller used — the RefSeq-accession translation described above is invisible to callers.

workers — chromosome-partitioned parallel fetch

300 real targets spread across 5 chromosomes, same file, same process:

workers=1: 0.37s, 300/300 rsIDs assigned
workers=8: 0.19s, 300/300 rsIDs assigned   (n_workers = min(8, 5 chroms) = 5 effective)

Both produce identical results — workers only changes how many pysam.VariantFile handles run concurrently (one per chromosome partition, mirroring TabixReader._query's pattern), not what gets matched. The gap above is small because these positions were already warm in page cache from prior exploration; the real-world win is much larger under cold, network-filesystem-latency-bound access. A production staged-annotation run (projects/vcfannotation/annotate_staged.py, chained with mcps/clinvar/ dbnsfp) against 29,269 rare variants took ~7–8 minutes with sequential, single-handle fetching, before DbSnpAnnotator had chromosome-partitioned workers support at all. After (a) restructuring the pipeline so dbSNP — which doesn't feed any CANDIDATES_FILTER rule — only runs on the small post-filter candidate set, and (b) adding this chromosome-partitioned workers support, the same sample's full three-stage pipeline (mcps → clinvar/dbnsfp → dbsnp) completed in ~15 seconds, byte-identical output.

Design notes & caveats

• One file for the whole genome, not per-chromosome. Every other tabix source iris reads (dbNSFP, CADD's TSV, MCPS, RGC-ME) is split into per-chromosome files opened on demand. dbSNP is a single 27.52 GiB VCF — DbSnpAnnotator opens it once per sequential call, or once per chromosome-partition worker thread when workers>1.
• Only the 25 primary NC_* contigs resolve. Records on the 673 NW_/NT_ scaffold/patch contigs are unreachable — see Contigs.
• FREQ needs manual reconstruction if you ever parse it directly with pysam — see the gotcha above. Not extracted into AlleleFrequency today.
• CLNSIG/CLNDN/etc. are populated but stale and numerically encoded differently from ClinVar's own vocabulary — use ClinvarAnnotator for clinical significance, not these fields. See above.
• COMMON is a derived boolean, not a frequency value. DbSnpAnnotator represents it as a VariantScore (name="dbsnp_common", value=1.0, only emitted when the flag is present — never 0.0) and merges dbSNP's exact >=1%/1000Genomes/founder-count criteria into VariantAnnotation.metadata via Metadata.merge() (non-destructive — won't clobber metadata from an earlier source in the pipeline), so the threshold definition travels with the data instead of living only in code comments.
• dbsnp_id is always overwritten on a hit, not just filled when missing — dbSNP157 is treated as the authoritative current rs# for that position, even if the input already carried a (possibly stale/merged) rsID from its source VCF's own ID column.

Reference

iris.adapters.driven.annotators.dbsnp

dbSNP rsID and COMMON-flag annotator via VCF random-access.

dbSNP ships as a single bgzipped, tabix-indexed VCF (e.g. GCF_000001405.40.gz for dbSNP157/GRCh38) whose contigs are NCBI RefSeq accessions (NC_000001.11, …) rather than UCSC/plain chromosome names. VcfReader._fetch resolves contigs via resolve_chromosome, which only toggles the chr prefix and can't bridge RefSeq accessions, so this annotator does its own contig resolution via resolve_human_chromosome before fetching.

Only the ID column (rs#) and the COMMON INFO flag are consumed. dbSNP also embeds a legacy CLN* block (CLNSIG etc.) using dbSNP's own numeric significance codes from an older ClinVar snapshot — still populated on a minority of records, but a stale, differently-encoded duplicate of what ClinvarAnnotator already provides from the authoritative ClinVar source, so it's deliberately not parsed here.

DbSnpAnnotator

Bases: VariantAnnotator

Annotates variants with dbSNP rsIDs and the COMMON flag.

Expects the standard dbSNP VCF shipped by NCBI (bgzipped, tabix-indexed), whose #CHROM values are RefSeq accessions such as NC_000001.11. version is stamped on the emitted COMMON score and should match the dbSNP build (e.g. "157").

COMMON is a boolean derived by dbSNP itself (>=1% MAF in >=1 1000Genomes population, >=2 founders — see _COMMON_CRITERIA), not a frequency value, so it's represented two ways: a VariantScore (name "dbsnp_common", value=1.0) for filtering via score_for(), and the criteria text merged into VariantAnnotation.metadata so the threshold travels with the data instead of living only in code comments.

workers class-attribute instance-attribute

workers = 1

Default parallel per-chromosome workers, used when annotate() omits an override.

from_registry classmethod

from_registry(source='dbsnp', registry=None, *, workers=1)

Build a DbSnpAnnotator resolving source from the iris DataRegistry.

Source code in src/iris/adapters/driven/annotators/dbsnp.py

@classmethod
def from_registry(
    cls, source: str = "dbsnp", registry=None, *, workers: int = 1
) -> "DbSnpAnnotator":
    """Build a DbSnpAnnotator resolving *source* from the iris DataRegistry."""
    from iris.config import get_registry

    r = registry or get_registry()
    return cls(path=r.resolve_path(source), workers=workers)

annotate

annotate(items, *, workers=None)

Annotate variants with dbSNP rsIDs and the COMMON flag.

PARAMETER	DESCRIPTION
items	Variants to annotate. TYPE: Sequence[VariantAnnotation]
workers	Parallel per-chromosome workers for this call; defaults to self.workers (set at construction, e.g. via from_registry(workers=...)) when omitted. TYPE: int | None DEFAULT: None

Source code in src/iris/adapters/driven/annotators/dbsnp.py

def annotate(
    self, items: Sequence[VariantAnnotation], *, workers: int | None = None
) -> list[VariantAnnotation]:
    """Annotate variants with dbSNP rsIDs and the COMMON flag.

    Args:
        items: Variants to annotate.
        workers: Parallel per-chromosome workers for this call; defaults
            to ``self.workers`` (set at construction, e.g. via
            ``from_registry(workers=...)``) when omitted.
    """
    targets: set[_Key] = set()
    for va in items:
        key = va.variant.as_key()
        if key is not None:
            targets.add(key)
    if not targets:
        return list(items)

    n_workers = self.workers if workers is None else workers
    hit_map: dict[_Key, tuple[str, bool]] = {}
    for record, key in self._fetch(Path(self.path), targets, workers=n_workers):
        rsid = record.id or ""
        if rsid:
            hit_map[key] = (rsid, "COMMON" in record.info)

    return [self._enrich(va, hit_map) for va in items]