10-MutationsBordeaux.Rmd

```{r setupmutationsfrench, 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)
```

# Mutations Bordeaux

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

```{r datamutationsbdx}
snv <- read_tsv("data/mutations/bordeaux/3P_mutations.tsv")
genome <- read_tsv("data/mutations/bordeaux/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 <- DBI::dbConnect(RSQLite::SQLite(), "data/mutations/bordeaux/strelka2_raw.sql") %>% tbl("mutations")
mutect2raw <- DBI::dbConnect(RSQLite::SQLite(), "data/mutations/bordeaux/mutect2_raw.sql") %>% tbl("mutations")
```

## Original mutations on 3P

I reported the original mutations (Tab. \@ref(tab:bdxTab)) from the [supplementary table 5](https://static-content.springer.com/esm/art%3A10.1038%2Fs41477-018-0172-3/MediaObjects/41477_2018_172_MOESM1_ESM.pdf) from @Plomion2018 and plotted them on the 3P genome (Fig. \@ref(fig:bdxFig)) after realignement from the haplome v2.3 to the public 3P genome.

```{r bdxTab}
kable(head(snv), caption = "SNVs in the 3P Oak. Rerported from @Plomion2018.") 
```

```{r bdxFig, fig.cap="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_text_repel(aes(y = POS/10^6, label = round(POS/10^6)), 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()) 
```

## Coverage

Coverage are around 160X, this is high but strangely lower than expected, besides all data have been used this time. (Fig. \@ref(fig:covbdx)).

```{r covbdx, fig.cap="Coverage for Bordeaux's libraries on chromosomes 1 to 4."}
cov <- list(L1 = "data/mutations/bordeaux/L1.10k.cov",
            L2 = "data/mutations/bordeaux/L2.10k.cov",
            L3 = "data/mutations/bordeaux/L3.10k.cov") %>% 
  lapply(read_tsv, c("CHROM", "start", "stop", "coverage")) %>% 
  bind_rows(.id = "branch") %>% 
  rowwise() %>% 
  mutate(POS = mean(c(start, stop))) %>% 
  filter(grepl("Chr", CHROM))
ggplot(filter(cov, CHROM %in% paste0("Qrob_Chr0", 1:4)), aes(POS/10^6, coverage)) +
  geom_line(col = "grey") +
  facet_grid(branch ~ CHROM, scales = "free_x") +
  ylim(0, 400) +
  xlab("Position (Mb)") +
  geom_hline(aes(yintercept = coverage), col = "red",
             data = group_by(cov, CHROM, branch) %>% 
               filter(CHROM %in% paste0("Qrob_Chr0", 1:4)) %>% 
               summarise(coverage = median(coverage))) 
```


```{r covbdx2, eval=F}
cov <- list(L1 = "data/mutations/bordeaux/L1.quantized.bed",
            L2 = "data/mutations/bordeaux/L2.quantized.bed",
            L3 = "data/mutations/bordeaux/L3.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 == "80:320")
```

## `Strelka2`

```{r, eval=F}
strelka2raw %>% 
  dplyr::select(CHROM, POS) %>% 
  collect() %>% 
  unique() %>% 
  nrow() 
# 3 056 425
```

`Strelka2` produced 3 millions of candidate mutations with unique position across the 3 branches.

### Overlap with original mutations

We found back all the 60 expected mutations (100%) (Tab. \@ref(tab:3PMutTabS)).

```{r 3PMutS}
# overlapSB <- lapply(1:nrow(snv), function(i)
#   strelka2raw %>%
#     filter(CHROM == local(snv[i,]$CHROM)) %>%
#     filter(POS == local(snv[i,]$POS)) %>%
#     collect()
# ) %>% bind_rows()
# write_tsv(overlapSB, file = "save/overlap_strelka2_bdx.tsv")
overlapSB <-  read_tsv(file = "save/overlap_strelka2_bdx.tsv")
overlapSBfiltered <- snv %>% 
  dplyr::select(CHROM, POS, origin) %>% 
  mutate(presence = recode(origin,
                           "3P – between XL1 and XL2" = "L2vsL1_L3vsL1",
                           "3P–between XL1 and XL2" = "L2vsL1_L3vsL1",
                           "3P – between XL2 and XL3" = "L3vsL1_L3vsL2",
                           "3P – between XL2 and L3" = "L3vsL1_L3vsL2",
                           "3P – L1 branch" = "L1vsL2_L1vsL3",
                           "3P – L2 branch" = "L2vsL1_L2vsL3")) %>% 
  separate_rows(presence) %>% 
  separate(presence, c("tumor", "normal"), "vs") %>% 
  left_join(dplyr::select(overlapSB, CHROM, POS, REF, ALT, tumor, normal, mutation_AF)) %>% 
  mutate(mutation_AF = ifelse(mutation_AF > 0.5, 1 - mutation_AF, mutation_AF)) %>% 
  group_by(CHROM, POS, origin, REF, ALT) %>% 
  summarise(mutation_AF = median(mutation_AF, na.rm = T)) %>% 
  filter(!is.nan(mutation_AF)) %>% 
  na.omit()
```

```{r 3PMutTabS}
overlapSB %>% 
  dplyr::select(CHROM, POS, tumor, normal, REF, ALT, mutation_altCountT1, mutation_refCountT1,
                normal_altCountT1, normal_refCountT1, mutation_AF) %>% 
  head() %>% 
  kable(caption = "Overlap between candidate mutations and original mutations.", 
        col.names = c("Chromosome", "Position", "Mutated", "Normal", "Ref", "Alt", 
                      "Mutated\nAltCount", "Mutated\nRefCount", "Normal\nAltCount", "Normal\nRefCount", "Allelic fraction"))
```

### Filtering

We filtered mutations with the original and a robust filters, adjusting the coverage to the coverage obtained:

Original:

* A read depth for the two sample between half and two times the mean coverage (`normal_DP <= 320, normal_DP >= 80, mutation_DP <= 320, mutation_DP >= 80`)
* 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`)
* An allelic frequency  inferior to 0.5 (`tumor_AF <= 0.5`)
* A coherent temporal pattern ("XL1-XL2", "XL2-L3", "XL1-L1", "XL2-L2")

Robust:

* A read depth for the two sample between half and two times the mean coverage (`normal_DP <= 320, normal_DP >= 80, mutation_DP <= 320, mutation_DP >= 80`)
* 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`)
* An allelic frequency  inferior to 0.5 (`tumor_AF <= 0.5`)
* `Strelka 2` automatic filtering (`Filter == PASS`)
* A coherent temporal pattern ("XL1-XL2", "XL2-L3", "XL1-L1", "XL2-L2")

