02-day5_variant_analysis.Rmd

---
title: "Day 5: Variant Analysis in R"
subtitle: "Applied Computational Genomics — SS26"
author: "Your Name"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 3
    theme: flatly
    highlight: tango
    df_print: paged
    code_folding: hide
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo    = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width  = 9,
  fig.height = 5
)
```

# Overview

Today we analyse the variants we called and annotated on **Day 3** (germline)
and **Day 4** (somatic) — entirely in R. No new bioinformatics tools.

```
Day 3 output  →  germline variants.tsv   →  R analysis
Day 4 output  →  somatic annotated .tsv  →  R analysis
                                                ↓
                                     Germline vs Somatic comparison
```

**Key question:** How does the variant allele frequency (VAF) distribution
differ between germline variants in a healthy individual and somatic mutations
in a tumour?

---

# Setup

## Packages

We only use **base R** and **ggplot2** today.

```{r packages, class.source="fold-show"}
# install.packages("ggplot2")   # run once if not installed
library(ggplot2)
```

## Input file paths

**Change these two paths to point to your own files.**
Everything else in this document will work automatically.

```{r paths, class.source="fold-show"}
# Germline variants TSV (you will generate this below)
GERMLINE_TSV <- "/omics/groups/OE0526/shared_data/HP-F9/03-vcfs/variants.tsv"

# Somatic variants TSV (generated at the end of Day 4)
SOMATIC_TSV  <- "/omics/groups/OE0526/shared_data/HP-F9/04-somatic_bams/results/somatic.filtered.annotated.tsv"
```

---

# Generating the input files

Run these commands **in the terminal** before opening R. The bash chunks below
are shown for reference and are not executed when you knit.

## Germline variants.tsv

We use `bcftools query` to extract the fields we need from the annotated VCF
into a plain TSV file that R can read directly.

```{bash generate-germline, eval=FALSE}
cd ~/day3_variants

# Write the header line
echo -e "CHROM\tPOS\tREF\tALT\tQUAL\tDP\tDP4\tGT\tgene\tbiotype\tConsequence\tamino_acid_change\tgnomad_popmax_af\tgnomad_filter" \
  > variants.tsv

# Extract one row per variant
bcftools query \
  -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%DP\t%DP4\t[%GT]\t%gene\t%biotype\t%Consequence\t%amino_acid_change\t%gnomad_popmax_af\t%gnomad_filter\n' \
  HG003_father.chr22.annotated.split.popfreq.vcf.gz >> variants.tsv

wc -l variants.tsv   # should be ~88,000 lines
```

**What each column means:**

| Column | Description |
|--------|-------------|
| CHROM / POS / REF / ALT | Genomic coordinates and alleles |
| QUAL | Variant quality score — higher means a more confident call |
| DP | Total sequencing depth at this position |
| DP4 | Read counts: ref_fwd, ref_rev, alt_fwd, alt_rev |
| GT | Genotype: `0/1` = heterozygous, `1/1` = homozygous alt |
| gene | Gene name (`.` = intergenic, not in any gene) |
| biotype | Gene type: `protein_coding`, `lncRNA`, pseudogene, etc. |
| Consequence | Functional impact: `missense`, `intron`, `stop_gained`, etc. |
| amino_acid_change | The protein-level change for coding variants |
| gnomad_popmax_af | Frequency in gnomAD across all populations. `-1` = absent from gnomAD |
| gnomad_filter | gnomAD quality flag |

## Somatic annotated TSV

This file was already generated at the end of Day 4 using `04-vcf_to_table.sh`
— no extra steps needed.

---

# Part 1: Germline variant analysis

## 1.1 Load and explore

Let's load the germline variants TSV. We have ~88,000 raw variant calls on
chromosome 22 from a healthy individual (HG003, the father in the
Genome in a Bottle trio).

```{r load-germline}
germ <- read.delim(GERMLINE_TSV, sep = "\t", header = TRUE,
                   stringsAsFactors = FALSE)
