A missing line in a VCF is not a clean result

A VCF from a 30x whole genome holds a few million records out of about three billion positions, because it lists differences only. A gVCF adds reference blocks, stretches where the caller looked and found the reference base, which lets a reader tell a position that matched from a position that was never read. A plain VCF cannot tell those apart: a missing line means either one. Genotyping chips differ again. They report every probe, reference or not, but only for the positions the chip was built to measure.

Silence has two meanings

When a plain VCF has no line at a position, the person may match the reference there, or the sequencer may never have covered it well enough to call. The file does not say which. For common positions the first reading is far more likely, so reports usually assume it. We do the same on a plain VCF, but the report labels those readings as inferred rather than read. It also cannot say "examined and clear" about a gene on a plain VCF, because nothing in the file shows the gene was covered. A gVCF can, because its reference blocks do.

The reference moved

No record means the person matches the reference base, but which reference? A file called against GRCh37 is silent where the person matches GRCh37. If the pipeline lifts that file to GRCh38 and then reads the silence against the GRCh38 base, it assumes a genotype the person was never tested for, wherever the two assemblies disagree.

We found this by comparing a whole genome with four chip exports from the same person. At rs4430796, a variant in one of our prostate-cancer scores, GRCh37 has G and GRCh38 has A, and the score counts copies of A. The person is G/G: all four chips agreed. On GRCh37 that is the reference, so the genome file had no record. After the lift, the pipeline read the missing record against GRCh38 and credited two copies of A, an allele the person does not carry. The chips were right and the whole genome was wrong.

Across the 920 variants in our five disease scores, the assemblies disagree at 15. We now pin the GRCh37 reference base of every scoring variant and read a GRCh37 file against it. Where that base cannot be pinned down, the variant is left out of the score and the coverage figure says so.

The existing checks missed it. The script that checks our reference bases asks whether the base is one of the two alleles the score names, and at rs4430796 the GRCh38 base is the effect allele, so it passed.

Depth is not a call

A gVCF block can say the position was sequenced 30 times and still carry the genotype ./., which means the caller looked and could not decide. We used to count any block with enough depth as covered. On one real 3.9 GB whole-genome gVCF that made about 99 percent of the polygenic-score positions look covered while not a single variant was observed, and every score came out at the very bottom of its range. We now count a position as covered only when it has enough depth and a called genotype.

The file that was never finished

That same file turned out to be more interesting than a failed sequencing run. Its header named GATK 3's CombineGVCFs, a step that merges per-chromosome gVCFs and deliberately leaves every genotype blank for a later step to fill in. The later step never ran. Every record read ./., but the genotype likelihoods, the PL field, were intact.

In the part of the file our panel covers, 2,379,287 of 2,417,662 records carried a PL that clearly favoured one genotype, including roughly 52,700 confident non-reference calls. We added a pass that turns a clearly favoured PL into a genotype, so a file in this state is read instead of reported as empty.

The gate that decides whether to run that pass had its own lesson. It first required half of the opening 2,000 records to have a clear PL. A whole genome starts at the end of chromosome 1, where coverage is thin, and a real file measured 49 percent there. The gate now counts clear PLs instead of taking a share.

What Aimosti would (and wouldn't) show you

Each card says how many of its positions were read from your file and how many were inferred as the reference base. On a GRCh37 file, a missing record is read against the GRCh37 reference, and a position where that base cannot be pinned down is left out of the score and counted as unread.

What we won't claim

We won't turn a missing record into a genotype without knowing which reference the file was called against, and we won't count a position as read because the file reported sequencing depth there. An empty section in a report is not a clean bill of health.

Bottom line. No record means "matches the reference the file was called against, if this position was read at all". Before a report turns that silence into a genotype, it has to know which reference, and whether the position was read.

Related: Coverage: what each file type can tell you. Restated from: The VCF 4.2 specification · GATK: GVCF, Genomic Variant Call Format · dbSNP: rs4430796 · Genome Reference Consortium: human assemblies.