How to Merge VCF Files

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-10

By the end of this guide you will be able to take two or more Variant Call Format (VCF) files and combine them correctly, whether that means adding new samples as genotype columns, stitching chromosome-by-chromosome shards into a single file, or building a consensus union from several different variant callers. You will also be able to prove the merge worked before you hand the output to anyone else.

You need five things on hand: a working bcftools installation (version 1.21 or newer is a safe baseline), htslib tools (bgzip and tabix), the reference FASTA used to call the variants plus its .fai index, the VCF files themselves, and a clear answer to one question. Are you combining different samples, or the same sample in different pieces? Almost every failed merge traces back to getting that answer wrong.

Optional extras depend on the job. Structural variant consensus needs SURVIVOR or an equivalent merger. Shard recombination inside a GATK pipeline can use GATK4 or Picard. Legacy scripts may expect VCFtools. Everything below uses free, open source software with commands you can copy directly.

Three Different Jobs People Call "Merging"

The word "merge" covers at least three operations in genomics, and they produce completely different files. Choosing the wrong one is the single most common source of broken downstream analysis. Establish which one you need before typing anything.

Job 1: Combining samples into one multi-sample file

You have sampleA.vcf.gz, sampleB.vcf.gz, and sampleC.vcf.gz, each called separately against the same reference. You want one file with three genotype columns so you can run a joint genotyping step, compute cohort allele frequencies, or run a case-control association test. The record count stays roughly at the union of variant positions across all three files. The sample count goes from one to three.

The standard tool is bcftools merge. GATK users sometimes reach for CombineGVCFs followed by GenomicsDBImport, but that path requires genomic VCF (gVCF) input with reference blocks, not ordinary VCFs.

Job 2: Concatenating shards of the same sample

Your pipeline split chromosome 1 into twenty chunks, called each chunk separately, and produced twenty small VCFs for one sample. Or you called each chromosome in parallel. Now you need one file for that sample covering the whole genome. The sample count stays at one. The record count becomes the sum of the inputs.

The standard tool is bcftools concat. bcftools merge is not a substitute, even though both accept a list of files.

Job 3: Building a consensus or union set across callers

You ran GATK HaplotypeCaller, DeepVariant, and Strelka on the same sample and you want one integrated call set. This is not a genotype merge. It is a reconciliation problem where different tools report overlapping but non-identical intervals with different quality semantics.

SV-MeCa, a structural variant meta-caller published in 2025, illustrates the pattern: it applies seven standalone SV callers and merges the resulting deletion and insertion calls into a union VCF using SURVIVOR, then extracts caller-specific quality metrics from each original VCF to feed an XGBoost classifier [1]. The union file is the scaffold. The per-caller evidence is what makes it useful.

Tools in this category include SURVIVOR, Jasmine, SVDB, and survivor-style mergers for small variants such as bcftools isec combined with custom logic. Comparative tools also exist that merge variants across samples for visualization, though a common limitation noted in that literature is that such tools merge records without automatically grouping variants across multiple samples for you [2].

What a VCF file actually contains

The VCF specification describes a text format for storing DNA polymorphism data such as SNPs, insertions, deletions, and structural variants together with rich annotations. Files are usually stored compressed and indexed for fast retrieval from position ranges [3]. Every VCF has three parts:

1. Meta-information lines beginning with ##, which declare the reference build, contig lengths, INFO and FORMAT field definitions, and any tool-specific metadata.
2. A column header line beginning with #CHROM, listing the eight fixed columns (CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO) followed by FORMAT and one genotype column per sample.
3. Data lines, one per variant record.

When you merge, the ## block and the #CHROM line are where most compatibility problems live. The data lines are comparatively forgiving.

Where VCF sits among bioinformatics file formats

Readers arriving from sequencing basics often ask about FASTQ versus FASTA. These are the upstream formats. A FASTQ file stores raw sequencing reads, each with a sequence and a per-base quality score encoded as ASCII characters. A FASTA file stores sequences only, with no quality information, and is used for reference genomes, protein sequences, and multiple alignments. VCF comes downstream of both: it holds the differences between reads and a reference, not the reads themselves. That distinction matters for merging because VCF merges are cheap bookkeeping operations on a few million lines, while read-level merges involve far larger files. The performance literature on sorted VCF merging exists precisely because a few million lines per file across thousands of files adds up [4].

Before You Run a Single Command: Five Checks

Skipping these checks is how people produce a merged file that looks fine and is silently wrong.

Check 1: Same reference build

Open the header of each file and read the ##reference line and any ##contig lines.

