GPTomics/bio-vcf-manipulation

Version Compatibility

Reference examples tested with: bcftools 1.19+, cyvcf2 0.30+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

The +fill-tags plugin ships with bcftools; --naive-force and -m snp-ins-del are recent additions -- confirm with bcftools concat --help / bcftools merge --help on the installed build.

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

VCF Manipulation

Combine, split, sort, intersect, and subset VCF/BCF files with bcftools.

"Combine, compare, or restructure my VCFs" -> Pick the operation from what changes (samples vs regions vs set membership), and normalize first so the same biological variant is recognized as the same row.

• CLI: bcftools merge (add samples), bcftools concat (add regions), bcftools isec (set operations), bcftools view (subset), bcftools sort / bcftools reheader (order and header fixes)

The governing principle: normalize BEFORE combining

Every combine operation here -- merge, concat -d dedup, isec, and downstream annotate -- keys on the (CHROM, POS, REF, ALT) tuple, and bcftools isec defaults to -c none (an ALT must match exactly to count as the same variant). An indel that is not left-aligned + parsimonious, an un-split multiallelic, or an un-decomposed MNP is a structurally valid VCF line carrying a different tuple for the same biological event. The combine then silently mis-joins: isec reports false discordance, merge emits duplicate rows and splits the allele frequency, dedup misses the duplicate. Nothing errors -- the counts are simply wrong.

Decision: normalize (left-align + parsimony against the SAME reference, split multiallelics) every input before any merge/concat-dedup/isec/annotate. This matters most for indels in homopolymers/STRs, where callers legitimately disagree on POS. It is safe to skip only when a single caller produced all inputs and no cross-file matching or database lookup follows. The canonical incantation is bcftools norm -m-any -f ref.fa; see variant-calling/variant-normalization for the full pipeline, MNP atomization, and the vt-vs-bcftools discordance -- do NOT re-derive that here, cross-reference it.

merge vs concat vs isec (choose by what differs)

| Operation | Inputs differ in | Produces | Requires index | Fails / misused when | |-----------|------------------|----------|----------------|----------------------| | bcftools merge | samples (same sites) | one multi-sample VCF (columns unioned) | yes | given the SAME sample split by region -> use concat; given single-sample VCFs and treated as joint genotyping -> fabricates 0/0 (see trap below) | | bcftools concat | regions (same samples, e.g. per-chromosome) | one VCF spanning all regions (rows appended) | only with -a | inputs from DIFFERENT samples -> use merge; inputs overlap without -a; --naive used when headers/sample-order differ | | bcftools isec | neither -- same cohort, compare membership | per-input private/shared partition dirs | yes | inputs not normalized to identical representation -> false discordance |

Common confusion: concat -a (--allow-overlaps) resolves duplicate records from the SAME sample across overlapping region files; it does NOT union genotypes across different samples -- that is merge. If bcftools reports a "different samples" error, the operation is inverted.

The merge trap: bcftools merge is not joint genotyping

When merging single-sample VCFs, a site called in sample A but simply absent from sample B's file is ambiguous: was B confidently homozygous reference there, or was B never covered/called? A project VCF cannot answer this. bcftools merge fills B's genotype with ./. (missing) by default, and -0/--missing-to-ref overrides it to 0/0 -- but both are guesses, because merge has no evidence for the unseen site. -0 therefore fabricates hom-ref genotypes and inflates the reference-allele count (see the vcf-basics ./.-is-not-0/0 distinction). Use -0 only when every input truly covered every site (e.g. gVCF-derived, or a shared target with confirmed coverage), never as a convenience to remove ./..

Decision: to build a multi-sample callset with correct hom-ref-vs-no-data resolution, joint-genotype gVCFs (GenomicsDBImport/CombineGVCFs -> GenotypeGVCFs), do not bcftools merge single-sample project VCFs -- see variant-calling/joint-calling. Reserve bcftools merge for combining already-jointly-genotyped cohorts, or samples that share a target with known coverage.

Merge also requires two harmonizations, or it silently drops or mis-collapses records:

• Consistent representation. All inputs must be normalized and split the same way first. If cohort A is split biallelic and cohort B keeps multiallelics, merge mis-collapses the shared site. Normalize all inputs identically (governing principle above).
• Matching ##contig headers and sample names. Merge unions sample columns; duplicate sample names abort unless --force-samples renames them, and mismatched contig naming (chr1 vs 1) prevents sites from aligning. Fix names/contigs with bcftools reheader first.

