Describe variant calling pipelines

HOW: the pipeline stages

Figure — Describe variant calling pipelines

The canonical DNA-seq pipeline flows: Raw reads → QC/trim → Align → Post-process → Call → Filter → Annotate.

1. Quality control & trimming

WHY: Sequencers append a per-base Phred quality score. Low-quality ends and adapter contamination create false variants downstream.
HOW: Assess quality with a QC-reporting tool (FastQC — it only produces reports, it does not trim). Then trim adapters and low-quality bases with dedicated trimming tools (Trimmomatic, fastp).

2. Alignment (mapping)

WHY: We must know where each read came from before we can compare it to the reference.
HOW: Align reads to the reference with a Burrows–Wheeler-based aligner (BWA-MEM, Bowtie2). Output = SAM/BAM file. Each read gets a MAPQ mapping-quality score (same Phred logic: probability the read is mis-placed).

3. Post-alignment processing

Sort & index BAM by coordinate.
Mark duplicates — WHY: PCR/optical duplicates are copies of the same original molecule; counting them as independent evidence fakes support for a variant. We flag them so they count once.
Base Quality Score Recalibration (BQSR) — WHY: sequencers systematically mis-estimate $Q$ ; recalibration corrects the numbers so the probability model is honest.

4. Variant calling proper — the Bayesian core

Formula Genotype likelihood via Bayes — from scratch

Goal: find the most probable genotype $G$ (e.g. AA, AG, GG) given observed reads $D$ . Step 1 — Bayes' theorem (WHY: we can only compute $P(\text{data}|\text{genotype})$ from the error model, but we want $P(\text{genotype}|\text{data})$ ): $P(G\mid D)=\frac{P(D\mid G)\,P(G)}{P(D)}$

Why this step? $P(D\mid G)$ is computable from base error rates; the prior $P(G)$ comes from population allele frequencies. Step 2 — likelihood of the data. Reads are (assumed) independent, so multiply per-read probabilities: $P(D\mid G)=\prod_{i=1}^{n} P(d_i\mid G)$
Why? Independence turns a hard joint probability into a simple product. Step 3 — single read under a heterozygote (genotype has two alleles $b_1,b_2$ , read shows base $d_i$ with error prob $\epsilon_i=10^{-Q_i/10}$ ): $P(d_i\mid G=b_1b_2)=\tfrac12 P(d_i\mid b_1)+\tfrac12 P(d_i\mid b_2)$ where $P(d_i\mid b)=1-\epsilon_i$ if $d_i=b$ , else $\epsilon_i/3$ .
Why the ½? Each of the two alleles is equally likely to be the template for a given read.
Why $\epsilon/3$ ? If it's an error, it could be any of the 3 other bases, split evenly. Step 5 — decide. Pick $G$ maximising $P(G\mid D)$ ; report confidence as a Phred-scaled QUAL score. Tools: GATK HaplotypeCaller, FreeBayes, bcftools, DeepVariant (CNN-based).

Output = a VCF (Variant Call Format) file: one line per variant with CHROM, POS, REF, ALT, QUAL, and per-sample genotypes.

5. Filtering