We obtained 460 candidates  (Fig \@ref(fig:mutFilteredSB)) for the original filter,
And in a robust dataset of 108 mutations (Fig \@ref(fig:mutRobustSB)).

```{r mutStrelka2Bdx, eval=F}
mcov <- 160
strelka2raw %>% 
  filter(normal_DP <= mcov*2, normal_DP >= mcov/2) %>%
  filter(mutation_DP <= mcov*2, mutation_DP >= mcov/2) %>% 
  filter(normal_altCountT1 == 0, normal_altCountT2 == 0) %>% 
  filter(mutation_altCountT1 >= 10, mutation_altCountT2 >= 10) %>% 
  mutate(comparison = paste0(tumor, "vs", normal)) %>% 
  collect() %>% 
  group_by(CHROM, POS) %>% 
  mutate(presence = paste(sort(comparison), collapse = "_")) %>% 
  mutate(sector = recode(presence,
                         "L2vsL1_L3vsL1" = "XL1-XL2",
                         "L3vsL1_L3vsL2" = "XL2-L3",
                         "L1vsL2_L1vsL3" = "XL1-L1",
                         "L2vsL1_L2vsL3" = "XL2-L2"
                         )) %>% 
  mutate(sector = ifelse(!(sector %in% c("XL1-XL2", "XL2-L3", "XL1-L1", "XL2-L2")), NA, sector)) %>% 
  filter(!is.na(sector)) %>%
  group_by(CHROM, POS, sector, REF, ALT) %>% 
  # group_by(CHROM, POS, sector) %>% 
  summarise_at("mutation_AF", mean) %>% 
  write_tsv(file = "save/strelka2_bdx_filtered.tsv")
strelka2raw %>% 
  filter(normal_DP <= mcov*2, normal_DP >= mcov/2) %>%
  filter(mutation_DP <= mcov*2, mutation_DP >= mcov/2) %>% 
  filter(normal_altCountT1 == 0, normal_altCountT2 == 0) %>% 
  filter(mutation_altCountT1 >= 10, mutation_altCountT2 >= 10) %>% 
  filter(mutation_AF <= 0.5) %>%
  filter(FILTER == "PASS") %>% 
  mutate(comparison = paste0(tumor, "vs", normal)) %>% 
  collect() %>% 
  group_by(CHROM, POS) %>% 
  mutate(presence = paste(sort(comparison), collapse = "_")) %>% 
  mutate(sector = recode(presence,
                         "L2vsL1_L3vsL1" = "XL1-XL2",
                         "L3vsL1_L3vsL2" = "XL2-L3",
                         "L1vsL2_L1vsL3" = "XL1-L1",
                         "L2vsL1_L2vsL3" = "XL2-L2"
                         )) %>% 
  mutate(sector = ifelse(!(sector %in% c("XL1-XL2", "XL2-L3", "XL1-L1", "XL2-L2")), NA, sector)) %>% 
  filter(!is.na(sector)) %>% 
  group_by(CHROM, POS, sector, REF, ALT) %>% 
  # group_by(CHROM, POS, sector) %>% 
  summarise_at("mutation_AF", mean) %>% 
  write_tsv(file = "save/strelka2_bdx_robust.tsv")
```

```{r mutFilteredSB, fig.cap="Mutations retained after original filtering: allele frequency (A) and positions on the 3P genome (B)."}
strelka2filtered <- read_tsv("save/strelka2_bdx_filtered.tsv") 
g1 <- ggplot(strelka2filtered, aes(mutation_AF, fill = sector)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(strelka2filtered))) +
  xlab("Allele frequency") +
  scale_fill_manual(values = c("darkblue", "plum", "darkgreen", "yellow2"))
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 = sector), size = 1,
             data = filter(strelka2filtered, grepl("Chr", CHROM))) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  scale_color_manual(guide = "none", values = c("darkblue", "plum", "darkgreen", "yellow2"))
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))   
```