bcftools view -h sampleA.vcf.gz | grep -E '^##(reference|contig)' | head -5

Expected output:

##reference=file:///refs/GRCh38.fa
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>

If one file says GRCh37 and another says GRCh38, stop. Positions are not comparable between builds. You need to lift over one set, and lift-over itself can drop or duplicate variants at boundaries.

Check 2: Contig naming convention

chr1 and 1 are different strings to a computer. A file using Ensembl-style names and a file using UCSC-style names will merge into a file with two separate contig namespaces, and downstream tools will process only one of them.

bcftools query -f '%CHROM\n' sampleA.vcf.gz | sort -u | head -5
bcftools query -f '%CHROM\n' sampleB.vcf.gz | sort -u | head -5

Both commands should print the same style. If they do not, rename contigs in the smaller file with bcftools annotate, rename-chrs rename.txt, where rename.txt is a two-column mapping file.

Check 3: Unique sample names

bcftools query -l sampleA.vcf.gz
bcftools query -l sampleB.vcf.gz

Expected output for a single-sample file:

SAMPLE_A

If both files report an identical sample name, bcftools will refuse the merge. This happens constantly when the same pipeline template was reused and nobody edited the sample column. Fix the names before merging rather than forcing the merge.

Check 4: Files compressed with bgzip and indexed

bgzip produces BGZF, a block-compressed variant of gzip that supports random access. Plain gzip output cannot be indexed and cannot be read by bcftools in region mode. tabix -p vcf builds the .tbi index.

for f in *.vcf
do
  bgzip -k "$f"
  tabix -p vcf "${f}.gz"
done
ls -1 *.vcf.gz.tbi

Check 5: Sorted, normalized, and left-aligned

Deletions and insertions can be represented at several equivalent positions. Two files can describe the same 4 bp deletion at positions 1000 and 1004. Without normalization, they merge as two separate records and your downstream allele frequencies are wrong.

Step-by-Step: How to Merge VCF Files

The walkthrough below takes three single-sample VCFs through to a validated multi-sample file, then covers concatenation and consensus union. Each step lists the command, the expected output, and the failure mode.

Step 1: Install and stage the tools

conda create -n vcfmerge -c bioconda -c conda-forge bcftools=1.21 htslib tabix survivor
conda activate vcfmerge
bcftools --version

Expected output:

bcftools 1.21
Using htslib 1.21
Copyright (C) 2024 [Genome Research](/blog/guides/genome-research) Ltd.

If bcftools, version prints nothing, the environment did not activate. If it prints a version below 1.10, several options used below will not exist.

Step 2: Compress and index each input

for f in sampleA.vcf sampleB.vcf sampleC.vcf
do
  bgzip -k "$f"
  tabix -p vcf "${f}.gz"
done

Expected result: no terminal output, and six new files (sampleA.vcf.gz, sampleA.vcf.gz.tbi, and so on). A successful tabix run is silent. That silence is the success signal.

Common error at this step:

[E::hts_open_format] Failed to open file "sampleA.vcf.gz" : not compressed with bgzip

The file was compressed with ordinary gzip. Decompress with gunzip and recompress with bgzip.

Another common error:

[E::idx_find_and_load] Could not retrieve index file for 'sampleB.vcf.gz'

The .tbi file is missing or stale. Rebuild it with tabix -f -p vcf sampleB.vcf.gz. A stale index happens when you regenerate the VCF but keep an old index beside it.

Step 3: Normalize and left-align every file

for s in A B C
do
  bcftools norm -f GRCh38.fa -m -any -Oz -o "sample${s}.norm.vcf.gz" "sample${s}.vcf.gz"
  tabix -p vcf "sample${s}.norm.vcf.gz"
done

Expected output for each file:

Lines   total/split/realigned/skipped:  412300/1180/940/0

Read those four numbers. Total is the input record count. Split counts multiallelic records broken into single-ALT records by -m -any. Realigned counts indels shifted to their leftmost equivalent position. Skipped counts records that could not be normalized, usually because the reference allele does not match.

The -f GRCh38.fa argument is what enables left alignment and reference checking. The FASTA must have a .fai index, created with samtools faidx GRCh38.fa.

If you see a reference mismatch error:

[E::vcf_parse] Contig 'chr1' is not defined in the header. (Quick workaround: index the file with tabix.)

the VCF header is missing contig declarations, or its contig names disagree with the FASTA. Check that the FASTA and the VCF use the same naming convention. The -check-ref option controls how strict the check is, with w warning and x excluding mismatching records.

Step 4: Sort each file

Normalization can move records leftward, which breaks coordinate order.