bcftools merge (combine different samples)

# Union samples across per-sample (already joint-genotyped or shared-target) VCFs
bcftools merge -l files.txt -Oz -o cohort.vcf.gz    # -l: one VCF path per line
bcftools index cohort.vcf.gz

• -m, --merge controls multiallelic collapse at shared sites (default both): -m none keeps a SNP and an indel at one POS as separate records; -m snps/-m indels restrict which types collapse. Leave the default unless a downstream tool needs types kept apart.
• --force-samples disambiguates colliding sample names; -r chr:beg-end restricts to a region (inputs must be indexed).

bcftools concat (stitch regions for the same samples)

# Genome-wide file from per-chromosome calls (same samples, disjoint regions)
bcftools concat chr{1..22}.vcf.gz chrX.vcf.gz -Oz -o genome.vcf.gz

• -a, --allow-overlaps is needed when region files overlap (e.g. windowed calling); pair with -d/--rm-dups <snps|indels|both|all|exact> to output a duplicate once. -a requires indexed inputs.
• -n, --naive concatenates BCF/VCF blocks WITHOUT recompression -- very fast for a large per-chromosome set, but it does only a header-compatibility check and requires identical headers and identical sample order across all files; it cannot reorder or reconcile anything. --naive-force skips even the header check and will silently produce a corrupt file if headers differ -- avoid it unless the files were provably produced identically.
• concat does NOT sort across file boundaries; overlapping unsorted inputs need -a, and the final file may still need bcftools sort.

bcftools sort (order by CHROM then POS)

bcftools sort -T /scratch/tmp -m 4G input.vcf.gz -Oz -o sorted.vcf.gz   # -T tempdir, -m spill threshold for large files

An unsorted VCF breaks everything downstream: tabix/bcftools index require coordinate-sorted input to build the index, and merge/isec/view -r all depend on that index for random access. Sort after any operation that can leave records out of order (naive concat of misordered files, some reheader edits). -T/-m bound memory for genome-scale files.

bcftools isec (set operations on call sets)

# Normalize BOTH first (governing principle), then partition
bcftools norm -m-any -f ref.fa gatk.vcf.gz     -Oz -o gatk.norm.vcf.gz
bcftools norm -m-any -f ref.fa freebayes.vcf.gz -Oz -o fb.norm.vcf.gz
bcftools isec -p comparison -Oz gatk.norm.vcf.gz fb.norm.vcf.gz

-p dir writes the four-way partition (-Oz to compress):

| File | Contents | |------|----------| | 0000.vcf[.gz] | private to file 1 | | 0001.vcf[.gz] | private to file 2 | | 0002.vcf[.gz] | shared, file-1 records (file-1 INFO/FORMAT) | | 0003.vcf[.gz] | shared, file-2 records (file-2 INFO/FORMAT) |

0002 and 0003 are the SAME sites with each file's own annotations -- pick by which annotations are needed downstream. Select membership instead of the full partition with -n and route records with -w (1-based file indices):

| Flag | Meaning | |------|---------| | -n=2 -w1 | present in exactly 2 files, output file-1 records | | -n+2 -w1 | present in >=2 files | | -n~10 -w1 | present in file1 but NOT file2 (boolean mask) | | -C | complement: positions only in file1, missing in the rest |

-c, --collapse sets what counts as "the same record"; the default none demands an exact REF+ALT match (why normalization is mandatory first), whereas -c all matches on position alone and ignores ALT -- rarely what a caller comparison wants.

Subsetting samples and regions (bcftools view)

bcftools view -s sample1,sample2 input.vcf.gz -Oz -o subset.vcf.gz   # -s ^s3 to EXCLUDE; -S file for a list
bcftools view -r chr1:1e6-2e6      input.vcf.gz -Oz -o region.vcf.gz  # -R file.bed for many regions

Two nuances that bite:

• -r/-R (regions) vs -t/-T (targets). -r/-R use the index to JUMP to regions (fast, require an index) and consider both POS and an indel's end; -t/-T STREAM the whole file filtering on POS (no index needed, slower). With -R, overlapping regions in the BED can emit a record MORE THAN ONCE and out of order -- deduplicate/sort after, or use non-overlapping regions.
• Stale INFO counts after subsetting. Dropping samples makes INFO AC/AN/AF wrong. bcftools view -s updates AC/AN by default (unless -I/--no-update), but recompute the full tag set explicitly: bcftools +fill-tags subset.vcf.gz -Oz -o out.vcf.gz -- -t AC,AN,AF.