```{r mutRobustSB, fig.cap="Mutations retained after robust filtering: allele frequency (A) and positions on the 3P genome (B)."}
strelka2robust <- read_tsv("save/strelka2_bdx_robust.tsv") 
g1 <- ggplot(strelka2robust, aes(mutation_AF, fill = sector)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(strelka2robust))) +
  xlab("Allele frequency") +
  scale_fill_manual(values = c("darkblue", "plum", "darkgreen", "yellow2"))
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 = sector), size = 1,
             data = filter(strelka2robust, grepl("Chr", CHROM))) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  scale_color_manual(guide = "none", values = c("darkblue", "plum", "darkgreen", "yellow2"))
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))   
```

## `Mutect2`

```{r, eval=F}
mutect2raw %>% 
  dplyr::select(CHROM, POS) %>% 
  collect() %>% 
  unique() %>% 
  nrow()
```

`Mutect2` produced 4.6 millions of candidate mutations with unique position across the 3 branches.

### Overlap with original mutations

We tried to find back original mutations to have a look to their metrics.
We found back 57 out of the 60 expected mutations (95%) (Tab. \@ref(tab:3PMutTabM).
which are exactly interestingly the same as `Strelka2`.

```{r 3PMutM}
# overlapSM <- lapply(1:nrow(snv), function(i)
#   mutect2raw %>%
#     filter(CHROM == local(snv[i,]$CHROM)) %>%
#     filter(POS == local(snv[i,]$POS)) %>%
#     collect()
# ) %>% bind_rows()
# write_tsv(overlapSM, file = "save/overlap_mutect2_bdx.tsv")
overlapSM <- read_tsv("save/overlap_mutect2_bdx.tsv")
overlapSMfiltered <- snv %>% 
  dplyr::select(CHROM, POS, origin) %>% 
  mutate(presence = recode(origin,
                           "3P – between XL1 and XL2" = "L2vsL1_L3vsL1",
                           "3P–between XL1 and XL2" = "L2vsL1_L3vsL1",
                           "3P – between XL2 and XL3" = "L3vsL1_L3vsL2",
                           "3P – between XL2 and L3" = "L3vsL1_L3vsL2",
                           "3P – L1 branch" = "L1vsL2_L1vsL3",
                           "3P – L2 branch" = "L2vsL1_L2vsL3")) %>% 
  separate_rows(presence) %>% 
  separate(presence, c("tumor", "normal"), "vs") %>% 
  left_join(dplyr::select(overlapSM, CHROM, POS, tumor, normal, tumor_AF)) %>% 
  mutate(mutation_AF = ifelse(tumor_AF > 0.5, 1 - tumor_AF, tumor_AF)) %>% 
  group_by(CHROM, POS, origin) %>% 
  summarise(tumor_AF = median(tumor_AF, na.rm = T)) %>% 
  filter(!is.na(tumor_AF))
```

```{r 3PMutTabM}
overlapSM %>% 
  dplyr::select(CHROM, POS, tumor, normal, REF, ALT, tumor_altCount, tumor_refCount,
                normal_altCount, normal_refCount, tumor_AF) %>% 
  head() %>% 
  kable(caption = "Overlap between candidate mutations and original mutations.", 
        col.names = c("Chromosome", "Position", "Mutated", "Normal", "Ref", "Alt", 
                      "Mutated\nAltCount", "Mutated\nRefCount", "Normal\nAltCount", "Normal\nRefCount", "Allelic fraction"))
```

### Filtering

We filtered mutations with the original filter, **adjusting the coverage to preliminary low coverage obtained**:

Original:

* A read depth for the two sample between half and two times the mean coverage (`normal_DP <= 320, normal_DP >= 80, mutation_DP <= 320, mutation_DP >= 80`)
* 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`)
* An allelic frequency  inferior to 0.5 (`tumor_AF <= 0.5`)
* A coherent temporal pattern ("XL1-XL2", "XL2-L3", "XL1-L1", "XL2-L2")

We obtained 948 candidates  (Fig \@ref(fig:mutFilteredFB)) for the original filter.

```{r mutMutect2Bdx, eval=F}
mcov <- 160
mutect2raw %>% 
  filter(normal_DP <= mcov*2, normal_DP >= mcov/2) %>% 
  filter(tumor_DP <= mcov*2, tumor_DP >= mcov/2) %>% 
  filter(normal_altCount == 0) %>% 
  filter(tumor_altCount >= 10) %>% 
  filter(tumor_AF <= 0.5) %>% 
  mutate(comparison = paste0(tumor, "vs", normal)) %>% 
  collect() %>% 
  group_by(CHROM, POS) %>% 
  mutate(presence = paste(sort(comparison), collapse = "_")) %>% 
  mutate(sector = recode(presence,
                         "L2vsL1_L3vsL1" = "XL1-XL2",
                         "L3vsL1_L3vsL2" = "XL2-L3",
                         "L1vsL2_L1vsL3" = "XL1-L1",
                         "L2vsL1_L2vsL3" = "XL2-L2"
                         )) %>% 
  mutate(sector = ifelse(!(sector %in% c("XL1-XL2", "XL2-L3", "XL1-L1", "XL2-L2")), NA, sector)) %>% 
  filter(!is.na(sector)) %>% 
  # group_by(CHROM, POS, sector) %>% 
  mutate(REF = stringr::str_sub(REF, -1, -1), ALT = stringr::str_sub(ALT, -1, -1)) %>% 
  group_by(CHROM, POS, sector, REF, ALT) %>% 
  summarise_at("tumor_AF", mean) %>% 
  write_tsv(file = "save/mutect2_bdx_filtered.tsv")
```

```{r mutFilteredFB, fig.cap="Mutations retained after original filtering: allele frequency (A) and positions on the 3P genome (B)."}
mutect2filtered <- read_tsv("save/mutect2_bdx_filtered.tsv") 
g1 <- ggplot(mutect2filtered, aes(tumor_AF, fill = sector)) +
  geom_histogram() +
  ggtitle("", paste("N =", nrow(mutect2filtered))) +
  xlab("Allele frequency") +
  scale_fill_manual(values = c("darkblue", "plum", "darkgreen", "yellow2"))
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 = sector), size = 1,
             data = filter(mutect2filtered, grepl("Chr", CHROM))) +
  coord_flip() +
  ylab("Position (Mb)") +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), axis.ticks.y = element_blank()) +
  scale_color_manual(guide = "none", values = c("darkblue", "plum", "darkgreen", "yellow2"))
