CompareVCF

Script: comparevcf.sh Package: var2 Class: CompareVCF.java

Performs set operations on VCF files: union, intersection, and subtraction. Supports BedMask region filtering to restrict comparisons to benchmark regions, left-alignment normalization for cross-caller concordance, and variant splitting to decompose complex alleles into comparable forms.

Basic Usage

comparevcf.sh in=<file,file,...> out=<file>

CompareVCF requires at least two input VCF files. By default it performs subtraction: variants present in the first file but absent from all subsequent files are written to the output. Files may be gzip-compressed.

Cross-Caller Concordance

For accurate comparison across different variant callers, enable both normalization and splitting before comparing. Different callers often represent the same variant differently; normalization ensures consistent representation:

comparevcf.sh in=caller1.vcf,caller2.vcf out=unique.vcf \
    ref=reference.fa normalize=t splitalleles=t canonize=t

Parameters

Parameters are organized by function. At least two input files are required; all other parameters are optional.

I/O Parameters

in=<file>: Input VCF files; must be at least 2 files, comma-separated (e.g. in=a.vcf,b.vcf,c.vcf). Files may be gzip-compressed. The order matters for subtraction mode: the first file is the minuend and all others are subtracted from it.
out=<file>: Output file for the result of the set operation. Written in VCF format, including the header from the first input file.
ref=<file>: Reference FASTA file. Required when normalize=t (left-alignment needs the reference sequence). Optional otherwise.
shist=<file>: (scorehist) Output file for a variant quality score histogram. Useful for analyzing the quality distribution of retained variants.
overwrite=f: (ow) Set to false to abort rather than overwrite an existing output file. Default is to overwrite.
bgzip=f: Use bgzip for gzip compression. bgzip-compressed files can be indexed with tabix for random access.

Mode Parameters (choose one only)

subtract=t: (difference, minus, diff) Subtract all other files from the first file. This is the default mode. Output contains variants present in the first input file but absent from every subsequent file. Useful for finding variants unique to one caller or sample.
union=f: (plus) Make a union of all files. Output contains all unique variants from all input files combined. Useful for building a comprehensive variant set from multiple callers.
intersection=f: (shared) Make an intersection of all files. Output contains only variants present in every input file. Useful for finding high-confidence variants confirmed by all callers.

Region Filtering Parameters

bed=<file>: (bedfile) BED file defining genomic intervals to restrict the comparison. Only variants inside these intervals are included in the comparison. Commonly used with GIAB high-confidence benchmark region sets (e.g. HG001_GRCh38_benchmark.bed) to ensure evaluation is performed only in reliably callable regions.
invertbed=f: (excludebed, bedexclude, bedinvert) Invert the BED filter: compare only variants OUTSIDE the specified intervals. Default false (keep variants inside the intervals).

Normalization and Canonicalization

normalize=f: (leftalign, norm) Left-align indels using the reference sequence. This is the standard normalization for cross-caller concordance: a deletion can be represented at multiple positions when located in a repeat region, and different callers may choose different positions. Left-alignment forces all equivalent representations to the leftmost position. Requires ref=.
canonize=t: (canonicalize, canonicize, trimtocanonical) Trim variants down to a minimal canonical representation by removing redundant reference bases from insertions and deletions. For example, "ATCG/A" and "TCG/" represent the same deletion; canonization produces the same minimal form from both. Enabled by default.

Variant Splitting Parameters

splitalleles=f: Split multi-allelic lines into individual lines. A line with ALT "A,T" becomes two separate variant lines, one for each alternate allele. Enables accurate comparison when one caller outputs multi-allelic records and another outputs biallelic records for the same variants.
splitsubs=f: (splitsnps) Split multi-base substitutions into individual SNPs. A MNP (multi-nucleotide polymorphism) such as "ATG/CGT" is decomposed into three individual substitutions. Enables comparison between callers that differ in MNP handling.
splitcomplex=f: Split complex variants (those involving both insertions and deletions) into their component parts. Useful when callers represent complex rearrangements differently.
sass=f: (split) Convenience alias: equivalent to setting splitalleles=t and splitsubs=t simultaneously.
splitall=f: (sascsss) Convenience alias: equivalent to setting splitalleles=t, splitsubs=t, and splitcomplex=t simultaneously. Splits all variant types into their simplest forms.

