domains / science / vcf-inspection / vcf-basics

vcf-basics

Skill card

combined

The two copies cover the same core task (viewing and querying VCF/BCF files with bcftools and cyvcf2) but the GPTomics copy adds substantial, non-trivial interpretive constraints that the FreedomIntelligence copy lacks: a governing principle about field level and Number semantics, a detailed QUAL vs GQ vs PL/GL comparison, the header contract, and specific traps (missing vs hom-ref, symbolic alleles, gVCF). The FreedomIntelligence copy provides a more complete, step-by-step command reference for bcftools view/query, format conversion, indexing, and cyvcf2 code patterns that the GPTomics copy truncates or omits. Crowning one would hide the other's real value.

Repositories

Other copies

Combined skill

---
name: bio-vcf-basics
description: View, query, and interpret VCF/BCF variant files with bcftools and cyvcf2. Use when inspecting variants, extracting fields, converting formats, or correctly interpreting field semantics (QUAL vs GQ vs PL, AD vs DP, GT phasing, Number A/R/G, symbolic alleles, gVCF).
tool_type: cli
primary_tool: bcftools
---

# VCF/BCF Basics

## Version Compatibility

Reference examples tested with: bcftools 1.19+, cyvcf2 0.30+, numpy 1.26+.

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

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

## The Governing Principle

A VCF field is only meaningful once its LEVEL and its `Number` are known. QUAL is a site-level property; GQ, PL, AD, DP, GT are per-sample. QUAL and GQ answer different questions and are NOT interchangeable. A field's header `Number` (A/R/G/.) dictates how many values it carries and how it must be re-subset after a multiallelic split. Several encodings are load-bearing traps: `.` (missing) is never `0/0` (hom-ref); a bare `*` ALT is a spanning-deletion placeholder, not an allele; a gVCF `<NON_REF>` record is a reference-confidence intermediate, not a filtered call. Read the header, read the Number, read the level -- a structurally valid VCF read at the wrong level silently produces wrong numbers with no error.

## Format Overview

| Format | Description | Use Case |
|--------|-------------|----------|
| VCF | Text format, human-readable | Debugging, small files |
| VCF.gz | Compressed VCF (bgzip) | Standard distribution |
| BCF | Binary VCF | Fast processing, large files |

## VCF Format Structure

```
##fileformat=VCFv4.2
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SAMPLE1
chr1    1000    rs123   A       G       30      PASS    DP=50   GT:DP   0/1:25
```

### Header Lines (##)
- `##fileformat` - VCF version
- `##INFO` / `##FORMAT` - INFO / FORMAT field definitions (ID, Number, Type)
- `##FILTER` - Filter definitions
- `##contig` - Reference contigs (required for region indexing and contig order)
- `##reference` - Reference genome

### The Header Contract

Every INFO/FORMAT tag used in the body MUST be declared in a `##INFO`/`##FORMAT` line giving its ID, Number, and Type; parsers (bcftools, cyvcf2, pysam) read these declarations to know how many values a field holds and how to type it. An out-of-sync header -- a tag used but not declared, or declared with the wrong Number/Type -- silently breaks parsing: a `Number=1` declaration over data that holds a vector, or a missing `##contig`, makes tools drop, mistype, or mis-subset values with NO error thrown. After any hand-edit or annotation that adds a field, update the header to match (`bcftools +fill-tags` and `bcftools annotate` manage this automatically).

### Data Columns

