Variant calls

The reading beside a variant call set: what the records add up to, how each figure was derived, and what about this file the numbers cannot be taken to mean.

The table itself — the records, the filters, the per-sample genotype columns — is described under the variant table. This page is about what the panel says the call set is.

What this computesCopy Anchor Link

Nothing here changes the file, and nothing here is a second counting: every figure is computed once, from the records actually in hand, and the table, the status line and this panel all read that one computation. The counts match what bcftools stats reports for the same file, which is the standard they are checked against.

The summary, section by sectionCopy Anchor Link

• Records by class — how many sites carry an SNV, an indel, an MNV, and how many carry something else (a symbolic allele such as a deletion called by length, or a breakend). A site counts in every class its alleles carry, so a site where one alternate allele is a substitution and another is an insertion is counted once under each. That is the convention bcftools uses, and it is why the class counts can sum to more than the number of records.
• Reference only — rows asserting that a position, or a block of positions, was examined and no variant was found. gVCFs and some callers emit them in quantity. They are records, and they are not variants. The name is deliberately not "no alternate allele": most such rows do state an alternate, they state one that means not a variant — a bare ., or one of the placeholders X, <X>, <*> and <NON_REF> that a caller writes when it wants a slot for "anything else that might be here". All five are read as reference-only, as is an alternate spelled identically to the reference, of any length. That is the rule htslib applies, so a gVCF's reference blocks are not mistaken for symbolic variant calls.
• Multiallelic sites — sites carrying more than one alternate allele. Reported separately for all classes and for SNVs alone, because the SNV figure uses an exact type match: a site mixing a substitution with an insertion is multiallelic but is not a multiallelic SNV site.
• Transitions and transversions — a transition is a purine-to-purine or pyrimidine-to-pyrimidine change (A↔G, C↔T); a transversion is any of the other eight. Both are counted twice over: once across every alternate allele, and once counting only each site's first alternate. The two arms differ exactly when the call set has multiallelic sites, and both are shown with their own counts so either can be checked rather than taken on the ratio alone.
• Ts/Tv — the ratio of the two, reported as undefined rather than as zero or infinity whenever it cannot be taken: a ratio with a zero denominator is a measurement nobody made. Two different situations produce that, and the panel says which. Transitions but no transversions is one. No SNV allele at all is the other, and it is ordinary in a locus read or in a gVCF, where most records state a reference block rather than a variant.
• Ts/Tv confidence interval — an exact (Clopper–Pearson) 95% interval on the fraction of SNV alleles that are transitions, carried onto the ratio scale, with the number of alleles it rests on.
• The twelve directed substitutions — A>C through T>G, each with its count.
• Indel lengths — the alternate allele's length minus the reference's, against the number of sites at that length. Negative is a deletion.
• Per contig — records, substitutions and indels for each contig, named as the file names them.
• FILTER — every filter value seen, with its count. PASS and an unfiltered . are counted apart, because "passed the filters" and "was never filtered" are different claims about a record.
• QUAL — the minimum, maximum and median of the scores present, and how many records state none.
• Per sample — for each sample: homozygous-reference, homozygous-alternate and heterozygous genotype counts, its own transitions, transversions and indels, its average depth, its singletons, its haploid calls, and its missing calls. These match bcftools stats per sample, which is worth a sentence of its own: bcftools' documentation says the ref/het/hom counts cover SNV sites only, and measured against the tool itself they do not — a call at a site carrying no SNV allele, such as a gVCF's reference block, is counted too. On a gVCF this is the difference between a hom-ref count of one and a hom-ref count of every block in the file. Haploid calls are kept in their own columns rather than folded into missing data, because a haploid reference call and an absent call are different observations.

What the Ts/Tv ratio can and cannot tell youCopy Anchor Link

This is the one figure the panel interprets, and the interpretation is bounded by what the literature actually supports: a low ratio implicates false positives. Substitutions do not arise at random with respect to class — transitions are chemically favoured — so a real set of calls runs well above the ratio you would get by chance, and a set diluted with sequencing and alignment artefacts is pulled back towards it.

