11-MutationsAngela.Rmd

```{r setupmutationsangela, include=FALSE}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(knitr)
library(readr)
library(dplyr)
library(stringr)
library(tidyr)
library(tibble)
library(ggplot2)
library(circlize)
library(ComplexHeatmap)
library(ape)
library(ggtree)
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 renameAngela, eval=F, echo=F}
sampleRef <- read_delim("data/mutations/angela/prep/sampleRef.txt", delim = " ", 
                        col_names = c("idGeno", "id")) %>% 
  mutate(idGeno = gsub("DBR_", "", idGeno)) %>% 
  separate(id, c("sample", "leaf", "replicate"), remove = F)
libraries <- read_table("data/mutations/angela/prep/raw_libraries.txt", col_names = "files") %>% 
  separate(files, c("idGeno", "file"), "/")
ids <- read_tsv("data/mutations/angela/prep/idBordeaux.tsv") %>% 
  filter(EspeceBotanique == "Dicorynia guianensis") %>% 
  filter(Tissu != "graines") %>% 
  dplyr::select(id, Branch, Replicate, Lumiere) %>% 
  rename(sample = id, branchOld = Branch, light = Lumiere) %>% 
  mutate(branchOld = ifelse(is.na(branchOld), Replicate, branchOld)) %>% 
  dplyr::select(-Replicate) %>% 
  mutate(light = recode(light, "soleil" = "L", "ombre" = "S")) %>% 
    mutate(light = ifelse(is.na(light), "", light)) %>% 
  mutate(branch = recode(branchOld, "C1" = "T1", "C Sud" = "T3", "C Nord Est" = "T2", "C Nord Ouest" = "T4", 
                         "B31" = "A", "B32" = "B", "B2" = "C", "B4" = "D", "B1" = "E")) 
ids %>% 
  full_join(sampleRef) %>% 
  full_join(libraries) %>% 
  rename(sampleOld = sample, fileOld = file, idOld = id) %>% 
  mutate(leaf = ifelse(light == "", "", leaf)) %>% 
  mutate(sample = paste0(branch, light, leaf)) %>% 
  group_by(sample) %>% 
  mutate(library = paste0(sample, "_", 1:2)) %>% 
  mutate(file = paste0(library, ".fq.gz")) %>% 
  dplyr::select(sample, branch, light, leaf, library, file, branchOld, sampleOld, idOld, idGeno, fileOld) %>% 
  write_tsv("data/mutations/angela/prep/samples.tsv")

read_tsv("data/mutations/angela/prep/samples.tsv") %>% 
  mutate(command = paste("mv", fileOld, file)) %>% 
  dplyr::select(command) %>% 
  write_tsv("data/mutations/angela/prep/rename.sh", col_names = F)

read_tsv("data/mutations/angela/prep/samples.tsv") %>% 
  mutate(command = paste("touch", file)) %>% 
  dplyr::select(command) %>% 
  write_tsv("data/mutations/angela/prep/dummy.sh", col_names = F)
```

# Mutations Angela