| Column | Description |
|--------|-------------|
| CHROM | Chromosome |
| POS | 1-based position of the first base in REF (contrast BED's 0-based half-open) |
| ID | Variant identifier (e.g., rs number) or `.` if novel |
| REF | Reference allele (matches the reference exactly) |
| ALT | Alternate allele(s), comma-separated. `*` = allele missing due to an overlapping deletion at this site |
| QUAL | Phred-scaled quality of the ALT assertion, `-10*log10 P(no variant)`; higher = more confident a variant exists (site-level, NOT per-sample) |
| FILTER | PASS or semicolon-separated filter names. `.` means filters were not applied |
| INFO | Semicolon-separated key=value pairs (site-level annotations) |
| FORMAT | Colon-separated format keys defining per-sample field order |
| SAMPLE | Colon-separated values matching FORMAT order |

## Critical Field Interpretation

### QUAL vs GQ vs PL/GL: three different confidences

| Field | Level | Scale | Question answered |
|-------|-------|-------|-------------------|
| QUAL (col 6) | Site | Phred: `-10*log10 P(no variant)` | "Is there ANY variant at this site?" |
| GQ (FORMAT) | Genotype | Phred, capped at 99 | "Is THIS sample's assigned genotype correct?" |
| PL (FORMAT) | Genotype | Phred, rebased to min=0 | Relative likelihood of every possible genotype |
| GL (FORMAT) | Genotype | log10, `<=0`, raw | Same info as PL, unscaled (`PL = -10*GL`, rebased) |

QUAL is computed once across all samples and SCALES with total depth, so a high-coverage artifact can carry a large QUAL -- hence QD (QUAL normalized by depth) is preferred for filtering. GQ is per-sample and does not scale with cohort size. They are NOT interchangeable: QUAL can be high while an individual genotype is uncertain (low GQ), and a sample can have a confident genotype (high GQ) at a site with only moderate QUAL. Filter site-level junk on QUAL/QD; no-call untrustworthy genotypes on GQ.

### PL/GL, and how GQ is derived

PL holds phred-scaled genotype likelihoods, rebased so the CALLED (most likely) genotype is exactly 0 and every other value is its phred penalty relative to that called genotype. GQ is derived from PL as `GQ = 99 - min(99, PL[1] - PL[0])` (or the second-smallest PL minus the smallest, capped at 99).

### AD vs DP and Allele Balance

- `DP` is total read depth at the site (all alleles).
- `AD` is allele-specific depth (e.g., `AD=25,5` means 25 reads support REF, 5 support ALT).
- Allele balance = AD[ALT] / DP. Extreme imbalance (e.g., 1/30) may indicate a technical artifact or a true low-frequency variant.

### GT Phasing, Ploidy, and PS

- `0/0` = Homozygous reference
- `0/1` = Heterozygous (unphased)
- `1/1` = Homozygous alternate
- `1/2` = Heterozygous (two different alts, unphased)
- `./.` = Missing (NOT the same as `0/0`)
- `0|1` = Phased heterozygous (alleles on known haplotypes)
- `PS` (Phase Set) in FORMAT groups genotypes into haplotypes; the same PS value means the genotypes are phased together.

### Symbolic Alleles and Spanning Deletions

- `<DEL>`, `<INS>`, `<DUP>`, `<INV>` are symbolic alleles; their actual sequence is in the `END` INFO field (for deletions) or the ALT string.
- A bare `*` in the ALT column is a placeholder for a spanning deletion; the actual deletion is described by the `END` field.
- `<NON_REF>` in a gVCF file is a reference-confidence record, not a variant call.

## bcftools view

### View VCF
```bash
bcftools view input.vcf.gz | head
```

### View Header Only
```bash
bcftools view -h input.vcf.gz
```

### View Without Header
```bash
bcftools view -H input.vcf.gz | head
```

### View Specific Region
```bash
bcftools view input.vcf.gz chr1:1000000-2000000
```

### View Specific Samples
```bash
bcftools view -s sample1,sample2 input.vcf.gz
```

### Exclude Samples
```bash
bcftools view -s ^sample3 input.vcf.gz
```

## bcftools query

Extract specific fields in custom format.

### Basic Query
```bash
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' input.vcf.gz
```

### Query with INFO Fields
```bash
bcftools query -f '%CHROM\t%POS\t%INFO/DP\t%INFO/AF\n' input.vcf.gz
```

### Query with Sample Fields
```bash
bcftools query -f '%CHROM\t%POS[\t%GT]\n' input.vcf.gz
```

### Query Specific Samples
```bash
bcftools query -f '%CHROM\t%POS[\t%SAMPLE=%GT]\n' -s sample1,sample2 input.vcf.gz
```

### Include Header
```bash
bcftools query -H -f '%CHROM\t%POS\t%REF\t%ALT\n' input.vcf.gz
```

### Common Format Specifiers

| Specifier | Description |
|-----------|-------------|
| `%CHROM` | Chromosome |
| `%POS` | Position |
| `%ID` | Variant ID |
| `%REF` | Reference allele |
| `%ALT` | Alternate allele |
| `%QUAL` | Quality score |
| `%FILTER` | Filter status |
| `%INFO/TAG` | INFO field value |
| `%TYPE` | Variant type (snp, indel, etc.) |
| `[%GT]` | Genotype (per sample) |
| `[%DP]` | Depth (per sample) |
| `[%SAMPLE]` | Sample name |
| `\n` | Newline |
| `\t` | Tab |

## Format Conversion

### VCF to BCF
```bash
bcftools view -Ob -o output.bcf input.vcf.gz
```

### BCF to VCF
```bash
bcftools view -Ov -o output.vcf input.bcf
```

### Compress VCF (bgzip)
```bash
bgzip input.vcf
# Creates input.vcf.gz
```

### Index VCF/BCF
```bash
bcftools index input.vcf.gz
# Creates input.vcf.gz.csi

bcftools index -t input.vcf.gz
# Creates input.vcf.gz.tbi (tabix index)
```

## Output Format Options

| Flag | Format |
|------|--------|
| `-Ov` | Uncompressed VCF |
| `-Oz` | Compressed VCF (bgzip) |
| `-Ou` | Uncompressed BCF |
| `-Ob` | Compressed BCF |

## Genotype Encoding

| Genotype | Meaning |
|----------|---------|
| `0/0` | Homozygous reference |
| `0/1` | Heterozygous |
| `1/1` | Homozygous alternate |
| `1/2` | Heterozygous (two different alts) |
| `./.` | Missing |
| `0\|1` | Phased heterozygous |

## cyvcf2 Python Alternative

### Open and Iterate
```python
from cyvcf2 import VCF

vcf = VCF('input.vcf.gz')
for variant in vcf:
    print(f'{variant.CHROM}:{variant.POS} {variant.REF}>{variant.ALT[0]}')
```

### Access Variant Properties
```python
from cyvcf2 import VCF

for variant in VCF('input.vcf.gz'):
    print(f'Chrom: {variant.CHROM}')
    print(f'Pos: {variant.POS}')
    print(f'ID: {variant.ID}')
    print(f'Ref: {variant.REF}')
    print(f'Alt: {variant.ALT}')  # List
    print(f'Qual: {variant.QUAL}')
    print(f'Filter: {variant.FILTER}')
    print(f'Type: {variant.var_type}')  # snp, indel, etc.
    break
```

### Access INFO Fields
```python
from cyvcf2 import VCF

for variant in VCF('input.vcf.gz'):
    dp = variant.INFO.get('DP')
    af = variant.INFO.get('AF')
    print(f'{variant.CHROM}:{variant.POS} DP={dp} AF={af}')
```

### Access Genotypes
```python
from cyvcf2 import VCF

vcf = VCF('input.vcf.gz')
samples = vcf.samples  # List of sample names

for variant in vcf:
    gts = variant.gt_types  # 0=HOM_REF, 1=HET, 2=UNKNOWN, 3=HOM_ALT
    for sample, gt in zip(samples, gts):
        gt_str = ['HOM_REF', 'HET', 'UNKNOWN', 'HOM_ALT'][gt]
        print(f'{sample}: {gt_str}')
    break
```

### Access Sample Fields
```python
from cyvcf2 import VCF

for variant in VCF('input.vcf.gz'):
    depths = variant.samples[0].DP  # Example: first sample's DP
    print(f'{variant.CHROM}:{variant.POS} DP={depths}')
```