← Về thư mục
📄 / / profiles / r-and-d / raw_evidence / OpenAI4S / skills / bioskills / bio-variant-calling-variant-annotation / references / usage-guide.md

Variant Annotation Usage Guide

Overview

Annotation attaches functional consequences, population frequencies, and pathogenicity evidence to variants. The organizing idea is that a consequence is not a property of the variant but of the tuple (variant, transcript model, engine, engine version, parameters): VEP, SnpEff, and ANNOVAR disagree, and so do RefSeq, Ensembl, and MANE. The job is not to find the one right tool but to PIN every axis (build, transcript set, engine+version, predictor, gnomAD version) and record it. This guide covers the bcftools annotate/csq mechanics plus the transcript, predictor, and frequency decisions that determine what an annotation means.

Prerequisites

Quick Start

Tell your AI agent what you want to do: - "Annotate my VCF with VEP on the MANE Select transcript" - "Add rsIDs from dbSNP and gnomAD frequencies to my VCF" - "Report the grpmax filtering allele frequency, not just global AF" - "Predict splice effects with SpliceAI delta scores" - "Add one calibrated missense predictor (REVEL) for PP3/BP4" - "Predict functional consequences of my variants using bcftools csq" - "Explain why VEP and SnpEff give different consequences for the same variant"

bcftools annotate Overview

bcftools annotate adds or removes annotations from VCF files. Common uses:

Adding Annotations from VCF

Add rsIDs from dbSNP

bcftools annotate -a dbsnp.vcf.gz -c ID input.vcf.gz -Oz -o with_rsids.vcf.gz
bcftools index with_rsids.vcf.gz

The -c ID tells bcftools to copy the ID column from the annotation file.

Add INFO Fields

# Add allele frequency
bcftools annotate -a gnomad.vcf.gz -c INFO/AF input.vcf.gz -Oz -o with_af.vcf.gz

# Add multiple fields
bcftools annotate -a gnomad.vcf.gz -c INFO/AF,INFO/AC,INFO/AN input.vcf.gz -Oz -o annotated.vcf.gz

Add All INFO Fields

bcftools annotate -a database.vcf.gz -c INFO input.vcf.gz -Oz -o annotated.vcf.gz

Column Specification Options

Option Description
ID Variant identifier
QUAL Quality score
FILTER Filter status
INFO All INFO fields
INFO/TAG Specific INFO field
+INFO/TAG Append to existing field
-INFO/TAG Remove field

Adding Annotations from BED/TAB Files

BED File Annotation

# regions.bed.gz has 4 columns: CHROM, START, END, REGION_NAME

# First add header definition
echo '##INFO=<ID=REGION,Number=1,Type=String,Description="Genomic region">' > header.txt

bcftools annotate \
    -a regions.bed.gz \
    -c CHROM,FROM,TO,INFO/REGION \
    -h header.txt \
    input.vcf.gz -Oz -o annotated.vcf.gz

TAB File Annotation

# scores.tab.gz has columns: CHROM, POS, SCORE

echo '##INFO=<ID=SCORE,Number=1,Type=Float,Description="Conservation score">' > header.txt

bcftools annotate \
    -a scores.tab.gz \
    -c CHROM,POS,INFO/SCORE \
    -h header.txt \
    input.vcf.gz -Oz -o annotated.vcf.gz

Removing Annotations

Remove Specific Fields

# Remove INFO fields
bcftools annotate -x INFO/DP,INFO/MQ input.vcf.gz -Oz -o clean.vcf.gz

# Remove FORMAT fields
bcftools annotate -x FORMAT/AD,FORMAT/PL input.vcf.gz -Oz -o clean.vcf.gz

# Remove ID column
bcftools annotate -x ID input.vcf.gz -Oz -o no_ids.vcf.gz

Remove All of a Type

# Remove all INFO fields
bcftools annotate -x INFO input.vcf.gz -Oz -o minimal.vcf.gz

# Remove all FORMAT fields except GT
bcftools annotate -x FORMAT input.vcf.gz -Oz -o gt_only.vcf.gz

Keep Only Specific Fields

# Remove all INFO except DP and AF
bcftools annotate -x ^INFO/DP,INFO/AF input.vcf.gz -Oz -o minimal.vcf.gz

Setting Variant IDs

Create IDs from Fields

bcftools annotate --set-id '%CHROM\_%POS\_%REF\_%ALT' input.vcf.gz -Oz -o with_ids.vcf.gz

