InvertVCF

Basic Usage

invertvcf.sh in=<mutations.vcf> out=<truth.vcf>

The input must be a VCF produced by mutate.sh. The output is an inverted VCF where all coordinates and alleles have been transformed into the coordinate space of the mutant genome, suitable for comparison with variant calls made against that mutant genome.

Variant Calling Accuracy Workflow

InvertVCF is designed for a specific benchmarking workflow using real reads. The complete pipeline is:

1. Start with a real genome and real reads from that genome (e.g., from a sequencing run or a public dataset).
2. Mutate the genome with mutate.sh to produce a mutant genome and a VCF describing the mutations.
3. Map the real reads to the mutant genome with bbmap.sh or a similar aligner.
4. Call variants with callvariants.sh, comparing the real reads against the mutant genome.
5. Invert the mutation VCF with invertvcf.sh to produce a ground-truth VCF in mutant-genome space.
6. Grade the variant calls with gradevcf.sh by comparing the called variants against the ground truth.

This workflow is realistic because it uses actual biological reads — not simulated reads — mapped to a slightly different genome. The mutations are known exactly, so variant-calling accuracy can be measured rigorously.

Parameters

Examples

# Step 1: Mutate the real genome (5% divergence) to create a mutant and a mutation VCF
mutate.sh in=genome.fa out=mutant.fa vcf=mutations.vcf id=0.95

# Step 2: Map real reads to the mutant genome
bbmap.sh ref=mutant.fa in=real_reads.fq out=mapped.sam

# Step 3: Call variants against the mutant genome
callvariants.sh in=mapped.sam ref=mutant.fa out=called.vcf ploidy=1

# Step 4: Invert the mutation VCF to mutant-genome coordinates
invertvcf.sh in=mutations.vcf out=truth.vcf

# Step 5: Grade the variant calls against ground truth
gradevcf.sh in=called.vcf truth=truth.vcf ref=mutant.fa

# Mutate genome
mutate.sh in=genome.fa out=mutant.fa vcf=mutations.vcf id=0.98

# Map paired-end reads
bbmap.sh ref=mutant.fa in1=reads_R1.fq in2=reads_R2.fq out=mapped.sam

# Call variants (paired-end provides better sensitivity)
callvariants.sh in=mapped.sam ref=mutant.fa out=called.vcf ploidy=1

# Invert and grade
invertvcf.sh in=mutations.vcf out=truth.vcf
gradevcf.sh in=called.vcf truth=truth.vcf ref=mutant.fa

# Mutate with only SNPs to evaluate SNP-calling accuracy
mutate.sh in=genome.fa out=mutant_snps.fa vcf=snp_mutations.vcf snprate=0.01 indelrate=0

# Mutate with only indels to evaluate indel-calling accuracy
mutate.sh in=genome.fa out=mutant_indels.fa vcf=indel_mutations.vcf snprate=0 indelrate=0.001

# Invert each mutation VCF separately
invertvcf.sh in=snp_mutations.vcf out=snp_truth.vcf
invertvcf.sh in=indel_mutations.vcf out=indel_truth.vcf

Algorithm Details

The Coordinate Inversion Problem

When mutate.sh inserts or deletes bases, it shifts the positions of all downstream variants. A deletion at position 100 that removes 5 bases means all variants originally recorded at positions 105, 110, 200, etc. are now at positions 100, 105, 195, etc. in the mutant genome. InvertVCF tracks this cumulative shift to correctly transform every variant's position.

Cumulative Indel Shift Tracking

InvertVCF processes the mutation VCF in two passes:

1. First pass: For each scaffold, the tool computes the total net length change from all insertions and deletions on that scaffold. This net shift per scaffold is used to update contig length annotations in the VCF header.
2. Second pass: The tool walks through variants in order, maintaining a running cumShift value. Before each variant is written, its position is adjusted by cumShift. After processing each variant, cumShift is updated by the delta of that variant (altLen - refLen). The cumulative shift resets to zero at the start of each new scaffold.

Allele Swapping

For each variant, the REF and ALT alleles are swapped: the original ALT (what was inserted into the mutant) becomes the new REF (what the mutant genome actually has), and the original REF (what the real genome has) becomes the new ALT (the variant being called). This correctly describes the variant from the perspective of someone who has the mutant genome as their reference and observes the original sequence in reads.

INFO Field Updates

The INFO field generated by mutate.sh contains three fields that InvertVCF modifies:

• STA=: Start coordinate of the variant. Adjusted by cumShift at the time of processing (before the current variant's delta is applied).
• STO=: Stop coordinate of the variant. Also adjusted by the same cumShift.
• TYP=: Variant type. INS is flipped to DEL, and DEL is flipped to INS. SNP and other types are left unchanged. This is correct: a base that was inserted to create the mutant looks like a deletion when you observe the real-genome reads mapped against the mutant.

Header Processing

The VCF header is updated in two ways:

• ##contig lines: The length= field for each scaffold is adjusted by the net indel shift on that scaffold, so the header accurately describes the mutant genome's scaffold lengths.
• ##Program line: A new ##InvertedBy=InvertVCF line is appended after the original program line, providing a record that this VCF has been inverted.

Ordering Requirement

InvertVCF assumes that variants within each scaffold are listed in positional order in the input VCF, as mutate.sh always produces them. Variants from different scaffolds may appear in any order, but mixing positions within a scaffold would produce incorrect cumulative shift calculations.

Related Tools

• mutate.sh: Creates the mutant genome and mutation VCF that invertvcf.sh inverts. The VCF format produced by mutate.sh is the expected input format for this tool.
• callvariants.sh: Calls variants by comparing mapped reads against a reference. Produces the VCF that is graded against the inverted ground truth.
• gradevcf.sh: Compares a called VCF against a truth VCF to compute sensitivity and precision. Used as the final step in the benchmarking workflow.
• bbmap.sh: Maps reads to a reference genome. Used in the workflow to map real reads to the mutant genome before variant calling.

Support

For questions and support:

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