For scale, the framework paper behind the widely used GATK calling pipeline (DePristo et al., Nature Genetics 43

–498, 2011) reports a recalibrated whole-genome human call set at 2.15 for known and 2.05 for novel variants, describing those as consistent with genome-wide expectations; call sets the same paper judges to be dominated by false positives come in at 1.2, 0.7 and 0.30.

Two cautions, and the panel states both rather than leaving them to be assumed:

• Whether that comparison applies at all depends on what was sequenced. A VCF does not record whether it came from a whole genome, a capture panel, an amplicon or a non-human organism, and the expected ratio differs between them — coding regions in particular run higher. The panel therefore never tells you which kind of call set you have, and you should not read its comparison as one.
• It is a property of the set, not of a record. A high overall ratio says nothing about whether any individual call is right.
• It needs enough SNV alleles to mean anything. The panel compares the ratio with the figures above only when its whole 95% interval falls in one of the ranges they describe. A single SNV allele reads 0.00 with an interval of 0.00 to 39.00, and 106 alleles at 1.36 still span 0.91 to 2.04, so for sets like these the panel says the counts are too few to compare rather than calling the set suspect. It flags a low ratio only when the interval's upper end sits below 1.5.

What the panel refuses to let a number meanCopy Anchor Link

Each of these appears beside the figures it affects, naming what it is about and why:

• A region read is the region's. Reading a locus out of an indexed file gives counts for that window. They are never presented as the file's, and the panel always says which of the two you are looking at. A record whose span reaches into the window is in it, even when its own position is before the start: a gVCF reference block beginning at 4,370 and ending at 4,383 is returned for a window starting at 4,380, because the question a locus asks is which calls cover these bases, not which calls begin inside them.
• A file that declares no ##fileformat line is read, and said so. Every version of the VCF spec asks for that line first, and htslib does not open a file without it, so bcftools and the tools built on it refuse such a file outright with unknown file type. This reader is more permissive: a #CHROM line names the columns on its own and the records parse. The panel names the file and the consequence, because a set that reads here and is refused at the next step is worth knowing about before that step.
• A file too large to read whole opens at its header. Then every count is of nothing until you read a locus, and the panel says so rather than showing zeroes.
• A read that stopped at a cap is a prefix, not the file.
• Genotypes not parsed means the per-sample section is empty because of a size limit, not because the file lacks genotypes.
• A genotype naming an allele its own record never declares is reported with the sample, the position and the offending allele. bcftools stats stops on such a file rather than reporting per-sample numbers for it, so where this appears the per-sample counts have nothing to be checked against — which is a stronger warning than the count being unusual.

Reading it in orderCopy Anchor Link

1. Open the call set. The conclusion at the top of the panel is one sentence: how many records, over what, dominated by which class, at what Ts/Tv, in how many samples.
2. Read the chain under it. Each link is one step of the derivation, and expands to say how that figure was arrived at — which counting convention, over which records, and with what left out.
3. Read the flags. If any figure cannot be taken at face value, it is here, with its reason.
4. Open the counts for the full tables, including both Ts/Tv arms with their own counts.
5. Narrow and repeat. Filtering to passing records, or to one contig, is a different measurement — and comparing the two is often the useful thing rather than either alone.

Doing this without the panelCopy Anchor Link

The same three operations are available wherever you work: reading a call set and querying it by locus, class, quality floor, filter status, sample or INFO key; summarising it exactly as above; and applying one called alternate allele to a sequence. That last one checks the reference allele against the sequence first and refuses on a mismatch, naming both, because applying a call from a different assembly produces a construct that looks correct and is not.

Not every called allele is applicable, and the ones that are not are refused by name. A symbolic allele (<DEL>, <DUP>, <INV>, <CNV> and the rest) and a breakend (written with a square bracket, such as G]chr2:321681]) name a rearrangement rather than the bases that replace the reference, so there is nothing to substitute: the record says a segment was lost here, not these letters go here. The * allele is refused for a different reason — it marks a site overlapped by a deletion called at another record, so the edit that allele stands for belongs to that other record and is applied from there. In each case the panel says which of the three it is rather than failing or, worse, writing the allele's own text into the sequence.