Read a variant call set at a locus

The route: open the call set → jump to the region → read the records → check the census against what you expected → apply a change to the sequence, or hand the region to the alignment.

Every number here came from running the reader on the real files in data/samples/VCF/, which the Variant calls (VCF) and BAM & SAM alignments starters open.

Step 1 — Open the call setCopy Anchor Link

Tool: Variant Calls. Where: open a .vcf or .vcf.gz from the lobby or File → Import.

Two shapes of file behave differently, and the difference is the whole point of this workflow:

The indexed path needs both files. A .vcf.gz or .bcf without its index has no random access — the lobby's BAM starter requires the same pairing for the same reason.

BCF is not a lesser citizen here. A binary call set answers a locus query exactly as a bgzipped text one does, and the numbers are identical because both end at the same reader: the same 109 records read 19 for contigA:3,000-5,000 whether you open the .bcf or the .vcf.gz. The difference you will notice is that the .bcf is smaller and needs no text parsing at all.

One thing that is specific to BCF and worth knowing: a BCF's index carries no contig names. The contigs come from the file's own header dictionary, which is why the region box offers exactly what that header declares — including a contig no record in the file actually uses.

When the file is too large to read wholeCopy Anchor Link

There is a ceiling on reading a call set in one pass — 128 MB of uncompressed text, which a whole-genome multi-sample set passes easily. What happens when you cross it depends entirely on whether the file is indexed, and it is the clearest reason to index before you need to:

• Indexed. The file opens. Nothing is parsed: you get its real header, and the panel says so in as many words, naming the size and the first contig to try. The region box stops being an optional narrowing and becomes the way in, and every query after that reads only the blocks it needs. The status line reads not read whole beside the sample count the header declares, so a census can never be mistaken for the file's.
• Not indexed. It is refused, and the message names the fix that matches what you are holding: tabix -p vcf if the file is already bgzipped, bgzip first if it is plain text. Subsetting with bcftools view -r remains the other route.

A 140 MB call set open at its header. Nothing has been read, and nothing pretends otherwise. The contig filter reads All contigs (0/1) — none of the file's contigs read, one known to the index — rather than a bare zero that would look like a call set with no contigs in it. The position filter, which narrows the records already in hand, is greyed out and says nothing read yet, because there is nothing yet to narrow.

Reading the size costs nothing, because it is not read by decompressing the file. Every block in a bgzipped stream records its own uncompressed length in its header, so the total is a walk over those headers — under a millisecond on a file that takes a quarter of a second to decompress.

Importing such a file costs nothing either, and this is newer than the reading path. Working out what a file is used to mean decompressing it, which on a 140 MB call set meant fifteen seconds and several gigabytes of memory before anything was shown — for a file the reader was then going to open at its header without decompressing at all. The header alone is enough to say what the file is, and it sits in the first few kilobytes of the stream, so that is all that is read: the size, the sample count, the contigs the header declares and the VCF version, in a few milliseconds. No record count is offered, because counting records means reading them, and a zero would be a made-up number about a file full of calls.

Once you name a locus, the file behaves like any other call set — a real table, a real census, the same hand-offs:

The same file after one region query: 19 of 19 records over 1,094 samples, Ts/Tv 0.78, 89.5% PASS, read straight out of the compressed blocks. The contig filter has moved to All contigs (1/1) — one contig read of the one the index holds — and the position filter is live again now that there are records for it to narrow.

Check: the header banner shows the fileformat and the sample count. 1000genomes.vcf reports VCFv4.1, 1,233 samples, 2 records; simple.vcf reports VCFv4.3, 3 samples, 5 records. A sample count of zero on a file you expected genotypes from means the FORMAT columns are not there.

Step 2 — Jump to the locusCopy Anchor Link

Type a region in the usual contig:start-end form. It is parsed and echoed back, so a readout can never claim a region it did not read:

20:1230000-1240000   →   20:1,230,000-1,240,000

A locus read out of a binary call set. The readout says region contigA:3,000-5,000 beside the counts — so a Ts/Tv computed over one window can never be misread as the whole file's.

On data/samples/VCF/indexed/volvox.test.vcf.gz:

query contigA:1-5,000     →  19 records
query contigA:1,000-2,000 →   0 records
query ctgA:1-1,000        →  REFUSED
     "ctgA" is not in this file's index. It has contigA.

The refusal, and what makes it useful: it names the contigs the file does have, and the table below still holds all 109 records. You are never left staring at an empty table wondering whether the region is empty or your query was wrong.