# How big is it?
dim(germ)
```

What does the data look like? Let's check column names and data types.

```{r explore-germline}
str(germ)
```

Let's get a quick feel for the numeric columns — what ranges do quality,
depth, and population frequency take?

```{r summary-germline}
summary(germ[, c("QUAL", "DP", "gnomad_popmax_af")])
```

> Notice that `gnomad_popmax_af` has a minimum of `-1`. This is a sentinel
> value meaning the variant was **not observed in gnomAD at all** — it is so
> rare it has never been seen in any of the ~76,000 genomes in the database.
> We will treat `-1` as "ultra-rare" and **keep** these variants in the rare
> filter step.

Since we used a chr22-restricted VCF, all variants should be on chromosome 22.
Let's confirm.

```{r chr-table}
sort(table(germ$CHROM), decreasing = TRUE)
```

A healthy diploid person carries two copies of every chromosome. Variants can
be **heterozygous** (one mutant copy, `GT = 0/1`) or **homozygous** (both
copies mutant, `GT = 1/1`). We expect far more heterozygous variants — most
rare mutations are only on one chromosome.

---

> **Exercise 1.1:** How many variants are heterozygous (`0/1`) vs homozygous
> (`1/1`)? Use `table()`.
>
> **Hint:** `table(germ$GT)`

```{r gt-table, eval=FALSE}
# --- Solution ---
table(germ$GT)
```

---

What consequence types does bcftools-csq annotate? This tells us the
**functional impact** of each variant on the nearest transcript.

```{r consequence-table}
sort(table(germ$Consequence), decreasing = TRUE)
```

> The vast majority of variants are intronic or intergenic — they fall inside
> gene bodies but not in exons, or entirely between genes. Only a small
> fraction affect protein-coding sequence. This reflects the actual composition
> of the human genome: exons make up only ~2% of it.

---

## 1.2 Compute Variant Allele Frequency (VAF)

The TSV does not have a VAF column. We need to calculate it ourselves from
the **DP4** field, which stores the raw read counts broken down by allele
and strand direction:

$$VAF = \frac{alt\_fwd + alt\_rev}{ref\_fwd + ref\_rev + alt\_fwd + alt\_rev}$$

The DP4 column looks like `"92,51,56,13"` — four comma-separated numbers per
variant. We can parse all rows at once using `read.csv()` with the `text`
argument, which treats each string as a line of CSV input.

```{r compute-vaf}
# Parse DP4 into four named columns — one row per variant
dp4 <- read.csv(text = germ$DP4, header = FALSE,
                col.names = c("ref_fwd", "ref_rev", "alt_fwd", "alt_rev"))
head(dp4)

# VAF = alt reads divided by total reads
germ$VAF <- (dp4$alt_fwd + dp4$alt_rev) /
            (dp4$ref_fwd + dp4$ref_rev + dp4$alt_fwd + dp4$alt_rev)

summary(germ$VAF)
```

> Some variants have `NaN` (Not a Number) VAF — these are positions where
> DP4 sums to zero, meaning no reads supported either allele in the DP4
> breakdown. We will drop these when plotting.

---

> **Exercise 1.2:** Does the VAF make biological sense? For a heterozygous
> variant (`0/1`), we expect ~50% of reads to carry the alt allele (VAF ≈ 0.5).
> For homozygous (`1/1`) we expect ~100% (VAF ≈ 1.0). Use `tapply()` to
> check the median VAF per genotype.
>
> **Hint:** `tapply(germ$VAF, germ$GT, median, na.rm = TRUE)`

```{r vaf-by-gt, eval=FALSE}
# --- Solution ---
tapply(germ$VAF, germ$GT, median, na.rm = TRUE)
```

> If the medians are close to 0.5 and 1.0, our calculation is correct.
> Variants that fall far outside these peaks usually have very low depth
> (DP), so a single extra read shifts the ratio dramatically.

---

## 1.3 Variant filtering funnel

In clinical genomics, raw variant calls are filtered step by step to
prioritise variants likely to have a real functional impact. Each step has a
clear biological justification — let's apply them one by one and watch the
numbers drop.

```{r filtering-funnel}
# Step 0 — start with all raw variants
v0 <- germ

