09-MutationsSwiss.Rmd

```{r setupmutationsswiss, include=FALSE}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(knitr)
library(tidyverse)
theme_set(bayesplot::theme_default())
opts_chunk$set(
  echo = F, message = F, warning = F, fig.height = 6, fig.width = 8,
  cache = T, cache.lazy = F)
```

```{r datamutationsswiss}
snv <- read_tsv("data/mutations/swiss/napoleon_mutations.tsv")
genome <- read_tsv("data/mutations/swiss/Qrob_PM1N.fa.fai", 
         col_names = c("CHROM", "length", "bytesindex", "basesperline", "bytesperline")) %>% 
  dplyr::select(CHROM, length) %>% 
  mutate(start = 0, stop = length/10^6)
```

```{r, cache=FALSE}
strelka2raw <- src_sqlite("data/mutations/swiss/strelka2_raw.sql") %>% tbl("mutations")
gatkraw <- src_sqlite("data/mutations/swiss/gatk_raw.sql") %>% tbl("mutations")
```

# Mutations Swiss

This chapter describes the reanalyses of data from @Schmid-Siegert2017 currently done in the [`swiss` branch of the `detectMutations` repository](https://github.com/sylvainschmitt/detectMutations/tree/swiss).

## Mutations from Schmid-Siegert on 3P

I reported the mutations (Tab. \@ref(tab:napoTab)) from the [supplementary table 2](https://static-content.springer.com/esm/art%3A10.1038%2Fs41477-017-0066-9/MediaObjects/41477_2017_66_MOESM1_ESM.pdf) from @Schmid-Siegert2017, and aligned them back on the 3P genome (Fig. \@ref(fig:napoFig)). 
I found back only 14 of the original 17 mutations from Napoleon

```{r napoTab}
kable(snv, caption = "SNVs in the Napoleon Oak. Rerported from @Schmid-Siegert2017.")
```

```{r napoFig, fig.cap="Napoleon's original mutations on the 3P genome."}
ggplot(genome, aes(x = CHROM, xend = CHROM)) +
  geom_segment(aes(y = start, yend = stop), size = 3, col = "grey") +
  geom_point(aes(y = POS/10^6), col = "red", data = snv) +
  ggrepel::geom_label_repel(aes(y = POS/10^6, label = Mutation), col = "red", data = snv) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), 
        axis.ticks.y = element_blank())
```

```{r covswiss2, eval=F}
cov <- list(lower = "data/mutations/swiss/lower.quantized.bed",
     upper = "data/mutations/swiss/upper.quantized.bed") %>% 
  lapply(vroom::vroom, col_names = c("chrom", "start", "end", "coverage")) %>% 
  bind_rows(.id = "branch") %>% 
  mutate(length = end - start) %>% 
  group_by(branch, coverage) %>% 
  summarise(length = sum(length)) %>% 
  group_by(branch) %>% 
  mutate(total = sum(length)) %>% 
  ungroup() %>% 
  mutate(pct = length/total)
cov %>% 
  filter(coverage == "30:120")
```

## `Strelka2`

`Strelka2` produced `r round(collect(tally(strelka2raw))$n/10^6, 1)` millions of candidate mutations.

### Overlap with mutations from Schmid-Siegert

I tried to find back Napoleon's original mutations to have a look to their metrics.
I found back only 12 out of the 14 expected mutations (86%) (Tab. \@ref(tab:napoMutTab) and Fig. \@ref(fig:napoMutFig)).
**Beware, `Strelka2` is detecting putative mutations in the normal sample !**
I looked at different metrics for each  (Fig. \@ref(fig:overlapMetrics)):

* `mutation_DP` and `normal_DP` are the read depth for the two sample, and shows as expected values between half and two times the mean coverage (60X)
* `normal_altCountT1` is the number of alternate allele count in the normal sample, should be 0, but is equal to 3 and 4 (9% of reads) for two SNVs
* `mutation_altCountT1` is the number of alternate allele count in the mutated sample, should be not too low, and is most the time above 5

**The main conclusion is that the mutations detected by @Schmid-Siegert2017 have not always no reads in the "normal" sample and that they show a wide variation of allelic frequency.**

```{r napoMut}
# overlapS <- lapply(1:nrow(snv), function(i)
#   strelka2raw %>%
#     filter(CHROM == local(snv[i,]$CHROM)) %>%
#     filter(POS == local(snv[i,]$POS)) %>%
#     collect() %>%
#     mutate(Mutation  = snv[i,]$Mutation)
# ) %>% bind_rows()
# write_tsv(overlapS, file = "save/overlap_strelka2.tsv")
overlapS <- read_tsv("save/overlap_strelka2.tsv")
overlapSfiltered <- overlapS %>% 
  group_by(Mutation) %>% 
  filter(normal_altCountT1 == min(normal_altCountT1))
```

```{r napoMutTab}
overlapS %>% 
  dplyr::select(Mutation, tumor, normal, REF, ALT, mutation_altCountT2, mutation_refCountT2,
                normal_altCountT1, normal_refCountT1, mutation_AF) %>% 
  kable(caption = "Overlap between candidate mutations and Napoleon's original mutations.", 
        col.names = c("Mutation", "Mutated", "Normal", "Ref", "Alt", 
                      "Mutated\nAltCount", "Mutated\nRefCount", "Normal\nAltCount", "Normal\nRefCount", "Allelic fraction"))
```

```{r napoMutFig, fig.cap="Overlap between candidate mutations and Napoleon's original mutations: allele frequency (A) and positions on the 3P genome (B)."}
g1 <- ggplot(overlapSfiltered, aes(mutation_AF)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(overlapSfiltered))) +
  xlab("Allele frequency")
g2 <- ggplot(genome, aes(x = CHROM, xend = CHROM)) +
  geom_segment(aes(y = start, yend = stop), size = 3, col = "lightgrey") +
  geom_point(aes(y = POS/10^6, col = mutation_AF), data = overlapSfiltered, size = 1) +
  ggrepel::geom_label_repel(aes(y = POS/10^6, label = Mutation), col = "red", 
                            data = dplyr::select(overlapSfiltered, CHROM, POS, Mutation) %>% unique()) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  viridis::scale_color_viridis(option = "inferno")
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))
```

```{r overlapMetrics, fig.cap="Evaluation of the overlap between candidate mutations and Napoleon's original mutations."}
overlapSfiltered %>% 
  group_by(Mutation) %>% 
  filter(normal_altCountT1 == min(normal_altCountT1)) %>% 
  dplyr::select(Mutation, FILTER, mutation_DP, normal_DP,
                mutation_altCountT1, normal_altCountT1, QSS) %>% 
  reshape2::melt(c("Mutation", "FILTER")) %>% 
  ggplot(aes(value, fill = FILTER)) +
  geom_histogram() +
  facet_wrap(~ variable, scales = "free") +
  theme(legend.position = c(0.8, 0.2)) +
  geom_vline(aes(xintercept = value, fill = NA), size = 1.1, col = "black", linetype = "dashed",
             data = data.frame(
               value = c(30-0.5, 120+0.5, 30-0.5, 120+0.5, 10-0.5, 0+0.5, 20),
               variable = c("mutation_DP", "mutation_DP", 
                            "normal_DP", "normal_DP", 
                            "mutation_altCountT1", "normal_altCountT1", "QSS")
             ))
```

### Filtering

We filtered mutations with following filters:

* A read depth for the two sample between half and two times the mean coverage (`normal_DP <= 120, normal_DP >= 30, mutation_DP <= 120, mutation_DP >= 30`)
* A null number of alternate allele count in the normal sample (`normal_altCount == 0`)
* A minimum of 10 alternate allele count in the mutated sample (`mutation_altCount >= 10`)

We obtained 223 candidates  (Fig \@ref(fig:mutFilteredS)).
We then used the suggested automatic filter of `Strelka2`,
resulting in a robust dataset of 87 mutations (Fig \@ref(fig:mutRobustS)).

```{r mutS}
mutFilteredS <- strelka2raw %>% 
  filter(normal_DP <= 120, normal_DP >= 30) %>% 
  filter(mutation_DP <= 120, mutation_DP >= 30) %>% 
  filter(normal_altCountT1 == 0, normal_altCountT2 == 0) %>% 
  filter(mutation_altCountT1 >= 10, mutation_altCountT2 >= 10) %>% 
  collect() 
mutRobustS <- filter(mutFilteredS, FILTER == "PASS")  
```

```{r mutFilteredS, fig.cap="Mutations retained after original filtering: allele frequency (A) and positions on the 3P genome (B)."}
g1 <- ggplot(mutFilteredS, aes(mutation_AF)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(mutFilteredS))) +
  xlab("Allele frequency")
g2 <- ggplot(genome, aes(x = CHROM, xend = CHROM)) +
  geom_segment(aes(y = start, yend = stop), size = 3, col = "lightgrey") +
  geom_point(aes(y = POS/10^6, col = mutation_AF), data = mutFilteredS, size = 1) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  viridis::scale_color_viridis(option = "inferno")
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))  
```

```{r mutRobustS, fig.cap="Mutations retained after robust filtering: allele frequency (A) and positions on the 3P genome (B)."}
g1 <- ggplot(mutRobustS, aes(mutation_AF)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(mutRobustS))) +
  xlab("Allele frequency")
g2 <- ggplot(genome, aes(x = CHROM, xend = CHROM)) +
  geom_segment(aes(y = start, yend = stop), size = 3, col = "lightgrey") +
  geom_point(aes(y = POS/10^6, col = mutation_AF), data = mutRobustS, size = 1) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  viridis::scale_color_viridis(option = "inferno")
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))  
```

## `GATK`

`GATK` produced `r round(collect(tally(gatkraw))$n/10^6, 1)` millions of candidates!

### Overlap with mutations from Schmid-Siegert

I tried to find back Napoleon's original mutations to have a look to their metrics.
I found back only 12 out of the 14 expected mutations (86%) (Tab. \@ref(tab:napoMutTab2) and Fig. \@ref(fig:napoMutFig2)).
I looked at different metrics for each  (Fig. \@ref(fig:overlapMetrics)):

* `mutation_DP` and `normal_DP` are the read depth for the two sample, and shows as expected values between half and two times the mean coverage (60X)
* `normal_altCountT1` is the number of alternate allele count in the normal sample, should be and is 0
* `mutation_altCountT1` is the number of alternate allele count in the mutated sample, should be not too low, and is most the time above 10

**The main conclusion is that the mutations detected by @Schmid-Siegert2017 have no reads in the "normal" sample using `GATK` with hard filtering which probably already removed low-DP copies in the normal sample, while `Strelka2` detect them.**


```{r napoMut2}
# overlapG <- lapply(1:nrow(snv), function(i)
#   gatkraw %>%
#     filter(CHROM == local(snv[i,]$CHROM)) %>%
#     filter(POS == local(snv[i,]$POS)) %>%
#     collect() %>%
#     mutate(Mutation  = snv[i,]$Mutation)
# ) %>% bind_rows()
# write_tsv(overlapG, file = "save/overlap_gatk.tsv")
overlapG <- read_tsv("save/overlap_gatk.tsv")
```

```{r napoMutTab2}
overlapG %>%
  dplyr::select(Mutation, tumor, normal, tumor_altCount, tumor_refCount, normal_altCount, normal_refCount, tumor_AF) %>%
  kable(caption = "Overlap between candidate mutations and Napoleon's original mutations.",
        col.names = c("Mutation", "Mutated", "Normal",
                      "Mutated\nAltCount", "Mutated\nRefCount", "Normal\nAltCount", "Normal\nRefCount", "Allelic fraction"))
```

```{r napoMutFig2, fig.cap="Overlap between candidate mutations and Napoleon's original mutations: allele frequency (A) and positions on the 3P genome (B)."}
g1 <- ggplot(overlapG, aes(tumor_AF)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(overlapG))) +
  xlab("Allele frequency")
g2 <- ggplot(genome, aes(x = CHROM, xend = CHROM)) +
  geom_segment(aes(y = start, yend = stop), size = 3, col = "lightgrey") +
  geom_point(aes(y = POS/10^6, col = tumor_AF), data = overlapG, size = 1) +
  ggrepel::geom_label_repel(aes(y = POS/10^6, label = Mutation), col = "red",
                            data = dplyr::select(overlapG, CHROM, POS, Mutation) %>% unique()) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  viridis::scale_color_viridis(option = "inferno")
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))
```

```{r overlapMetrics2, fig.cap="Evaluation of the overlap between candidate mutations and Napoleon's original mutations."}
overlapG %>%
  group_by(Mutation) %>%
  dplyr::select(Mutation, tumor_DP, normal_DP, tumor_altCount, normal_altCount) %>%
  reshape2::melt(c("Mutation")) %>%
  ggplot(aes(value)) +
  geom_histogram() +
  facet_wrap(~ variable, scales = "free") +
  geom_vline(aes(xintercept = value), size = 1.1, col = "black", linetype = "dashed",
             data = data.frame(
               value = c(30-0.5, 120+0.5, 30-0.5, 120+0.5, 10-0.5, 0+0.5),
               variable = c("tumor_DP", "tumor_DP",
                            "normal_DP", "normal_DP",
                            "tumor_altCount", "normal_altCount")
             ))
```

### Filtering

We filtered mutations with following filters:

* A read depth for the tumor sample between half and two times the mean coverage (`tumor_DP <= 120, tumor_DP >= 30`)
* A null number of alternate allele count in the normal sample (`normal_altCount == 0`)
* A minimum of 10 alternate allele count in the mutated sample (`tumor_altCount >= 10`)
* An allelic frequency  inferior to 0.5 (`tumor_AF <= 0.5`)

We obtained 510 611 candidates (Fig \@ref(fig:mutFilteredG)).
We then looked for the overlap between `GATK` candidates and the suggested automatic filter of `Strelka2`,
resulting in a robust dataset of 47 mutations (Fig \@ref(fig:mutRobustG)).

```{r mutG}
mutFilteredG <- gatkraw %>% 
  filter(tumor_DP <= 120, tumor_DP >= 30) %>% 
  filter(normal_altCount == 0) %>% 
  filter(tumor_altCount >= 10) %>% 
  filter(tumor_AF <= 0.5) %>% 
  collect() 
mutRobustG <- left_join(mutFilteredG,
          dplyr::select(mutRobustS, CHROM, POS) %>%
            mutate(Robust = 1)) %>%
  mutate(Robust = as.factor(ifelse(is.na(Robust), 0, 1))) %>% 
  filter(Robust == 1) 
```

```{r mutFilteredG, fig.cap="Mutations retained after original filtering: allele frequency (A) and positions on the 3P genome (B)."}
g1 <- ggplot(mutFilteredG, aes(tumor_AF)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(mutFilteredG))) +
  xlab("Allele frequency")
g2 <- ggplot(genome, aes(x = CHROM, xend = CHROM)) +
  geom_segment(aes(y = start, yend = stop), size = 3, col = "lightgrey") +
  geom_point(aes(y = POS/10^6, col = tumor_AF), data = filter(mutFilteredG, grepl("Chr", CHROM)), size = 1) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  viridis::scale_color_viridis(option = "inferno")
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B")) 
```

```{r mutRobustG, fig.cap="Mutations retained after robust filtering: allele frequency (A) and positions on the 3P genome (B)."}
g1 <- ggplot(mutRobustG, aes(tumor_AF)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(mutRobustG))) +
  xlab("Allele frequency")
g2 <- ggplot(genome, aes(x = CHROM, xend = CHROM)) +
  geom_segment(aes(y = start, yend = stop), size = 3, col = "lightgrey") +
  geom_point(aes(y = POS/10^6, col = tumor_AF), data = mutRobustG, size = 1) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  viridis::scale_color_viridis(option = "inferno")
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))
```

## Conclusion

```{r datasetsswiss}
datasets <- list(
  "Schmid" = overlapSfiltered %>% 
    ungroup() %>% 
    mutate(SNV = paste0(CHROM, "_", "Pos", as.integer(POS))) %>%
    dplyr::rename(AF = mutation_AF) %>%
    mutate(type = paste0(REF, "->", ALT)) %>% 
    dplyr::select(CHROM, POS, SNV, AF, type),
  "GATK" = mutFilteredG %>% 
    mutate(SNV = paste0(CHROM, "_", "Pos", as.integer(POS))) %>%
    dplyr::rename(AF = tumor_AF) %>%
    dplyr::select(CHROM, POS, SNV, AF),
  "Strelka2" = mutFilteredS %>% 
    mutate(SNV = paste0(CHROM, "_", "Pos", as.integer(POS))) %>%
    dplyr::rename(AF = mutation_AF) %>%
    mutate(type = paste0(REF, "->", ALT)) %>% 
    dplyr::select(CHROM, POS, SNV, AF, type),
  "Robust" = mutRobustS %>% 
    mutate(SNV = paste0(CHROM, "_", "Pos", as.integer(POS))) %>%
    dplyr::rename(AF = mutation_AF) %>%
     mutate(type = paste0(REF, "->", ALT)) %>% 
    dplyr::select(CHROM, POS, SNV, AF, type)
)
```

```{r datasetsswissTab}
lapply(datasets, nrow) %>%
  bind_rows(.id = "dataset") %>%
  reshape2::melt(variable.name = "Dataset", value.name = "Number of candidates") %>%
  mutate(`Estimated total` = round(`Number of candidates`*100/72)) %>% 
  kable(caption = "Size of the different datasets.", format.args = list(big.mark = " ")) 
```

```{r datasetsswissAF, fig.cap="Allele frequency for the different datasets."}
bind_rows(datasets, .id = "dataset") %>%
  ggplot(aes(AF, col = dataset), fill = NA) +
  geom_density(size = 1.3) +
  scale_color_discrete("Dataset") +
  xlab("Allelic frequency")  
```

```{r datasetsswissVenn2}
ggvenn::ggvenn(lapply(datasets, `[[`, "SNV"), show_percentage = F)  
```


```{r swissMutType}
bind_rows(datasets, .id = "dataset") %>% 
  na.omit() %>% 
  group_by(dataset, type) %>% 
  summarise(N = n()) %>% 
  group_by(dataset) %>% 
  mutate(P = N / sum(N) * 100) %>% 
  ggplot(aes(type, P, fill = dataset)) +
  geom_col(position = "dodge") +
  coord_flip() +
  ylab("Percentage") + xlab(" ") 
```

```{r fig4, eval=F}
bind_rows(datasets, .id = "dataset") %>%
  write_tsv("save/candidates_swiss.tsv")
```