That refusal is the design, and it is worth understanding. This file's header declares ##contig=<ID=20> and ##contig=<ID=ctgA>, while every record in it is on contigA. Declaring contigs in a VCF header is optional and unenforced, so real headers are routinely stale, partial, or copied from another call set — and an early version of this reader shipped a check that compared the header's declarations against the index and warned when they disagreed. It fired on the very first real file it met. It was deleted, and a test pins the finding so nobody re-adds it.

What protects you instead is that nothing can be silently wrong: a query for a contig the index does not have is refused, and the refusal lists the contigs the index actually has — which is the one piece of information you need to spot a mismatched index, without anyone guessing on your behalf.

An empty region is likewise an empty region, not an error. contigA:1000-2000 genuinely contains no calls in this file.

Check: the contig list offered by the picker comes from the index, not from the header. Those are different lists and the index is the one that can answer a query.

Step 3 — Read the recordsCopy Anchor Link

Each record carries what the spec says it carries, and the reader keeps the distinctions the spec makes:

contigA:3000  G>T          SNP    qual 100  filters ["PASS"]     end 3000
contigA:3105  G>C          SNP    qual 100  filters ["q10"]      end 3105
contigA:3152  A><DEL:ME>   other  qual 100  filters ["q10","s50"] end 3152
contigA:3171  G><INS>      other  qual 100  filters ["PASS"]     end 3171

Three things there are easy to get wrong and are not:

• end is not pos. A symbolic allele with INFO/END spans its real interval; a plain SNV ends where it starts. A <DEL> drawn as one base would be a lie on every canvas downstream.
• An empty filters list means PASS was written. A missing one means the column was . — "not filtered" and "filtered and passed" are different statements, and only one of them is a claim by the caller.
• Variant typing follows htslib, per ALT allele plus an OR'd mask for the record, so the labels match what bcftools would tell you. A higher-level library with its own var_type definition would have made every count on this page disagree with the tool you check against.

Step 4 — Read the census, and let it catch your mistakesCopy Anchor Link

The statistics are the same sections bcftools stats reports, computed over what you loaded:

The census over the window from Step 2, not over the file: 19 records, 16 SNVs, 3 symbolic, Ts 7 against Tv 9. Every section here is one bcftools stats reports, and each is checked against real bcftools output file by file.

Plus, per file: the 12 directed substitution types, indel lengths by size (simple.vcf: one −2 and one +1), per-contig record counts, every FILTER value seen with its tally (simple.vcf: PASS ×4, q10 ×1), QUAL min/median/max (3 / 47 / 67), and per-sample genotype counts.

Ts/Tv is the number that catches a bad call set. A real human SNP set sits near 2.0 genome-wide and higher in exons; a value near 1.0 is what random noise looks like. On these tiny sample files it is 1.00, and that is honest — two records cannot support a transition/transversion ratio. The statistic is reported twice, over all ALTs and over first-ALTs only, because multi-allelic sites make those differ.

The per-sample table gives every sample its own row of the same counts:

Per-sample counts for 1,094 samples over the same window. Hom-ref, het and hom-alt count SNV sites only, which is what bcftools reports and a caveat most reimplementations drop silently.

An empty per-sample table always has a stated reason rather than looking like a file with no samples.

One refusal you may meet here, and it is a real one. If a genotype names an allele the record never declares — a 0/3 where the record has one ALT — that call is left out of these counts and named for you, with its sample and position. It is not a display quirk: bcftools refuses to produce a per-sample section for such a file at all, so those counts would have nothing to check them against. The site counts on the Summary tab are unaffected and stay comparable.

Step 5 — Do something with the regionCopy Anchor Link

Three hand-offs, all from the same loaded records:

• Apply a variant to the sequence. With the reference open in the editor, applying a call edits the sequence and reports what it did — which is how you get from "there is a C→T at 3105" to a molecule you can then digest, translate or design primers against.
• Filter and export. Narrow by type, filter status, quality or region, then write back a VCF, a BED of the positions, or a TSV for a spreadsheet. The VCF writer round-trips the header.
• Look at the reads underneath. If you have the alignment, Alignment (BAM) opens the same locus as a pileup — which is the only way to tell a real heterozygote from a strand-biased artefact.

Check before you apply anything: the record's end, and whether the ALT is symbolic. A <DEL> or <INS> names a change without spelling it out, and applying one to a sequence is not the same operation as substituting a base.

What you did, and where each piece livesCopy Anchor Link

If it does not workCopy Anchor Link

RelatedCopy Anchor Link

The tools, in depth: Variant Calls · Alignment (BAM) · Table for the exported TSV

What follows: applying a call and then running the changed sequence through Properties or Protein Analysis — the point of the workbench being one editor is that a variant becomes a molecule you can keep working on.

The canvases this draws on: the variant table · the pileup · the sequence map.