cowplot::plot_grid(g1, g2, nrow = 2, labels = c("A", "B"))  
```

## Conclusion

```{r datasetsbdx}
datasets <- list(
  "Plomion" = overlapSBfiltered %>% 
    mutate(SNV = paste0(CHROM, "_", "Pos", as.integer(POS))) %>%
    dplyr::rename(AF = mutation_AF) %>%
    mutate(type = paste0(REF, "->", ALT)) %>% 
    ungroup() %>% 
    dplyr::select(CHROM, POS, SNV, AF, type),
  "Mutect2" = mutect2filtered %>% 
    mutate(SNV = paste0(CHROM, "_", "Pos", as.integer(POS))) %>%
    dplyr::rename(AF = tumor_AF) %>%
    mutate(type = paste0(REF, "->", ALT)) %>% 
    dplyr::select(CHROM, POS, SNV, AF, type),
  "Strelka2" = strelka2filtered %>% 
    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" = strelka2robust %>% 
    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 datasetsbdxTab}
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/84)) %>% 
  kable(caption = "Size of the different datasets.", format.args = list(big.mark = " "))    
```

```{r datasetsbdxAF, 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 datasetsbdxVenn2}
ggvenn::ggvenn(lapply(datasets, `[[`, "SNV"), show_percentage = F)
```

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

## Swiss & Bordeaux

```{r candidatesAll, eval=F}
list("Schmid-Siegert et al., (2017)" = "save/candidates_swiss.tsv", "Plomion et al. (2018)" = "save/candidates_bdx.tsv") %>% 
  lapply(read_tsv) %>% 
  bind_rows(.id = "study") %>% 
  vroom::vroom_write("save/candidates_all.tsv")
```

```{r candidatesAllSpectra, eval=F}
mutations <- read_tsv("save/candidates_all.tsv")
ref <- Biostrings::readDNAStringSet("data/mutations/bordeaux/Qrob_PM1N.fa")
ref <- data.frame(CHROM = names(ref), seq = as.character(ref, use.names = F)) %>% 
  as.tbl()
mutations <- mutations %>% 
  left_join(ref) %>% 
  mutate(POS5 = stringr::str_sub(seq, POS-1, POS-1)) %>%
  mutate(POS3 = stringr::str_sub(seq, POS+1, POS+1)) %>%
  dplyr::select(-seq)
write_tsv(mutations, "save/candidates_all.tsv")
```

```{r candidatesAllannot, eval=F}
library(genomeIntervals)
mutations <- read_tsv("save/candidates_all.tsv")
mutationsI <- new(
  "Genome_intervals_stranded",
  as.matrix(mutations[c("POS", "POS")]),
  closed = T,
  annotation = data.frame(
    seq_name = mutations$CHROM,
    inter_base = F,
    strand = factor("+", levels=c("+", "-") )
  )
)
genes <- genomeIntervals::readGff3("data/mutations/bordeaux/Qrob_PM1N_genes_20161004.gff")
te <- genomeIntervals::readGff3("data/mutations/bordeaux/Qrob_PM1N_refTEs.gff")

mutations$mRNA <- as.numeric(unlist(lapply(interval_overlap(mutationsI, genes[genes$type == "mRNA"]),
                                           function(x) !length(x) == 0)))
mutations$CDS <- as.numeric(unlist(lapply(interval_overlap(mutationsI, genes[genes$type == "CDS"]),
                                           function(x) !length(x) == 0)))
