Combine, split, sort, intersect, and subset VCF/BCF files with bcftools merge, concat, isec, sort, view, and reheader. Use when merging different samples into a cohort VCF, concatenating per-chromosome or per-region call sets for the same samples, intersecting or complementing call sets from different callers, subsetting samples/regions, harmonizing sample names and
Combine, split, sort, intersect, and subset VCF/BCF files with bcftools merge, concat, isec, sort, view, and reheader. Use when merging different samples into a cohort VCF, concatenating per-chromosome or per-region call sets for the same samples, intersecting or complementing call sets from different callers, subsetting samples/regions, harmonizing sample names and
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.
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)
-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.
-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.
-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)