Quality Filtering Parameters

minscore=-99999: (minqual, minq) Minimum variant quality score (QUAL field in VCF). Variants with quality below this threshold are excluded before set operations. Default of -99999 passes all variants. Set to e.g. minscore=20 to keep only variants with QUAL ≥ 20.
lines=<number>: Maximum number of input lines to process. Set to -1 for no limit (default). Useful for quick testing with large files.

Miscellaneous Processing Parameters

addsamples=t: Include sample columns from all input files in the output lines. When enabled, sample information from all input files is preserved in the output VCF.
verbose=f: Enable verbose output for debugging. Produces additional diagnostic information to stderr.

Java Parameters

-Xmx: Set Java's memory usage, overriding autodetection. -Xmx20g specifies 20 GB of RAM; -Xmx200m specifies 200 MB. CompareVCF loads all variants into memory, so large variant sets may require explicit memory allocation. The default allocation is 4 GB.
-eoom: Exit the process if an out-of-memory exception occurs. Requires Java 8u92+.
-da: Disable assertions for slightly better performance in production environments.

Examples

Basic Subtraction (Default Mode)

comparevcf.sh in=caller1.vcf,caller2.vcf out=unique_to_caller1.vcf

Find variants called by caller1 but not by caller2. The first input file is the one being filtered; subsequent files define what to subtract.

Union Operation

comparevcf.sh in=caller1.vcf,caller2.vcf,caller3.vcf out=all_variants.vcf union=t

Combine all unique variants from three callers into a single comprehensive variant set.

Intersection Operation

comparevcf.sh in=caller1.vcf,caller2.vcf,caller3.vcf out=consensus.vcf intersection=t

Keep only variants called by all three callers. These high-confidence variants are confirmed by independent algorithms.

BedMask: Restrict to Benchmark Regions

comparevcf.sh in=calls.vcf,truth.vcf out=false_positives.vcf \
    bed=HG001_GRCh38_benchmark.bed

Compare calls to a truth set, but only within GIAB high-confidence benchmark regions. Variants outside the BED intervals are excluded from the comparison entirely, preventing penalization for calls in difficult regions.

BedMask: Exclude Repetitive Regions

comparevcf.sh in=caller1.vcf,caller2.vcf out=outside_repeats.vcf \
    union=t bed=repeats.bed invertbed=t

Make a union of variants from both callers, keeping only variants that fall outside known repetitive regions (invertbed=t reverses the filter).

Normalization for Cross-Caller Concordance

comparevcf.sh in=gatk.vcf,deepvariant.vcf out=gatk_unique.vcf \
    ref=GRCh38.fa normalize=t splitalleles=t canonize=t

Find variants called by GATK but not DeepVariant, with full normalization to ensure equivalent variants are recognized as identical regardless of caller-specific representation choices. normalize=t requires ref=.

Benchmark Evaluation with Normalization and BedMask

# False positives: in calls but not in truth, within benchmark regions
comparevcf.sh in=calls.vcf,truth.vcf out=false_positives.vcf \
    ref=GRCh38.fa normalize=t splitalleles=t canonize=t \
    bed=HG001_GRCh38_benchmark.bed

# False negatives: in truth but not in calls, within benchmark regions
comparevcf.sh in=truth.vcf,calls.vcf out=false_negatives.vcf \
    ref=GRCh38.fa normalize=t splitalleles=t canonize=t \
    bed=HG001_GRCh38_benchmark.bed

Standard workflow for evaluating a variant caller against a GIAB truth set. Both normalization and BedMask filtering are applied for a fair comparison. Use GradeVCF for a more complete precision/recall analysis.

Quality Filtering with Splitting

comparevcf.sh in=raw1.vcf,raw2.vcf out=high_quality_union.vcf \
    union=t minscore=30 splitall=t canonize=t

Create a union of high-quality variants (QUAL ≥ 30), fully splitting and canonizing all variants before comparison.

Score Histogram Output

comparevcf.sh in=variants1.vcf,variants2.vcf out=difference.vcf \
    shist=quality_distribution.txt

Perform subtraction and simultaneously generate a histogram of quality scores for the retained variants.

Algorithm Details

Set Operation Implementation