mutations$five_prime_UTR <- as.numeric(unlist(lapply(interval_overlap(mutationsI, genes[genes$type == "five_prime_UTR"]),
                                           function(x) !length(x) == 0)))
mutations$three_prime_UTR <- as.numeric(unlist(lapply(interval_overlap(mutationsI, genes[genes$type == "three_prime_UTR"]),
                                           function(x) !length(x) == 0)))
mutations$te <- as.numeric(unlist(lapply(interval_overlap(mutationsI, te), function(x) !length(x) == 0)))

mutations <- mutations %>% 
  filter(grepl("Chr", CHROM)) %>% 
  filter(dataset %in% c("Robust")) %>% 
  rowwise() %>% 
  mutate(exon = as.numeric(sum(c(CDS, five_prime_UTR, three_prime_UTR)) > 0)) %>% 
  mutate(annotation1 = ifelse(mRNA == 1, "genic", "intergenic")) %>% 
  mutate(annotation2 = ifelse(annotation1 == "genic" & exon == 1, "exon", "intron")) %>% 
  mutate(annotation2 = ifelse(annotation1 == "intergenic" & te == 1, "TE", annotation2)) %>% 
  mutate(annotation2 = ifelse(annotation1 == "intergenic" & te == 0, "non-TE", annotation2)) %>% 
  group_by(study, CHROM, annotation1, annotation2) %>% 
  summarise(N = n()) %>% 
  group_by(study, CHROM) %>% 
  mutate(P = N / sum(N) * 100)

write_tsv(mutations, "save/mut_anot_genom.tsv")

mutations
all_mut <- mutations %>% 
  dplyr::select(study) %>% 
  unique() %>% 
  mutate(CHROM = list(mutations %>% 
                        dplyr::select(CHROM) %>% 
                        unique() %>% 
                        unlist())) %>% 
  unnest(CHROM) %>% 
  mutate(annotation1 = list(mutations %>% 
                              dplyr::select(annotation1) %>% 
                              unique() %>% 
                              unlist())) %>%
  unnest(annotation1) %>% 
  mutate(annotation2 = ifelse(annotation1 == "genic",
                              list(mutations %>% 
                                     filter(annotation1 == "genic") %>% 
                                     dplyr::select(annotation2) %>% 
                                     unique() %>% 
                                     unlist()),
                              list(mutations %>% 
                                     filter(annotation1 != "genic") %>% 
                                     dplyr::select(annotation2) %>% 
                                     unique() %>% 
                                     unlist())
                              )) %>% 
  unnest(annotation2) %>% unique()

mutations %>% 
  full_join(all_mut) %>% 
  mutate_at(c("N", "P"), funs(ifelse(is.na(.), 0, .))) %>% 
  write_tsv("save/mut_anot_genom0.tsv")


read_tsv("save/mut_anot_genom0.tsv") %>% 
  ggplot(aes(paste(annotation1, annotation2), P,  fill = study)) +
  geom_boxplot() +
  coord_flip()


genic <- interval_union(genes[genes$type == "mRNA"])
te <- interval_union(te)
genic_te <- interval_intersection(te, genic)
exon <- interval_union( genes[genes$type %in% c("CDS", "five_prime_UTR", "three_prime_UTR")])
genome <- Biostrings::readDNAStringSet("data/mutations/bordeaux/Qrob_PM1N.fa")
genome_width <- data.frame(CHROM = names(genome), total = width(genome)) %>% 
  filter(grepl("Chr", CHROM))
exon_width <- data.frame(CHROM = exon@annotation$seq_name, exon = width(exon)) %>% 
  filter(grepl("Chr", CHROM)) %>% 
  group_by(CHROM) %>% 
  summarise(exon = sum(exon))
genic_width <- data.frame(CHROM = genic@annotation$seq_name, genic = width(genic)) %>% 
  filter(grepl("Chr", CHROM)) %>% 
  group_by(CHROM) %>% 
  summarise(genic = sum(genic))
te_width <- data.frame(CHROM = te@annotation$seq_name, te = width(te)) %>% 
  filter(grepl("Chr", CHROM)) %>% 
  group_by(CHROM) %>% 
  summarise(te = sum(te))
genic_te_width<- data.frame(CHROM = genic_te@annotation$seq_name, genic_te = width(genic_te)) %>% 
  filter(grepl("Chr", CHROM)) %>% 
  group_by(CHROM) %>% 
  summarise(genic_te = sum(genic_te))
all_widths <- genome_width %>% 
  left_join(exon_width) %>% 
  left_join(genic_width) %>% 
  left_join(te_width) %>% 
  left_join(genic_te_width) %>% 
  mutate(intron = genic - exon) %>% 
  mutate(te_intergenic = te - genic_te) %>% 
  mutate(nonte_intergenic = total - exon - intron - te_intergenic)
all_widths %>% 
  dplyr::rename(exon_genic = exon, intron_genic = intron) %>% 
  dplyr::select(CHROM, exon_genic, intron_genic, te_intergenic, nonte_intergenic) %>% 
  reshape2::melt(id.var = "CHROM", variable.name = "annotation", value.name = "N") %>% 
  separate(annotation, c("annotation2", "annotation1")) %>% 
  mutate(annotation2 = recode(annotation2, "te" = "TE", "nonte" = "non-TE")) %>% 
  mutate(study = "Genome") %>% 
  group_by(study, CHROM) %>% 
  mutate(P = N / sum(N) * 100) %>% 
  write_tsv("save/genome_anot.tsv")
