GradeVCF

Script: gradevcf.sh Package: var2 Class: GradeVCF.java

Grades a VCF against a truth set and produces concordance metrics and cumulative histograms. Given a candidate VCF (from CallVariants or filtervcf) and a benchmark truth VCF (e.g. GIAB), GradeVCF counts true positives (TP), false positives (FP), and false negatives (FN), then reports precision, recall, F1, FPR, and FNR at a configurable operating point. Optionally applies a neural network scorer and writes cumulative TP/FP/FN curves suitable for ROC or precision-recall plotting. Replaces a multi-step filtervcf + comparevcf pipeline with a single fast tool.

Basic Usage

gradevcf.sh in=calls.vcf truth=giab.vcf

Both in= and truth= are required. The input VCF is typically produced by callvariants.sh or filtervcf.sh clearfilters. The truth VCF should be a normalized benchmark such as a GIAB high-confidence call set.

Counting Contract

GradeVCF follows a marking contract to avoid double-counting compound sites. The truth set is loaded into a hash map using the same variant identity (equals/hashCode) as CompareVCF. Each truth variant is marked at most once — when a call matches it, the maximum matching-call score is recorded. TP equals the number of distinct marked truth variants; FN equals unmarked truth variants; FP equals call variants that matched no truth entry. This means a multi-allelic call that splits into two identical components counts as a single TP, not two.

Memory and Threading

GradeVCF uses 42% of available memory by default (4 GB minimum, overridable with -Xmx) and caps worker threads at 8. The truth set is loaded once, then variant lines are processed in parallel across threads. Thread results are merged by keeping the maximum matching-call score per truth variant before computing final metrics.

Parameters

I/O Parameters

in=<file>: Input candidate VCF. Required. Typically the output of callvariants.sh or filtervcf.sh clearfilters. Must contain ##contig lines (or ref= must be provided) so that the scaffold map can be populated.
truth=<file>: Truth VCF. Required. A benchmark call set such as a GIAB high-confidence VCF. Should be normalized and, for best results, restricted to the same high-confidence regions as the BED file.
out=<file>: Optional. Write variants that pass the operating-point threshold to this file. In QUAL mode the threshold is minscore; in NN mode it is netcutoff. The VCF header from the input file is preserved.
ref=<file>: Reference FASTA. Optional unless net= is used (required when using a neural network, because the network recomputes contigEndDist and homopolymerCount features from reference bases). Also required when normalize=t is set.
overwrite=t: Grant permission to overwrite existing output files. Set to false to abort rather than overwrite.

Region Parameters

bed=<file>: Restrict evaluation to variants inside these BED intervals (e.g. a GIAB high-confidence benchmark region set). Truth variants outside the BED are excluded from the truth map; calls outside the BED are excluded from FP counting. The number of loaded intervals and scaffolds is reported to stderr.
bedweight=0.0: Concordance weight for variants outside the high-confidence BED region. Default 0 reproduces the old hard BED restriction: off-region variants contribute nothing. Setting bedweight=0.5 makes off-region variants count half. Setting bedweight=1.0 disables region restriction entirely while still accepting a BED file.
invertbed=f: When true, the high-confidence region becomes the complement of the BED (i.e. grade only variants OUTSIDE the BED intervals). Requires bed= to be set.

Quality Cutoff Parameters

minscore=0.0: Reject variants with QUAL below this Phred-scaled threshold. This is the operating-point threshold when no neural network is used. Variants below this value are still scored and counted for histograms unless a varFilter pre-filter also rejects them. Aliases: minqual, minphred.

Neural Network Parameters

net=<file>: Neural network file (.bbnet). When provided, GradeVCF scores each variant through the network and uses the NN score as the operating-point axis instead of QUAL. Requires ref= (the network recomputes contigEndDist and homopolymerCount from reference bases; without ref these features are silently corrupted).
netmode=: Feature vector mode controlling which input features are assembled for the network. Recognized values: ump45, elba, lawrence, donovan. Must match the mode used during network training.
netcutoff=0.5: NN output threshold. Variants with NN score below this are excluded from the operating-point TP/FP/FN summary and, if out= is set, from the output VCF. Use netcutoff=auto to read the cutoff that was embedded in the .bbnet file at training time. Alias: cutoff.
includescore=f: Include the composite QUAL score as an additional feature in the NN input vector.

