A hands-on session built on one real whole-genome HIV-1 deep-sequencing sample — 9,131 positions, median depth 4,002×.
Partway through, we will produce a confident, FDR-significant result. It will be wrong. Catching it is the point.
→ next ← back F fullscreen
A per-position pileup: for every base of the HIV-1 genome, how many reads showed A, C, G and T.
That per-base count is the whole game. A consensus FASTA collapses each position to a single letter and throws the population away. These counts let us see underneath the consensus.
Everything today — types, diversity, burden, minority variants — is computed from these four columns.
| Column | Meaning |
|---|---|
| Nucleotide pos | position, 1–9181 |
| ref | reference base |
| strain | this sample's consensus |
| depth (HQ) | reads passing filter |
| A/C/G/T count | the population |
| indel count | reads supporting indel |
A change in sequence relative to a reference. Three things must be pinned down before the word means anything at all.
Reference genome? The patient's own baseline? Population consensus? "Mutation" is always a comparison — name the thing you compared to.
A virus is a population of genomes. A change in 100% of them and in 2% of them are both real — and mean completely different things.
Sequencers make errors. A 0.3% signal may be a mutation, or may be noise. This question will bite us later.
Our file offers two comparisons: ref vs strain → consensus mutations (Sanger sees these). strain vs base counts → minority variants (deep sequencing only).
This sample is almost certainly a non-subtype-B virus measured against a subtype-B reference (HXB2 / NC_001802.1).
Those 1,287 differences are subtype differences — evolution that happened decades before this infection existed. They did not "arise" in this patient.
A difference from the reference is not a mutation that arose in your sample.
The reference is a historical accident — one man's virus, sequenced in 1983. It is not a wild-type ideal, and it is not "normal."
This single confusion produces more bad clinical genomics than any statistical error.
| Type | What happens | In this file | Count |
|---|---|---|---|
| Substitution (SNV/SNP) | one base swapped for another | ref & strain both 1 base | 9,005 |
| Insertion | bases added | strain longer than ref | 9 |
| Deletion | bases removed | ref longer than strain | 17 |
An indel sized in multiples of 3 adds or removes whole codons — the protein is shorter or longer, but still readable.
Any other size shifts the reading frame: every downstream codon is garbled and the protein is usually destroyed. Same "one mutation" — wildly different consequence.
| Frame effect | n |
|---|---|
| In-frame (×3) | 21 |
| Frameshift | 5 |
Mostly in-frame — consistent with real viral variation. Random error has no reason to respect codon boundaries.
Transition (Ti) — purine↔purine A↔G, or pyrimidine↔pyrimidine C↔T. Chemically similar bases; easy swaps.
Transversion (Tv) — purine↔pyrimidine. A bigger chemical change.
There are 8 possible transversions but only 4 transitions. So pure chance predicts Ti/Tv = 0.5.
Real data runs far higher — transitions are mechanistically favoured (deamination of methyl-C → T) and more often synonymous, so selection removes them less often.
These consensus differences carry real evolutionary signal.
Remember 0.5. It is the cheapest QC metric you own, and in twenty minutes it will save us from publishing an artefact.
Needs the reading frame — can't be read off a pileup alone. But it's what clinicians actually care about.
| Class | Effect on protein | Example |
|---|---|---|
| Synonymous (silent) | codon changes, amino acid does not | GGA→GGG, both Gly |
| Missense | one amino acid replaced | HIV RT K103N |
| Nonsense | codon → STOP, protein truncated | usually lethal to the virus |
| Frameshift | indel not ×3, everything downstream garbled | usually knockout |
| Non-coding | outside genes; may affect expression | LTR promoter |
Ratio of nonsynonymous to synonymous change.
< 1 purifying · ≈ 1 neutral · > 1 positive (diversifying)
dN/dS > 1 is the signature of immune escape and drug pressure.
Germline (inherited) vs somatic (acquired — the basis of cancer genomics).
For a virus: transmitted vs acquired — which, for HIV resistance, decides whether a first-line regimen was ever going to work.
Minority calls, Shannon entropy, nucleotide diversity, burden — all of it is arithmetic on this ratio.
Sanity check: our computed majority base agrees with the file's own consensus call at 100.00% of sites. The arithmetic is sound.
Depth first, always. We keep only sites with ≥100× coverage: 8,924 of 9,005 (99.1%). You cannot call a variant you did not sequence deeply enough to see.
| Mutation frequency | Mutation rate | |
|---|---|---|
| Question | How many mutants are here now? | How fast do new ones appear? |
| Units | a proportion (dimensionless) | per base per replication |
| Needs | one sample | ≥2 timepoints, or known generations |
| From this file? | Yes | No |
Frequency is a snapshot — the product of mutation and selection and drift and elapsed time. A high-frequency variant may be common because it arose often, or because it was favoured. One timepoint cannot tell you which.
That second number is the whole ballgame. Resistance is not a question of whether it is generated — it is generated constantly. It is a question of whether the drug gives it a reason to expand.
Divergence asks: how far from the reference?
Diversity asks: how variable is the population within itself?
A sample can sit far from the reference yet be perfectly uniform, or sit right on it and be a seething swarm. This is the quasispecies concept: an infection is not a virus — it is a cloud of related genomes.
0 = invariant · 1.386 = all four bases equally likely
The probability that two randomly drawn reads differ at a site.
The n/(n−1) correction matters. Without it, low-coverage sites look artificially less diverse — purely because you sampled fewer reads. That is a sampling artefact masquerading as biology.
HIV-1 is a masterclass in genomic compression: genes overlap, in different reading frames. vpu sits inside env; tat/rev exon 2 sit inside env.
A naive case_when() assigns each site to one gene and silently discards every overlap — 685 sites here. Use a long-format join so a site can belong to several genes.
env sits high — the surface protein, under relentless antibody pressure, and tolerant of change. pol — the drug target, structurally constrained — sits lower. Diversity maps onto selection.
Sanger sequencing sees ~20% and above. Anything rarer is invisible — buried in the chromatogram noise. Deep sequencing, at thousands of reads per site, can in principle reach 1% or lower.
Variants in that gap — present, real, sub-consensus — are minority variants. In HIV they matter enormously:
We did it properly. If the per-base error rate is ε, then under the null "this is just error" the alt count is Binomial(n, ε). So we tested every site against ε = 0.2%, and corrected for multiple testing with Benjamini–Hochberg.
A clean, defensible, fully corrected result. We could write this up this afternoon.
So let us try to destroy it.
The next slide is the most important one in this workshop.
This is damning. At 0.1% the ratio is 0.09 — not merely random (0.5), but far below it. Random error would give 0.5. Something worse than random is happening: a systematic, direction-specific artefact.
At our 1% threshold, Ti/Tv sits at exactly 0.50. Roughly half of those 451 "significant" variants are noise.
Look at what dominates low VAF: A>C, T>G, G>T, C>A — all transversions.
And notice they form reverse-complement pairs (A>C ≡ T>G; G>T ≡ C>A). That strand symmetry is the fingerprint of a chemical/technical artefact — not biology.
Above 2% the spectrum flips to G>A, A>G, C>T — transitions. Real biology.
At ~4,000× depth, a 1% variant means ~40 supporting reads. Against a null of ε = 0.2%, forty reads is overwhelmingly significant — p is astronomically small.
So the test passes every single call. It never had a chance to reject anything.
The binomial test assumes error is uniform at 0.2% everywhere. Real error is context- and site-specific — at some sites it runs at 2%. The model cannot see that, so it confidently certifies artefacts.
Statistical significance is not biological reality.
A p-value only ever tests the model you handed it. With enough depth you can get vanishingly small p-values for pure garbage.
Deep sequencing does not protect you from this — it guarantees it. Depth makes everything significant.
Filter out the artefact signature (A>C, T>G, G>T, C>A) and raise the threshold to 2% — justified empirically, from Ti/Tv and the spectrum, not by convention.
Ti/Tv jumps from 0.50 to 5.25. We discarded 83% of our variants, and what remains finally looks like biology.
| VAF threshold | n (naive) | Ti/Tv (naive) | n (filtered) | Ti/Tv (filtered) |
|---|---|---|---|---|
| ≥ 1% | 451 | 0.50 | 177 | 5.81 |
| ≥ 2% | 129 | 0.95 | 75 | 5.25 |
| ≥ 5% | 20 | 0.82 | 11 | 4.50 |
If our surviving calls are real, they should carry biological signal. HIV has a famous one.
Host APOBEC3G deaminates cytosine during reverse transcription, producing G>A hypermutation. The virus's vif gene exists precisely to destroy APOBEC3G — that is its whole job.
G>A climbs from 1.6% to 20% as VAF rises. Sequencing error has no reason whatsoever to do that. A specific, known, host-driven biological process does.
This is independent confirmation that the surviving calls are real — and that the low-VAF tier is not.
Single-variant tests work for common variants. For rare ones they fail — and the reason is purely statistical. If a variant appears in 3 people out of 5,000, nothing on Earth has power to associate it with anything.
The burden test's move is collapsing: stop testing variants, start testing units — a gene, a region, a domain.
Instead of "does variant X associate with disease?", ask "do carriers of any rare variant in gene G differ from non-carriers?"
You no longer know which variant matters — but you can now detect that the gene matters at all.
| Test | Collapsing rule | Best when |
|---|---|---|
| CAST | binary: carries ≥1 rare variant | rare variants, same direction |
| Burden / MB | count, weighted by rarity | effects same direction |
| SKAT | variance-component | mixed directions |
| SKAT-O | data-driven blend | you don't know — the safe default |
The plain burden test's fatal flaw: it assumes every variant in the unit pushes the same way. If half are harmful and half protective, they cancel — and it sees nothing. That is why SKAT exists.
A classical burden test compares cases vs controls across many individuals. We have one sample. We cannot run one — and you should be suspicious of anyone who claims otherwise.
The same machinery, in its within-sample form: collapse variants into genes, test each against a background rate. This is exactly how regional mutation-enrichment scans work in cancer genomics — and it drills the identical logic.
Minority variants fall uniformly across the genome, at the genome-wide rate. Under that null, the count in a gene is Binomial(gene length, background rate). Test it. Correct it. Report the q-value.
Nine genes = nine tests. At α = 0.05 you expect ~0.45 false positives by chance. Always report the q-value.
fold 1.78 → 1.44 · q 0.0005 → 0.65
Its "enrichment" was driven by artefact-signature calls — 77% of vif's 1% calls were A>C/T>G/G>T/C>A. A confident, FDR-corrected, entirely false result.
fold 1.33 → 1.74 · q 0.0007 → 0.011
It got stronger when the noise was removed — the hallmark of real signal, and exactly what HIV biology predicts for the antibody-facing envelope protein.
Our q-value was 0.0005 and the result was still garbage.
Multiple-testing correction did not save us. It was never designed to. FDR controls for chance — it has nothing at all to say about systematic bias.
No amount of statistics rescues contaminated input. Only the biology told us the truth:
Every one of those checks is a domain check, not a statistical one. They work because we know what HIV does — and what a sequencer does when it is wrong.
Know your data's failure modes, or your p-values will lie to you with perfect confidence.
State: relative to what, at what frequency, with what confidence. 14% of this genome differs from HXB2 because it is a different subtype — not because the patient mutated.
One timepoint gives a snapshot, never a speed.
Independent axes — and diversity maps onto selection (env high, pol constrained).
451 variants at 1%; 75 at 2% with filtering. Neither is "the truth." Report both numbers, or report nothing.
0.5 = noise. If your minority calls sit at 0.5, they are noise — no matter how small the p-value.
We got q = 0.0005 for vif and it was pure artefact. Only the biology told us the truth.
Consensus Ti/Tv was 1.73. Before running the minority analysis, predict whether minority Ti/Tv will be higher or lower — and why. Then check.
Compute diversity in 100 bp windows across env only. Can you locate V1–V5 as peaks? What does it mean clinically if a patient's env diversity is very low?
We justified 2% empirically. Now defend it to a clinician. What must you state alongside "we found N minority variants" for the claim to be reproducible?
In run_burden(), make half the variants protective. Re-run. Watch power collapse. Then explain in one sentence why SKAT exists.
Everything on these slides is reproduced in mutations_workshop.Rmd — every number is computed from the CSV, not typed in. Change params at the top and the whole document re-derives itself.
Your job is to know when it's lying.
A · C · G · T