```

```{r fig4Amutswbdx}
study_cols <- RColorBrewer::brewer.pal(2, "Dark2") # 1 - Schmid-Siegert, 2 - Plomion
nucleotide_cols <- c('"#0E8F19"', '"#9C0C22"', '"#120b96"', '"#000000"') # 1 - A green, 2 - T red, 3 - C blue, 4 - G black
g0 <- lapply(list(swiss = "save/candidates_swiss.tsv",
            bordeaux = "save/candidates_bdx.tsv"), read_tsv) %>% 
  bind_rows(.id = "study") %>% 
  mutate(study = recode(study, "bordeaux" = "Plomion et al. (2018)", "swiss" = "Schmid-Siegert et al., (2017)")) %>% 
  mutate(study = factor(study, levels = c("Schmid-Siegert et al., (2017)", "Plomion et al. (2018)"))) %>% 
  mutate(dataset = recode(dataset, "Plomion" = "Original studies", "Schmid" = "Original studies",
                          "Robust" = "Strelka2 EVS", "Strelka2" = "Strelka2 All")) %>% 
  mutate(dataset = factor(dataset, levels = c("Original studies", "GATK", "Mutect2",
                                              "Strelka2 All", "Strelka2 EVS"))) %>% 
  mutate(study = recode(study,
                        "Schmid-Siegert et al., (2017)" = paste0("<span style = 'color: " , study_cols[1],
                                                                 ";'> Schmid-Siegert et al., (2017) </span>"), 
                        "Plomion et al. (2018)" = paste0("<span style = 'color: ", study_cols[2], 
                                                         ";'> Plomion et al. (2018) </span>")
                        )) %>% 
  ggplot(aes(AF, col = dataset), fill = NA) +
    geom_text(aes(y = y, col = dataset, label = label), show_guide  = FALSE,
             data.frame(study = rep(c("Schmid-Siegert et al., (2017)", "Plomion et al. (2018)"), each = 4),
                        dataset = c("Original studies", "GATK", "Strelka2 All", "Strelka2 EVS",
                                    "Original studies", "Mutect2", "Strelka2 All", "Strelka2 EVS"),
                        label = c("N=12", "N=510 611", "N=51", "N=41",
                                  "N=60", "N=464", "N=121", "N=87"),
                        AF = c(0.4, 0.5, 0.15, 0.15,
                               0.45, 0.2, 0.25, 0.35),
                        y = c(10, 5, 8, 9,
                              4, 10, 7, 5)) %>%
               mutate(study = recode(study,
                        "Schmid-Siegert et al., (2017)" = paste0("<span style = 'color: " , study_cols[1],
                                                                 ";'> Schmid-Siegert et al., (2017) </span>"), 
                        "Plomion et al. (2018)" = paste0("<span style = 'color: ", study_cols[2], 
                                                         ";'> Plomion et al. (2018) </span>")
                        )) %>% 
               mutate(dataset = factor(dataset, levels = c("Original studies", "GATK", "Mutect2",
                                                           "Strelka2 All", "Strelka2 EVS")))) +
  geom_density(size = 1.3) +
  scale_color_brewer("", palette = "Set1") +
  xlab("Allelic frequency")  + ylab("Frequency") +
  facet_wrap(~study, scales = "free_y", nrow = 1) +
  theme(legend.position = "bottom") +
  ggtitle("A.") +
  theme(strip.text = ggtext::element_markdown()) +
  guides(col=guide_legend(nrow=2,byrow=TRUE))