Header harmonization (bcftools reheader)

printf 'old_name\tnew_name\n' > rename.txt
bcftools reheader -s rename.txt input.vcf.gz -o renamed.vcf.gz   # -s renames samples only, no record rewrite

reheader rewrites only the header (fast, no record pass): -s maps sample names, -h swaps in a whole new header, -f ref.fa.fai fixes ##contig lines to match a reference. Harmonize sample names and contigs BEFORE merge so columns and sites align.

Structural variants merge differently -- do NOT use bcftools merge

For SVs (<DEL>/<DUP>/<INV>/BND), "the same event" is fuzzy: breakpoints disagree by CIPOS/CIEND margins, so tuple-exact bcftools operations treat one deletion called by two tools as two variants. SV merging needs coordinate-and-size (ideally sequence) aware tools -- Truvari, SURVIVOR, or Jasmine -- whose distance/size parameters ARE the result. Use bcftools here only for small variants; route SV consensus to variant-calling/structural-variant-calling.

Quick Reference

| Task | Command | |------|---------| | Union samples | bcftools merge -l files.txt -Oz -o cohort.vcf.gz | | Stitch regions | bcftools concat chr{1..22}.vcf.gz -Oz -o genome.vcf.gz | | Fast stitch (identical headers) | bcftools concat --naive chr*.bcf -Ob -o all.bcf | | Sort | bcftools sort -T tmp input.vcf -Oz -o sorted.vcf.gz | | Compare callers | bcftools isec -p dir a.norm.vcf.gz b.norm.vcf.gz | | Shared only | bcftools isec -n=2 -w1 a.vcf.gz b.vcf.gz -Oz -o shared.vcf.gz | | Subset samples | bcftools view -s s1,s2 in.vcf.gz -Oz -o out.vcf.gz | | Recompute AC/AN/AF | bcftools +fill-tags in.vcf.gz -- -t AC,AN,AF | | Rename samples | bcftools reheader -s names.txt in.vcf.gz |

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | different samples on concat | merge/concat inverted (different samples given to concat) | Use merge for samples, concat for regions | | False discordance in isec | inputs not normalized to one representation | bcftools norm -m-any -f ref.fa both first (see variant-normalization) | | Duplicate rows / split AF after merge | inputs represented inconsistently, or un-normalized indels | Normalize + split all inputs identically before merge | | Fabricated 0/0 genotypes, inflated ref-allele count | -0/--missing-to-ref on single-sample merge (not joint genotyping) | Drop -0; joint-genotype gVCFs instead (joint-calling) | | not sorted / index build fails | unsorted records | bcftools sort then re-index | | --naive output corrupt | headers or sample order differ across inputs | Reheader to a common header, or drop --naive | | Records duplicated / out of order after -R | overlapping regions in the BED | Use non-overlapping regions, then sort/dedup | | Sample-name conflict aborts merge | duplicate sample names across files | --force-samples, or reheader -s first | | Stale AF after subsetting samples | INFO not fully recomputed | bcftools +fill-tags -- -t AC,AN,AF |

Related Skills

• variant-calling/variant-normalization - Normalize (left-align, split, atomize) before any merge/isec -- the load-bearing prerequisite
• variant-calling/joint-calling - Joint-genotype gVCFs instead of merging single-sample VCFs (correct hom-ref vs no-data)
• variant-calling/vcf-basics - VCF fields, the ./.-is-not-0/0 distinction, sample/region query
• variant-calling/structural-variant-calling - SV merging by breakpoint fuzz (Truvari/SURVIVOR/Jasmine), not bcftools
• variant-calling/filtering-best-practices - Filter call sets before combining
• variant-calling/vcf-statistics - Sanity-check Ti/Tv and counts after manipulation
• variant-calling/variant-calling - Upstream variant discovery that produces input VCFs

References

• Danecek P, Bonfield JK, Liddle J, et al. Twelve years of SAMtools and BCFtools. GigaScience. 2021;10(2):giab008. doi:10.1093/gigascience/giab008 (bcftools merge/concat/isec/norm/view/reheader reference implementation)
• Tan A, Abecasis GR, Kang HM. Unified representation of genetic variants. Bioinformatics. 2015;31(13):2202-2204. doi:10.1093/bioinformatics/btv112 (why normalization before tuple-keyed merge/isec is mandatory)