CompareVCF uses Java HashSet data structures to perform set operations. Each VCF file is parsed into VCFLine objects that serve as the atomic units for comparison. The three operations are:

Difference (subtract): Load the first file into a set; for each subsequent file, remove any variants that appear in it. What remains is unique to the first file.
Union: Iteratively add variants from each file to a single shared HashSet. Duplicate variants (same position and alleles) are automatically collapsed.
Intersection: Load the first file into a set; for each subsequent file, use HashSet.retainAll() to keep only variants present in both. The result is variants common to all files.

Variant Processing Pipeline

For each input file, variants pass through this pipeline before entering the set operation:

VCF Parsing: Each file is parsed by VCFFile, which extracts the header, sample names, and all variant lines into a map keyed by VCFLine objects.
Variant Splitting: If any splitting flag is enabled (splitalleles, splitsubs, splitcomplex), multi-allelic and complex variants are decomposed into simpler forms using VCFLine.split().
Left-Alignment: If normalize=t, each variant is left-aligned against the reference sequence using VCFLine.leftAlign(). This shifts indels to the leftmost equivalent position in repeat regions, ensuring callers that chose different positions for the same indel are recognized as identical.
BedMask Filtering: If bed= is specified, variants outside the defined intervals are discarded (or kept if invertbed=t). The BedMask.contains() check is O(log n) per variant using an interval tree over the loaded BED regions.
Quality Filtering: Variants with QUAL below minscore are discarded.
Canonicalization: If canonize=t, variants are trimmed to their minimal representation before being added to the HashSet. This step runs globally via VCFLine.TRIM_TO_CANONICAL, ensuring consistent representation during hash and equality comparisons.

Normalization and Left-Alignment

Left-alignment is the standard normalization method for VCF indel comparison. The same insertion or deletion in a repeat region can be validly placed at multiple genomic positions. Without normalization, two callers reporting the same biological indel at different positions within the repeat will appear as discordant calls. Left-alignment forces all equivalent representations to the leftmost position, resolving this ambiguity. The reference FASTA is required because the shifting algorithm needs to compare allele sequence to the adjacent reference bases.

BedMask Region Filtering

BedMask filtering is applied per-variant during file loading, before any set operations. When a BED file is provided, only variants whose position falls within a BED interval are included in the comparison. This is critical for fair benchmark evaluation: GIAB high-confidence regions exclude low-complexity and structurally variable regions where no caller is expected to perform reliably. Restricting comparison to these regions avoids penalizing callers for variants in difficult regions where the truth is uncertain.

At startup, the number of loaded BED intervals and scaffolds is printed to stderr for verification.

Variant Equality and Hashing

VCFLine objects implement equals() and hashCode() based on chromosome, position, reference allele, and alternate allele. Canonicalization ensures that variants with redundant leading or trailing reference bases hash and compare identically to their minimal forms. Splitting ensures that multi-allelic records and MNPs are decomposed to the same granularity across callers before comparison.

Memory and Performance

CompareVCF loads all variants from all input files into memory. Memory usage scales with the total number of unique variants across all files. For typical human WGS variant call sets (a few million variants per file), the default 4 GB allocation is sufficient. Very large sets or many input files may require increasing memory with -Xmx. HashSet operations (add, contains, remove) are O(1) average case, making set operations fast even for large variant counts.

Output and Header Management

The VCF header from the first input file is written to the output. Sample names from all input files are collected and can be included in output lines when addsamples=t. Output variants are sorted by genomic position before writing.

Related Tools

CallVariants — BBTools' variant caller. Produces the VCF files that CompareVCF operates on.
FilterVCF — Filters a single VCF file by score, type, allele frequency, or other criteria. Use FilterVCF to preprocess individual files before CompareVCF comparison.
GradeVCF — Evaluates a VCF call set against a truth set and computes precision and recall statistics. For complete benchmark analysis, GradeVCF is more informative than CompareVCF subtraction.
InvertVCF — Produces a VCF containing variants not present in the input. Conceptually related to CompareVCF subtraction, but operates on a single file against all possible variants.
LabelVCF — Annotates variants in one VCF based on presence or absence in another. Similar to CompareVCF but preserves the full original file with added labels rather than filtering.

Support

For questions and support:

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