g0
```

```{r fig4Bmutswbdx}
study_cols <- RColorBrewer::brewer.pal(2, "Dark2") # 1 - Schmid-Siegert, 2 - Plomion
nucleotide_cols <- c('"#0E8F19"', '"#9C0C22"', '"#120b96"', '"#000000"') # 1 - A green, 2 - T red, 3 - C blue, 4 - G black
g1 <- bind_rows(vroom::vroom("save/mut_anot_genom0.tsv"),
                vroom::vroom("save/genome_anot.tsv")) %>% 
  mutate(study = recode(study, "Genome" = "Genomic expectation")) %>% 
  mutate(study = factor(study, levels = c("Genomic expectation",
                                          "Schmid-Siegert et al., (2017)", 
                                          "Plomion et al. (2018)"))) %>% 
  ggplot(aes(x = study, y = P, fill = study)) +
  geom_col(aes(y = m), position = "dodge",
           data = bind_rows(vroom::vroom("save/mut_anot_genom0.tsv"),
                            vroom::vroom("save/genome_anot.tsv")) %>%
             mutate(study = recode(study, "Genome" = "Genomic expectation")) %>%
             mutate(study = factor(study, levels = c("Genomic expectation",
                                                     "Schmid-Siegert et al., (2017)", 
                                                     "Plomion et al. (2018)"))) %>% 
             group_by(study, annotation1, annotation2) %>%
             summarise(m = mean(P), sd = sd(P)) %>%
             mutate(l = m - sd, h = m +sd) %>%
             mutate(l = ifelse(l < 0, 0, l))) +
  geom_linerange(aes(ymin = l, ymax = h),
                 data = bind_rows(vroom::vroom("save/mut_anot_genom0.tsv"),
                                  vroom::vroom("save/genome_anot.tsv")) %>%
                   mutate(study = recode(study, "Genome" = "Genomic expectation")) %>%
                   mutate(study = factor(study, levels = c("Genomic expectation",
                                                           "Schmid-Siegert et al., (2017)", 
                                                           "Plomion et al. (2018)"))) %>% 
                   group_by(study, annotation1, annotation2) %>%
                   summarise(m = mean(P), sd = sd(P), P = mean(P)) %>%
                   mutate(l = m - sd, h = m +sd) %>%
                   mutate(l = ifelse(l < 0, 0, l))) +
  ylab("Percentage") + xlab(" ") +
  facet_wrap(~ annotation1 + annotation2, nrow = 1) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(), legend.position = "bottom") +
  scale_fill_manual("", values = c("grey", study_cols[1], study_cols[2])) +
  ggtitle("B.")  +
  guides(fill = guide_legend(nrow=2,byrow=TRUE)) +
  ggpubr::stat_compare_means(comparisons = list(c("Genomic expectation", "Plomion et al. (2018)"),
                                                c("Genomic expectation", "Schmid-Siegert et al., (2017)")),
                             label = "p.signif", method = "t.test", label.y = c(80, 90)
  ) +
  ylim(0, 100)
g1
```

```{r Fig4Cbdxswiss}
study_cols <- RColorBrewer::brewer.pal(2, "Dark2") # 1 - Schmid-Siegert, 2 - Plomion
nucleotide_cols <- c('"#0E8F19"', '"#9C0C22"', '"#120b96"', '"#000000"') # 1 - A green, 2 - T red, 3 - C blue, 4 - G black

data0 <- read_tsv("save/candidates_all.tsv") %>% 
  mutate(study = factor(study, levels = c("Schmid-Siegert et al., (2017)", "Plomion et al. (2018)"))) %>% 
  mutate(dataset = recode(dataset, "Schmid" = "Original studies",
                          "Plomion" = "Original studies", 
                          "Strelka2" = "Strelka2 All", "Robust" = "Strelka2 EVS")) %>% 
  mutate(dataset = factor(dataset, levels = c("Original studies", "GATK", "Mutect2",
                                              "Strelka2 All", "Strelka2 EVS"))) %>% 
  na.omit() %>% 
  mutate(strand = recode(type,
                         "G->T" = "-",
                         "C->A" = "+",
                         "G->C" = "-",
                         "C->G" = "+",
                         "G->A" = "-",
                         "C->T" = "+",
                         "A->T" = "-",
                         "T->A" = "+",
                         "A->G" = "-",
                         "T->C" = "+",
                         "A->C" = "-",
                         "T->G" = "+"
                       )) %>% 
  mutate(base = type) %>% 
  mutate(type = recode(type, 
                       "G->T" = "C->A",
                       "G->C" = "C->G",
                       "G->A" = "C->T",
                       "A->T" = "T->A",
                       "A->G" = "T->C",
                       "A->C" = "T->G"
  )) %>% 
  mutate(POS5r = recode(POS5, "A" = "T", "T" = "A", "C" = "G", "G" = "C")) %>% 
  mutate(POS3r = recode(POS3, "A" = "T", "T" = "A", "C" = "G", "G" = "C")) %>% 
  mutate(POS5f = ifelse(strand == "+", POS5, POS3r)) %>% 
  mutate(POS3f = ifelse(strand == "+", POS3, POS5r))