This chapter describes the analyses of Angela mutations currently done in the the [`angela` branch of the `detectMutations` repository](https://github.com/sylvainschmitt/detectMutations/tree/angela).

## Quality check

The results of multiQC are available for the haplotypes HS1 [here](data/mutations/angela/qc/multiqc_report_HS1.html) and HS2 [here](data/mutations/angela/qc/multiqc_report_HS2.html). They globally all gave a green light for subsequent analyses.

## Genome

We compared the two haplophases and chose HS1 for contiguity and length (Fig. \@ref(fig:AngelaGenomFig) and Fig. \@ref(fig:AngelaBusco)) as well as mean coverage and homogeneity (Fig. \@ref(fig:circosHS1g)).

```{r AngelaGenomeData}
hs1 <- read_tsv("data/mutations/angela/genome/Dgu_HS1_HYBRID_SCAFFOLD.fa.fai", 
         col_names = c("name", "length", "bytesindex", "basesperline", "bytesperline")) %>% 
  dplyr::select(name, length) %>% 
  mutate(name = gsub("Super-Scaffold_", "", name)) %>% 
  mutate(start = 0, end = length) %>% 
  dplyr::select(-length) %>%  
  mutate(trsc = letters[1:nrow(.)], exon = LETTERS[1:nrow(.)])
hs2 <- read_tsv("data/mutations/angela/genome/Dgu_HS2_HYBRID_SCAFFOLD.fa.fai", 
         col_names = c("name", "length", "bytesindex", "basesperline", "bytesperline")) %>% 
  dplyr::select(name, length) %>% 
  mutate(name = gsub("Super-Scaffold_", "", name)) %>% 
  mutate(start = 0, end = length) %>% 
  dplyr::select(-length) %>%  
  mutate(trsc = letters[1:nrow(.)], exon = LETTERS[1:nrow(.)])
gc_hs1 <- vroom::vroom("data/mutations/angela/genome/Dgu_HS1_HYBRID_SCAFFOLD.gc", skip = 1,
             col_names = c("chrom", "start", "stop", "AT", "GC", "A", "C", "G", "T", "N", "other", "length")) %>% 
  mutate(pos = start + 500) %>% 
  mutate(chrom = gsub("Super-Scaffold_", "", chrom))
gc_hs2 <- vroom::vroom("data/mutations/angela/genome/Dgu_HS2_HYBRID_SCAFFOLD.gc", skip = 1,
             col_names = c("chrom", "start", "stop", "AT", "GC", "A", "C", "G", "T", "N", "other", "length")) %>% 
  mutate(pos = start + 500) %>% 
  mutate(chrom = gsub("Super-Scaffold_", "", chrom))
cov_hs1 <- vroom::vroom("data/mutations/angela/genome/DS2_HS1.regions.bed", col_names = c("chrom", "start", "stop", "coverage")) %>% 
  mutate(chrom = gsub("Super-Scaffold_", "", chrom))
cov_hs2 <- vroom::vroom("data/mutations/angela/genome/DS2_HS2.regions.bed", col_names = c("chrom", "start", "stop", "coverage")) %>% 
  mutate(chrom = gsub("Super-Scaffold_", "", chrom))
busco_hs1 <- read_tsv("data/mutations/angela/genome/short_summary.specific.viridiplantae_odb10.Dgu_HS1_HYBRID_SCAFFOLD_busco.txt",
         skip = 10, n_max = 4, col_names = c("X1", "N", "Type")) %>% 
  dplyr::select(-X1)
busco_hs2 <- read_tsv("data/mutations/angela/genome/short_summary.specific.viridiplantae_odb10.Dgu_HS2_HYBRID_SCAFFOLD_busco.txt",
         skip = 10, n_max = 4, col_names = c("X1", "N", "Type")) %>% 
  dplyr::select(-X1)
covs_files <- list.files("data/mutations/angela/genome/dists/", full.names = T)
names(covs_files) <- list.files("data/mutations/angela/genome/dists/")
covs <- lapply(covs_files, read_tsv, col_names = c("chromosome", "coverage", "proportion"), cols(
  chromosome = col_character(),
  coverage = col_double(),
  proportion = col_double()
)) %>% 
  bind_rows(.id = "file") %>% 
  mutate(library = gsub(".mosdepth.global.dist.txt", "", file)) %>% 
  filter(chromosome == "total")
cov_lims <- covs %>% 
  filter(proportion == 0.05 | proportion == 0.95 | proportion == 0.5) %>% 
  group_by(library, proportion) %>% 
  summarise(coverage = mean(coverage)) %>% 
  mutate(proportion = recode(as.character(proportion), "0.05" = "high_cov", "0.95" = "low_cov", "0.5" = "med_cov")) %>% 
  reshape2::dcast(library ~ proportion)
```

```{r AngelaGenomFig, fig.cap="Anchored scaffolds size-distribution for haplotypes HS1 and HS2 in Angela's genome."}
cowplot::plot_grid(
  plotlist = 
    lapply(list(HS1 = hs1, HS2 = hs2), function(genome)
      ggplot(genome, aes(reorder(name, end), end/10^6, fill = end > 2*10^6)) +
        geom_col() +
        coord_flip() +
        ylab("Length (mb)") + 
        xlab("Anchored super-scaffold") + 
        scale_fill_discrete(">2Mb")),
  labels = c("HS1", "HS2")
)
```

```{r AngelaBusco, fig.cap="BUSCO results for haplotypes HS1 and HS2."}
list(HS1 = busco_hs1, HS2 = busco_hs2) %>% 
  bind_rows(.id = "haploytpe") %>% 
  mutate(Type = factor(Type, levels= c("Complete and single-copy BUSCOs (S)",
                                       'Complete and duplicated BUSCOs (D)',
                                       "Fragmented BUSCOs (F)",
                                       "Missing BUSCOs (M)"))) %>% 
  ggplot(aes(haploytpe, N, fill = Type)) +
  geom_col() +
  geom_text(aes(y = N, label = N), col = "white", position = position_stack(vjust = .5)) +
  scale_y_sqrt() +
  theme(axis.title.x = element_blank())
```

```{r AngelaCovs, fig.cap="Coverage distribution. Distribution of the number of locations in the reference genome with a given depth of coverage."}
ggplot(covs, aes(coverage, proportion*100, col = library)) +
  geom_line() +
  xlim(0, 500) +
  geom_hline(yintercept = 5, col = "red", linetype = "dashed") +
  geom_hline(yintercept = 95, col = "red", linetype = "dashed") +
  scale_color_discrete("") + 
  xlab("Coverage (X)") +
  ylab("% of base in the genome with at least X reads")
```

```{r circosHS1, eval=F}
hs1_f <- filter(hs1, end > 2*10^6)
gc_hs1_f <- filter(gc_hs1, chrom %in% hs1_f$name)
cov_hs1_f <- filter(cov_hs1, chrom %in% hs1_f$name) %>% 
  mutate(coverage = ifelse(coverage > 500, 500, coverage))
hs2_f <- filter(hs2, end > 2*10^6)
gc_hs2_f <- filter(gc_hs2, chrom %in% hs2_f$name)
cov_hs2_f <- filter(cov_hs2, chrom %in% hs2_f$name) %>% 
  mutate(coverage = ifelse(coverage > 500, 500, coverage))
png("save/AngelaCircosHS1CovGC.png", width = 1000, height = 1000)
circos.initializeWithIdeogram(as.data.frame(hs1_f))
circos.genomicTrack(cov_hs1_f, 
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "blue", ...)
})

circos.genomicTrack(dplyr::select(gc_hs1_f, chrom, start, stop, GC), 
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "red", ...)
})
text(0, 0, "HS1", cex = 1.5)
legend("bottomleft", pch = 16, col = c("blue", "red"), legend = c("coverage", "GC"), bty = "n")
dev.off()
png("save/AngelaCircosHS2CovGC.png", width = 1000, height = 1000)
circos.initializeWithIdeogram(as.data.frame(hs2_f))
circos.genomicTrack(cov_hs2_f,
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "blue", ...)
})

circos.genomicTrack(dplyr::select(gc_hs2, chrom, start, stop, GC),
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "red", ...)
})
text(0, 0, "HS2", cex = 1.5)
legend("bottomleft", pch = 16, col = c("blue", "red"), legend = c("coverage", "GC"), bty = "n")
dev.off()
```

```{r circosHS1g, fig.cap="Genome coverage and GC-content on a 10-kb windows for the HS1 haplotype (scaffolds > 2-Mb)."}
include_graphics("save/AngelaCircosHS1CovGC.png")
```

```{r circosHS2g, fig.cap="Genome coverage and GC-content on a 10-kb windows for the HS2 haplotype (scaffolds > 2-Mb)."}
include_graphics("save/AngelaCircosHS2CovGC.png")
```

```{r}
rm(busco_hs1, busco_hs2, cov_hs1, cov_hs2, covs, gc_hs1, gc_hs2, covs_files, hs2)
```

### Annotation

```{r, eval=F}
annotation <- DBI::dbConnect(RSQLite::SQLite(), dbname = "data/mutations/angela/annotation/angela.annot.sqlite")
orfs <- tbl(annotation, "ORF") %>% 
  mutate(start = str_sub(peptide, 1, 1)) %>% 
  mutate(stop = str_sub(peptide, -1, -1)) %>% 
  mutate(length = rend - lend) %>% 
  collect()
DBI::dbDisconnect(annotation) ; rm(annotation)
n <- orfs %>% 
  filter(start == "M", stop == "*", length >= 100) %>% 
  nrow()
orfs %>% 
  dplyr::select(orf_id, start, stop, length) %>% 
  reshape2::melt(c("orf_id", "length")) %>% 
  ggplot(aes(value)) +
  geom_bar() +
  facet_grid((length > 100) ~ variable) +
  ggtitle(paste("Angela:", n, "ORFs"), "Start M, Stop *, length >= 100")
rm(orfs, n)
```

```{bash, eval=F}
sqlite3 angela.annot.sqlite
sqlite> .headers on
sqlite> .mode csv
sqlite> .output angela.orf.csv
sqlite> SELECT * FROM ORF;
sqlite> .quit
```

```{r, circosAnnot, eval=F}
min_length <- 20*10^6
window <- 100*10^3
hs1_f <- filter(hs1, end > min_length)
orfs <- vroom::vroom("data/mutations/angela/annotation/angela.orf.csv") %>% 
  mutate(start = str_sub(peptide, 1, 1)) %>% 
  mutate(stop = str_sub(peptide, -1, -1)) %>% 
  mutate(length = rend - lend) %>% 
  filter(start == "M", stop == "*", length >= 100) %>% 
  separate(transcript_id, c("gene_id", "isoform"), remove = F)
gene_f <- vroom::vroom("data/mutations/angela/annotation/braker.gtf",
                   col_names = c("scf", "source", "type", "start", "end", 
                                 "score", "strand", "phase", "attributes")) %>% 
   filter(type  == "gene") %>% 
  filter(attributes %in% unique(orfs$gene_id)) %>% 
  mutate(name = gsub("Super-Scaffold_", "", scf)) %>% 
  filter(name %in% hs1_f$name) %>% 
  mutate(interval = cut(start, breaks = seq(0, max(start)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(gene = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T) %>% 
  ungroup()
te_f <- read_tsv("data/mutations/angela/annotation/Dgu_HS1_HYBRID_SCAFFOLD.fa.out.gff", skip = 3, 
               col_names = c("scf", "source", "type", "start", "end", "score", "strand", "phase", "attributes")) %>% 
  mutate(name = gsub("Super-Scaffold_", "", scf)) %>% 
  filter(name %in% hs1_f$name) %>% 
  mutate(interval = cut(start, breaks = seq(0, max(start)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(TE = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T)
png("save/AngelaCircosAnnot.png", width = 1000, height = 1000, bg = NA)
circos.par(gap.degree = c(rep(1, nrow(hs1_f)-2), 5, rep(1, 1)),
           track.height = 0.2,
           points.overflow.warning = F)
circos.initializeWithIdeogram(as.data.frame(hs1_f),
                              labels.cex = 4, axis.labels.cex = 1)
circos.genomicTrack(gene_f, track.height = 0.2,
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "#7070ff", border = "#7070ff", area = T, ...)
})
circos.yaxis(side = "left", at = c(5, 10, 15), labels.cex = 1)
circos.genomicTrack(te_f, track.height = 0.2,
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "#66ccff", border = "#66ccff", area = T, ...)
})
circos.yaxis(side = "left", at = c(200,400), labels.cex = 1)
dev.off()
circos.clear()
```

```{r}
include_graphics("save/AngelaCircosAnnot.png")
```

### Mutations

```{r}
min_length <- 20*10^6
window <- 100*10^3
hs1_f <- filter(hs1, end > min_length)
mutation_f <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  dplyr::select(SNV, CHROM, POS) %>% 
  unique() %>% 
  mutate(name = gsub("Super-Scaffold_", "", CHROM)) %>% 
  filter(name %in% hs1_f$name) %>% 
  mutate(interval = cut(POS, breaks = seq(0, max(POS)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(mutations = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T) %>% 
  ungroup()
candidates <- bind_rows(
  read_tsv("data/mutations/fruits/angela_fruits_candidate_mutations.tsv") %>% 
    mutate(mutation = 1:n(), tree = "Angela"),
  read_tsv("data/mutations/fruits/sixto_fruits_candidate_mutations.tsv") %>% 
    mutate(mutation = 125:(124+n()), tree = "Sixto")
) %>% 
  dplyr::select(CHROM, POS, REF, ALT, af, branch, mutation, tree) %>% 
  mutate(CHROM = as.numeric(gsub("Super-Scaffold_", "", CHROM)))
transmitted_f <- readxl::read_xlsx("data/mutations/fruits/fruit_mutations_igv.xlsx", "classification") %>% 
  mutate(mutation = as.numeric(gsub("SNV", "", SNV))) %>% 
  dplyr::select(mutation, Type, SNV) %>% 
  left_join(candidates) %>% 
  mutate(transmitted = ifelse(mutation %in% c(6, 13, 31, 54, 57, 107, 128, 132, 151, 160), "yes", "no")) %>% 
  filter(tree == "Angela", transmitted == "yes") %>% 
  mutate(name = as.character(CHROM), start = POS, stop = POS, snv = SNV) %>% 
  dplyr::select(name, start, stop, snv)
orfs <- vroom::vroom("data/mutations/angela/annotation/angela.orf.csv") %>% 
  mutate(start = str_sub(peptide, 1, 1)) %>% 
  mutate(stop = str_sub(peptide, -1, -1)) %>% 
  mutate(length = rend - lend) %>% 
  filter(start == "M", stop == "*", length >= 100) %>% 
  separate(transcript_id, c("gene_id", "isoform"), remove = F)
gene_f <- vroom::vroom("data/mutations/angela/annotation/braker.gtf",
                   col_names = c("scf", "source", "type", "start", "end", 
                                 "score", "strand", "phase", "attributes")) %>% 
   filter(type  == "gene") %>% 
  filter(attributes %in% unique(orfs$gene_id)) %>% 
  mutate(name = gsub("Super-Scaffold_", "", scf)) %>% 
  filter(name %in% hs1_f$name) %>% 
  mutate(interval = cut(start, breaks = seq(0, max(start)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(gene = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T) %>% 
  ungroup()
te_f <- read_tsv("data/mutations/angela/annotation/Dgu_HS1_HYBRID_SCAFFOLD.fa.out.gff", skip = 3, 
               col_names = c("scf", "source", "type", "start", "end", "score", "strand", "phase", "attributes")) %>% 
  mutate(name = gsub("Super-Scaffold_", "", scf)) %>% 
  filter(name %in% hs1_f$name) %>% 
  mutate(interval = cut(start, breaks = seq(0, max(start)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(TE = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T)
```

```{r, eval=FALSE}
png("save/AngelaCircosmutations.png", width = 1000, height = 1000, bg = NA)
circos.par(gap.degree = c(rep(1, nrow(hs1_f)-2), 5, rep(1, 1)),
           track.height = 0.2,
           points.overflow.warning = F)
circos.initializeWithIdeogram(as.data.frame(hs1_f),
                              labels.cex = 4, axis.labels.cex = 1)
circos.genomicTrack(gene_f, track.height = 0.1,
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "#7070ff", border = "#7070ff", area = T, ...)
})
circos.yaxis(side = "left", at = c(5, 10, 15), labels.cex = 1)
circos.genomicTrack(te_f, track.height = 0.1,
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "#66ccff", border = "#66ccff", area = T, ...)
})
circos.yaxis(side = "left", at = c(200,400), labels.cex = 1)
circos.genomicTrack(mutation_f,
    panel.fun = function(region, value, ...) {
        circos.genomicLines(region, value, col = "#ffc000", border = "#ffc000", area = T, ...)
})
circos.yaxis(side = "left", labels.cex = 1)
circos.genomicLabels(transmitted_f, labels.column = 4, side = "inside", 
                     col = "#cc0000", cex = 1.5)
dev.off()
circos.clear()
```

```{r}
include_graphics("save/AngelaCircosmutations.png")
```

```{r, eval=F}
library(chromoMap)
min_length <- 10*10^6
window <- 100*10^3
hs1_f <- filter(hs1, end > min_length)
mutation_f <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  dplyr::select(SNV, CHROM, POS) %>% 
  unique() %>% 
  mutate(name = gsub("Super-Scaffold_", "", CHROM)) %>% 
  filter(name %in% hs1_f$name) %>% 
  mutate(interval = cut(POS, breaks = seq(0, max(POS)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(mutations = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T) %>% 
  ungroup()
gene_f <- read_tsv("data/mutations/angela/annotation/braker.gtf",
                   col_names = c("scf", "source", "type", "start", "end", "score", "strand", "phase", "attributes")) %>% 
  mutate(name = gsub("Super-Scaffold_", "", scf)) %>% 
  filter(name %in% hs1_f$name) %>% 
  filter(type  == "gene") %>% 
  mutate(interval = cut(start, breaks = seq(0, max(start)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(gene = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T) %>% 
  ungroup()
t1 <- hs1_f %>%  
  mutate(name = paste0("Sc", name)) %>% 
  mutate(start = as.numeric(start+1),
         end = as.numeric(end+1)) %>%
  dplyr::select(name, start, end) %>% 
  as.data.frame()
t2 <- mutation_f %>%   
  # left_join(gene_f) %>% 
  mutate(annot = paste0("A", 1:n())) %>% 
  mutate(name=paste0("Sc", name)) %>% 
    mutate(start = as.numeric(start+1),
         end = as.numeric(stop+1),
         mutations = as.numeric(mutations)) %>%
  dplyr::select(annot, name, start, end, mutations) %>% 
  as.data.frame()
chromoMap(list(t1), list(t2), data_based_color_map = T,
          data_type = "numeric", legend = T)
```

```{r}
# library(RIdeogram)
```

### Comp

```{r}
library(ivs)
window <- 10^5
windows <- hs1 %>% 
  filter(end > window) %>% 
  rowwise() %>% 
  mutate(window_end = list(seq(from = (0+window), to = end, by = window))) %>% 
  unnest_longer(window_end) %>% 
  mutate(window_start = lag(window_end)) %>% 
  mutate(window_start = ifelse(is.na(window_start), 0, window_start)) %>% 
  filter(window_end - window_start == window) %>% 
  mutate(window  = ivs::iv(window_start, window_end)) %>% 
  dplyr::select(name, window)
genes <- vroom::vroom("data/mutations/angela/annotation/braker.gtf",
                   col_names = c("scf", "source", "type", "start", "end", 
                                 "score", "strand", "phase", "attributes")) %>% 
   filter(type  == "gene") %>% 
  filter(attributes %in% unique(orfs$gene_id)) %>% 
  mutate(name = gsub("Super-Scaffold_", "", scf)) %>% 
  group_by(name) %>% 
  mutate(gene = ivs::iv(start, end)) %>% 
  dplyr::select(name, gene)
gene_pct <- lapply(unique(windows$name), function(n){
  locations <- iv_locate_overlaps(
    filter(windows, name == n)$window,
    filter(genes, name == n)$gene,
    type = "contains", 
    no_match = "drop"
  )
  iv_align(filter(windows, name == n)$window, filter(genes, name == n)$gene,
           locations = locations) %>% 
    group_by(needles) %>% 
    reframe(range = iv_groups(haystack)) %>% 
    mutate(length = iv_end(range) - iv_start(range)) %>% 
    group_by(needles) %>% 
    summarise(length = sum(length)) %>% 
    mutate(percentage_gene = length/window) %>% 
    mutate(start = iv_start(needles), stop = iv_end(needles)) %>% 
    mutate(name = n) %>% 
    dplyr::select(name, start, stop, percentage_gene)
}) %>% bind_rows()
mutations <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  dplyr::select(SNV, CHROM, POS) %>% 
  unique() %>% 
  mutate(name = gsub("Super-Scaffold_", "", CHROM)) %>% 
  mutate(interval = cut(POS, breaks = seq(0, max(POS)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(mutations = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T) %>% 
  ungroup()
all <- gene_pct %>% 
  left_join(mutations) %>% 
  mutate(mutations = ifelse(is.na(mutations), 0, mutations)) %>% 
  mutate(percentage_gene_interval = cut(percentage_gene, 
                                        breaks = seq(0, 1, by = 0.1)))
ggplot(all, aes(percentage_gene, mutations)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, col = "red") +
  ggtitle("Angela", "Windows of 10^5")
glm(mutations ~ percentage_gene, all, family = "poisson") %>% 
  summary()
```

```{r}
library(ivs)
window <- 10^5
windows <- hs1 %>% 
  filter(end > window) %>% 
  rowwise() %>% 
  mutate(window_end = list(seq(from = (0+window), to = end, by = window))) %>% 
  unnest_longer(window_end) %>% 
  mutate(window_start = lag(window_end)) %>% 
  mutate(window_start = ifelse(is.na(window_start), 0, window_start)) %>% 
  filter(window_end - window_start == window) %>% 
  mutate(window  = ivs::iv(window_start, window_end)) %>% 
  dplyr::select(name, window)
tes <- read_tsv("data/mutations/angela/annotation/Dgu_HS1_HYBRID_SCAFFOLD.fa.out.gff", skip = 3, 
               col_names = c("scf", "source", "type", "start", "end", "score", "strand", "phase", "attributes")) %>% 
  mutate(name = gsub("Super-Scaffold_", "", scf)) %>% 
  group_by(name) %>% 
  mutate(te = ivs::iv(start, end)) %>% 
  dplyr::select(name, te)
tes_pct <- lapply(unique(windows$name), function(n){
  locations <- iv_locate_overlaps(
    filter(windows, name == n)$window,
    filter(tes, name == n)$te,
    type = "contains", 
    no_match = "drop"
  )
  iv_align(filter(windows, name == n)$window, filter(tes, name == n)$te,
           locations = locations) %>% 
    group_by(needles) %>% 
    reframe(range = iv_groups(haystack)) %>% 
    mutate(length = iv_end(range) - iv_start(range)) %>% 
    group_by(needles) %>% 
    summarise(length = sum(length)) %>% 
    mutate(percentage_te = length/window) %>% 
    mutate(start = iv_start(needles), stop = iv_end(needles)) %>% 
    mutate(name = n) %>% 
    dplyr::select(name, start, stop, percentage_te)
}) %>% bind_rows()
mutations <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  dplyr::select(SNV, CHROM, POS) %>% 
  unique() %>% 
  mutate(name = gsub("Super-Scaffold_", "", CHROM)) %>% 
  mutate(interval = cut(POS, breaks = seq(0, max(POS)+window, by = window))) %>% 
  group_by(name, interval) %>% 
  summarise(mutations = n()) %>% 
  mutate(interval = gsub("(", "", interval, fixed = T)) %>% 
  mutate(interval = gsub("]", "", interval, fixed = T)) %>% 
  separate(interval, c("start", "stop"), ",", convert = T) %>% 
  ungroup()
all <- tes_pct %>% 
  left_join(mutations) %>% 
  mutate(mutations = ifelse(is.na(mutations), 0, mutations))
ggplot(all, aes(percentage_te, mutations)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, col = "red") +
  ggtitle("Angela", "Windows of 10^5")
glm(mutations ~ percentage_te, all, family = "poisson") %>% 
  summary()
```


## Heterozygosity

We used k-mer analyses (21-mers) using `jellyfish` and the `GenomeScope` to estimate genetic diversity $\pi=0.901%$ (Fig. \@ref(fig:genomescopeAngela)).
**Should I also estimate non-synonymous diversity with ORF annotation using `transdecoder` and `SNPdat` (https://phdpages.netlify.app/6b189d4142e8224bff99e28abd11cbfdd50c51b1/symcapture/annotation.html#synonymy)?**
We detected and filtered SNPs as follow:

* Raw=5M : raw result from `GATK HaplotypeCaller + GenomicsDBImport + GenotypeGVCFs`
* Biallelic=5.3M : biallelic sites with `bcftools`
* SNP=4.7M : SNPs with `GATK`
* Filtered=3.85M : SNPs with QUAL < 30, QD < 2, FS > 60, SOR > 3 using `GATK`
* Non-missing=3.84M : SNPs in all genotypes and individuals with `plink`
* Shared=3.84M : shared SNPs by at least 32 out of 33 individuals with  `bcftools`

```{r genomescopeAngela, fig.cap="GenomeScope Profile. Full results here: http://genomescope.org/analysis.php?code=ayEfgONEP1cFWw2yjhfb."}
knitr::include_graphics("data/mutations/angela/hz/genomescope.png")
```

## Mutations filters

```{r, cache=FALSE}
trunk <- DBI::dbConnect(RSQLite::SQLite(), "data/mutations/angela/raw/cambium_nonhz_mutations.sql") %>% tbl("mutations")
leaf <- DBI::dbConnect(RSQLite::SQLite(), "data/mutations/angela/raw/leaf_nonhz_mutations.sql") %>% tbl("mutations")
```

We removed all candidate mutations present in the heterozygous sites and in the raw candidate mutations from the cambium.
We kept only mutations with no copy of the mutated allele in the normal sample (NAC, `normal_altCount == 0`).
We tested independently the effect of 4 filters and look their effect on (1) the percentage kept, (2) the distribution of allelic frequencies, and (3) the overlap between filters:

1. **MAC**: a minimum of 5 copies of the mutated allele in the mutated sample (`mutation_altCount >= 5`)
1. **DP**: a read depth for the two sample between the 5th quantile and the 95th quantile of the coverage of the corresponding library (`normal_DP <= high_cov, normal_DP >= low_cov, mutation_DP <= high_cov, mutation_DP >= low_cov`)
1. **BIO**: the mutation is present in at least two biological replicates (2 leaves from the crown)
1. **EVS**: `Strelka 2` automatic filtering based on the empiric variant score (`Filter == "PASS"`)

**MAC** and **EVS** are the most stringent filters (Fig. \@ref(fig:pctmutfilter), 
and mostly **MAC** filter changes the allele frequencies distribution (Fig. \@ref(fig:afmutfilter).
Most filters share mutations, except **EVS** and **BIO** that rejected individually a lot mutations shared by the two others (Fig. \@ref(fig:vennmutfilter).
And all filters except **MAC** are not sensitive to the library coverage (Fig. \@ref(fig:nmutcovfilter).
We will thus use two filters for next steps:

1. **base**: 
    * **NAC**: no copy of the mutated allele in the normal sample (`normal_altCountT1 == 0`)
    * **MAC**: a minimum of 5 copies of the mutated allele in the mutated sample (`mutation_altCount >= 5`)
    * **DP**: a read depth for the two sample between the 10th quantile and the 90th quantile of the coverage of the mutated library (`normal_DP <= high_cov, normal_DP >= low_cov, mutation_DP <= high_cov, mutation_DP >= low_cov`)
    * **BIO**: the mutation is present in at least two biological replicates (2 leaves from the crown)
1. **evs**: 
    * **base**
    * **EVS**: `Strelka 2` automatic filtering based on the empiric variant score (`Filter == "PASS"`)

```{r, eval=F}
library(dtplyr)
# pot trunk
trunk_pot_mut <- trunk %>% 
  mutate(SNV = paste0(CHROM, "#", as.integer(POS))) %>% 
  dplyr::select(SNV) %>% collect() %>% unlist()
# all leaves
leaf_mutations <- leaf %>% 
  filter(normal_altCountT1 == 0, normal_altCountT2 == 0) %>% 
  mutate(SNV = paste0(CHROM, "#", as.integer(POS))) %>% 
   dplyr::select(tumor, normal,SNV, normal_altCountT1, normal_altCountT2, normal_DP, mutation_DP, 
                 mutation_altCountT1, mutation_altCountT2, FILTER, mutation_AF) %>% 
  collect() %>% 
  filter(!(SNV %in% trunk_pot_mut))
# define filters
leaf_mutations <- leaf_mutations %>% 
  left_join(rename(cov_lims, tumor = library, mutation_high_cov = high_cov, mutation_low_cov = low_cov, mutation_med_cov = med_cov)) %>% 
  left_join(rename(cov_lims, normal = library, normal_high_cov = high_cov, normal_low_cov = low_cov, normal_med_cov = med_cov)) %>% 
  mutate(DP = (normal_DP <= normal_high_cov &
                 normal_DP >= normal_low_cov &
                 mutation_DP <= mutation_high_cov &
                 mutation_DP >= mutation_low_cov)) %>% 
  mutate(MAC = (mutation_altCountT1 >= 10 & mutation_altCountT1 >= 10)) %>% 
  mutate(EVS = (FILTER == "PASS")) %>% 
  group_by(SNV) %>% 
  mutate(replicate = n()) %>% 
  ungroup() %>% 
  mutate(BIO = replicate >= 2) %>% 
  dplyr::select(tumor, SNV, mutation_AF, normal_DP, mutation_DP, replicate, DP, MAC, EVS, BIO) %>% 
  reshape2::melt(c("tumor", "SNV", "mutation_AF", "normal_DP", "mutation_DP", "replicate"), variable.name = "filter")
vroom::vroom_write(leaf_mutations, "data/mutations/angela/filters/filters_mutations.stats")
# percentanges
lazy_dt(leaf_mutations) %>% 
  group_by(tumor, filter, value) %>% 
  summarise(N = n())  %>% 
  as_tibble() %>% 
  reshape2::dcast(tumor + filter ~ value, value.var = 'N') %>% 
  mutate(percentage = `TRUE`/(`TRUE`+`FALSE`)*100) %>% 
  vroom::vroom_write("data/mutations/angela/filters/filters_mutations.pct")
# AF
g <- ggplot(leaf_mutations, aes(mutation_AF, col = value)) +
  geom_density() +
  scale_x_sqrt() +
  facet_wrap(~ filter) +
  xlab("Mutation allelic frequency") +
  theme(legend.position = "bottom")
ggsave(g, filename = "data/mutations/angela/filters/filters_mutations_af.png", width = 4, height = 4, bg = "white")
# sharing
t_venn <- lazy_dt(leaf_mutations) %>% 
  filter(filter != "BASE") %>% 
  filter(value == TRUE) %>%
  dplyr::select(-mutation_AF, -tumor, -value) %>% 
  unique() %>% 
  as_tibble()
g <- ggvenn::ggvenn(lapply(split(t_venn, t_venn$filter), `[[`, "SNV"), show_percentage = F)
ggsave(g, filename = "data/mutations/angela/filters/filters_mutations_venn.png", bg = "white")
rm(t_venn, g)
```

```{r pctmutfilter, fig.cap="Percentage of kept mutations per filter."}
vroom::vroom("data/mutations/angela/filters/filters_mutations.pct") %>% 
  ggplot(aes(filter, percentage)) +
  geom_boxplot() +
  geom_jitter(aes(col = tumor)) +
  scale_y_log10() +
  scale_color_discrete("") +
  xlab("Filter") + ylab("Percentage of mutations kept")     
```

```{r afmutfilter, fig.cap="Alleles frequencies of filtered mutations.", fig.width=4, fig.height=4}
knitr::include_graphics("data/mutations/angela/filters/filters_mutations_af.png")  
```

```{r vennmutfilter, fig.cap="Filtered mutations sharing across filters.", fig.width=4, fig.height=4}
knitr::include_graphics("data/mutations/angela/filters/filters_mutations_venn.png") 
```

```{r nmutcovfilter, fig.cap="Link between the number of detected mutations and the coverage across samples depending on the filter used."}
vroom::vroom("data/mutations/angela/filters/filters_mutations.pct")  %>%  
  left_join(dplyr::rename(cov_lims, tumor = library)) %>% 
  group_by(tumor, filter) %>% 
  ggplot(aes(med_cov, `TRUE`, col = filter, label = tumor)) +
  geom_point() +
  ggrepel::geom_text_repel() +
  geom_smooth(se = F, method = "lm") +
  scale_y_log10() +
  xlab("Coverage") + ylab("Number of mutations") 
```

## Leaf mutations

```{r, eval=F}
library(dtplyr)
trunk_pot_mut <- trunk %>% 
  mutate(SNV = paste0(CHROM, "#", as.integer(POS))) %>% 
  dplyr::select(SNV) %>% collect() %>% unlist()
leaf_mutations <- leaf %>% 
  filter(!(SNV %in% trunk_pot_mut)) %>% 
  mutate(SNV = paste0(CHROM, "#", as.integer(POS))) %>% 
  filter(normal_altCountT1 == 0, normal_altCountT2 == 0) %>% 
  filter(mutation_altCountT1 >= 5, mutation_altCountT2 >= 5) %>% 
  collect() %>% 
  lazy_dt() %>% 
  left_join(rename(cov_lims, tumor = library, mutation_high_cov = high_cov, mutation_low_cov = low_cov, mutation_med_cov = med_cov)) %>% 
  left_join(rename(cov_lims, normal = library, normal_high_cov = high_cov, normal_low_cov = low_cov, normal_med_cov = med_cov)) %>% 
  filter(normal_DP <= normal_high_cov & normal_DP >= normal_low_cov & mutation_DP <= mutation_high_cov & mutation_DP >= mutation_low_cov) %>%
  group_by(SNV) %>% 
  mutate(replicate = n()) %>% 
  filter(replicate >= 2) %>% 
  ungroup() %>% 
  mutate(Filter = ifelse(FILTER == "PASS", "evs", "base")) %>% 
  dplyr::select(-ChromKey, -NT, -ID, -QUAL, -FILTER, -SOMATIC, -PNOISE, -PNOISE2, -normal, -caller, 
                -mutation_high_cov, -mutation_low_cov, -mutation_med_cov, -normal_high_cov, -normal_low_cov, -normal_med_cov) %>% 
  group_by(tumor, SNV, Filter, CHROM, POS, REF, ALT, SGT) %>% 
  summarise_all(median) %>% 
  mutate(branch = str_sub(tumor, 1, 1), tip = str_sub(tumor, 2, 2), leaf = str_sub(tumor, 3, 3))
vroom::vroom_write(as_tibble(leaf_mutations), "data/mutations/angela/filters/leaf_mutations.tsv")
rm(leaf_mutations, trunk_pot_mut) ; invisible(gc())
```

Only 30 to 89 mutations passed evs filtering across samples (Fig. \@ref(fig:nmutTrunkFiltTab)).
Most of mutations are low frequency (Fig. \@ref(fig:afLeafFiltFig)).
Most of mutations are not shared by biological replicates (Fig. \@ref(fig:vennLeafFiltevsFig) and Fig. \@ref(fig:vennLeafFiltbaseFig)).

```{r nmutLeafFiltTab, fig.cap="Number of filtered mutations per leaf."}
vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  group_by(tumor, branch, tip, leaf, Filter) %>% 
  summarise(N = n()) %>% 
  # filter(Filter == "evs") %>% summary()
  ggplot(aes(branch, N, fill = tip)) +
  geom_boxplot() +
  facet_wrap(~ Filter, scales = "free_x") +
  coord_flip() +
  ylab("Number of mutations") + xlab("Branch")   
```

```{r afLeafFiltFig, fig.cap="Alleles-frequencies of filtered mutations in leaf samples."}
vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  ggplot(aes(mutation_AF, col = tumor), fill = NA) +
  geom_density(size = 1.3) +
  scale_color_discrete("Sample") +
  xlab("Allelic frequency") +
  facet_wrap(~ Filter)     
```

```{r vennLeafFiltevsFig, fig.cap="Filtered mutations sharing across leaf samples within tips for evs filtering.", fig.height=8}
t <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% filter(Filter == "evs")
g1 <- lapply(LETTERS[1:5], function(B)
  ggvenn::ggvenn(lapply(split(filter(t, branch == B, tip == "S"),
                              filter(t, branch == B, tip == "S")$tumor),
                        `[[`, "SNV"), show_percentage = F)
)
g2 <- lapply(LETTERS[1:5], function(B)
  ggvenn::ggvenn(lapply(split(filter(t, branch == B, tip == "L"),
                              filter(t, branch == B, tip == "L")$tumor),
                        `[[`, "SNV"), show_percentage = F)
)
cowplot::plot_grid(plotlist = c(g1,g2))     
```

```{r vennLeafFiltbaseFig, fig.cap="Filtered mutations sharing across leaf samples within tips for base filtering.", fig.height=8}
t <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv")
g1 <- lapply(LETTERS[1:5], function(B)
  ggvenn::ggvenn(lapply(split(filter(t, branch == B, tip == "S"),
                              filter(t, branch == B, tip == "S")$tumor),
                        `[[`, "SNV"), show_percentage = F)
)
g2 <- lapply(LETTERS[1:5], function(B)
  ggvenn::ggvenn(lapply(split(filter(t, branch == B, tip == "L"),
                              filter(t, branch == B, tip == "L")$tumor),
                        `[[`, "SNV"), show_percentage = F)
)
cowplot::plot_grid(plotlist = c(g1,g2))     
```

## Mutations architecture

```{r archiMutationsEVSFig}
archi_mutations <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  filter(Filter == "evs") %>% 
  mutate(tip = paste0(branch, tip)) %>% 
  group_by(SNV, CHROM, POS, REF, ALT, branch, tip) %>% 
  summarise(mutation_AF = mean(mutation_AF)) %>% 
  group_by(SNV, CHROM, POS, REF, ALT) %>% 
  summarise(tip = paste(tip, collapse = " "), mutation_AF = mean(mutation_AF), n_tips = n(),
             n_branches = length(unique(branch)), branch = paste(sort(unique(branch)), collapse = " ")) %>% 
  ungroup() %>% 
  arrange(desc(n_tips), desc(mutation_AF)) %>% 
  mutate(origin1 = ifelse(n_tips == 1, tip, "branch")) %>% 
  mutate(origin1 = ifelse(origin1 == "branch" & n_branches == 1, branch, origin1)) %>% 
  mutate(origin1 = ifelse(origin1 == "branch" & branch == "A B", "AB", origin1)) %>% 
  mutate(origin1 = ifelse(origin1 == "branch" & branch == "C D", "CD", origin1)) %>% 
  mutate(origin1 = ifelse(origin1 == "branch", "X", origin1)) %>% 
  mutate(origin2 = ifelse(n_tips == 1, "tip", "branch")) %>% 
  mutate(origin2 = ifelse(origin1 %in% c("CD", "AB"), "carpenter", origin2)) %>% 
  mutate(origin2 = ifelse(origin1 == "X", "crown", origin2))
angela <- 
  "((((AS:3.83,AL:6.04)A:6.44,(BS:3.34,BL:10.64)B:1.6)AB:6.5,((CS:3.8,CL:5.13)C:14.45,(DS:3.5,DL:7.76)D:6.2)CD:1.7,(ES:0,EL:2.2)E:3.7)X:26.5)T:1.3;"
angela <- phylobase::phylo4d(read.tree(text = angela))
angela@data <- data.frame(label = angela@label) %>% 
  left_join(archi_mutations %>% 
              group_by(origin1) %>% 
              summarise(mutations = n()) %>% 
              rename(label = origin1)) %>% 
  mutate(mutations = ifelse(is.na(mutations), 0, mutations)) %>% 
  dplyr::select(-label)
ggtree(angela, right = T) +
  geom_point2(aes(size = mutations, col = mutations)) +
  geom_tiplab(offset = 1) +
  geom_nodelab(nudge_x = 0.8) +
  geom_label(aes(label = mutations, fill = mutations), vjust=1.5, col = "white") +
  theme_tree() +
  theme(legend.position = "right") +
  xlim(20, 50) +
  viridis::scale_color_viridis(limits = c(1, 260), trans = scales::log10_trans()) +
  viridis::scale_fill_viridis(limits = c(1, 260), trans = scales::log10_trans()) +
  scale_size_continuous(guide = "none", limits = c(1, 260)) +
  coord_flip() +
  ggtitle("Mutations origin in the crown", "At least in 2-leaves with evs filtering.")
```


```{r archiMutationsEVSFigAF}
vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  left_join(dplyr::select(archi_mutations, SNV, origin2)) %>% 
  group_by(origin2, SNV) %>% 
  filter(!is.na(origin2)) %>% 
  ggplot(aes(x = SNV, y = mutation_AF, color = SNV)) +
  geom_boxplot() +
  scale_color_discrete(guide = "none") +
  coord_flip() +
  ylab("Allelic frequencies (AF)") + xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.line.y = element_blank()) +
  theme(legend.position = "bottom", legend.title = element_blank()) +
  facet_wrap(~origin2, nrow = 1, scales = "free_y")
```


```{r archiMutationsBaseFig}
archi_mutations <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  mutate(tip = paste0(branch, tip)) %>% 
  group_by(SNV, CHROM, POS, REF, ALT, branch, tip) %>% 
  summarise(mutation_AF = mean(mutation_AF)) %>% 
  group_by(SNV, CHROM, POS, REF, ALT) %>% 
  summarise(tip = paste(tip, collapse = " "), mutation_AF = mean(mutation_AF), n_tips = n(),
             n_branches = length(unique(branch)), branch = paste(sort(unique(branch)), collapse = " ")) %>% 
  ungroup() %>% 
  arrange(desc(n_tips), desc(mutation_AF)) %>% 
  mutate(origin1 = ifelse(n_tips == 1, tip, "branch")) %>% 
  mutate(origin1 = ifelse(origin1 == "branch" & n_branches == 1, branch, origin1)) %>% 
  mutate(origin1 = ifelse(origin1 == "branch" & branch == "A B", "AB", origin1)) %>% 
  mutate(origin1 = ifelse(origin1 == "branch" & branch == "C D", "CD", origin1)) %>% 
  mutate(origin1 = ifelse(origin1 == "branch", "X", origin1)) %>% 
  mutate(origin2 = ifelse(n_tips == 1, "tip", "branch")) %>% 
  mutate(origin2 = ifelse(origin1 %in% c("CD", "AB"), "carpenter", origin2)) %>% 
  mutate(origin2 = ifelse(origin1 == "X", "crown", origin2))
angela <- 
  "((((AS:3.83,AL:6.04)A:6.44,(BS:3.34,BL:10.64)B:1.6)AB:6.5,((CS:3.8,CL:5.13)C:14.45,(DS:3.5,DL:7.76)D:6.2)CD:1.7,(ES:0,EL:2.2)E:3.7)X:26.5)T:1.3;"
angela <- phylobase::phylo4d(read.tree(text = angela))
angela@data <- data.frame(label = angela@label) %>% 
  left_join(archi_mutations %>% 
              group_by(origin1) %>% 
              summarise(mutations = n()) %>% 
              rename(label = origin1)) %>% 
  mutate(mutations = ifelse(is.na(mutations), 0, mutations)) %>% 
  dplyr::select(-label)
ggtree(angela, right = T) +
  geom_point2(aes(size = mutations, col = mutations)) +
  geom_tiplab(offset = 1) +
  geom_nodelab(nudge_x = 0.8) +
  geom_label(aes(label = mutations, fill = mutations), vjust=1.5, col = "white") +
  theme_tree() +
  theme(legend.position = "right") +
  xlim(20, 50) +
  viridis::scale_color_viridis(limits = c(1, 10849), trans = scales::log10_trans()) +
  viridis::scale_fill_viridis(limits = c(1, 10849), trans = scales::log10_trans()) +
  scale_size_continuous(guide = "none", limits = c(1, 10849)) +
  coord_flip() +
  ggtitle("Mutations origin in the crown", "At least in 2-leaves with base filtering.")
```

```{r archiMutationsEVSFig2AF}
vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  left_join(dplyr::select(archi_mutations, SNV, origin2)) %>% 
  filter(!is.na(origin2)) %>%
  sample_n(10^3) %>% 
  group_by(origin2, SNV) %>% 
  
  ggplot(aes(x = SNV, y = mutation_AF, color = SNV)) +
  geom_boxplot() +
  scale_color_discrete(guide = "none") +
  coord_flip() +
  ylab("Allelic frequencies (AF)") + xlab("") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.line.y = element_blank()) +
  theme(legend.position = "bottom", legend.title = element_blank()) +
  facet_wrap(~origin2, nrow = 1, scales = "free_y")
```

```{r archiMutationsEVSFig2AF2}
vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  left_join(dplyr::select(archi_mutations, SNV, origin2)) %>% 
  filter(!is.na(origin2)) %>% 
  ggplot(aes(x = mutation_AF, color = origin2)) +
  geom_density() + scale_x_log10() +
  xlab("Allelic frequencies (AF)")
```

## Phylogeny

```{r tipsvcfs, eval=F}
library(vcfR)
mutations <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  mutate(tumor = paste0(branch, tip)) %>% 
  select(tumor, SNV, CHROM, POS, REF, ALT, Filter) %>% 
  unique()
# vcf header
ref_mut_vcf <- read.vcfR("data/mutations/angela/raw/AL1_vs_T2.raw.vcf")
header <- ref_mut_vcf@meta
header[64] <- "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">"
header[65] <- "##FILTER=<ID=base,Description=\"Filter based on minimum allele count for mutation (>10) and normal sample (=0), sequencing depth (5th-95th), and mutations sharing across leaves (at least 2 in the crown).\">"
header[66] <- "##FILTER=<ID=evs,Description=\"Filter based on minimum allele count for mutation (>10) and normal sample (=0), sequencing depth (5th-95th), and mutations sharing across leaves (at least 2 in the crown).\">"
header <- header[-c(67:88)]
# filters
tsv2vcf <- function(data){
  fix <- data %>% 
    mutate(QUAL = NA, INFO=NA) %>% 
    dplyr::rename(ID = SNV, FILTER = Filter) %>% 
    dplyr::select(CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO) %>% 
    arrange(CHROM, POS) %>% 
    mutate(POS = as.character(POS))
  # head(as.matrix(fix))
  gt <- dplyr::select(data, SNV) %>% 
    mutate(tumor = list(unique(data$tumor))) %>% 
    unnest() %>% 
    left_join(select(data, SNV, tumor, Filter)) %>% 
    mutate(GT = ifelse(is.na(Filter), "0/0", "1/1")) %>% 
    mutate(FORMAT = "GT") %>% 
    dplyr::select(SNV, FORMAT, GT, tumor) %>% 
    unique() %>% 
    reshape2::dcast(SNV + FORMAT ~ tumor, value.var = "GT") %>% 
    dplyr::rename(ID = SNV) %>% 
    mutate(R = "0/0")
  gt <- left_join(dplyr::select(fix, ID), gt) %>% 
    dplyr::select(-ID)
  # head(as.matrix(gt))
  mutation_vcf <- ref_mut_vcf
  mutation_vcf@meta <- header
  mutation_vcf@fix <- as.matrix(fix)
  mutation_vcf@gt <- as.matrix(gt)
  return(mutation_vcf)
}
write.vcf(tsv2vcf(mutations), "data/mutations/angela/trees/leaf_mutations.vcf.gz")
# wget https://raw.githubusercontent.com/edgardomortiz/vcf2phylip/master/vcf2phylip.py
# for file in $(ls *.vcf.gz); do python3 vcf2phylip.py -i $file -n ; done
files <- list.files(path = "data/mutations/angela/trees/", pattern = ".nexus")
invisible(lapply(files, function(f){
  mutations <- phangorn::read.phyDat(file.path("data/mutations/angela/trees", f), format = "nexus")
  mutations <- ape::as.DNAbin(mutations)
  ips::write.fas(mutations, file.path("data/mutations/angela/trees", gsub("nexus", "fa", f)))
}))
# for file in $(ls *.min4.fa); do iqtree2 -s $file -m TIM2+I+G -B 1000 -redo -o R ; done
```

```{r angleatipsphylo, fig.width=12}
files <- list.files(path = "data/mutations/angela/trees/", pattern = ".min4.fa.treefile")
phylos <- lapply(files, function(f) read.tree(file.path("data/mutations/angela/trees", f)))
names(phylos) <- gsub(".min4.fa.treefile", "", files)
cowplot::plot_grid(plotlist = lapply(phylos, function(p) 
  ggtree(phylobase::phylo4d(p)) + geom_tippoint(aes(col = substr(label, 1, 1)), size = 4) + geom_tiplab()
), labels = names(phylos)) 
```

```{r rev2phyloangela, eval=F}
library(vcfR)
library(ape)
library(ggtree)
mutations <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  arrange(branch, tip, SNV) %>% 
  group_by(branch, tip, SNV) %>% 
  filter(n() == 3) %>% 
  ungroup() %>% 
  mutate(tumor = paste0(branch, tip)) %>% 
  select(tumor, SNV, CHROM, POS, REF, ALT, Filter) %>% 
  unique()
# vcf header
ref_mut_vcf <- read.vcfR("data/mutations/angela/raw/AL1_vs_T2.raw.vcf")
header <- ref_mut_vcf@meta
header[64] <- "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">"
header[65] <- "##FILTER=<ID=base,Description=\"Filter based on minimum allele count for mutation (>10) and normal sample (=0), sequencing depth (5th-95th), and mutations sharing across leaves (at least 2 in the crown).\">"
header[66] <- "##FILTER=<ID=evs,Description=\"Filter based on minimum allele count for mutation (>10) and normal sample (=0), sequencing depth (5th-95th), and mutations sharing across leaves (at least 2 in the crown).\">"
header <- header[-c(67:88)]
# filters
tsv2vcf <- function(data){
  fix <- data %>% 
    mutate(QUAL = NA, INFO=NA) %>% 
    dplyr::rename(ID = SNV, FILTER = Filter) %>% 
    dplyr::select(CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO) %>% 
    arrange(CHROM, POS) %>% 
    mutate(POS = as.character(POS))
  # head(as.matrix(fix))
  gt <- dplyr::select(data, SNV) %>% 
    mutate(tumor = list(unique(data$tumor))) %>% 
    unnest() %>% 
    left_join(select(data, SNV, tumor, Filter)) %>% 
    mutate(GT = ifelse(is.na(Filter), "0/0", "1/1")) %>% 
    mutate(FORMAT = "GT") %>% 
    dplyr::select(SNV, FORMAT, GT, tumor) %>% 
    unique() %>% 
    reshape2::dcast(SNV + FORMAT ~ tumor, value.var = "GT") %>% 
    dplyr::rename(ID = SNV) %>% 
    mutate(R = "0/0")
  gt <- left_join(dplyr::select(fix, ID), gt) %>% 
    dplyr::select(-ID)
  # head(as.matrix(gt))
  mutation_vcf <- ref_mut_vcf
  mutation_vcf@meta <- header
  mutation_vcf@fix <- as.matrix(fix)
  mutation_vcf@gt <- as.matrix(gt)
  return(mutation_vcf)
}
write.vcf(tsv2vcf(mutations), "data/mutations/angela/rev2//leaf_mutations.vcf.gz")
# wget https://raw.githubusercontent.com/edgardomortiz/vcf2phylip/master/vcf2phylip.py
# for file in $(ls *.vcf.gz); do python3 vcf2phylip.py -i $file -n ; done
files <- list.files(path = "data/mutations/angela/rev2/", pattern = ".nexus")
invisible(lapply(files, function(f){
  mutations <- phangorn::read.phyDat(file.path("data/mutations/angela/rev2", f), format = "nexus")
  mutations <- ape::as.DNAbin(mutations)
  ips::write.fas(mutations, file.path("data/mutations/angela/rev2", gsub("nexus", "fa", f)))
}))
```

```{r angleatipsphylo2, fig.width=12}
files <- list.files(path = "data/mutations/angela/trees/", pattern = ".min4.fa.treefile")
phylos <- lapply(files, function(f) read.tree(file.path("data/mutations/angela/trees", f)))
names(phylos) <- gsub(".min4.fa.treefile", "", files)
cowplot::plot_grid(plotlist = lapply(phylos, function(p) 
  ggtree(phylobase::phylo4d(p)) + geom_tippoint(aes(col = substr(label, 1, 1)), size = 4) + geom_tiplab()
), labels = names(phylos)) 
```

## Light & annotation

```{r, eval=F}
library(Biostrings)
library(genomeIntervals)
genome <- readDNAStringSet("data/mutations/angela/genome/Dgu_HS1_HYBRID_SCAFFOLD.fa")
genome_index <- read_tsv("data/mutations/angela/genome/Dgu_HS1_HYBRID_SCAFFOLD.fa.fai",
                         col_names = c("CHROM", "length", "bytesindex", "basesperline", "bytesperline")) %>% 
  dplyr::select(CHROM, length) %>% 
  mutate(start = 0, end = length) %>% 
  dplyr::select(-length)
mutations <- vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  group_by(SNV, CHROM, POS, REF, ALT) %>% 
  unique() %>% 
  ungroup()
seqs <- lapply(1:nrow(mutations), function(i)
  as.character(unlist(extractAt(genome[mutations$CHROM[i]], IRanges(mutations$POS[i]-1, width=3)))))
t <- mutations %>% 
  mutate(spectra = unlist(seqs)) %>% 
  mutate(REF5 = str_sub(spectra, 1, 1), REFspectra = str_sub(spectra, 2, 2), REF3 = str_sub(spectra, 3, 3)) %>% 
  dplyr::select(-REFspectra) %>% 
  mutate(type = paste0(REF, "->", ALT)) %>% 
  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(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(POS5 = REF5, POS3 = REF3) %>% 
  mutate(POS5r = recode(POS5, "A" = "T", "T" = "A", "C" = "G", "G" = "C")) %>%
  mutate(REFr = recode(REF, "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(REFf = ifelse(strand == "+", REF, REFr)) %>%
  mutate(POS3f = ifelse(strand == "+", POS3, POS5r)) %>% 
  mutate(spectra = paste0(POS5f, REFf, POS3f)) %>% 
  dplyr::select(CHROM, POS, SNV, REF, ALT, tumor, branch, tip, mutation_AF, Filter, type, spectra)
write_tsv(t, "data/mutations/angela/filters/leaf_mutations_spectra.tsv") 
```

```{r, fig.cap="Number of mutations per filter and light condition across branches."}
 vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
   group_by(branch, tip, Filter) %>%
  summarise(N = n()) %>%
  reshape2::dcast(branch + tip ~ Filter, value.var = "N") %>%
  mutate_all(funs(ifelse(is.na(.), 0, .))) %>% 
  mutate(base = base + evs) %>% 
  reshape2::melt(c("branch", "tip"), variable.name = "filter", value.name = "mutations") %>% 
  mutate(condition = recode(tip, "L" = "Light", "S" = "Shadow")) %>% 
  ggplot(aes(filter, mutations, fill = condition)) +
  geom_boxplot() +
  scale_y_sqrt() +
  xlab("") + ylab("Number of mutations") + scale_color_discrete("")  +
  facet_wrap(~ filter, scales = "free") +
  ggpubr::stat_compare_means(label = "p.signif", method = "t.test")
```

```{r}
vroom::vroom("data/mutations/angela/filters/leaf_mutations_spectra.tsv") %>% 
  filter(!(Filter %in% c("base1", "evs1"))) %>% 
  mutate(condition = recode(tip, "L" = "Light", "S" = "Shadow")) %>% 
  mutate(sample = paste0(branch, tip)) %>% 
  group_by(type, sample, condition) %>% 
  summarise(N = n()) %>% 
  group_by(sample) %>% 
  ggplot(aes(type, N, fill = condition)) +
  geom_boxplot(position = "dodge") + 
  theme(axis.text.x = element_blank(), axis.title.x = element_blank(),
        axis.line.x = element_blank(), axis.ticks.x = element_blank(),
        legend.position = "bottom") +
  facet_wrap(~ type, nrow = 1, scales = "free_x") +
  xlab("Percentage of mutations (%) per sample") +
  ggpubr::stat_compare_means(label = "p.signif", method = "t.test") +
  ggtitle("Mutation type", "At least in 2-leaves with evs filtering.") +
  ylab("Number of mutations")
```

```{r, fig.width=12}
vroom::vroom("data/mutations/angela/filters/leaf_mutations_spectra.tsv") %>% 
    filter(!(Filter %in% c("base1", "evs1"))) %>% 
  mutate(condition = recode(tip, "L" = "Light", "S" = "Shadow")) %>% 
  mutate(sample = paste0(branch, tip)) %>% 
  group_by(spectra) %>% 
  group_by(type, spectra, sample, condition) %>% 
  summarise(N = n()) %>% 
  ggplot(aes(spectra, N, fill = condition)) +
  geom_boxplot(position = "dodge") + 
  facet_wrap(~ type, scales = "free_x") +
  theme(axis.text.x = element_text(angle = 90), legend.position = "bottom", axis.title.x = element_blank()) +
  theme(panel.spacing = unit(0, "lines"), panel.border = element_rect(size = 0.25, color = "black", fill = NA)) +
  scale_y_log10() +
  ylab("Number of mutations") + ggtitle("Mutation spectra", "At least in 2-leaves with evs filtering.")
```

## Rate

```{r}
vroom::vroom("data/mutations/angela/filters/leaf_mutations.tsv") %>% 
  mutate(Filter = recode(Filter, "evs" = "evs base")) %>% 
  separate_rows(Filter) %>% 
  group_by(Filter) %>% 
  arrange(desc(mutation_AF)) %>% 
  mutate(mutations = 1) %>% 
  mutate(cummulated_mutations = cumsum(mutations)) %>% 
  ggplot(aes(mutation_AF, cummulated_mutations, col = Filter)) +
  geom_line() +
  scale_y_log10() + scale_x_reverse() +
  xlab("Allelic frequency") + ylab("Cummulated mutations with decreasing AF")
```

## Mutations for Fruits

```{r fruitMutations}
fruit_mutations <- vroom::vroom("data/mutations/angela/filters/filtered_tip_mutations_spectra.tsv") %>% 
  filter(branch %in% c("A", "B", "C", "D")) %>% 
  filter(Filter == "evs") %>% 
  group_by(SNV, CHROM, POS, REF, ALT) %>% 
  summarise(af = mean(mutation_AF), replicate = n(), branch = paste(sort(unique(branch)), collapse = ", ")) %>% 
  ungroup() %>% 
  filter(replicate > 2) %>% 
  arrange(desc(replicate), desc(af)) %>% 
  mutate(rank = 1:n())    
# write_tsv(fruit_mutations, "data/mutations/angela/filters/angela_fruits_candidate_mutations.tsv")
```

Selection of mutations for transmission in fruits.

```{r fruitMutationsTab}
fruit_mutations %>% 
  group_by(replicate) %>% 
  summarise(N = n(), minAF = min(af), meanAF = mean(af), maxAF = max(af)) %>% 
  ungroup() %>% 
  arrange(desc(replicate)) %>% 
  mutate_all(funs(ifelse(is.na(.), 0, .))) %>% 
  # dplyr::select(-n_branches) %>% 
  kable(caption = "Number of mutations per branch.", digits = 2)   
```

```{r fruitMutationsFig}
ggplot(filter(hs1, name %in% gsub("Super-Scaffold_", "", unique(fruit_mutations$CHROM))),
       aes(x = name, xend = name)) +
  geom_segment(aes(y = start/10^6, yend = end/10^6), size = 3, col = "grey") +
  coord_flip() +
  ylab("Position (Mb)") +
  geom_point(aes(y = POS/10^6, col = af, size = replicate), 
             data = mutate(fruit_mutations, name = gsub("Super-Scaffold_", "", CHROM))) +
  theme(axis.line.y = element_blank(), axis.title.y = element_blank(), 
        axis.ticks.y = element_blank()) +
  viridis::scale_color_viridis("Allele\nfrequency", trans = scales::log_trans()) +
  scale_size_continuous("Number\nof branches") +
  ggtitle(paste(nrow(fruit_mutations), "candidate mutations for transmission to fruits."))    
```

```{r fruitMutationsTab2}
fruit_mutations %>% 
  dplyr::select(-SNV) %>% 
  mutate(af = round(af, 2)) %>% 
  dplyr::select(CHROM, POS, branch, replicate, af, rank) %>% 
  DT::datatable(colnames = c("Super scaffold", "position", "branch(es)", "replicate", "allele frequency", "rank")) 
```

## Genes

```{r, eval=F}
genes <- read_tsv("data/mutations/angela/annotation/braker.gtf",
                   col_names = c("CHROM", "source", "type", "start", "end", "score", "strand", "phase", "attributes")) %>% 
  filter(type == "transcript")
t <- archi_mutations %>%
  split(.$SNV) %>%
  lapply(function(x) left_join(x, genes, by = "CHROM") %>%
           filter(start <= POS, end >= POS)) %>%
  bind_rows(.id = "SNV")
nrow(t)/nrow(archi_mutations)*100 # 2067 in trsc
annotation <- DBI::dbConnect(RSQLite::SQLite(), dbname = "data/mutations/angela/annotation/angela.annot.sqlite")
orfs <-  tbl(annotation, "ORF") %>%
  mutate(start_aa = str_sub(peptide, 1, 1)) %>%
  mutate(stop_aa = str_sub(peptide, -1, -1)) %>%
  mutate(length = rend - lend) %>%
  collect()
t <- t %>%
  left_join(orfs, by = c("attributes" = "transcript_id")) %>%
  filter(!is.na(orf_id))
rm(orfs)
DBI::dbDisconnect(annotation) ; rm(annotation)
write_tsv(t, "save/tmp.tsv")
read_tsv("save/tmp.tsv") %>% 
  mutate_at(c("lend", "rend"), as.numeric) %>% 
  mutate(posTrsc = POS - start + 1) %>% 
  filter(posTrsc >= lend, posTrsc <=   rend) %>% # 567 in CDS
  dplyr::select(attributes, posTrsc, ALT) %>% 
  write_tsv("data/mutations/angela/annotation/SNPsOnTrsc.tsv", col_names = F)
```

```{bash, eval=F, echo=T}
cd data/mutations/angela/annotation
~/Tools/gffread/gffread trsc.fa.transdecoder.gff3 -T -o trsc.fa.transdecoder.gtf
perl ~/Tools/SNPdat_package_v1.0.5/SNPdat_v1.0.5.pl \
  -i SNPsOnTrsc.tsv \
  -g trsc.fa.transdecoder.gtf \
  -f trsc.fa \
  -s synonymy.summary \
  -o synonymy.output
```

```{r, eval=F}
read_tsv("save/tmp.tsv") %>% 
  mutate_at(c("lend", "rend"), as.numeric) %>% 
  mutate(posTrsc = POS - start + 1) %>% 
  filter(posTrsc >= lend, posTrsc <=   rend) %>% 
  mutate(`Chromosome Number` = toupper(attributes), `SNP Position` = posTrsc) %>% 
  dplyr::select(SNV, `Chromosome Number`, `SNP Position`) %>% 
  left_join(read_tsv("data/mutations/angela/annotation/synonymy.output", col_names = T)) %>% 
  group_by(SNV, synonymous) %>% 
  summarise(N = n()) %>% 
  reshape2::dcast(SNV ~ synonymous) %>% 
  filter(!is.na(N)) %>% 
  nrow()
```

```{r goeAngela}
go <- read_tsv("data/mutations/angela/annotation/trinotate_go_annotations.txt", col_names = c("gene", "id")) %>% 
  separate_rows(id, sep = ",") %>% 
  left_join(src_sqlite("data/mutations/angela/annotation/angela.annot.sqlite") %>% 
              tbl("go") %>% 
              collect())
genes <- read_tsv("save/tmp.tsv") %>% 
  mutate_at(c("lend", "rend"), as.numeric) %>% 
  mutate(posTrsc = POS - start + 1) %>% 
  filter(posTrsc >= lend, posTrsc <=   rend) %>% 
  mutate(`Chromosome Number` = toupper(attributes), `SNP Position` = posTrsc) %>% 
  dplyr::select(SNV, `Chromosome Number`, `SNP Position`) %>% 
  left_join(read_tsv("data/mutations/angela/annotation/synonymy.output", col_names = T)) %>% 
  filter(synonymous == "N") %>% 
  separate(`Chromosome Number`, c("gene", "trsc_id"), remove = F) %>% 
  dplyr::select(gene) %>% 
  unique() %>% 
  mutate(gene = tolower(gene))
goe <- clusterProfiler::enricher(gene = unique(genes$gene),
                          pvalueCutoff = 0.2, 
                          pAdjustMethod = "BH",
                          universe = unique(go$gene),
                          minGSSize = 10, 
                          maxGSSize = 500, 
                          qvalueCutoff = 0.2, 
                          TERM2GENE = dplyr::select(go, id, gene),
                          TERM2NAME = dplyr::select(go, id, name))
```

```{r, fig.cap="Mutations types."}
data.frame(type = c("All", "Within transcripts with genes", "Within CDS", "Non-synonymous"),
           count = c(15468, 2069, 567, 314)) %>% 
  ggplot(aes(type, count)) +
  geom_col() +
  coord_flip() +
  scale_y_log10() +
  ggtitle("Type of mutations for Angela") +
  xlab("") +
  geom_text(aes(y = 2, label = format(count, big.mark = " ")), col = "white")
```

```{r, fig.cap="Mutation allelic fraction per type."}
# genes <- genomeIntervals::readGff3("data/mutations/angela/annotation/trsc.fa.transdecoder.gff3")
# genes <- genomeIntervals::interval_union(genes)
# sum(genomeIntervals::width(genes))/sum(hs1$end)*100
data <- archi_mutations %>% 
  mutate(all = 1) %>% 
  left_join(read_tsv("save/tmp.tsv") %>% 
              dplyr::select(SNV) %>% 
              mutate(transcript = 1)) %>% 
  left_join(read_tsv("save/tmp.tsv") %>% 
              mutate_at(c("lend", "rend"), as.numeric) %>% 
              mutate(posTrsc = POS - start + 1) %>% 
              filter(posTrsc >= lend, posTrsc <=   rend) %>% 
              dplyr::select(SNV) %>% 
              mutate(CDS = 1)) %>% 
  left_join(read_tsv("save/tmp.tsv") %>% 
              mutate_at(c("lend", "rend"), as.numeric) %>% 
              mutate(posTrsc = POS - start + 1) %>% 
              filter(posTrsc >= lend, posTrsc <=   rend) %>% 
              mutate(`Chromosome Number` = toupper(attributes), `SNP Position` = posTrsc) %>% 
              dplyr::select(SNV, `Chromosome Number`, `SNP Position`) %>% 
              left_join(read_tsv("data/mutations/angela/annotation/synonymy.output", col_names = T)) %>% 
              group_by(SNV, synonymous) %>% 
              summarise(N = n()) %>% 
              reshape2::dcast(SNV ~ synonymous) %>% 
              filter(!is.na(N)) %>% 
              dplyr::select(SNV) %>% 
              mutate(nonsynonymous = 1)) %>% 
  left_join(filter(goe@result, qvalue < 0.05) %>% 
              dplyr::select(geneID) %>% 
              separate_rows(geneID) %>% 
              unique() %>% 
              left_join(read_tsv("save/tmp.tsv") %>% 
                          dplyr::select(SNV, attributes) %>% 
                          separate(attributes, c("geneID", "transcript"))) %>% 
              dplyr::select(SNV) %>% 
              mutate(goenriched = 1))
g1 <- data %>% 
  filter(mutation_AF > 0.01) %>% 
  dplyr::select(SNV, mutation_AF, CDS) %>% 
  unique() %>% 
  mutate(type = ifelse(is.na(CDS), "intergenic", "CDS")) %>% 
  ggplot(aes(type, mutation_AF)) +
  geom_violin() +
  ggpubr::stat_compare_means(method = "t.test") +
  ylab("Allelic fraction") +
  ggtitle("All mutations AF > 0.01", "ratio:0.034 (<expected:0.139)") +
  geom_point(col = "red", size = 2,
             data = data.frame(mutation_AF = c(median(filter(data, is.na(CDS))$mutation_AF),
                                               median(filter(data, CDS == 1)$mutation_AF)),
                               type = c("intergenic", "CDS")))
g2 <- filter(data, CDS == 1) %>% 
  filter(mutation_AF > 0.01) %>% 
  dplyr::select(SNV, mutation_AF, nonsynonymous) %>% 
  mutate(synonymy = ifelse(is.na(nonsynonymous), "S", "NS")) %>% 
  dplyr::select(SNV, synonymy, mutation_AF) %>% 
  unique() %>% 
  ggplot(aes(synonymy, mutation_AF)) +
  geom_violin() +
  ggpubr::stat_compare_means(method = "t.test") +
  ylab("Allelic fraction") +
  ggtitle("Mutations in CDS, AF > 0.01", "dN/dS=1.67") +
  geom_point(col = "red", size = 2,
             data = data.frame(mutation_AF = c(median(filter(data, CDS == 1, nonsynonymous == 1)$mutation_AF),
                                               median(filter(data, CDS == 1, is.na(nonsynonymous))$mutation_AF)),
                                            synonymy = c("NS", "S")))
cowplot::plot_grid(g2, g1)
```

```{r, fig.cap="Synonymy with allelic fraction."}
data <- archi_mutations %>% 
  mutate(all = 1) %>% 
  left_join(read_tsv("save/tmp.tsv") %>% 
              dplyr::select(SNV) %>% 
              mutate(transcript = 1)) %>% 
  left_join(read_tsv("save/tmp.tsv") %>% 
              mutate_at(c("lend", "rend"), as.numeric) %>% 
              mutate(posTrsc = POS - start + 1) %>% 
              filter(posTrsc >= lend, posTrsc <=   rend) %>% 
              dplyr::select(SNV) %>% 
              mutate(CDS = 1)) %>% 
  left_join(read_tsv("save/tmp.tsv") %>% 
              mutate_at(c("lend", "rend"), as.numeric) %>% 
              mutate(posTrsc = POS - start + 1) %>% 
              filter(posTrsc >= lend, posTrsc <=   rend) %>% 
              mutate(`Chromosome Number` = toupper(attributes), `SNP Position` = posTrsc) %>% 
              dplyr::select(SNV, `Chromosome Number`, `SNP Position`) %>% 
              left_join(read_tsv("data/mutations/angela/annotation/synonymy.output", col_names = T)) %>% 
              group_by(SNV, synonymous) %>% 
              summarise(N = n()) %>% 
              reshape2::dcast(SNV ~ synonymous) %>% 
              filter(!is.na(N)) %>% 
              dplyr::select(SNV) %>% 
              mutate(nonsynonymous = 1)) %>% 
  left_join(filter(goe@result, qvalue < 0.05) %>% 
              dplyr::select(geneID) %>% 
              separate_rows(geneID) %>% 
              unique() %>% 
              left_join(read_tsv("save/tmp.tsv") %>% 
                          dplyr::select(SNV, attributes) %>% 
                          separate(attributes, c("geneID", "transcript"))) %>% 
              dplyr::select(SNV) %>% 
              mutate(goenriched = 1))
write_tsv(data, "save/mutations_angela_sixto.tsv")
filter(data, CDS == 1) %>% 
  filter(mutation_AF > 0.01) %>% 
  dplyr::select(SNV, mutation_AF, nonsynonymous) %>% 
  mutate(synonymy = ifelse(is.na(nonsynonymous), "S", "NS")) %>% 
  dplyr::select(SNV, synonymy, mutation_AF) %>% 
  unique() %>% 
  mutate(AF = cut(mutation_AF, seq(0, 0.5, by = 0.005))) %>% 
  group_by(AF, synonymy) %>% 
  summarise(N = n()) %>% 
  reshape2::dcast(AF ~ synonymy) %>% 
  mutate_at(c("NS", "S"), funs(ifelse(is.na(.),0,.))) %>% 
  mutate(dNds = NS/S, N = NS+S) %>% 
  separate(AF, c("X1", "X2"), ",") %>% 
  mutate(X1 = as.numeric(gsub("(", "", X1, fixed = T))) %>% 
  mutate(X2 = as.numeric(gsub("]", "", X2, fixed = T))) %>% 
  mutate(AF = (X1+X2)/2) %>% 
  ggplot(aes(AF, dNds, size = N)) +
  geom_smooth(method = "lm", size = 1) +
  geom_point() +
  geom_hline(yintercept = 1, linetype = "dashed")
```

```{r, fig.cap="GO enrichment for Angela non-synonymous mutations."}
goe@result %>% 
  filter(qvalue < 0.05) %>% 
  na.omit() %>% 
  ggplot(aes(reorder(Description, qvalue), Count, fill = -log(qvalue))) +
  geom_bar(stat = "identity") +
  coord_flip() +
  viridis::scale_fill_viridis() +
  theme(axis.title = element_blank()) +
  ggtitle("GO enrichment for Angela non-synonymous mutations")
```

```{r, fig.cap="Enrichment map of GO enrichment for Angela non-synonymous mutations.", eval=F}
enrichplot::emapplot(goe) +
  ggtitle("GO enrichment for Angela non-synonymous mutations")
```