for s in A B C
do
  bcftools sort -Oz -o "sample${s}.sorted.vcf.gz" "sample${s}.norm.vcf.gz"
  tabix -p vcf "sample${s}.sorted.vcf.gz"
done

Expected output:

Writing to /tmp/bcftools.XXXXXX
Merging 2 temporary files
Cleaning
Done

The temporary directory is controlled by -T. On a machine with limited /tmp space, point -T at a scratch volume with room for roughly the uncompressed file size.

An unsorted input produces this error downstream:

[E::hts_idx_push] Chromosome blocks not continuous

That message means the indexer encountered a position out of order. Sort the file and rebuild the index.

Step 5: Merge the samples

bcftools merge -Oz --threads 8 -o cohort.merged.vcf.gz \
  sampleA.sorted.vcf.gz sampleB.sorted.vcf.gz sampleC.sorted.vcf.gz
tabix -p vcf cohort.merged.vcf.gz

Expected result: a compressed VCF with three genotype columns and a record count close to, but generally larger than, the largest input. Confirm the sample columns:

bcftools query -l cohort.merged.vcf.gz

Expected output:

SAMPLE_A
SAMPLE_B
SAMPLE_C

The parameters that matter in bcftools merge

The -m option controls what happens when two records in the same file share a position. The default is -m both, which merges SNPs and indels at the same position separately. -m none leaves overlapping records as separate lines and is the safest choice when you do not want any within-file collapsing. -m id merges only records with matching ID values.

The -0 flag (long form missing-to-ref) treats absent genotypes in a sample as homozygous reference rather than missing. This is appropriate when every sample was called over the same confident regions and absence genuinely implies no variant. It is inappropriate when coverage varies, because low-coverage regions will be silently coded as reference.

The -R option restricts the merge to regions listed in a BED or tab-delimited file, which is useful for targeted panels.

The -l option reads the input file list from a text file, one path per line. Use it when the argument list would be too long for the shell.

Duplicate sample names produce:

[E::bcf_sr_add_reader] Error: Duplicate sample names (SAMPLE_A), use --force-samples to proceed

The force-samples flag renames the second occurrence by appending an index to the sample name so the merge can continue. That is a rescue, not a fix. Rename the samples properly with bcftools reheader -s names.txt unless you genuinely intend to merge two libraries from the same individual as separate columns.

A FORMAT column mismatch produces:

[E::vcf_parse_format] Number of columns at chr1:12345 does not match the number of samples (4 vs 3)

This means a data line has the wrong number of tab-separated fields for the declared sample count. Usually it comes from hand-editing a file in a spreadsheet. Validate the offending file with bcftools view on that region and repair the source.

Step 6: Concatenate shards of one sample

When the inputs are pieces of the same sample, concat is the correct command and merge is wrong.

bcftools concat -a -d exact -Oz --threads 8 -o sampleA.genome.vcf.gz \
  chr1.vcf.gz chr2.vcf.gz chr3.vcf.gz
tabix -p vcf sampleA.genome.vcf.gz

The -a flag allows overlaps, which is required before duplicate removal can be requested. The -d exact flag removes records that appear identically in more than one input, which happens when chunk boundaries overlap by one base. Without -a, an overlapping input triggers an error, and without -d you get duplicated records at the seams.

For non-overlapping, already-sorted, already-bgzipped inputs, the -naive option skips validation and runs much faster. Use it only when you are certain about the inputs.

Record count after concatenation should equal the sum of the input record counts, minus any exact duplicates removed.

Step 7: Build a consensus union across callers

For multi-caller reconciliation on structural variants, SURVIVOR is the common choice. Write the VCF paths to a list file, one per line, then merge.

ls caller1.vcf caller2.vcf caller3.vcf > survivor_list.txt
SURVIVOR merge survivor_list.txt 1000 1 1 1 0 50 consensus.vcf

The positional arguments after the list file are: maximum distance in base pairs between calls considered the same event (1000 here), minimum number of callers supporting a record (1 here), minimum size threshold for a call to enter the merge, whether genotypes are taken into account, whether distances are estimated, and the minimum SV size in base pairs (50 here). The final argument is the output path.

Setting the minimum support to 1 produces a union file containing every call from every caller. Raising it to 2 or 3 produces an intersection-leaning file. Neither is more correct. The choice depends on whether your downstream filter can absorb false positives or must avoid false negatives.

The union file alone does not resolve which calls are real. That requires the per-caller quality evidence, which is why meta-callers like SV-MeCa re-extract quality metrics from the original VCFs after the SURVIVOR merge rather than trusting the union file on its own [1].

