← Về thư mục
📄 / / profiles / r-and-d / raw_evidence / OpenAI4S / skills / bioskills / bio-variant-calling-vcf-basics / SKILL.md

name: bio-vcf-basics description: View, query, and interpret VCF/BCF variant files with bcftools and cyvcf2. Use when inspecting variants, extracting fields with query format strings, converting VCF/BCF, or correctly reading a field -- QUAL (site) vs GQ (genotype) vs PL/GL likelihoods, AD vs DP and allele balance, GT phasing/ploidy/PS and missing-vs-hom-ref, INFO/FORMAT Number A/R/G semantics, symbolic alleles (, , spanning *) and END, or telling a raw gVCF apart from a filtered callset. origin: openai4s category: bioskills/variant-calling metadata: tool_type: cli primary_tool: bcftools third_party: name: GPTomics/bioSkills repository: https://github.com/GPTomics/bioSkills commit: d91ed3d563019e649dc854c56ccd62551359488a license: MIT

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.

VCF/BCF Basics

"Show me and extract fields from this VCF" -> Parse the VCF/BCF format, then view, subset, or pull specific columns into a flat table. - CLI: bcftools view / bcftools query -f - Python: cyvcf2.VCF (iterate records with attribute access)

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. And 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 (##)

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

What each field actually measures -- and what it does not -- drives every filtering and interpretation decision. QUAL, GQ, and PL answer three DIFFERENT questions.

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 call. For a biallelic diploid site PL is ordered [PL(0/0), PL(0/1), PL(1/1)] -- the index of the 0 IS the genotype the caller assigned. GL is the same information as raw log10 likelihoods (<=0, larger is better). GQ = the difference between the two SMALLEST PL values, i.e. phred confidence in the call versus the next-best genotype; GQ=0 means the top two genotypes are tied (uninformative), GQ is capped at 99 by convention.

For a site with n alleles, diploid genotype j/k (j<=k) sits at PL index k*(k+1)/2 + j (this is the Number=G ordering). Getting this index formula wrong is the classic bug when re-parsing PL after a multiallelic split -- the vector must be re-subset by the formula, never sliced positionally.

AD vs DP, and allele balance

AD (FORMAT, Number=R) is per-allele read depth [ref_depth, alt1_depth, ...], REF first. DP is total depth. sum(AD) is often LESS than DP -- expected, not an error: - DP counts all reads spanning the position, including uninformative reads (low base quality, ambiguous alignment, filtered reads). - AD counts only reads that confidently support a specific allele. - INFO/DP (site-level, summed across samples) differs from FORMAT/DP (per-sample).

Allele balance for a het is DERIVED from AD (GATK does not emit it directly): AB = alt_AD / (ref_AD + alt_AD). A true het sits near 0.5; hets far from 0.5 (e.g. <0.2 or >0.8) suggest a mapping artifact, CNV, or contamination.

INFO/FORMAT Number semantics (A / R / G / .)

Every ##INFO/##FORMAT header declares a Number telling a parser how many values a field holds AND how to re-subset it when a multiallelic record is split:

Number One value per Examples On multiallelic split
A ALT allele AF, AC take the k-th element for the k-th ALT
R allele incl. REF AD REF value first, then per-ALT (off-by-one vs A)
G genotype PL, GL re-subset via the k*(k+1)/2+j index formula
. variable/unknown -- parser CANNOT auto-subset; carried whole onto every record
0 flag (presence only) -- --

Load-bearing for correctness: bcftools norm -m- uses these codes to reapportion fields on split. A field mis-declared Number=. when it is really A keeps its full multiallelic vector on every split record, so downstream tools read the WRONG allele's value with no error. See variant-calling/variant-normalization for split/join reapportionment.

Key INFO Annotations for Filtering

Annotation Meaning What It Detects
QD QUAL / allele depth Low values suggest variant quality not supported by reads
FS Fisher strand bias (phred-scaled) Variant reads predominantly on one strand (artifact)
SOR Strand odds ratio Same as FS but handles high-depth sites better
MQ Root mean square mapping quality Low values indicate reads map ambiguously (paralogous regions)
MQRankSum MQ difference: ref vs alt reads Very negative = alt reads map much worse than ref (suspicious)
ReadPosRankSum Read position: ref vs alt reads Very negative = variant only at read ends (misalignment artifact)

Genotype Encoding

Genotype Meaning
0/0 Homozygous reference (confidently called ref)
0/1 Heterozygous
1/1 Homozygous alternate
1/2 Heterozygous for two different ALT alleles (compound het at multiallelic site)
./. Missing genotype (no confident call)
0\|1 Phased heterozygous (allele before \| is on haplotype 1)

Phased vs Unphased

Phase Sets (PS)

A | is only meaningful WITHIN a phase set. The FORMAT/PS tag (an integer, usually the POS of the block's first variant) groups variants phased relative to EACH OTHER; 0|1 in two different PS blocks are not guaranteed to lie on the same physical haplotype. Read-backed phasers (WhatsHap) and trio phasing emit PS. A | with no consistent PS across records carries no global phase -- a subtle trap when merging phased VCFs.

Ploidy and Missing vs Hom-Ref

Ploidy is read from the NUMBER of alleles in GT: 0/1 diploid, 0 haploid (chrY, chrM, male chrX outside the PAR), 0/1/1 triploid. Per-region ploidy (PAR, chrX in males, mito) must match the sample karyotype.

. (missing) is NOT reference. ./. = no-call (genotype could not be determined, usually low depth); 0/0 = confidently called homozygous reference. Treating ./. as 0/0 inflates the reference-allele count and biases allele frequencies, missingness, and burden tests. This is load-bearing: never impute ./. as reference. In a gVCF, the ABSENCE of a record also does not mean reference -- see the reference-confidence model below.

Multiallelic Genotypes

At multiallelic sites (e.g. ALT = G,T), allele indices reference the comma-separated ALT list: 0=REF, 1=first ALT, 2=second ALT. 1/2 means one copy of each ALT. Splitting multiallelics into biallelic records with bcftools norm -m- converts 1/2 into two 0/1 records, losing compound-heterozygosity information -- see variant-calling/variant-normalization for caveats.

Symbolic Alleles, END, and Spanning Deletions

Not every ALT spells out a sequence. Symbolic alleles are angle-bracketed placeholders for events whose sequence is not given inline:

ALT Meaning
<DEL> <DUP> <INS> <INV> <CNV> Structural-variant classes (sequence not spelled out)
<NON_REF> gVCF: "any allele not yet observed" (reference-confidence model)
<*> Same role as <NON_REF> in some callers' gVCF/mpileup output
* (bare) Spanning deletion: allele MISSING because an upstream deletion on ANOTHER line overlaps this position

Two parsing traps: - INFO/END gives the end coordinate of a symbolic/large event. A tool that infers a record's span from len(REF) is WRONG for symbolic alleles -- it must read END. gVCF reference blocks also use END to mark the last position of the band. - The bare * ALT is interpretable only relative to the overlapping deletion on another record; it is not a real alternate allele here. Splitting/subsetting can strand a * from the deletion it references (see variant-calling/variant-normalization).

gVCF and the <NON_REF> Reference-Confidence Model

A gVCF (GATK HaplotypeCaller -ERC GVCF) is fundamentally different from a filtered callset: it emits a record for EVERY position or block, not just variant sites. Non-variant stretches are compressed into END-delimited blocks (bands) grouped by GQ, so a gVCF is not one line per base.

bcftools view

Goal: View, subset, and convert VCF/BCF files from the command line.

Approach: Use bcftools view with flags for header control, region selection, sample extraction, and format conversion.

bcftools view input.vcf.gz | head           # full records
bcftools view -h input.vcf.gz               # header only
bcftools view -H input.vcf.gz | head        # skip header
bcftools view input.vcf.gz chr1:1000000-2000000   # region (needs index)
bcftools view -s sample1,sample2 input.vcf.gz      # keep samples
bcftools view -s ^sample3 input.vcf.gz             # exclude samples

bcftools query

Goal: Extract specific fields from a VCF in a custom tabular format.

Approach: Use bcftools query -f with format specifiers for CHROM, POS, INFO, and FORMAT fields. Square brackets [...] loop over samples.

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\n' input.vcf.gz
bcftools query -f '%CHROM\t%POS\t%INFO/DP\t%INFO/AF\n' input.vcf.gz
bcftools query -f '%CHROM\t%POS[\t%GT]\n' input.vcf.gz              # per-sample GT
bcftools query -f '%CHROM\t%POS[\t%SAMPLE=%GT]\n' -s sample1 input.vcf.gz
bcftools query -H -f '%CHROM\t%POS\t%REF\t%ALT\n' input.vcf.gz     # column header

Common Format Specifiers

Specifier Description
%CHROM %POS %ID Position fields
%REF %ALT Alleles
%QUAL %FILTER Site quality / filter status
%INFO/TAG INFO field value
%TYPE Variant type (snp, indel, etc.)
[%GT] [%DP] [%AD] [%GQ] Per-sample FORMAT fields (loop in [...])
[%SAMPLE] Sample name
\n \t Newline / tab

Format Conversion and Indexing

Goal: Convert between VCF, compressed VCF, and BCF, and index for region queries.

Approach: Use bcftools view output flags (-Ov/-Oz/-Ou/-Ob), then bgzip + index.

bcftools view -Ob -o output.bcf input.vcf.gz   # VCF -> BCF
bcftools view -Ov -o output.vcf input.bcf      # BCF -> VCF
bgzip input.vcf                                 # -> input.vcf.gz (bgzip, NOT gzip)
bcftools index input.vcf.gz                     # -> .csi index
bcftools index -t input.vcf.gz                  # -> .tbi (tabix) index

Output Format Options

Flag Format
-Ov Uncompressed VCF
-Oz Compressed VCF (bgzip)
-Ou Uncompressed BCF (fast piping)
-Ob Compressed BCF

BCF is the binary encoding of VCF: faster to parse and smaller for large callsets. Region queries (chr1:1-1000) require a bgzipped+indexed VCF or a BCF -- plain .gz (gzip) is not seekable and fails.

cyvcf2 Python Alternative

Goal: Read, query, and write VCF files programmatically in Python.

Approach: Use cyvcf2's VCF reader to iterate variants with attribute access to fields, and Writer to emit filtered output.

"Parse this VCF in Python" -> Open with cyvcf2 and iterate variant records.

Open, Iterate, and Access Fields

from cyvcf2 import VCF

vcf = VCF('input.vcf.gz')
for variant in vcf:
    # ALT is a list; QUAL is site-level and may be None
    print(variant.CHROM, variant.POS, variant.REF, variant.ALT)
    print(variant.ID, variant.QUAL, variant.FILTER, variant.var_type)
    dp = variant.INFO.get('DP')   # INFO field, None if absent
    af = variant.INFO.get('AF')
    break
vcf.close()

Access Genotypes and Per-Sample Fields

from cyvcf2 import VCF

vcf = VCF('input.vcf.gz')
samples = vcf.samples
for variant in vcf:
    # gt_types: 0=HOM_REF, 1=HET, 2=UNKNOWN(missing), 3=HOM_ALT
    gts = variant.gt_types
    depths = variant.format('DP')   # numpy array, one row per sample
    gqs = variant.format('GQ')      # per-sample genotype quality
    ad = variant.format('AD')       # per-allele depth, Number=R
    print(dict(zip(samples, gts)))
    break
vcf.close()

Note: cyvcf2 codes missing genotypes as gt_types == 2 (UNKNOWN) -- treat that as no-call, never as HOM_REF.

Fetch Region and Read the Header

from cyvcf2 import VCF

vcf = VCF('input.vcf.gz')
print(vcf.samples, vcf.seqnames)      # sample names, contig names
for info in vcf.header_iter():
    if info['HeaderType'] == 'INFO':
        print(info['ID'], info['Description'])
for variant in vcf('chr1:1000000-2000000'):   # requires an index
    print(variant.CHROM, variant.POS)

Write Filtered VCF

from cyvcf2 import VCF, Writer

vcf = VCF('input.vcf.gz')
writer = Writer('output.vcf', vcf)   # inherit the input header
for variant in vcf:
    if variant.QUAL is not None and variant.QUAL > 30:   # QUAL is site-level
        writer.write_record(variant)
writer.close()
vcf.close()

Quick Reference

Task bcftools cyvcf2
View VCF bcftools view file.vcf.gz VCF('file.vcf.gz')
View header bcftools view -h file.vcf.gz vcf.header_iter()
Get region bcftools view file.vcf.gz chr1:1-1000 vcf('chr1:1-1000')
Query fields bcftools query -f '%CHROM\t%POS\n' Loop with properties
Count variants bcftools view -H file.vcf.gz \| wc -l sum(1 for _ in vcf)
VCF to BCF bcftools view -Ob -o out.bcf in.vcf.gz Use Writer

Common Errors

Error Cause Solution
no BGZF EOF marker Not bgzipped (plain gzip) Recompress with bgzip, not gzip
index required / region query fails Missing index Run bcftools index (-t for tabix)
sample not found Wrong sample name Check with bcftools query -l
INFO/FORMAT field missing or mistyped Header out of sync with body Fix ##INFO/##FORMAT Number/Type; use bcftools +fill-tags
Every hom-alt or missing site vanishes on filter Treated ././. as failing or as ref Missing != hom-ref; make missing pass, never impute 0/0
Wrong allele's AF/AD after split Number=. field not re-subset Declare the true Number (A/R/G) so bcftools reapportions

References