Example output: chr1_12345_A_G

Only Set Missing IDs

# '+' prefix means only set if ID is currently '.'
bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%ALT' input.vcf.gz -Oz -o with_ids.vcf.gz

ID Format Specifiers

Specifier Value
%CHROM Chromosome
%POS Position
%REF Reference allele
%ALT Alternate allele
%FIRST_ALT First alternate allele

Chromosome Renaming

UCSC to Ensembl Style

# rename.txt: old_name<tab>new_name
cat > rename.txt << EOF
chr1    1
chr2    2
chr3    3
chrX    X
chrY    Y
chrM    MT
EOF

bcftools annotate --rename-chrs rename.txt input.vcf.gz -Oz -o renamed.vcf.gz

Ensembl to UCSC Style

cat > rename.txt << EOF
1   chr1
2   chr2
3   chr3
X   chrX
Y   chrY
MT  chrM
EOF

bcftools annotate --rename-chrs rename.txt input.vcf.gz -Oz -o renamed.vcf.gz

Consequence Prediction with bcftools csq

Requirements

  1. Reference FASTA
  2. GFF3 gene annotation (from Ensembl, NCBI, etc.)

Basic Usage

bcftools csq -f reference.fa -g genes.gff3.gz input.vcf.gz -Oz -o consequences.vcf.gz

Output Format

Adds BCSQ INFO field:

BCSQ=missense|GENE1|ENST00000123456|protein_coding|+|123A>V|c.368C>T

Fields are pipe-separated: 1. Consequence type 2. Gene name 3. Transcript ID 4. Biotype 5. Strand 6. Amino acid change 7. DNA change

Common Consequences

Consequence Impact
synonymous Low - same amino acid
missense Moderate - different amino acid
stop_gained High - premature stop
frameshift High - reading frame shift
splice_donor High - affects splicing
splice_acceptor High - affects splicing
intron Low - intronic
intergenic Low - between genes

Query Consequences

# List all consequences
bcftools query -f '%CHROM\t%POS\t%INFO/BCSQ\n' consequences.vcf.gz | head

# Extract just consequence and gene
bcftools query -f '%INFO/BCSQ\n' consequences.vcf.gz | cut -d'|' -f1,2 | sort | uniq -c

Filter by Consequence

# Keep high-impact variants
bcftools view -i 'INFO/BCSQ~"stop_gained" || INFO/BCSQ~"frameshift"' \
    consequences.vcf.gz -Oz -o high_impact.vcf.gz

Annotation Database Setup

dbSNP

# Download
wget https://ftp.ncbi.nih.gov/snp/latest_release/VCF/GCF_000001405.40.gz
wget https://ftp.ncbi.nih.gov/snp/latest_release/VCF/GCF_000001405.40.gz.tbi

# Annotate
bcftools annotate -a GCF_000001405.40.gz -c ID input.vcf.gz -Oz -o with_rsids.vcf.gz

gnomAD

gnomAD files are large. Download only needed chromosomes:

# Download chromosome 1
wget https://gnomad-public-us-east-1.s3.amazonaws.com/release/4.0/vcf/genomes/gnomad.genomes.v4.0.sites.chr1.vcf.bgz
wget https://gnomad-public-us-east-1.s3.amazonaws.com/release/4.0/vcf/genomes/gnomad.genomes.v4.0.sites.chr1.vcf.bgz.tbi

# Annotate. v4 renamed popmax -> grpmax; the filtering-AF fields are
# fafmax_faf95_max / fafmax_faf95_max_joint (prefer these over global AF, see below).
bcftools annotate -a gnomad.genomes.v4.0.sites.chr1.vcf.bgz \
    -c INFO/AF,INFO/AF_grpmax,INFO/fafmax_faf95_max \
    -r chr1 \
    input.vcf.gz -Oz -o with_gnomad.vcf.gz

ClinVar

# Download
wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz
wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz.tbi

# Annotate
bcftools annotate -a clinvar.vcf.gz \
    -c INFO/CLNSIG,INFO/CLNDN \
    input.vcf.gz -Oz -o with_clinvar.vcf.gz

Complete Annotation Workflows

Basic Pipeline

# 1. Add rsIDs
bcftools annotate -a dbsnp.vcf.gz -c ID input.vcf.gz | \
# 2. Add population frequencies
bcftools annotate -a gnomad.vcf.gz -c INFO/AF | \
# 3. Output
bgzip -c > annotated.vcf.gz