Step 8: Validate the output

bcftools stats cohort.merged.vcf.gz > cohort.merged.stats
grep '^SN' cohort.merged.stats | head -8

Expected output shape:

SN	0	number of samples:	3
SN	0	number of records:	4620341
SN	0	number of no-ALTs:	0
SN	0	number of SNPs:	4100221
SN	0	number of MNPs:	0
SN	0	number of indels:	520120
SN	0	number of others:	0
SN	0	number of multiallelic sites:	38210

Three checks follow. First, the sample count must equal the number of input samples. Second, the record count must be plausible (see the arithmetic below). Third, the transition-to-transversion ratio reported in the TSTV block should resemble the input files. A large shift points to a problem.

Step 9: Recompute cohort-level annotations

Per-sample INFO fields do not survive a sample merge. Allele count, allele number, and allele frequency must be recomputed across the new sample set.

bcftools +fill-tags cohort.merged.vcf.gz -Oz -o cohort.merged.tags.vcf.gz -- -t AC,AN,AF
tabix -p vcf cohort.merged.tags.vcf.gz

Expected: the AC, AN, and AF INFO fields in the output now reflect all three samples rather than whichever file happened to carry them. Without this step, downstream filters that use AF or AC silently use the wrong values.

Table: Key Options and When to Use Them

Checking That the Merge Is Right: A Worked Example

Suppose sampleA.sorted.vcf.gz holds 4,120,000 records and sampleB.sorted.vcf.gz holds 4,050,000 records, both called over the same whole genome. A merged file containing 8,170,000 records is wrong. That number is the sum, which means the two files were concatenated rather than merged, or that the merge wrote every record twice.

The correct merged count is the size of the union of variant positions. In practice, two well-called human genomes share the large majority of common variant sites, so a merge of two such files typically lands somewhere between the larger input and roughly 1.5 times the larger input, not at the sum. The exact number depends on how many rare, sample-specific variants each file carries.

For a more direct structural check, compare the input and output at a fixed region:

bcftools view -H -r chr1:1000000-1010000 cohort.merged.vcf.gz | head -3

Each line should have three genotype columns, and the ALT alleles should be consistent across samples at shared positions. If a position is present in the inputs but missing from the output, investigate before proceeding.

For set-level comparison between two files you expect to be similar, bcftools isec writes four output files: records private to the first input, records private to the second, records shared, and records in the union. Counts from those files tell you exactly how much the two call sets diverged.

For any comparison of this kind, bcftools stats is the fastest sanity check, and VCFtools remains a usable alternative for validation, merging, and comparison of VCF files [3].

Scaling Up: Hundreds or Thousands of VCFs

A single-machine bcftools merge handles a few hundred files comfortably. Beyond that, throughput becomes the constraint, and the constraint is not where most people assume.

The -threads option in bcftools accelerates compression of the output. It does not parallelize the work done on variant records, which stays single-core. That is the specific problem pbcftools was built to address: it splits the genome into chunks, runs an ordinary bcftools command on each chunk in parallel, and reassembles the outputs. On 1000 Genomes Phase 3 data, merging 100 VCF files ran 19.2 times faster with 32 cores, and parallel output was identical to serial output across every command tested [5].

Two other patterns are worth knowing. First, region-chunked merging followed by concatenation parallelizes cleanly on any cluster scheduler and does not require a special wrapper. Second, at genuinely large scale the bottleneck shifts to input and output, and distributed systems help: custom divide-and-conquer schemas on Hadoop MapReduce, HBase, and Spark have been designed specifically for sorted merging of massive VCF collections, splitting the job into ordered phases so that no stage becomes a bottleneck [4].

There is also the option of not merging at all. If you control the calling pipeline, generating the VCF deterministically in one pass avoids the merge. elPrep 5 executes the full GATK Best Practices variant calling pipeline and produces identical BAM and VCF output to GATK4 while running 8 to 16 times faster on whole-exome and whole-genome data on the same hardware, by parallelizing and merging the execution of pipeline steps [6]. That is a pre-call optimization, not a post-call merge, but for teams regenerating calls frequently it can remove the merge step entirely.

Common Mistakes and Limitations

Merging files from different genome builds

Positions from GRCh37 and GRCh38 are not interchangeable. A merge across builds produces a file that parses cleanly and is biologically meaningless. Lift over first, then merge.

Mixing chr1 and 1

bcftools will happily merge a file using chr1 with a file using 1 into a two-namespace output. Most downstream tools then process only the namespace they recognize and quietly ignore the rest. Standardize contig names before merging.