data <- data0 %>% 
  filter(dataset  == "Strelka2 EVS") %>% 
  group_by(POS5f, POS3f, study, dataset, type) %>% 
  summarise(N = n()) %>% 
  group_by(study, dataset) %>% 
  mutate(P = N / sum(N) * 100) %>% 
  mutate(POS3f = factor(POS3f, levels = c("A", "C", "G", "T"))) %>% 
  mutate(POS5f = factor(POS5f, levels = c("T", "G", "C", "A"))) %>% 
  separate(type, c("REF", "ALT"), remove = F) %>% 
  full_join(data.frame(type = c("C->A", "C->G", "C->T", "T->A", "T->C", "T->G")) %>%
              mutate(POS5f = list(c("A", "C", "T", "G"))) %>%
              unnest(POS5f) %>%
              mutate(POS3f = list(c("A", "C", "T", "G"))) %>%
              unnest(POS3f) %>%
              separate(type, c("REF", "ALT"), remove = F) %>%
              mutate(study = list(c("Schmid-Siegert et al., (2017)", "Plomion et al. (2018)"))) %>%
              unnest(study)) %>%
    mutate(POS5f = paste("<span style = 'color: ",
                       recode(POS5f, "A" = nucleotide_cols[1], "T" = nucleotide_cols[2],
                              "C" = nucleotide_cols[3], "G" = nucleotide_cols[4]),
                       ";'>",
                       POS5f,
                       "</span>", sep = "")) %>%
  mutate(REF = paste("<span style = 'color: ",
                       recode(REF, "A" = nucleotide_cols[1], "T" = nucleotide_cols[2],
                              "C" = nucleotide_cols[3], "G" = nucleotide_cols[4]),
                       ";'>",
                       REF,
                       "</span>", sep = "")) %>%
  mutate(POS3f = paste("<span style = 'color: ",
                        recode(POS3f, "A" = nucleotide_cols[1], "T" = nucleotide_cols[2],
                              "C" = nucleotide_cols[3], "G" = nucleotide_cols[4]),
                       ";'>",
                       POS3f,
                       "</span>", sep = "")) %>%
    mutate(Spectrum = paste0(POS5f, REF, POS3f)) %>% 
  mutate(P = ifelse(is.na(P), 0.01, P)) %>% 
  ungroup() %>% 
  dplyr::select(study, type, Spectrum, P) %>% 
  mutate(x = paste0(type, Spectrum)) %>% 
  bind_rows(data0 %>% 
              dplyr::select(study, type) %>% 
              unique() %>% 
              mutate(x = paste0(type, "Z"), P = NA)) %>% 
  mutate(study = factor(study, levels = c("Schmid-Siegert et al., (2017)", 
                                          "Plomion et al. (2018)")))

data_label <- data %>% 
  group_by(x, Spectrum, study) %>% 
  summarise(P = max(P)) %>% 
  filter(study == "Plomion et al. (2018)") %>% 
  ungroup() %>% 
  mutate(id = 1:n())
number_of_bar <- nrow(data_label)
angle <- 90 - 360 * (data_label$id-0.5) /number_of_bar     # I substract 0.5 because the letter must have the angle of the center of the bars. Not extreme right(1) or extreme left (0)
data_label$hjust <- ifelse( angle < -90, 1, 0)
data_label$angle <- ifelse(angle < -90, angle+180, angle)

data_sgt <- data %>% 
   filter(study == "Plomion et al. (2018)") %>% 
  dplyr::select(type, x) %>% 
  unique() %>% 
  group_by(type) %>% 
  arrange(x) %>% 
  summarise(start = x[1], stop = x[17], mid = x[9]) %>% 
  mutate(class = ifelse(type %in% c("C->T", "T->C"), "transition", "transversion")) %>% 
  mutate(type = recode(type, 
                       "C->A" = sprintf("C \u2192 A"),
                       "C->G" = sprintf("C \u2192 G"),
                       "C->T" = sprintf("C \u2192 T"),
                       "T->A" = sprintf("T \u2192 A"),
                       "T->C" = sprintf("T \u2192 C"),
                       "T->G" = sprintf("T \u2192 G")
  ))

g2 <- ggplot(data,aes(x, P, fill = study)) +
  geom_bar(stat="identity", position = "dodge") +
  ylim(-20,NA) +
  theme_minimal() +
  coord_polar(start = 0) +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    plot.margin = unit(c(-5, -2, -2, -3), "cm"),
    plot.title = element_text(debug = F, vjust = -40, hjust = 0.25),
    panel.background = element_rect(fill = NA, color = NA)
  ) +
  ggtitle("C.") +
  scale_x_discrete(breaks=data$x, labels=data$Spectrum) +
  scale_fill_manual(guide = "none", values = c(study_cols[1], study_cols[2])) +
  ggtext::geom_richtext(data = data_label,
                        aes(x=x, y=-1, label=Spectrum),
                        fill = NA, label.color = NA,
                        size=2, angle= data_label$angle, inherit.aes = FALSE) +
  geom_segment(data = data_sgt,
               aes(x = start, y = -3, xend = stop, yend = -3, 
                   size = class),
               colour = "black", alpha=0.8, inherit.aes = FALSE ) +
  geom_text(data = data_sgt, 
            aes(x = mid, y = -8, label = type), 
            colour = "black", alpha=0.8, size=4, fontface="bold", inherit.aes = FALSE) +
  annotate("text", 
           x = rep(data$x[grepl("Z", data$x)], each = 3), 
           y = rep(c(0, 5, 10), sum(grepl("Z", data$x))),
           label = rep(c("0", "5", "10") , sum(grepl("Z", data$x))),
           color="grey", size=2 , angle=0, fontface="bold", hjust=1) +
  annotate("text", 
           x = max(data$x), 
           y = -17,
           label = "Pearson’r=0.49",
           color="black", size=5 , angle=0, fontface="bold") +
  annotate("text", 
           x = max(data$x), 
           y = -20,
           label = "p<1%*%10^-5",
           color="black", size=5 , angle=0, fontface="bold", parse=TRUE) +
  scale_size_manual(guide = "none", values = c(2, 0.6))

# g2
```

```{r Fig4,  fig.height=8, fig.width=12, eval=F}
gridExtra::grid.arrange(
  grobs = list(g0,  g1, g2),
  layout_matrix = rbind(c(1, 3),
                        c(2, 3)),
  widths = c(1.5,2),
  labels = c("A.", "B.", "C.")
)
```