Bioinformatics File Formats Explained: FASTA, FASTQ, SAM/BAM, VCF, GTF and BED

Chapter 4 of 5: Variants and annotation – VCF, BCF, GFF3, GTF and BED

VCF: the variant call format

A VCF file has three parts: meta-information lines starting with ##, one header line starting with #CHROM, and data lines, one per variant site. Let’s write a small VCF by hand. Its four variants are real events in our toy reads: r001 carries a 4-base insertion and a 1-base deletion, read x3 a 4-base insertion, and read x1 a T→A mismatch. The quality values are made up for practice (illustration, not a real call set).

VCF columns must be separated by tabs, and tabs are awkward to type. So we type spaces and let sed convert them on every line that does not start with ## (the meta-lines contain spaces inside descriptions that must stay).

cat << 'EOF' | sed -E '/^##/!s/ +/\t/g' > toy.vcf
##fileformat=VCFv4.2
##contig=<ID=ref,length=45>
##contig=<ID=ref2,length=40>
##FILTER=<ID=LowQual,Description="Quality below 20">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total read 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 toy
ref 14 . A AAGAG 30 PASS DP=2 GT:DP 0/1:2
ref 18 . AG A 25 PASS DP=2 GT:DP 0/1:2
ref2 14 . C CAAAT 40 PASS DP=4 GT:DP 0/1:4
ref2 17 . T A 12 LowQual DP=5 GT:DP 0/1:5
EOF
cat -A toy.vcf | sed -n '8,9p'

• cat << 'EOF' … EOF: a “heredoc”: everything you type until the line EOF is sent into the pipe.
• sed -E '/^##/!s/ +/\t/g': on lines not (!) starting with ##, replace each run of spaces with one tab.
• cat -A: shows hidden characters. Tabs appear as ^I and line ends as $.

Why does the deletion start at 18, not 19?

In r001, the deleted base is the G at position 19. But an allele in VCF may never be empty, so a deletion is written with the base before it: REF AG (positions 18–19), ALT A, POS 18. Insertions work the same way: A → AAGAG at 14 means “after the A at 14, insert AGAG”. This padding base is required by the specification for simple indels.

Genotypes

Step 9: Compress, index and query the VCF

bgzip toy.vcf
tabix -p vcf toy.vcf.gz
bcftools view -H toy.vcf.gz
bcftools view -H -r ref2 toy.vcf.gz
bcftools view -H -f PASS toy.vcf.gz | wc -l

• bgzip: compresses with BGZF and replaces toy.vcf with toy.vcf.gz. Use bgzip, not gzip.
• tabix -p vcf: builds toy.vcf.gz.tbi (bcftools index -t does the same).
• bcftools view -H: print data lines without the header.
• -r ref2: only one region (needs the index). -f PASS: only sites that passed all filters.

bcftools query turns a VCF into a clean table, which is the easiest way to get variants into R, Python or Excel:

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%FILTER[\t%GT]\n' toy.vcf.gz
bcftools norm -f toy.fa -c e -Ou toy.vcf.gz > /dev/null && echo "REF alleles match toy.fa"

• -f '…': the output template. %CHROM, %POS … are site fields; anything inside [ ] repeats for every sample.
• bcftools norm -f toy.fa -c e: checks every REF allele against the reference and stops with an error (e) on a mismatch. -Ou writes uncompressed BCF, which we throw away.

The REF check is the best single test that a VCF belongs to a reference. A VCF from GRCh37 checked against GRCh38 fails on thousands of sites.

# Left-align an indel (simplified idea, pseudo-code)
while REF and ALT end with the same base:
    remove that last base from both
    if REF or ALT becomes empty:
        prepend the reference base before POS to both; POS -= 1

GFF3: gene models with a hierarchy

GFF3 describes features on a sequence: genes, mRNAs, exons, CDS, UTRs. It has nine tab-separated columns, and its 9th column links features into a tree with ID and Parent. Let’s write a tiny gene, toyA, on ref: two exons, and a coding sequence (CDS) that starts with the ATG at positions 4–6 and ends with the TAG stop at 26–28. (A made-up gene for practice, not real biology.)