bcftools index annotated.vcf.gz

Clinical Variant Analysis

# Normalize
bcftools norm -f reference.fa -m-any input.vcf.gz | \
# Add clinical significance
bcftools annotate -a clinvar.vcf.gz -c INFO/CLNSIG,INFO/CLNDN | \
# Predict consequences
bcftools csq -f reference.fa -g genes.gff3.gz | \
# Filter pathogenic
bcftools view -i 'INFO/CLNSIG~"Pathogenic"' -Oz -o pathogenic.vcf.gz

Rare Variant Analysis

# Annotate with population frequency
bcftools annotate -a gnomad.vcf.gz -c INFO/AF input.vcf.gz | \
# Filter rare variants (AF < 1%)
bcftools filter -i 'INFO/AF<0.01 || INFO/AF="."' | \
# Predict consequences
bcftools csq -f reference.fa -g genes.gff3.gz -Oz -o rare_consequences.vcf.gz

Parsing SnpEff ANN in Python

VEP CSQ parsing is shown in SKILL.md. SnpEff writes an ANN field with a fixed pipe order:

from cyvcf2 import VCF

ann_fields = ['Allele', 'Annotation', 'Impact', 'Gene_Name', 'Gene_ID',
              'Feature_Type', 'Feature_ID', 'Transcript_BioType', 'Rank',
              'HGVS_c', 'HGVS_p', 'cDNA_pos', 'CDS_pos', 'Protein_pos', 'Distance']

for variant in VCF('snpeff_output.vcf'):
    ann = variant.INFO.get('ANN')
    if not ann:
        continue
    for block in ann.split(','):
        parsed = dict(zip(ann_fields, block.split('|')[:len(ann_fields)]))
        # Impact is a SnpEff triage bucket, not ACMG evidence; check the SO Annotation term
        if parsed['Impact'] == 'HIGH':
            print(variant.CHROM, variant.POS, parsed['Gene_Name'], parsed['Annotation'])

Deciding transcript, predictor, and frequency metric

These three decisions determine what the annotation means; SKILL.md covers the reasoning, summarized here as prompts the agent resolves:

Troubleshooting

"no such tag"

The annotation file doesn't have the requested field:

# Check available fields
bcftools view -h annotation.vcf.gz | grep "^##INFO"

"chromosome not found"

Chromosome names don't match between files:

# Check chromosomes in each file
bcftools view -H input.vcf.gz | cut -f1 | sort -u | head
bcftools view -H annotation.vcf.gz | cut -f1 | sort -u | head

# Rename if needed
bcftools annotate --rename-chrs rename.txt input.vcf.gz

"index required"

Annotation file needs an index:

tabix -p vcf annotation.vcf.gz

Slow Annotation

For large files, restrict to regions:

bcftools annotate -a database.vcf.gz -c INFO/AF \
    -r chr1,chr2,chr3 \
    input.vcf.gz -Oz -o annotated.vcf.gz

Example Prompts

Basic annotation

"Add rsIDs from dbSNP and gnomAD allele frequencies to my VCF"

"Predict functional consequences of my variants using bcftools csq"

Clinical-grade annotation

"Annotate with VEP on GRCh38 reporting the MANE Select transcript, and explain any coordinate change from our old Ensembl-canonical pipeline"

"Why do VEP and SnpEff report different consequences for this indel, and which should I trust for HGVS?"

"Add one calibrated missense predictor for PP3/BP4 and tell me the calibrated REVEL thresholds"

"Score splice effects with SpliceAI and tell me whether the window is wide enough for a deep-intronic variant"

"Report the grpmax filtering allele frequency for this variant and whether presence in gnomAD rules out pathogenicity for a recessive disorder"

What the Agent Will Do

  1. Normalize (left-align + split multiallelics) before annotation, and rely on the engine for 3'-shifted HGVS rather than deriving it from POS
  2. Confirm genome build and chromosome naming match across the VCF and every annotation source
  3. Annotate on MANE Select (with MANE Plus Clinical where assigned), not worst-consequence, pinning engine+version+transcript set
  4. Attach one calibrated predictor per evidence type and the grpmax filtering AF, avoiding correlated-predictor stacks and global AF cutoffs
  5. Validate output fields, report annotation hit rate, and record the pinned axes for reproducibility

Tips