# Step 1 — quality filter
#   QUAL > 20 : the probability of a wrong call is < 1% (phred scale)
#   DP   > 10 : we need at least 10 reads to trust the call
v1 <- v0[v0$QUAL > 20 & v0$DP > 10, ]

# Step 2 — keep only protein-coding genes
#   Variants in lncRNAs, pseudogenes, etc. are very rarely clinically relevant
v2 <- v1[v1$biotype == "protein_coding", ]

# Step 3 — keep only exonic variants
#   Intronic and UTR variants generally do not change the protein sequence
remove_csq <- c("intron", "3_prime_utr", "5_prime_utr")
v3 <- v2[!v2$Consequence %in% remove_csq, ]

# Step 4 — keep only rare variants
#   gnomad_popmax_af < 0.01  : present in < 1% of the general population
#   gnomad_popmax_af == -1   : never observed in gnomAD (ultra-rare)
#   Common variants are almost certainly benign polymorphisms — not disease-causing
v4 <- v3[v3$gnomad_popmax_af < 0.01 | v3$gnomad_popmax_af == -1, ]

cat(sprintf("All variants              : %6d\n", nrow(v0)))
cat(sprintf("Quality  (QUAL>20, DP>10) : %6d\n", nrow(v1)))
cat(sprintf("Protein-coding genes      : %6d\n", nrow(v2)))
cat(sprintf("Exonic only               : %6d\n", nrow(v3)))
cat(sprintf("Rare  (gnomAD AF < 1%%)   : %6d\n", nrow(v4)))
```

Let's visualise the funnel. Because the range spans from ~88,000 to ~49,
we need a **log10 y-axis** to make all steps visible.

```{r plot-funnel}
funnel_n <- c(
  "All variants"   = nrow(v0),
  "Quality filter" = nrow(v1),
  "Protein-coding" = nrow(v2),
  "Exonic only"    = nrow(v3),
  "Rare (AF<1%)"   = nrow(v4)
)

bp <- barplot(
  funnel_n,
  main      = "Germline variant filtering funnel — HG003 chr22",
  ylab      = "Number of variants (log10 scale)",
  col       = "#005B9F",
  border    = "white",
  log       = "y",
  las       = 2,
  cex.names = 0.85
)
text(x = bp, y = funnel_n, labels = funnel_n, pos = 3, cex = 0.85)
```

> **Discussion:** Which single step removes the most variants? What does that
> tell you about the composition of the human genome?
>
> The answer reveals something important: even within protein-coding genes,
> the vast majority of variants are **intronic** — they lie within gene bodies
> but between exons, and are spliced out before translation.

---

## 1.4 Consequence distribution

Now let's look at the **types** of variants we have in protein-coding genes.
We use `v2` (after quality + coding filter) so we can see the full spectrum
including intronic variants within genes.

```{r plot-consequence-all}
csq_df <- as.data.frame(sort(table(v2$Consequence), decreasing = TRUE))
colnames(csq_df) <- c("consequence", "n")

ggplot(csq_df, aes(x = reorder(consequence, n), y = n)) +
  geom_col(fill = "#005B9F") +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(title    = "Variant consequences — all protein-coding variants",
       subtitle = "HG003 chr22",
       x = NULL, y = "Count") +
  theme_bw(base_size = 13)
```

After all four filters, how many candidate functional variants remain in `v4`?

```{r plot-consequence-rare}
csq_rare_df <- as.data.frame(sort(table(v4$Consequence), decreasing = TRUE))
colnames(csq_rare_df) <- c("consequence", "n")