Using merge where concat belongs

If all inputs come from the same sample, bcftools merge treats them as separate samples and inflates the sample count. The output looks structured but the genotypes are scattered across columns that should not exist. Check the sample count in the output every time.

Treating a merged file as a consensus call set

bcftools merge combines records. It does not adjudicate disagreements. If sampleA calls a variant and sampleB does not, the merged file shows one genotype and one missing genotype. It does not tell you which is right, and it does not re-call anything. A merged file inherits every error from every input. Discordance resolution requires evidence, which is why meta-callers extract per-caller quality metrics into a trained model rather than relying on the union file alone [1].

Skipping left alignment

Equivalent indel representations merge as separate records, which double-counts variants and distorts allele frequency calculations. Normalize every input with bcftools norm -f before merging, every time.

Assuming the threads option scales the whole run

Adding -threads 32 to a merge that is bottlenecked on record processing buys almost nothing. Measure first, then parallelize across regions or use a wrapper that does it for you [5].

Losing the provenance of FILTER and QUAL

After merging, a PASS in the merged file may come from one caller and a failing filter from another. If you discard the original files, you cannot tell which. Keep the pre-merge files, or add a source field to INFO during the merge.

Forgetting that INFO fields collide

Two files that both define an INFO field named DP may mean different things by it. bcftools combines fields it recognizes and warns on conflicts, but semantic conflicts between similarly named fields from different pipelines are invisible. Read the ##INFO header lines of the inputs before merging.

Overlooking privacy when aggregating across institutions

Aggregating VCFs means aggregating identifiable genotype data. Tools have been built specifically to make this safer: GrabBlur aggregates and shares NGS-derived single nucleotide variant data in a public database while keeping individual samples unidentifiable, and it combines phenotype and VCF data through a local interface or command line [7]. If your merge crosses an institutional boundary, address consent and de-identification before you run the command.

Limitation: merging is not joint calling

The strongest multi-sample analyses come from joint genotyping, where all samples are called together so that genotype likelihoods are evaluated across the cohort. Merging post-hoc VCFs cannot recover that information. A merged VCF from single-sample calls is a reasonable approximation for many uses and a poor substitute for others.

Limitation: no tool fixes bad inputs

Every step above assumes reasonably well-formed inputs. A merge of two broken files produces a broken file with more lines.

Frequently Asked Questions

What is the difference between bcftools merge and bcftools concat?

bcftools merge combines different samples into one file with multiple genotype columns, while bcftools concat joins pieces of the same sample into one file with a single genotype column. Use merge when the sample names differ, and concat when the sample names are identical and the files cover different regions.

Can I merge VCF files called against different genome builds?

No. Positions are not comparable between builds, so you must lift one file over to the other build and validate the result before merging.

Why does bcftools refuse to merge my files with a duplicate sample error?

Both files declare the same sample name in their header. Rename one sample with bcftools reheader -s and re-index before merging.

Why is my merged VCF smaller than expected?

Common causes are a wrong genome build causing position mismatches, differing contig names, unnormalized indels collapsing incorrectly, or a region restriction left in place. Compare bcftools stats output between inputs and output to find the discrepancy.

Should I normalize VCF files before merging them?

Yes. Normalization splits multiallelic records and left-aligns indels so that equivalent variants share a representation, which is required for a correct merge.

How many VCF files can bcftools merge at once?

No fixed limit exists in the tool, but single-core record processing becomes the bottleneck at scale. Chunk the genome by region and merge in parallel, or use a parallel wrapper, when the input count reaches the hundreds [5].

What is the difference between FASTQ and FASTA files?

A FASTQ file stores sequencing reads with per-base quality scores, while a FASTA file stores sequences without quality information. VCF files are produced after reads are aligned and variants are called, so they are downstream of both.

How do I merge structural variant calls from multiple callers?

Use a consensus merger such as SURVIVOR to build a union file, then bring in per-caller quality evidence to rank the consensus calls rather than trusting the union alone [1].

Sources

1. SV-MeCa: an XGBoost-based meta-caller approach for structural variant calling from short-read data.
2. ShinyVar: a web-based application for comparative Influenza variant analysis supporting structure-guided approaches to vaccine and antiviral drug design.
3. The variant call format and VCFtools.
4. Optimized distributed systems achieve significant performance improvement on sorted merging of massive VCF files.
5. pbcftools: parallel execution of bcftools for large variant call sets.
6. Multithreaded variant calling in elPrep 5.
7. GrabBlur--a framework to facilitate the secure exchange of whole-exome and -genome SNV data using VCF files.