Histogram Output Parameters

histnn=<file>: Write a cumulative NN score histogram. The file contains two sections: (1) raw per-bin true and false counts across 101 bins (0.00–1.00 in steps of 0.01); (2) a cumulative threshold/TP/FP/FN table sweeping from high score to low. Suitable for plotting precision-recall or ROC curves. Alias: nnhist.
histqual=<file>: Write a cumulative QUAL histogram. The file contains two sections: (1) raw per-bin true and false counts across 241 bins (0.00–60.00 in steps of 0.25); (2) a cumulative threshold/TP/FP/FN table. Useful for tuning QUAL cutoffs. Alias: qualhist.

Variant Normalization Parameters

normalize=f: Left-align and normalize indels before matching, using the loaded reference sequence. Requires ref=. Applied to both the truth set and the candidate calls. Aliases: leftalign, norm.
splitalleles=f: Split multi-allelic VCF records into individual biallelic records before matching. Applied to both truth and calls.
splitsubs=f: Split multi-nucleotide polymorphisms (MNPs) into individual SNPs before matching. Applied to both truth and calls. Alias: splitsnps.
split=f: Convenience flag that sets splitalleles=t, splitsubs=t, and normalize=t simultaneously. Equivalent to sass=t. Requires ref=.

Variant-Quality Pre-Filtering Parameters

These parameters apply a VarFilter pre-filter before grading. Variants failing the filter are ignored entirely — they neither mark truth variants nor count as FP. If no filter parameters are set, no pre-filtering is applied (equivalent to clearfilters).

minreads=0: Ignore variants seen in fewer than this many reads.
minqualitymax=0: Ignore variants with max base quality below this threshold.
minedistmax=0: Ignore variants with max distance from read ends below this threshold.
minmapqmax=0: Ignore variants with max mapping quality below this threshold.
minidmax=0: Ignore variants with max read identity below this threshold.
minscore=0.0: Ignore variants with Phred-scaled composite score below this threshold.
clearfilters: Reset all variant filter thresholds to zero, disabling pre-filtering. Use this when the input VCF already contains only the variants to evaluate (e.g. output of filtervcf clearfilters).

Java Parameters

-Xmx: Set Java's memory usage, overriding autodetection. -Xmx20g will specify 20 GB of RAM, and -Xmx200m will specify 200 MB. The default is 42% of physical memory (minimum 4 GB). The max is typically 85% of physical memory.
-eoom: Cause the process to exit if an out-of-memory exception occurs. Requires Java 8u92+.
-da: Disable assertions.

Examples

Basic Grading Against GIAB

# Grade calls against a GIAB truth set restricted to high-confidence regions
gradevcf.sh in=calls.vcf truth=HG001_GRCh38_GIAB_highconf.vcf bed=HG001_GRCh38_GIAB_highconf.bed

Loads all truth variants inside the BED, matches call variants against them, and prints TP/FP/FN counts with precision, recall, F1, FPR, and FNR to stderr. No files need to be written for this summary.

Grading with Normalization and Output VCF

# Normalize indels, split multi-allelics, write passing calls to a file
gradevcf.sh in=calls.vcf truth=truth.vcf ref=GRCh38.fa \
            bed=highconf.bed split=t minscore=20 out=passing.vcf

split=t sets splitalleles=t, splitsubs=t, and normalize=t simultaneously. Variants with QUAL < 20 are excluded from the output VCF and from the operating-point summary. The histogram still reflects all variants above score 0.

Neural Network Scoring with Histogram Output

# Score with a trained .bbnet model and write precision-recall curve data
gradevcf.sh in=calls.vcf truth=truth.vcf ref=GRCh38.fa \
            net=callvars.bbnet netcutoff=auto \
            bed=highconf.bed \
            histnn=nn_curve.tsv histqual=qual_curve.tsv

