Before you start
- You've worked through Read Alignment and the samtools survival kit. This lesson stitches them together.
- A free Google account to run the notebook. No installs on your machine.
- Any term new? The Glossary has it.
Learning objectives
By the end of this lesson you will be able to: take raw sequencing reads through alignment to a list of variants (a VCF), mark duplicates, and judge a called variant by its depth and quality.
What "calling a variant" means
Every person's genome differs from the reference at millions of spots. A variant is one of those differences: a single letter changed (a SNP), or a few letters inserted or deleted (an indel). Variant calling is the process of reading the stacked-up reads at each position and deciding, statistically, "the reference says A here, but the reads clearly say G, so this sample carries a variant."
The skill is knowing the difference between a real variant and a sequencing error. One read disagreeing with the reference is probably noise. Twenty reads all agreeing on the same change is a variant. A variant caller does exactly this reasoning, position by position, across the whole genome.
Decode the jargonSNP and indel
A SNP (single-nucleotide polymorphism) is a one-letter difference, like an A where the reference has a G. An indel is a small insertion or deletion of bases. Together they're the most common kinds of variants, and the ones you'll call in this lesson.
The whole pipeline on one screen
You already know the first four boxes. This lesson adds the last one, bcftools, and connects them all into a single runnable pipeline.
▶ Run the full pipeline in Google Colab
The notebook builds a reference, plants known variants in a sample, then recovers them, so you can check the answer.
1Install and prepare
We need three tools: bwa to align, samtools to sort and index, and bcftools to call variants.
!apt-get -qq install -y bwa samtools bcftools
The notebook then creates a reference and a sample genome that's identical except for a few planted SNPs at positions we record. Reads are generated from the sample, so when we call variants we should recover exactly those planted changes, a built-in answer key.
2Align, sort, index
This is everything from the last two lessons, compressed into the standard three lines. Notice how each tool hands off to the next.
!bwa index reference.fasta !bwa mem reference.fasta sample.fastq | samtools sort -o sample.sorted.bam !samtools markdup sample.sorted.bam sample.markdup.bam # flag PCR/optical duplicates !samtools index sample.markdup.bam
Why mark duplicates?
PCR and optical duplicates are copies of the same original fragment. Left in, they inflate the apparent support for a position and cause false variant calls, so real pipelines flag them before calling. On paired-end data you run samtools fixmate -m (on a name-sorted BAM) before the coordinate sort so markdup has the mate information it needs.
3Call variants with bcftools
Calling is two steps piped together. bcftools mpileup walks the genome and, at every position, gathers the pile of read bases sitting there. bcftools call then judges each pile and decides whether it's a variant.
!bcftools mpileup -f reference.fasta sample.markdup.bam | \ bcftools call -mv -Oz -o calls.vcf.gz !bcftools norm -f reference.fasta calls.vcf.gz -Oz -o calls.norm.vcf.gz # left-align indels !bcftools index calls.norm.vcf.gz
The extra bcftools norm step left-aligns and normalises indels so the same insertion or deletion is always written the same way, which matters the moment you compare or annotate variants on real data.
The flags read as: -m use the modern multiallelic caller, -v output only variant sites (skip positions that match the reference), and -Oz write a compressed VCF.
Decode the jargonPileup
A pileup is the column of read bases stacked over a single genome position, like reading straight down through all the reads that cover that spot. The caller looks at each column and asks: do the reads agree with the reference, or do enough of them disagree to call a variant?
4Read your VCF
VCF (Variant Call Format) is the standard text file listing variants. After header lines (starting with #), each row is one variant.
!bcftools view calls.vcf.gz | grep -v '^##' | head
| Column | Meaning |
|---|---|
CHROM | Which chromosome the variant is on. |
POS | The position of the variant. |
REF | The base in the reference genome. |
ALT | The alternative base seen in your sample, the change. |
QUAL | How confident the caller is that this variant is real. |
INFO | Extra details, such as depth (DP) at this site. |
Compare the POS and ALT columns to the variants the notebook planted. If they match, you've just validated your own pipeline end to end.
⚠️ Not every call is real
Variant callers make mistakes, especially at low coverage or in repetitive regions. Real projects filter the VCF by quality and depth (for example, dropping calls with QUAL below a threshold or very low DP) before trusting it. Raw calls are a starting point, not a final answer.
🏆 Your Track 2 deliverable
Turn the notebook into something you can show. Build a documented variant-calling mini-pipeline:
- Run the full FASTQ → VCF pipeline end to end in the notebook.
- Add a markdown cell at the top explaining, in your own words, what each stage does and why.
- Report two numbers: your mapping rate (from
flagstat) and how many variants you called. - Confirm the called variants match the planted ones, and write one sentence on what would change with real data.
Save it to your GitHub and you have a portfolio piece that shows you can run the core genomics workflow.
🔶 Level up: real data and GATK
For human germline work, the field standard is the GATK Best Practices pipeline, which adds duplicate marking, base-quality recalibration, and joint genotyping. bcftools (what you used here) is lighter, faster, and excellent for many projects, especially microbial genomes and quick analyses. Same concepts, more rigor. Try the pipeline next on a small real dataset, like an E. coli run from the ENA.
Now do it on real data: E. coli from the ENA
The toy run above is honest about being a rehearsal. Here is the same pipeline on a real sequencing run of E. coli, downloaded from the European Nucleotide Archive (ENA). Everything runs in Colab; the download and alignment take a few minutes, not seconds, because this is real data. The reads are a real isolate from the Lenski long-term evolution experiment (run SRR2584866), and the reference is the standard E. coli K-12 genome.
!apt-get -qq install -y bwa samtools bcftools fastp # 1. Reference genome (E. coli K-12 MG1655) !wget -q https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz !gunzip -f GCF_000005845.2_ASM584v2_genomic.fna.gz && mv GCF_000005845.2_ASM584v2_genomic.fna ecoli.fasta # 2. A real paired-end sequencing run from the ENA !wget -q https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR258/006/SRR2584866/SRR2584866_1.fastq.gz !wget -q https://ftp.sra.ebi.ac.uk/vol1/fastq/SRR258/006/SRR2584866/SRR2584866_2.fastq.gz # 3. QC and trim (paired-end: fastp auto-detects adapters from the read overlap) !fastp -i SRR2584866_1.fastq.gz -I SRR2584866_2.fastq.gz -o r1.fq.gz -O r2.fq.gz -h fastp.html # 4. Align, then fixmate -> sort -> markdup (the real paired-end dedup flow) !bwa index ecoli.fasta !bwa mem ecoli.fasta r1.fq.gz r2.fq.gz | samtools sort -n -o namesorted.bam !samtools fixmate -m namesorted.bam fixmate.bam !samtools sort -o sorted.bam fixmate.bam !samtools markdup sorted.bam markdup.bam && samtools index markdup.bam # 5. How well did the reads ACTUALLY map? (this number will not be 100%) !samtools flagstat markdup.bam # 6. Call and normalise !bcftools mpileup -f ecoli.fasta markdup.bam | bcftools call -mv --ploidy 1 -Oz -o ecoli.vcf.gz !bcftools norm -f ecoli.fasta ecoli.vcf.gz -Oz -o ecoli.norm.vcf.gz && bcftools index ecoli.norm.vcf.gz !bcftools stats ecoli.norm.vcf.gz | grep "number of records:"
What real data breaks (that the toy run hid)
Look at flagstat: the mapping rate is high but not 100%, because real reads include contamination, low-quality ends, and repeats that map ambiguously. markdup now flags a real percentage of duplicates. There is no planted answer key; you judge each call by its depth (DP) and quality (QUAL) instead of checking it against a list. And you get hundreds of variants, some multi-allelic and some indels, which is exactly why the norm step matters. This is the gap between "ran the notebook" and "ran an analysis," and you just crossed it.
Check your understanding
bcftools mpileup and bcftools call each contribute?mpileup gathers the evidence: at each position it assembles the pile of read bases and their qualities. call makes the decision: it judges each pile statistically and outputs the positions it believes are genuine variants.