ggplot(csq_rare_df, aes(x = reorder(consequence, n), y = n,
                         fill = consequence)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(title    = "Rare exonic variant consequences",
       subtitle = "After all four filters — candidate functional variants",
       x = NULL, y = "Count") +
  theme_bw(base_size = 13)
```

> **Missense** variants change one amino acid to another — they may or may not
> affect protein function. **Loss-of-function (LoF)** variants (stop_gained,
> frameshift, splice_donor, splice_acceptor) are more likely to disrupt the
> protein completely.

How many of each do we have?

```{r lof-count}
lof_types  <- c("stop_gained", "frameshift", "splice_donor", "splice_acceptor")
n_missense <- sum(v4$Consequence == "missense", na.rm = TRUE)
n_lof      <- sum(v4$Consequence %in% lof_types, na.rm = TRUE)

cat(sprintf("Rare missense variants         : %d\n", n_missense))
cat(sprintf("Rare loss-of-function variants : %d\n", n_lof))
```

---

> **Exercise 1.3:** Which genes carry the rare loss-of-function variants?
> These are the most interesting candidates for disease-associated mutations.
>
> **Hint:** Use `subset(v4, Consequence %in% lof_types)` to filter, then
> `table()` on the `gene` column.

```{r lof-genes-exercise, eval=FALSE}
# --- Solution ---
lof_genes <- subset(v4, Consequence %in% lof_types)
sort(table(lof_genes$gene), decreasing = TRUE)
```

---

## 1.5 Ti/Tv ratio

Every single nucleotide variant (SNV) is either a **transition** or a
**transversion**:

- **Transition (Ti):** A↔G or C↔T — a change within the same class (both purines, or both pyrimidines). These are chemically more likely to occur.
- **Transversion (Tv):** all other changes — a purine mutates to a pyrimidine or vice versa.

There are 4 possible transitions and 8 possible transversions, so by random
chance alone the Ti/Tv ratio would be 0.5. In reality, biological processes
favour transitions:

- **Whole genome:** Ti/Tv ≈ 2.0
- **Exome:** Ti/Tv ≈ 2.8–3.0 (even more transitions in coding sequence)

A Ti/Tv ratio much below 2.0 is a red flag for technical artefacts or
contamination.

We define a helper function that classifies any SNV:

```{r titv-function}
classify_titv <- function(ref, alt) {
  ti_pairs <- c("A>G", "G>A", "C>T", "T>C")
  change   <- paste0(toupper(ref), ">", toupper(alt))
  ifelse(change %in% ti_pairs, "Transition", "Transversion")
}
```

Now apply it to the germline SNVs (single-base changes only — indels are excluded):

```{r titv-germline}
# Keep only SNVs: REF and ALT must each be exactly one base
snvs_germ <- germ[nchar(germ$REF) == 1 & nchar(germ$ALT) == 1, ]

snvs_germ$titv         <- classify_titv(snvs_germ$REF, snvs_germ$ALT)
snvs_germ$substitution <- paste0(toupper(snvs_germ$REF), ">",
                                  toupper(snvs_germ$ALT))

titv_counts <- table(snvs_germ$titv)
titv_ratio  <- titv_counts["Transition"] / titv_counts["Transversion"]
cat(sprintf("Germline Ti/Tv ratio: %.2f\n", titv_ratio))
```

Let's look at all 12 possible substitution types to see which are most common:

```{r plot-titv-germline}
sub_counts <- sort(table(snvs_germ$substitution), decreasing = TRUE)
sub_df <- data.frame(
  substitution = names(sub_counts),
  n            = as.numeric(sub_counts),
  type         = classify_titv(
    sub("(.)>.", "\\1", names(sub_counts)),
    sub(".>(.)", "\\1", names(sub_counts))
  )
)

ggplot(sub_df, aes(x = reorder(substitution, -n), y = n, fill = type)) +
  geom_col() +
  geom_text(aes(label = scales::comma(n)), vjust = -0.4, size = 3.5) +
  scale_fill_manual(
    values = c("Transition" = "#E74C3C", "Transversion" = "#3498DB")
  ) +
  labs(title    = sprintf("Germline substitution spectrum  —  Ti/Tv = %.2f",
                          titv_ratio),
       subtitle = "HG003 chr22 SNVs",
       x = "Substitution", y = "Count", fill = NULL) +
  theme_bw(base_size = 13)
```

> C→T and G→A transitions dominate — these are the most common germline
> mutations because cytosine spontaneously deaminates to uracil (which is
> read as thymine), producing C→T changes.

---

## 1.6 Top genes by variant burden

Which genes on chr22 carry the most germline variants in this individual?
High variant burden in a gene is not necessarily bad — it could reflect the
gene's length, its tolerance for variation, or a recently duplicated region.

```{r plot-top-genes-germline}
# Remove intergenic rows where gene is "." (not in any gene)
v2_genes    <- v2[v2$gene != "." & v2$gene != "", ]
gene_counts <- sort(table(v2_genes$gene), decreasing = TRUE)
top20       <- head(as.data.frame(gene_counts), 20)
colnames(top20) <- c("gene", "n")

ggplot(top20, aes(x = reorder(gene, n), y = n)) +
  geom_col(fill = "#005B9F") +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(title    = "Top 20 genes by germline variant count",
       subtitle = "Protein-coding variants, HG003 chr22",
       x = NULL, y = "Number of variants") +
  theme_bw(base_size = 13)
```

---

# Part 2: Somatic variant analysis

## 2.1 Load and explore

Now let's switch to the **somatic** calls — variants present in the tumour
but absent from the matched normal. These were called by Mutect2 in a
tumour-normal paired analysis.

```{r load-somatic}
# The somatic TSV has a trailing tab on every data row (a bcftools quirk),
# giving 20 data columns against 19 header columns.
# We provide explicit column names to handle this cleanly.
som_cols <- c("CHROM","POS","REF","ALT","FILTER","gene","transcript","biotype",
              "Consequence","amino_acid_change","dna_change","gnomad_popmax_af",
              "gnomad_filter","normal_GT","normal_AF","normal_DP",
              "tumor_GT","tumor_AF","tumor_DP")

som <- read.delim(SOMATIC_TSV, sep = "\t", header = FALSE, skip = 1,
                  col.names = c(som_cols, "extra"),
                  stringsAsFactors = FALSE)
som$extra <- NULL   # drop the trailing-tab artefact column

dim(som)
head(som)
```

> Unlike germline variants, somatic calls have a `FILTER` column with
> Mutect2's assessment of call quality. `PASS` = high confidence. Other labels
> (e.g. `weak_evidence`, `strand_bias`, `normal_artifact`) are reasons a call
> was flagged as potentially unreliable.

---

> **Exercise 2.1:** What FILTER labels did Mutect2 assign, and how many
> variants have each label? Sort from most to least common.
>
> **Hint:** `sort(table(som$FILTER), decreasing = TRUE)`

```{r somatic-filters, eval=FALSE}
# --- Solution ---
sort(table(som$FILTER), decreasing = TRUE)
```

---

## 2.2 Filtering somatic variants

First, keep only `PASS` variants — these are Mutect2's high-confidence calls.

```{r somatic-filter-pass}
som_pass <- som[som$FILTER == "PASS", ]
cat(sprintf("PASS variants : %d / %d total\n", nrow(som_pass), nrow(som)))
```

Next, remove common germline polymorphisms that may have slipped through the
normal sample filter. We use the same gnomAD logic as for germline.

```{r somatic-filter-rare}
som_rare <- som_pass[
  som_pass$gnomad_popmax_af < 0.01 | som_pass$gnomad_popmax_af == -1, ]
cat(sprintf("PASS + rare   : %d\n", nrow(som_rare)))
```

---

## 2.3 Somatic consequence distribution

What types of functional changes do the somatic mutations cause?

```{r plot-somatic-consequences}
som_csq_df <- as.data.frame(
  sort(table(som_rare$Consequence[som_rare$Consequence != "."]),
       decreasing = TRUE)
)
colnames(som_csq_df) <- c("consequence", "n")

ggplot(som_csq_df, aes(x = reorder(consequence, n), y = n,
                        fill = consequence)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(title    = "Somatic variant consequences  (PASS, rare)",
       subtitle = "Mutect2 tumour-normal calls",
       x = NULL, y = "Number of variants") +
  theme_bw(base_size = 13)
```

---

## 2.4 Ti/Tv ratio — somatic

We already defined `classify_titv()` in Section 1.5. Now let's apply it to
the somatic mutations. The Ti/Tv ratio in tumours can differ dramatically from
germline — different mutational processes leave characteristic signatures.

---

> **Exercise 2.2:** Compute the Ti/Tv ratio for somatic PASS SNVs.
>
> 1. Filter `som_pass` to SNVs only: `nchar(REF) == 1 & nchar(ALT) == 1`
> 2. Add a `titv` column using `classify_titv()`
> 3. Compute `table()` and divide Transitions by Transversions
>
> **Hint:** The steps are identical to Section 1.5 — just change the data frame name.

```{r titv-somatic, eval=FALSE}
# --- Solution ---
snvs_som <- som_pass[nchar(som_pass$REF) == 1 & nchar(som_pass$ALT) == 1, ]

snvs_som$titv         <- classify_titv(snvs_som$REF, snvs_som$ALT)
snvs_som$substitution <- paste0(toupper(snvs_som$REF), ">",
                                 toupper(snvs_som$ALT))

som_titv_counts <- table(snvs_som$titv)
som_titv_ratio  <- som_titv_counts["Transition"] / som_titv_counts["Transversion"]
cat(sprintf("Somatic Ti/Tv ratio: %.2f\n", som_titv_ratio))
```

Now let's plot the somatic substitution spectrum:

```{r plot-titv-somatic}
# Recompute snvs_som for the plot (in case Exercise 2.2 was not run)
snvs_som <- som_pass[nchar(som_pass$REF) == 1 & nchar(som_pass$ALT) == 1, ]
snvs_som$titv         <- classify_titv(snvs_som$REF, snvs_som$ALT)
snvs_som$substitution <- paste0(toupper(snvs_som$REF), ">", toupper(snvs_som$ALT))
som_titv_counts <- table(snvs_som$titv)
som_titv_ratio  <- som_titv_counts["Transition"] / som_titv_counts["Transversion"]

som_sub_counts <- sort(table(snvs_som$substitution), decreasing = TRUE)
som_sub_df <- data.frame(
  substitution = names(som_sub_counts),
  n            = as.numeric(som_sub_counts),
  type         = classify_titv(
    sub("(.)>.", "\\1", names(som_sub_counts)),
    sub(".>(.)", "\\1", names(som_sub_counts))
  )
)

ggplot(som_sub_df, aes(x = reorder(substitution, -n), y = n, fill = type)) +
  geom_col() +
  geom_text(aes(label = n), vjust = -0.4, size = 3.5) +
  scale_fill_manual(
    values = c("Transition" = "#E74C3C", "Transversion" = "#3498DB")
  ) +
  labs(title    = sprintf("Somatic substitution spectrum  —  Ti/Tv = %.2f",
                          som_titv_ratio),
       subtitle = "PASS somatic SNVs",
       x = "Substitution", y = "Count", fill = NULL) +
  theme_bw(base_size = 13)
```

> **Discussion:** The somatic Ti/Tv is much lower than the germline ratio.
> This tumour is dominated by **transversions** — particularly C→A changes.
> C→A (or G→T on the reverse strand) is the hallmark of **oxidative DNA
> damage** (8-oxoguanine). Compare this to the germline, where C→T
> transitions dominate due to spontaneous cytosine deamination.
> Different mutational processes leave different fingerprints.

---

## 2.5 Normal AF vs Tumour AF

One of the most powerful quality checks for somatic calls is to plot the
allele frequency in the normal sample against the allele frequency in the
tumour. Where do we expect real somatic mutations to fall?

- **True somatic mutations:** high tumour AF, near-zero normal AF — the tumour has the mutation, the normal does not
- **Germline variants leaking through:** both AFs similar (~0.5 or ~1.0) — present in both samples
- **Sequencing artefacts:** low AF in both — likely noise

```{r plot-normal-vs-tumor}
ggplot(som_pass, aes(x = normal_AF, y = tumor_AF)) +
  geom_point(alpha = 0.4, color = "#E74C3C", size = 1.8) +
  geom_abline(intercept = 0, slope = 1,
              linetype = "dashed", color = "grey40") +
  labs(title    = "Normal AF vs Tumour AF  —  PASS variants",
       subtitle = "Points near x = 0 are true somatic. Points on the diagonal may be germline.",
       x = "Normal allele frequency",
       y = "Tumour allele frequency") +
  theme_bw(base_size = 13)
```

> The dashed line is where tumour AF = normal AF (the diagonal). Points
> close to the left edge (normal AF ≈ 0) with varying tumour AF are the
> true somatic mutations we are interested in.

---

## 2.6 Tumour VAF by consequence

Are high-VAF somatic mutations enriched in certain consequence types?

A variant's VAF in the tumour reflects when the mutation arose during tumour
evolution:

- **High VAF (> 0.3):** likely a **clonal** mutation — present in all tumour cells, arose early
- **Low VAF (< 0.1):** likely a **subclonal** mutation — present in only a fraction of cells, arose late

```{r plot-vaf-by-consequence}
keep_csq <- c("missense", "synonymous", "intron", "stop_gained",
              "3_prime_utr", "5_prime_utr", "splice_region")
som_sub  <- som_pass[som_pass$Consequence %in% keep_csq, ]

ggplot(som_sub,
       aes(x = reorder(Consequence, tumor_AF, median), y = tumor_AF)) +
  geom_boxplot(fill = "#E74C3C", alpha = 0.6, outlier.size = 0.8) +
  geom_jitter(width = 0.2, alpha = 0.25, size = 0.8) +
  coord_flip() +
  labs(title    = "Tumour VAF by consequence",
       subtitle = "Ordered by median VAF — higher = more clonal",
       x = NULL, y = "Tumour VAF") +
  theme_bw(base_size = 13)
```

---

## 2.7 Top mutated genes — somatic

Which genes are most frequently mutated in this tumour?

---

> **Exercise 2.3:** Find the top 20 most mutated genes in `som_rare`.
> Filter out intergenic rows (`gene != "."`), count variants per gene with
> `table()`, and plot as a horizontal bar chart.
>
> **Hint:** The pattern is the same as Section 1.6 — just use `som_rare`
> instead of `v2_genes`, and colour bars `"#E74C3C"` to distinguish somatic
> from germline.

```{r plot-top-genes-somatic, eval=FALSE}
# --- Solution ---
som_genes       <- som_rare[som_rare$gene != "." & som_rare$gene != "", ]
som_gene_counts <- sort(table(som_genes$gene), decreasing = TRUE)
som_top20       <- head(as.data.frame(som_gene_counts), 20)
colnames(som_top20) <- c("gene", "n")

ggplot(som_top20, aes(x = reorder(gene, n), y = n)) +
  geom_col(fill = "#E74C3C") +
  geom_text(aes(label = n), hjust = -0.2, size = 3.5) +
  coord_flip() +
  labs(title    = "Top mutated genes — somatic (PASS, rare)",
       subtitle = "Mutect2 tumour-normal",
       x = NULL, y = "Number of somatic variants") +
  theme_bw(base_size = 13)
```

---

# Part 3: Germline vs Somatic VAF comparison

## 3.1 The key question

This is the conceptual centrepiece of the day. A germline variant is present
in **every cell** from conception — it was inherited. A somatic mutation arose
in **one cell** during the person's lifetime and expanded as that cell divided.

This single difference in origin leads to a completely different VAF signature:

| Variant type | Expected VAF | Reason |
|---|---|---|
| Germline heterozygous | ~0.50 | One of two alleles is mutant |
| Germline homozygous | ~1.00 | Both alleles are mutant |
| Somatic clonal | 0.2–0.5 | All tumour cells carry it, but normal cells do not |
| Somatic subclonal | 0.01–0.15 | Only a fraction of tumour cells carry it |

Let's stack both datasets into one data frame and plot them together.

```{r vaf-comparison-data}
germ_vaf <- data.frame(VAF  = germ$VAF[is.finite(germ$VAF)],
                       type = "Germline")
som_vaf  <- data.frame(VAF  = som_pass$tumor_AF,
                       type = "Somatic")
vaf_df   <- rbind(germ_vaf, som_vaf)
```

```{r plot-vaf-comparison, fig.height=7}
ggplot(vaf_df, aes(x = VAF, fill = type)) +
  geom_histogram(binwidth = 0.02, color = "white", linewidth = 0.15) +
  scale_fill_manual(
    values = c("Germline" = "#005B9F", "Somatic" = "#E74C3C")
  ) +
  facet_wrap(~type, ncol = 1, scales = "free_y") +
  labs(title    = "VAF distribution: Germline vs Somatic",
       subtitle = "Germline is bimodal (het = 0.5, hom = 1.0).  Somatic is low and broad.",
       x = "Variant Allele Frequency (VAF)", y = "Number of variants",
       fill = NULL) +
  theme_bw(base_size = 13) +
  theme(legend.position = "none",
        strip.text      = element_text(size = 12, face = "bold"))
```

> The germline plot shows two sharp peaks — this is the expected diploid
> signature. The somatic plot shows a broad, low distribution — mutations
> are present at variable frequencies because they arose at different times
> during tumour growth and are therefore present in different proportions
> of cells.

## 3.2 Overlay (density)

A density plot overlays both distributions on the same axis, making the
contrast even clearer.

```{r plot-vaf-density}
ggplot(vaf_df, aes(x = VAF, fill = type, color = type)) +
  geom_density(alpha = 0.4, linewidth = 0.8) +
  scale_fill_manual(values  = c("Germline" = "#005B9F", "Somatic" = "#E74C3C")) +
  scale_color_manual(values = c("Germline" = "#003F72", "Somatic" = "#C0392B")) +
  labs(title = "VAF density: Germline vs Somatic",
       x = "Variant Allele Frequency (VAF)", y = "Density",
       fill = NULL, color = NULL) +
  theme_bw(base_size = 13)
```

---

> **Exercise 3.1:** Redo the VAF comparison using **quality-filtered** germline
> variants (`v1`) instead of all raw variants. Does the bimodal shape at
> 0.5 and 1.0 become sharper and cleaner?
>
> **Hint:** Replace `germ$VAF` with `v1$VAF` when building `germ_vaf`.
> The rest of the code is identical.

```{r vaf-exercise, eval=FALSE}
# --- Solution ---
germ_vaf_filtered <- data.frame(VAF  = v1$VAF[is.finite(v1$VAF)],
                                type = "Germline")
vaf_df2 <- rbind(germ_vaf_filtered, som_vaf)

ggplot(vaf_df2, aes(x = VAF, fill = type)) +
  geom_histogram(binwidth = 0.02, color = "white", linewidth = 0.15) +
  scale_fill_manual(
    values = c("Germline" = "#005B9F", "Somatic" = "#E74C3C")
  ) +
  facet_wrap(~type, ncol = 1, scales = "free_y") +
  labs(title    = "VAF distribution: Quality-filtered Germline vs Somatic",
       x = "Variant Allele Frequency (VAF)", y = "Number of variants") +
  theme_bw(base_size = 13) +
  theme(legend.position = "none",
        strip.text = element_text(size = 12, face = "bold"))
```

---

# Session info

```{r session-info}
sessionInfo()
```