/bcftools-variant-manipulation

name: "bcftools-variant-manipulation"
description: "CLI for VCF/BCF: filter, merge, annotate, query, normalize, compute stats. Core post-variant-calling: quality filtering, multi-sample merging, rsID annotation, genotype extraction. Samtools companion in HTSlib. Use GATK for complex indel realignment during calling; use VCFtools for population genetics stats."
license: "MIT"

bcftools — VCF/BCF Variant Manipulation Toolkit

Overview

bcftools is the standard command-line toolkit for processing VCF (Variant Call Format) and BCF (Binary Call Format) files in the HTSlib ecosystem. It covers the complete post-variant-calling workflow: format conversion, quality filtering, variant normalization, multi-sample merging, annotation with external databases, genotype extraction, and QC statistics. bcftools uses streaming by design — most commands read from stdin and write to stdout, making it ideal for memory-efficient pipelines on large cohorts.

When to Use

• Filtering variants by quality (QUAL, DP, AF) after variant calling
• Merging VCF files from multiple samples into a joint call set
• Adding rsIDs or gene annotations to variant calls
• Extracting specific fields (genotypes, allele depths) as tabular output
• Normalizing indel representations and splitting multi-allelic records
• Calling variants from pileup output (mpileup + call)
• Computing per-sample and overall VCF QC statistics
• Use `GATK HaplotypeCaller` instead when calling variants with local realignment in human samples
• Use `VCFtools` instead for population genetics statistics (Fst, LD, Hardy-Weinberg)
• Use `bcftools` in the HTSlib pipeline; use `picard` for duplicate-marking and library metrics

Prerequisites

• **Installation**: bcftools 1.17+ (part of HTSlib suite with samtools)
• **Input requirements**: VCF or BGzipped+tabix-indexed VCF (`.vcf.gz + .vcf.gz.tbi`) for region queries
• **Companion tools**: `samtools` for BAM processing; `tabix` for VCF indexing

> **Check before installing**: The tool may already be available in the current environment (e.g., inside a `pixi` / `conda` env). Run `command -v bcftools` first and skip the install commands below if it returns a path. When running inside a pixi project, invoke the tool via `pixi run bcftools` rather than bare `bcftools`.

# Bioconda (recommended — installs HTSlib suite)
conda install -c bioconda bcftools

# Homebrew (macOS)
brew install bcftools

# Verify
bcftools --version | head -1
# bcftools 1.20

# Index a VCF for region queries
bcftools index -t variants.vcf.gz   # creates .tbi
bcftools index -c variants.vcf.gz   # creates .csi (for chromosomes > 512 Mb)

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: filterExpression
    kind: required
    source: user
    ask: "Which variants should survive - by depth, genotype quality, allele balance, or population frequency?"
    default: null

  - id: D2
    param: multiallelicHandling
    kind: required
    source: user
    ask: "Should multi-allelic records be split into one row per alternate allele before filtering and annotation?"
    default: "split - most annotation and comparison tools assume one alternate per row"

  - id: D3
    param: leftAlignment
    kind: required
    source: upstream
    depends_on: [D2]
    ask: "Should indels be left-aligned against the reference so the same variant is written identically across callers?"
    default: "left-aligned when a reference FASTA is available"

  - id: D4
    param: sampleSubset
    kind: optional
    source: user
    ask: "Should the output be restricted to particular samples?"
    default: "all samples"

  - id: D5
    param: regionSubset
    kind: optional
    source: user
    ask: "Should records be restricted to particular regions?"
    default: "whole file"

  - id: D6
    param: outputFormat
    kind: optional
    source: upstream
    ask: "Which container should the result be written in?"
    default: "bgzip-compressed VCF with an index"

  - id: D7
    param: threads
    kind: never_ask
    source: data
    reason: "Compression threads affect runtime only, not the records"
    default: "min(4, available_cores)"

D2 and D3 together decide whether two VCFs of the same sample can be compared at all. Unsplit multi-allelic rows and right-aligned indels produce apparent private variants that are the same variant written differently - and every downstream intersection then reports a difference that does not exist.

Quick Start

# Typical post-calling workflow: normalize → filter → annotate → extract
bcftools norm -d any -f reference.fa variants.vcf.gz \
  | bcftools filter -i 'QUAL>20 && DP>10' \
  | bcftools annotate -a dbSNP.vcf.gz -c ID \
  | bcftools view -O z -o final.vcf.gz

# Index the output
bcftools index -t final.vcf.gz

# Count variants at each stage
bcftools stats final.vcf.gz | grep "^SN"

Core API

Module 1: VCF/BCF I/O and Format Conversion

Convert between text VCF and binary BCF; compress and index for random access.

# VCF → compressed BCF (fastest format for piping)
bcftools view -O b -o variants.bcf variants.vcf

# BCF → VCF (for human-readable output)
bcftools view -O v -o variants.vcf variants.bcf

# VCF → bgzipped + indexed (standard archive format)
bcftools view -O z -W -o variants.vcf.gz variants.vcf
# -W automatically creates .tbi index after writing

# Extract specific samples
bcftools view -s sample1,sample2 -O z -o subset.vcf.gz variants.vcf.gz

# Exclude samples (prefix with ^)
bcftools view -s ^outlier_sample -O z -o cleaned.vcf.gz variants.vcf.gz

# Extract by region (fast; requires index)
bcftools view -r chr1:1000000-2000000 variants.vcf.gz -O v -o chr1_region.vcf

# Streaming pipeline: no intermediate files
samtools mpileup -Ou input.bam | bcftools call -m -Oz -o calls.vcf.gz

Module 2: Variant Filtering

Apply quality thresholds and FLAG-based filters to retain high