The netcutoff=auto flag reads the operating-point cutoff embedded in the .bbnet file. Both histogram files contain a raw per-bin section and a cumulative TP/FP/FN section. Plot the cumulative section to visualize the full precision-recall or ROC tradeoff across all thresholds.

Algorithm Details

Variant Matching

Variant matching uses the same identity contract (equals/hashCode on VCFLine) as CompareVCF. Two variants are considered identical if they have the same chromosome, position, reference allele, and alternate allele after optional normalization and splitting. This means that pre-normalizing the truth set with an external tool (e.g. bcftools norm) and setting normalize=t in GradeVCF should produce identical results as either alone, since both left-align indels in the same way.

Marking Contract and Double-Counting Prevention

The truth set is stored in a HashMap (key == value) so that a call variant can retrieve the canonical truth object via get(). When a call matches a truth variant, the truth object is marked by inserting it into a per-thread matched map with the maximum {NN score, QUAL score} of any matching call. After all threads finish, the per-thread maps are merged by taking the per-truth maximum. This design ensures:

Each truth variant is counted as TP at most once, even if multiple calls match it.
A compound site that splits into two identical components does not inflate TP counts.
The TP score recorded is the highest-confidence call that matched a given truth variant, giving a conservative operating-point curve.

Histogram Construction

Two parallel histograms are maintained: one indexed on NN score (101 bins, 0.00–1.00) and one on QUAL score (241 bins, 0.00–60.25). For each histogram:

trueCounts[bin]: Each distinct matched truth variant is binned once, at the bin of its maximum matching-call score. Built after thread merge, so no truth is counted twice.
falseCounts[bin]: Each unmatched call is binned at its own score. Summed directly across threads.

The cumulative section of each output file sweeps from high threshold to low. Cumulative FP sums false counts from high bins down; cumulative FN starts at the number of truth variants never matched by any call (callerMissed) and decreases as the threshold drops. At threshold 0 all matched truth variants become TP and FN equals callerMissed.

BED Region Weighting

With bedweight=0 (default), truth variants outside the BED are never added to the truth map, and off-region calls receive weight 0 (contribute nothing to FP counts or histograms). This exactly reproduces a hard-filter approach. When bedweight is set to a value between 0 and 1, off-region variants participate at fractional weight, allowing soft evaluation across confidence regions. Setting bedweight=1 makes the BED file act as a partitioning hint only, with no effect on counts.

Neural Network Scoring Requirements

When net= is used, the feature vector is assembled by VectorUMP45.makeVector, which recomputes two reference-derived features: contigEndDist (capped at nScan using scaffold length) and homopolymerCount (count of matching bases in a window). These features require the actual reference base sequence, not just scaffold lengths. For this reason, ref= is mandatory when net= is specified. Without ref=, these features would be silently corrupted rather than raising an error, producing incorrect NN scores.

A future improvement would read CED= and HMP= tags already written to the VCF INFO field by callvariants, eliminating the ref= requirement for post-hoc scoring entirely.

Operating-Point Summary

The summary printed to stderr at the end of each run reports metrics at a single operating point: the configured minscore (QUAL mode) or netcutoff (NN mode). These are point estimates. For the full precision-recall or ROC curve, write histnn= or histqual= and plot the cumulative section.

Reported metrics:

Precision = TP / (TP + FP)
Recall = TP / (TP + FN) = TP / truthSize
F1 = 2 × Precision × Recall / (Precision + Recall)
FPR = FP / (TP + FP)
FNR = FN / truthSize

Related Tools

callvariants.sh — Calls variants from a mapped BAM/SAM, producing the input VCF for GradeVCF.
filtervcf.sh — Applies quality filters to a VCF before grading. Use clearfilters to pass all variants for GradeVCF to evaluate.
comparevcf.html — Earlier per-file comparison tool. GradeVCF replaces the filtervcf+comparevcf two-step pipeline.
invertvcf.sh — Invert a VCF filter to switch passing and failing variants.

Support

For questions and support:

Email: bbushnell@lbl.gov
Documentation: bbmap.org