cat << 'EOF' | sed -E '/^#/!s/ +/\t/g' > toy.gff3
##gff-version 3
ref bhl gene 1 40 . + . ID=gene1;Name=toyA
ref bhl mRNA 1 40 . + . ID=tx1;Parent=gene1
ref bhl exon 1 15 . + . ID=exon1;Parent=tx1
ref bhl exon 20 40 . + . ID=exon2;Parent=tx1
ref bhl CDS 4 15 . + 0 ID=cds1;Parent=tx1
ref bhl CDS 20 28 . + 0 ID=cds1;Parent=tx1
EOF
grep -v "^#" toy.gff3 | cut -f 3 | sort | uniq -c

Read the tree from the bottom up: each exon and CDS has Parent=tx1, and tx1 has Parent=gene1. A gene with three transcripts has three mRNA lines, and an exon shared by two transcripts can list both: Parent=tx1,tx2. Both CDS lines share ID=cds1, which is how GFF3 says “these pieces form one CDS”.

# Phase of each CDS piece, in transcription order (pseudo-code)
phase = 0
for piece in cds_pieces_in_5prime_to_3prime_order:
    write(piece, phase)
    length = piece.end - piece.start + 1
    phase = (3 - ((length - phase) % 3)) % 3

# toyA: piece 1 is 4-15 (12 bp, phase 0): (3 - (12 % 3)) % 3 = 0, so piece 2 has phase 0

GTF: the format most RNA-seq tools want

GTF (also called GFF2.5) has the same first eight columns as GFF3, but a different 9th column: key "value"; pairs separated by spaces, and every line carries gene_id and transcript_id instead of an ID/Parent tree. An Ensembl GTF line looks like this:

1	havana	exon	3069168	3069296	.	+	.	gene_id "ENSG00000142611"; gene_version "17"; transcript_id "ENST00000511072"; transcript_version "5"; exon_number "1"; gene_name "PRDM16"; …

(The first exon line of the human GTF from Ensembl release 115; the attribute list is shortened here. Note the chromosome is 1, not chr1.)

Use the GTF and the FASTA from the same source and release. Convert between the two with gffread (gffread in.gff3 -T -o out.gtf) rather than by hand.

BED: intervals, 0-based

BED is the simplest interval format: chromosome, start, end, and up to nine optional columns. Only the first three are required. The specification allows spaces or tabs between fields, but tools such as bedtools need tabs.

ref	0	40	toyA	0	+	3	28	0	2	15,21	0,19

Step 10: Convert GFF3 to BED and prove it is right

Turn the exons and CDS of toy.gff3 into BED with awk. The only change is start − 1:

awk -v OFS='\t' '$3=="exon" {print $1, $4-1, $5, "exon", 0, $7}' toy.gff3 > exons.bed
awk -v OFS='\t' '$3=="CDS" {print $1, $4-1, $5, "CDS", 0, $7}' toy.gff3 > cds.bed
cat exons.bed cds.bed

• -v OFS='\t': print columns separated by tabs.
• $3=="exon": only lines whose 3rd column is exon.
• $4-1, $5: 1-based start minus one; the end stays the same.

Now check that the BED interval and the 1-based region give exactly the same bases:

bedtools getfasta -fi toy.fa -bed exons.bed
samtools faidx toy.fa ref:1-15
bedtools getfasta -fi toy.fa -bed cds.bed

Same 15 bases from ref:0-15 (BED, 0-based) and ref:1-15 (samtools, 1-based). Joined, the two CDS pieces read ATG TTA GAT AAG CTG TGC TAG: start codon, five codons, stop. The coordinates are right.

Step 11: Combine formats

Real questions need several formats at once. Which variants hit the coding sequence? How many reads overlap the exons?

bedtools intersect -a toy.vcf.gz -b cds.bed | cut -f 1-7
samtools view -c -L exons.bed toy.bam

• bedtools intersect -a VCF -b BED: prints the VCF lines that overlap any BED interval. bedtools handles the 1-based VCF and 0-based BED for you.
• samtools view -c -L exons.bed: counts alignments that overlap the BED regions.

Only the 4-base insertion at ref:14 falls in the CDS. Four bases is not a multiple of three, so in a real gene it would be a frameshift. The deletion at 18–19 sits in the intron (16–19), so it is not coding.

bedtools genomecov -bga -ibam toy.bam | sort -k1,1 -k2,2n > toy.bedgraph
cut -f 1,2 toy.fa.fai > toy.chrom.sizes
bedGraphToBigWig toy.bedgraph toy.chrom.sizes toy.bw