WHY: Even good callers make errors. Remove variants with low QUAL, low depth (DP), strand bias, etc. (Hard filters or GATK's VQSR machine-learning recalibration.)

6. Annotation

WHY: A position is meaningless until you know what gene/effect it hits.
HOW: Tools like ANNOVAR, VEP, SnpEff label variants (missense, nonsense, synonymous, coding vs intron), and cross-reference dbSNP, ClinVar, gnomAD.

Worked examples

Worked example Example 2 — Is this a variant?

Position has 20 reads, all $Q=30$ ( $\epsilon=0.001$ ): 20 show A, reference is A. Compare genotype AA vs AG.

$P(D\mid AA)=(1-\epsilon)^{20}=0.999^{20}\approx 0.980$ . Why? All reads match A, each correct with prob $1-\epsilon$ .
Per-read under het (read shows A): $\tfrac12(1-\epsilon)+\tfrac12(\epsilon/3)=\tfrac12(0.999)+\tfrac12(0.000333)\approx 0.49967\;(\approx 0.5)$ .
$P(D\mid AG)=(0.49967)^{20}\approx 9.4\times10^{-7}$ . Why? Under het, each read matching A has prob $\approx 0.5$ .
AA wins overwhelmingly → no variant. Correct: pure reference support ≠ variant.

Common mistakes

Recall Feynman: explain to a 12-year-old

You copied a huge book by tearing it into confetti and scanning each piece. To find typos, you first figure out where each scrap belongs in the original book (alignment). Then, at every spot, you look at all the scraps covering it. If most scraps agree the letter changed, it's a real typo (variant). If only one weird scrap disagrees, it's probably a smudge (sequencing error). You even throw away scraps that are just Xerox-copies of the same page so you don't get fooled into thinking a smudge is real. That's variant calling!

Flashcards

What does variant calling identify?

Positions where a sample's genome differs from the reference genome, classified as SNVs, indels, or structural variants.

Define the Phred score formula and its inverse.

Q=-10\log_{10}P

and

P=10^{-Q/10}

, where P is the base-call error probability.

What error probability corresponds to Q30?

10^{-3}=0.001

(1 in 1000).

Does FastQC perform trimming?

No — FastQC only generates quality-control reports; trimming is done by tools like Trimmomatic or fastp.

What is the standard output file of alignment?

SAM/BAM (BAM = binary compressed SAM).

What is MAPQ?

Mapping quality: Phred-scaled probability that a read is aligned to the wrong location.

Why mark PCR duplicates before calling?

Duplicates are copies of the same original molecule; treating them as independent evidence violates the independence assumption in the genotype likelihood product and inflates false confidence.

State Bayes' theorem for genotype calling.

P(G\mid D)=P(D\mid G)P(G)/P(D)

; choose G maximizing the posterior.

Why is

P(D\mid G)

a product over reads?

Reads are assumed independent, so joint probability = product of per-read probabilities.

Under a heterozygote AG, what is P(read shows A)?

\tfrac12(1-\epsilon)+\tfrac12(\epsilon/3)\approx 0.49967\approx 0.5

, where

\epsilon

is the base error probability.

What file format stores called variants?

VCF (Variant Call Format).

What does BQSR correct?

Systematic miscalibration of base quality scores by the sequencer.

Name the ordered pipeline stages.

QC/trim → Align → Post-process (sort, mark duplicates, BQSR) → Call → Filter → Annotate.

Why annotate variants?

To assign biological meaning (gene, effect: missense/nonsense/synonymous) and cross-reference databases like dbSNP/ClinVar/gnomAD.

Connections

Phred Quality Scores
Sequence Alignment (BWA / Burrows-Wheeler Transform)
Bayes' Theorem
VCF and BAM File Formats
SNPs and Genetic Variation
GATK Best Practices
Next-Generation Sequencing
Genome Annotation

Concept Map

Hinglish (regional understanding)

Intuition Hinglish mein samjho

Socho tumne ek reference genome ki ek "asli book" maan li, aur sequencer ne tumhare sample ko crore chote tukdo (reads) me tod diya. Variant calling ka matlab hai: har read ko wapas sahi jagah pe rakho (alignment), aur phir dekho ki kahan tumhara genome reference se alag hai — koi single base badla (SNP), koi chota insertion/deletion (indel), ya bada structural change. Yeh sirf "difference dhundhna" nahi hai; asli sawal hai — yeh difference real hai ya sequencing machine ki galti (noise)?

Isko solve karne ke liye hum probability use karte hain. Har base call ke saath ek Phred score $Q$ aata hai, jahan $Q=-10\log_{10}P$ ; yaani $Q30$ ka matlab error probability $0.001$ , yaani accha base. Ek position pe agar 10 reads me se 5 'A' aur 5 'G' dikha rahe hain, to hum Bayes theorem se compute karte hain: $P(genotype|data) \propto P(data|genotype)\times P(genotype)$ . Jo genotype (AA, AG, GG) sabse zyada probable ho wahi call karte hain. Isliye 5-5 wale case me AG heterozygous winner hota hai.

Ek chhoti si baat yaad rakho: QC step me FastQC sirf report banata hai (quality dikhata hai), trimming nahi karta — trimming ke liye Trimmomatic ya fastp use hote hain. Aur ek important trap: PCR duplicates. Agar 50 reads asal me ek hi molecule ki copies hain, to woh independent evidence nahi hain — isliye "mark duplicates" step lagta hai, taki fake confidence na aaye. Pipeline ke end me filtering aur annotation zaroori hai, kyunki VCF me jo variant hai woh sirf candidate hai, guaranteed sach nahi. Yeh pipeline medical genomics, cancer mutation detection, aur rare disease diagnosis me backbone hai — isliye ise samajhna high-yield hai.

WHAT is variant calling?

HOW: the pipeline stages

1. Quality control & trimming

2. Alignment (mapping)

3. Post-alignment processing

4. Variant calling proper — the Bayesian core

5. Filtering

6. Annotation

Worked examples

Common mistakes

Flashcards

Connections

Concept Map

Hinglish (regional understanding)

Test yourself — Bioinformatics & Computational Biology

Connections