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)

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.