Notebook_GrowthCoefficients_AgeAtOnset_PossibleAssociations.Rmd

---
title: "The longitudinal dynamics and natural history of clonal haematopoiesis"
subtitle: "Growth coefficient inference, age at onset estimation and possible associations"
author: 
    - José Guilherme de Almeida
    - Moritz Gerstung
output: 
  html_document:
    theme: lumen
    highlight: tango
    code_folding: hide
    df_print: paged
    toc: true
    toc_float: true
    toc_collapsed: true
    toc_depth: 4
    number_sections: true
--- 

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,dpi = 300)
```

```{r import_functions, include=FALSE}
source("Scripts/vaf_dynamics_functions.R")
model_file_name <- "models/model_ch.RDS"
model_id <- gsub('\\.','',str_match(model_file_name,'model_.*\\.'))
dir.create("figures/data_exploration/",showWarnings = F)
dir.create("data_output/",showWarnings = F)
set.seed(42)
```

# Description

We have a cohort with roughly $400$ individuals from Sardinia. These individuals were target-sequenced at a high depth (1000X), which enables us to detect clonal haematopoiesis with high precision, and at different time points, thus enabling the modelling of clonal dynamics. We can, from this and using some assumptions regarding clonal evolution which have validated using simulations, infer the age at onset for the clones captured by our model. 

# Data preparation

Data preparation here is concerned mostly with the following aspects:

* Selecting genes based on their significance for our analysis 
* Excluding individuals who have/had haematological malignancies
* Determining mutations which occur across more than 1 individual

```{r data_preparation}
source("Scripts/prepare_data.R")
```

# Data exploration

An initial inspection of the trajectories does not allow us to make any ambitious claims on their dynamics considering age as the only determinant of dynamics. 

Some conclusions can be drawn at this point:

* No apparent single dynamics underlies most VAF trajectories when considering the age of individuals, but this becomes neater when using relative timepoints (relative time instead of absolute time) - this points towards the necessity of an individual-based offset. This also makes the likeliness of other factors (i.e. other mutations) affecting the dynamics of VAF higher

```{r plot_vaf_trajectories_per_individual, fig.width=6, fig.height=5}
full_data %>% 
  subset(single_occurring == F & truncating == F) %>% 
  mutate(q005=qbeta(p=0.05,MUTcount_Xadj+1,TOTALcount+1-MUTcount_Xadj),
         q095=qbeta(p=0.95,MUTcount_Xadj+1,TOTALcount+1-MUTcount_Xadj)) %>% 
  select(TOTALcount,MUTcount_Xadj,SardID,Age,VAF,q005,q095,Phase,amino_acid_change,Gene) %>% 
  rowwise() %>% 
  ggplot(aes(x = Age,y = MUTcount_Xadj/(TOTALcount),ymin = q005,ymax = q095,color = as.factor(SardID))) + 
  geom_linerange(size = 0.25) + 
  geom_line(size = 0.25) + 
  geom_point(size = 0.5) + 
  facet_wrap(~ amino_acid_change,scales = 'free',ncol=6) + 
  theme_gerstung(base_size=6) + 
  scale_color_discrete(guide=FALSE) +
  ylab("VAF") + 
  ggtitle("Trajectory view for each site with more than two mutations (one color per individual)") + 
  ggsave(useDingbats=FALSE,"figures/data_exploration/vaf_trajectories_per_individual.pdf",width=6, height=5)
```

```{r plot_vaf_trajectories_per_site, fig.width=6, fig.height=6}
full_data %>% #subset(single_occurring == F) %>% 
  mutate(q005=qbeta(p=0.05,MUTcount_Xadj+1,TOTALcount+1-MUTcount_Xadj),
         q095=qbeta(p=0.95,MUTcount_Xadj+1,TOTALcount+1-MUTcount_Xadj)) %>% 
  select(TOTALcount,MUTcount_Xadj,SardID,Age,VAF,q005,q095,Phase,amino_acid_change) %>% 
  group_by(SardID) %>% 
  filter(any(VAF > 0.05)) %>%
  ggplot(aes(x = Age,y = MUTcount_Xadj/(TOTALcount),ymin = q005,ymax = q095,color = amino_acid_change)) + 
  geom_linerange(size = 0.25) + 
  geom_line(size = 0.25) + 
  geom_point(size = 0.5) + 
  facet_wrap(~ SardID,scales = 'free') + 
  theme_gerstung(base_size=6) + 
  scale_color_discrete(guide=FALSE) +
  ylab("VAF") + 
  ggtitle("Trajectory view for each individual (one color per mutation)") + 
  scale_x_continuous(n.breaks = 4,
                     breaks = function(x) floor(x[1]) + seq(1,x[2],by = round((x[2] - x[1])/5))) + 
  ggsave(useDingbats=FALSE,"figures/data_exploration/vaf_trajectories_per_site.pdf",width=6, height=6)
```

With this data we can also investigate the relevance and prevalence of co-occurring mutations. With these data, we find no evidence for statistically significant co-occurrence of mutations in either genes or sites.

```{r cooccurrence, fig.width=7,fig.height=5}
unpack_square_matrix <- function(square_matrix,V) {
  square_matrix <- as.data.frame(square_matrix)
  colnames(square_matrix) <- V
  square_matrix$labels_2 <- V
  square_matrix %>% 
    gather(key = "labels_1",value = "value",-labels_2) %>%
    return
}

tmp <- full_data %>%
  select(amino_acid_change,SardID,Gene,single_occurring) %>%
  mutate(Gene = str_match(Gene,"[A-Z0-9]+")) %>% 
  distinct()
genes <- names(gene_colours) %>% 
  str_match('[A-Z0-9]+') %>%
  unique %>%
  sort
usm <- unique(
  full_data$amino_acid_change[full_data$single_occurring == F])

co_occurrence_genes <- matrix(0,length(genes),length(genes))
co_occurrence_sites <- matrix(0,length(usm),length(usm))
for (individual in unique(full_data$SardID)) {
  individual_data <- tmp %>% 
    subset(SardID==individual) 
  individual_sites <- individual_data %>%
    subset(single_occurring == F) %>% 
    select(amino_acid_change) %>% 
    unlist() %>%
    unique
  individual_genes <- individual_data %>%
    select(Gene) %>% 
    unlist %>% 
    unique
  if (length(individual_sites) >= 2) {
    site_combinations <- t(combn(individual_sites,2))
    for (i in 1:nrow(site_combinations)) {
      x <- site_combinations[i,]
      co_occurrence_sites[usm == x[1],
                          usm == x[2]] <- co_occurrence_sites[usm == x[1],usm == x[2]] + 1
    }
  }
  if (length(individual_genes) >= 2) {
    gene_combinations <- t(combn(individual_genes,2))
    for (i in 1:nrow(gene_combinations)) {
      x <- gene_combinations[i,]
      co_occurrence_genes[genes == x[1],
                          genes == x[2]] <- co_occurrence_genes[genes == x[1],genes == x[2]] + 1
    }
  }
}

data_unique_gene <- full_data %>%
  select(SardID,Gene) %>% 
  distinct
chi_squared_pvalues <- combn(unique(str_match(full_data$Gene,'[A-Z0-9]+')),2) %>% 
  t %>%
  apply(1,function(x) {
    tmp <- data_unique_gene %>% 
      group_by(SardID) %>%
      summarise(
        HasA = and(sum(Gene %in% x[1])>0,sum(Gene %in% x[2])==0),
        HasB = and(sum(Gene %in% x[1])==0,sum(Gene %in% x[2])>0),
        HasNone = and(sum(Gene %in% x[1])==0,sum(Gene %in% x[2])==0),
        HasBoth = and(sum(Gene %in% x[1])>0,sum(Gene %in% x[2])>0)
      ) %>%
      gather("key","value",HasA,HasB,HasNone,HasBoth) %>%
      group_by(key) %>%
      summarise(N = sum(value)) 
    O <- fisher.test(matrix(c(tmp$N[4],tmp$N[2],tmp$N[1],tmp$N[3]),nrow=2)) 
    return(c(O$estimate[1],pValue=O$p.value,
             GeneA=x[1],GeneB=x[2],
             HasGeneA=tmp$N[1],HasGeneB=tmp$N[2],
             HasNone=tmp$N[4],HasBoth=tmp$N[3]))
  }) %>% 
  t %>% 
  as_tibble() %>%
  mutate(`odds ratio` = as.numeric(`odds ratio`),
         pValue = as.numeric(pValue))

chi_squared_pvalues$pValue_adj <- p.adjust(chi_squared_pvalues$pValue)
chi_squared_pvalues %>%
  subset(pValue_adj < 0.05 & `odds ratio` < Inf)

total_counts <- full_data %>%
  select(SardID,labels_1=Gene) %>%
  mutate(labels_1 = str_match(labels_1,"[A-Z0-9]+")) %>%
  distinct() %>% 
  group_by(labels_1) %>% 
  summarise(value = length(SardID)) %>%
  mutate(labels_2 = 'Total',labels_1 = labels_1,value = value)

co_occurrence_genes %>%
  unpack_square_matrix(genes) %>%
  merge(total_counts[,-3],by = 'labels_1') %>%
  merge(total_counts[,-3],by.x = 'labels_2',by.y = 'labels_1') %>%
  mutate(
    intersection = value.x,
    union = value.y + value - value.x
  ) %>%
  select(labels_1,labels_2,intersection,union) %>%
  mutate(IoU = intersection / union) %>% 
  mutate(IoU = ifelse(labels_1 == labels_2,1,IoU)) %>% 
  filter(IoU > 0) %>% 
  ggplot(aes(x = labels_1,y = labels_2,fill = IoU,label = round(IoU,2))) +
  geom_tile() + 
  geom_text(aes(colour = IoU),size=3) +
  theme_gerstung(base_size = 15) + 
  rotate_x_text() + 
  scale_fill_material(palette = "deep-purple",
                      name = "Jaccard\nIndex",
                      na.value="white") + 
  scale_colour_material(palette = "purple",
                        na.value="white",
                        trans = "reverse",
                        guide=F) + 
  xlab("") + 
  ylab("") +
  theme(panel.grid = element_blank()) + 
  scale_y_discrete() + 
  ggtitle("Jaccard index for different gene pairs")
```

We also assess the distribution of mutations for single genes in individuals. We can define, for a given gene, an outlierness for an individual characterised by having more mutations than those expected from prevalence alone. Remarkably, we find that the genes with the most outliers are *TET2* and *DNMT3A*, each with 3 individuals having more mutations than expected.

```{r recurrence_prevalence,warning=F,message=F,fig.height=2,fig.width=3.3}
N_ind <- length(unique(full_data$SardID))
NMutIndividual <- full_data %>%
  group_by(SardID) %>%
  summarise(N = length(unique(amino_acid_change)))

NMutGeneIndividual <- full_data %>%
  mutate(Gene = str_match(Gene,"[0-9A-Z]+")) %>%
  group_by(SardID,Gene) %>% 
  summarise(N = length(unique(amino_acid_change)))

PrevalenceGenes <- full_data %>%
  mutate(Gene = str_match(Gene,"[0-9A-Z]+")) %>%
  group_by(Gene) %>% 
  summarise(N = length(unique(paste(amino_acid_change,SardID)))) %>%
  ungroup %>%
  mutate(P = N / sum(N))

simulations <- list()
n_sim <- 100
for (i in 1:n_sim) {
  simulations[[i]] <- NMutIndividual %>%
    apply(1,
      function(x) {
        S <- rmultinom(1,x[2],PrevalenceGenes$P) %>%
          t 
         S_idx <- which(S > 0)
        return(data.frame(Gene = PrevalenceGenes$Gene[S_idx],Count = S[S_idx]))
        }
    ) %>% 
    do.call(what = rbind) %>% 
    as.tibble() %>%
    group_by(Gene,Count) %>% 
    summarise(N = length(Gene)) %>% 
    group_by(Gene) %>%
    mutate(Total = sum(N)) %>% 
    spread(key = Count,value = N,fill = 0) %>%
    #mutate(`0` = N_ind - Total) %>%
    gather(key = "Count",value = "N",-Gene,-Total) %>%
    select(-Total)
}

simulation_df <- simulations %>% 
  do.call(what = rbind)

ecdf_list <- list()
for (x in PrevalenceGenes$Gene) {
  S <- simulation_df %>% 
    subset(Gene == x) 
  ecdf_list[[x]] <- S %>% 
    apply(1,function(x) rep(x[2],x[3])) %>% 
    unlist %>%
    ecdf
}

expected_N <- simulation_df %>%
  group_by(Gene,Count) %>%
  summarise(ExpectedTotal = round(sum(N)/n_sim)) %>%
  select(Gene,N = Count,ExpectedTotal) %>%
  mutate(N = as.numeric(N),
         Gene = as.character(Gene))

NTests <- length(unique(paste(NMutGeneIndividual$Gene,NMutGeneIndividual$N)))
NMutGene <- NMutGeneIndividual %>% 
  group_by(Gene,N) %>%
  summarise(Total = length(unique(SardID))) 

r_div <- 0.3
PlotData <- NMutGeneIndividual %>% 
  rowwise %>% 
  mutate(Prob = 1 - ecdf_list[[Gene]](N)) %>% 
  group_by(Gene,N,Prob) %>%
  summarise(Total = length(unique(SardID))) %>% 
  ungroup() %>% 
  mutate(Significant = ifelse(Prob < 0.05 / NTests,"Significant","Non-significant"),
         N = as.numeric(N)) %>% 
  merge(expected_N,by = c("Gene","N"),all.x = T) %>%
  mutate(Excess = Total - ExpectedTotal) %>%
  mutate(Missing = Excess < 0) %>%
  mutate(Missing = ifelse(is.na(Missing),F,Missing)) %>% 
  mutate(Proportion = Total/ExpectedTotal) %>%
  mutate(ProportionA = ifelse(Proportion > 1,1,Proportion)) %>%
  mutate(ProportionB = 1 - ProportionA) %>%
  mutate(ProportionExcess = Proportion - 1) %>% 
  mutate(ProportionExcess = ifelse(ProportionExcess < 0,0,ProportionExcess)) %>% 
  mutate(ProportionExcess = ifelse(is.infinite(ProportionExcess),1,ProportionExcess)) %>%
  mutate(ProportionExecssInv = 1 - ProportionExcess) %>%
  mutate(ProportionExcessExtra = ifelse(ProportionExcess > 1,ProportionExcess - 1,0)) %>%
  mutate(ProportionExcessExtraInv = 1 - ProportionExcessExtra) %>%
  mutate(Radius = r_div/2 + r_div*(Total-min(Total))/(max(Total)-min(Total))) %>%
  mutate(RadiusAB = ifelse(Radius <= r_div/2,NA,Radius))

GeneOrder <- PlotData %>% 
  group_by(Gene) %>% 
  summarise(Total = sum(Total),
            MaxN = max(N)) %>% 
  arrange(-MaxN,-Total) %>% 
  select(Gene) %>% 
  unlist

ggplot() + 
  geom_scatterpie(data = PlotData,aes(x = as.numeric(factor(Gene,levels = GeneOrder)),
                                      y = N,group = paste(Gene,N),r = RadiusAB),
                  cols = c("ProportionA","ProportionB"),colour = NA) + 
  geom_scatterpie(data = PlotData,aes(x = as.numeric(factor(Gene,levels = GeneOrder)),
                                      y = N,group = paste(Gene,N),r = Radius),
                  cols = c("ProportionExcess","ProportionExecssInv"),colour = NA) + 
  geom_scatterpie(data = PlotData,aes(x = as.numeric(factor(Gene,levels = GeneOrder)),
                                      y = N,group = paste(Gene,N),r = Radius),
                  cols = c("ProportionExcessExtra","ProportionExcessExtraInv"),colour = NA) + 
  geom_scatterpie_legend(PlotData$Radius,x = 10,y = 7,n = 2) + 
  coord_equal() +
  scale_fill_manual(values = c(ProportionA="grey80",ProportionB=NA,
                               ProportionExcess="grey40",ProportionExcessInv=NA,
                               ProportionExcessExtra="black",ProportionExecssExtraInv=NA),guide = F) + 
  theme_gerstung(base_size = 6) + 
  scale_x_continuous(labels = GeneOrder,breaks = c(1:17),expand = c(0.01,0)) + 
  scale_y_continuous(breaks = c(1:10),expand = c(0.01,0)) + 
  xlab("") + 
  ylab("Number of unique\nmutations per individual") + 
  rotate_x_text() + 
  theme(axis.text.x = element_text(face = "italic")) + 
  ggsave(useDingbats=FALSE,"figures/data_exploration/PrevalenceByGenePie.pdf",height = 1.5,width = 2) 
```

We can also confirm a few previously stated conclusions regarding clonal haematopoiesis - that its magnitude, prevalence, and clonality increase with age, with gene specific signatures being relevant for prevalence (namely for *U2AF1*).

```{r summary_figure_first,fig.height=2.5,fig.width=2}
pd <- position_dodge(0.5)
colour_palette <- colorRampPalette(c("grey80","orange"))
load_data_keep() %>%
  mutate(DetectableAt0005 = VAF > 0.005) %>%
  mutate(AgeGroup = cut(Age,c(0,60,65,70,75,80,85,150),
                        c("<60","60-65","65-70","70-75","75-80","80-85",">85"))) %>% 
  group_by(SardID,AgeGroup) %>%
  summarise(N = length(unique(AAChange.refGene[Gene %in% load_included_genes()]))) %>%
  mutate(N = ifelse(N > 5,5,N)) %>%
  group_by(AgeGroup,N) %>%
  summarise(count = length(unique(SardID))) %>%
  group_by(AgeGroup) %>%
  mutate(freq = count / sum(count)) %>%
  ggplot(aes(x = AgeGroup,y = freq,fill = as.factor(N))) + 
  geom_bar(stat="identity",
           colour = "black",
           size = 0.1) + 
  scale_fill_manual(breaks = c(0,1,2,3,4,5),
                    values = colour_palette(6),
                    labels = c(0,1,2,3,4,"5+")) +
  theme_gerstung(base_size = 6) +
  scale_y_continuous(expand = c(0,0,0,0)) + 
  guides(fill=guide_legend(title="No. of\nmutations",ncol=1)) +
  xlab("Age") + 
  ylab("Prevalence") + 
  theme(legend.margin = margin(),
        legend.box.margin = margin(),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(0.2,"cm"),
        legend.key.size = unit(0.3,"cm")) + 
  rotate_x_text() +
  ggsave(useDingbats=FALSE,"figures/data_exploration/prevalence_through_time.pdf",width=1.8, height=2)

full_data %>% 
  group_by(SardID,Phase,Age) %>%
  summarise(VAF = mean(VAF)) %>%
  mutate(AgeGroup = cut(Age,c(0,60,65,70,75,80,85,90,150),
                        c("<60","60-65","65-70","70-75","75-80","80-85","85-90",">90"))) %>% 
  ggplot(aes(x = AgeGroup,y = VAF)) + 
  geom_jitter(height = 0,width = 0.2,
              alpha = 0.2,size = 0.3) + 
  geom_boxplot(aes(group = AgeGroup),alpha = 0.4,
               size = 0.3,
               outlier.alpha = 0) + 
  theme_gerstung(base_size = 6) + 
  scale_y_continuous(trans = 'log10') + 
  ylab("Average VAF") + 
  xlab("Age") +
  rotate_x_text() +
  ggsave(useDingbats=FALSE,"figures/data_exploration/mean_vaf_through_time.pdf",width=1.8, height=2)

gene_order <- load_data_keep() %>%
  mutate(Gene = str_match(Gene,"[0-9A-Z]+")) %>%
  mutate(Total = length(unique(SardID))) %>%
  subset(Gene %in% load_included_genes()) %>%
  group_by(Gene,Total) %>%
  summarise(N = length(unique(SardID))) %>%
  arrange(-N) %>% 
  select(Gene) %>%
  unlist

prevalence_gene <- load_data_keep() %>%
  mutate(Gene = str_match(Gene,"[0-9A-Z]+")) %>%
  mutate(Total = length(unique(SardID))) %>%
  subset(Gene %in% load_included_genes()) %>%
  mutate(Gene = factor(Gene,levels = gene_order)) %>% 
  group_by(Gene,Total) %>%
  summarise(N = length(unique(SardID))) %>% 
  ggplot(aes(x = Gene,y = N/Total,fill = Gene)) + 
  geom_bar(stat = "identity") + 
  scale_y_continuous(expand = c(0,0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  scale_fill_manual(values = gene_colours,guide=F) +
  theme_gerstung(base_size = 6) + 
  rotate_x_text() + 
  theme(axis.text.x = element_text(face = "italic")) + 
  ylab("Prevalence") + 
  xlab("")

prevalence_age_gene <- load_data_keep() %>%
  mutate(Gene = str_match(Gene,"[0-9A-Z]+")) %>%
  subset(Gene %in% load_included_genes()) %>%
  mutate(Gene = factor(Gene,levels = gene_order)) %>% 
  group_by(Gene) %>% 
  mutate(AgeGroup = cut(Age,c(0,70,80,90,150),
                        c("<70","70-80","80-90",">90"),include.lowest = T)) %>% 
  group_by(Gene,AgeGroup) %>%
  summarise(N = length(unique(SardID[VAF > 0.005]))) %>% 
  group_by(Gene) %>%
  mutate(Total = sum(N)) %>%
  ggplot(aes(x = Gene,y = N/Total,fill = as.factor(AgeGroup))) + 
  geom_bar(stat = "identity",
           colour = "black",size = 0.1) + 
  scale_y_continuous(expand = c(0,0),breaks = c(0,0.5,1.0)) + 
  scale_x_discrete(expand = c(0,0)) + 
  theme_gerstung(base_size = 6) + 
  rotate_x_text() + 
  theme(axis.text.x = element_text(face = "italic")) + 
  scale_fill_manual(
    values = colorRampPalette(c("yellow","lightgreen","royalblue"))(4),
    name = "Age"
  ) + 
  ylab("Relative prevalence") + 
  xlab("") 

prevalence_pie_age_gene_order <- c(
  "TET2","DNMT3A",
  "SF3B1", "SRSF2","SRSF2-P95H","U2AF1",
  "PPM1D","TP53",
  "ASXL1","CBL","CTCF",
  "BRCC3","JAK2","PTPN11",
  "GNB1","KRAS",
  "IDH1","IDH2")
prevalence_pie_age <- full_data %>% 
  mutate(Gene = str_match(Gene,"[A-Z0-9]+")) %>%
  mutate(Gene = ifelse(amino_acid_change == "SRSF2-P95H",
                       "SRSF2-P95H",
                       Gene)) %>% 
  group_by(SardID,amino_acid_change,Gene) %>% 
  summarise(Age = min(Age)) %>% 
  mutate(AgeGroup = cut(Age,c(0,60,80,200),c("<60","60-80",">80"))) %>% 
  group_by(AgeGroup,Gene) %>% 
  summarise(N = length(unique(SardID))) %>% 
  group_by(AgeGroup) %>% 
  mutate(Total = sum(N)) %>% 
  mutate(Gene = factor(Gene,levels = prevalence_pie_age_gene_order)) %>% 
  ggplot(aes(x = factor(1),y = N/Total,fill = Gene)) + 
  geom_col(size = 0.25,colour = "black") + 
  coord_polar("y") + 
  facet_wrap(~ AgeGroup) + 
  scale_fill_manual(values = c(gene_colours,`SRSF2-P95H` = "red"),name = NULL) + 
  theme_gerstung(base_size = 6) + 
  theme(axis.line = element_blank(),axis.text = element_blank(),
        legend.position = "bottom",axis.ticks = element_blank(),
        legend.text = element_text(face = "italic"),
        legend.key.size = unit(0.2,"cm")) + 
  xlab("") + 
  ylab("") + 
  guides(fill = guide_legend(nrow = 3,byrow = T))

plot_grid(prevalence_gene + theme(panel.background = element_blank()),
          prevalence_age_gene + theme(legend.position = "none"),
          get_legend(prevalence_age_gene + 
                       guides(fill = guide_legend(ncol = 2,byrow = F)) + 
                       theme(legend.margin = margin(),
                             legend.box.margin = margin(),
                             legend.key.size = unit(0.15,"cm"))),
          rel_heights = c(2,2,0.6),
          ncol = 1) + 
  ggsave(useDingbats=FALSE,"figures/data_exploration/prevalence_relative_prevalence.pdf",
         width = 1.7,height = 3)
```

```{r,fig.width=1.8,fig.height=1}
bootstrap_ci <- function(m,n,p=0.5,repeats=100) {
  hits <- c(rep(1,m),rep(0,n-m))
  trials <- rep(1,n)
  sapply(
    seq(1,repeats),
    function(x) {
      s <- sample(n,round(p*n),replace = F)
      return(sum(hits[s]))
    }
  )
}

tmp_data <- load_data_keep() %>%
  subset(VAF >= 0.005) %>% 
  group_by(SardID) %>%
  filter(Age == max(Age)) %>%
  mutate(AgeGroup = cut(Age,c(0,75,80,85,200),c("<75","75-80","80-85",">85"))) %>% 
  group_by(SardID,AgeGroup,Age,Gender) %>%
  summarise(TET2 = sum(grepl("TET2",unique(AAChange.refGene))),
            DNMT3A = sum(grepl("DNMT3A",unique(AAChange.refGene))),
            Splicing = sum(grepl("SF3B1|U2AF1|SRSF2",unique(AAChange.refGene))),
            TotalClones = length(unique(AAChange.refGene)),
            Total = 1)

tmp_data_splicing <- load_data_keep() %>%
  subset(VAF >= 0.005) %>% 
  group_by(SardID) %>%
  filter(Age == max(Age)) %>%
  mutate(AgeGroup = cut(Age,c(0,75,80,85,200),c("<75","75-80","80-85",">85"))) %>% 
  group_by(SardID,AgeGroup,Age,Gender) %>%
  summarise(SF3B1 = sum(grepl("SF3B1",unique(AAChange.refGene))),
            U2AF1 = sum(grepl("U2AF1",unique(AAChange.refGene))),
            SRSF2 = sum(grepl("SRSF2",unique(AAChange.refGene))),
            TotalClones = length(unique(AAChange.refGene)),
            Total = 1)

tmp_data_regression <- load_data_keep() %>%
  subset(VAF >= 0.005) %>% 
  group_by(SardID) %>% 
  mutate(Age = max(Age)) %>% 
  group_by(SardID,Age,Gender) %>%
  summarise(TET2 = any(grepl("TET2",unique(AAChange.refGene))),
            DNMT3A = any(grepl("DNMT3A",unique(AAChange.refGene))),
            Splicing = any(grepl("U2AF1|SRSF2|SF3B1",unique(AAChange.refGene))))

tmp_data_per_individual <- load_data_keep() %>%
  subset(VAF >= 0.005) %>% 
  group_by(SardID) %>% 
  mutate(Age = max(Age)) %>% 
  group_by(SardID,Age,Gender) %>%
  summarise(TET2 = sum(grepl("TET2",unique(AAChange.refGene))),
            DNMT3A = sum(grepl("DNMT3A",unique(AAChange.refGene))),
            Splicing = sum(grepl("U2AF1|SRSF2|SF3B1",unique(AAChange.refGene)))) %>%
  mutate(AgeGroup = cut(Age,c(0,75,80,85,200),c("<75","75-80","80-85",">85"))) %>%
  mutate(AgeGroup = factor(AgeGroup,levels = c("<75","75-80","80-85",">85"))) %>%
  group_by(AgeGroup) %>% 
  mutate(AgeGroup_ = sprintf("%s\n(N=%s)",AgeGroup,length(unique(SardID)))) %>%
  gather(key = "key",value = "value",-Age,-AgeGroup,-AgeGroup_,-SardID,-Gender) %>%
  mutate(key = factor(key, levels = c("DNMT3A","TET2","Splicing"),
                      labels = c("DNMT3A","TET2","Splicing"))) 

tmp_data_agg <- tmp_data %>% 
  group_by(AgeGroup) %>% 
  summarise(TotalTET2 = sum(TET2),
            TotalDNMT3A = sum(DNMT3A),
            TotalSplicing = sum(Splicing),
            TotalClones = sum(TotalClones),
            Total = length(unique(paste(SardID,Age)))) %>%
  gather(key = "key",value = "value",-AgeGroup,-Total,-TotalClones) %>%
  rowwise() %>%
  mutate(p05 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.05),
         p95 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.95)) %>% 
  mutate(key = factor(key, levels = c("TotalDNMT3A","TotalTET2","TotalSplicing"),
                      labels = c("DNMT3A","TET2","Splicing"))) %>% 
  mutate(AgeGroup = factor(AgeGroup,levels = c("<75","75-80","80-85",">85"))) %>%
  mutate(AgeGroup_ = sprintf("%s\n(N=%s)",AgeGroup,Total)) 

tmp_data_splicing_agg <- tmp_data_splicing %>% 
  group_by(AgeGroup) %>% 
  summarise(TotalSF3B1 = sum(SF3B1),
            TotalSRSF2 = sum(SRSF2),
            TotalU2AF1 = sum(U2AF1),
            TotalClones = sum(TotalClones),
            Total = length(unique(paste(SardID,Age)))) %>%
  gather(key = "key",value = "value",-AgeGroup,-Total,-TotalClones) %>%
  rowwise() %>%
  mutate(p05 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.05),
         p95 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.95)) %>% 
  mutate(key = factor(key, levels = c("TotalSF3B1","TotalSRSF2","TotalU2AF1"),
                      labels = c("SF3B1","SRSF2","U2AF1"))) %>% 
  mutate(AgeGroup = factor(AgeGroup,levels = c("<75","75-80","80-85",">85"))) %>%
  mutate(AgeGroup_ = sprintf("%s\n(N=%s)",AgeGroup,Total)) 

lvls <- tmp_data_agg %>%
  select(AgeGroup,Total) %>%
  distinct %>%
  arrange(AgeGroup) %>%
  ungroup %>%
  transmute(AgeGroup_ = sprintf("%s\n(N=%s)",AgeGroup,Total)) %>%
  select(AgeGroup_) %>%
  unlist

mm <- tmp_data_agg$p95 %>%
  max

tmp_data_agg %>% 
  ggplot(aes(x = factor(AgeGroup_,levels = lvls),y = value/Total,fill = key,group = key)) + 
  geom_bar(stat = "identity",position = position_dodge(width = 0.9)) + 
  geom_linerange(aes(ymin = p05,ymax = p95),position = position_dodge(width = 0.9),
                 size = 0.25) +
  geom_point(data = tmp_data_per_individual,aes(y = value/10),size = 0.1,
             alpha = 0.5,
             position = position_jitterdodge(dodge.width = 0.9,jitter.width = 0.2,
                                             jitter.height = 0.02)) +
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = c(gene_colours,Splicing="grey50"),name = NULL,guide = F) + 
  theme(axis.title.x = element_blank(),
        legend.key.size = unit(0.2,"cm"),
        legend.position = "bottom",
        legend.box.spacing = unit(0,"cm")) + 
  scale_y_continuous(expand = c(0,0),limits = c(0,mm),
                     sec.axis = sec_axis(trans = ~ .*10,name = "Mutations per individual (.)")) + 
  ylab("Average number of\nclones per individual") + 
  ggsave(useDingbats=FALSE,sprintf("figures/data_exploration/dnmt3a_tet2_splicing_prevalence.pdf",model_id),
         height = 1,width = 1.8)

tmp_data_splicing_agg %>% 
  ggplot(aes(x = factor(AgeGroup_,levels = lvls),y = value/Total,fill = key,group = key)) + 
  geom_bar(stat = "identity") + 
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = c(gene_colours,Splicing="grey50"),name = NULL,guide = F) + 
  theme(axis.title.x = element_blank(),
        legend.key.size = unit(0.2,"cm"),
        legend.position = "bottom",
        legend.box.spacing = unit(0,"cm")) + 
  scale_y_continuous(expand = c(0,0),limits = c(0,mm)) + 
  ylab("Average number of\nclones per individual") + 
  ggsave(useDingbats=FALSE,sprintf("figures/data_exploration/splicing_prevalence.pdf",model_id),
         height = 1,width = 1.8)

glm(TET2 ~ Age + Gender,data = tmp_data_regression,family = stats::binomial()) %>%
  summary

glm(DNMT3A ~ Age + Gender,data = tmp_data_regression,family = stats::binomial()) %>%
  summary

glm(Splicing ~ Age + Gender,data = tmp_data_regression,family = stats::binomial()) %>%
  summary

rbind(tmp_data_agg,tmp_data_splicing_agg) %>%
  select(-AgeGroup_,TotalIndividuals = Total) %>%
  write.csv("data_output/figure_1d.csv")
```

Same as before but considering only *TET2* and *DNMT3A*.

```{r,fig.height=2,fig.width=2}
load_data_keep() %>%
    group_by(SardID) %>%
    filter(Age == max(Age)) %>%
    mutate(AgeGroup = cut(Age,c(0,75,80,85,200),c("<75","75-80","80-85",">85"))) %>% 
    group_by(SardID,AgeGroup,Age,Gender) %>%
    summarise(TET2 = sum(ifelse(VAF >= 0.005,grepl("TET2",unique(AAChange.refGene)),F)),
              DNMT3A = sum(ifelse(VAF >= 0.005,grepl("DNMT3A",unique(AAChange.refGene)),F)),
              Splicing = sum(ifelse(VAF >= 0.005,grepl("SF3B1|U2AF1|SRSF2",unique(AAChange.refGene)),F)),
              TotalClones = length(unique(AAChange.refGene)),
              Total = 1) %>%
    group_by(AgeGroup) %>% 
    summarise(TotalTET2 = sum(TET2,na.rm=T),
              TotalDNMT3A = sum(DNMT3A,na.rm=T),
              TotalSplicing = sum(Splicing,na.rm=T),
              TotalClones = sum(TotalClones),
              Total = length(unique(paste(SardID,Age)))) %>%
    gather(key = "key",value = "value",-AgeGroup,-Total,-TotalClones) %>%
    rowwise() %>%
    mutate(p05 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.05),
           p95 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.95)) %>% 
    mutate(key = factor(key, levels = c("TotalDNMT3A","TotalTET2","TotalSplicing"),
                        labels = c("DNMT3A","TET2","Splicing"))) %>% 
    mutate(AgeGroup = factor(AgeGroup,levels = c("<75","75-80","80-85",">85"))) %>%
    mutate(AgeGroup_ = sprintf("%s\n(N=%s)",AgeGroup,Total)) %>% 
    subset(key != 'Splicing') %>%
    group_by(AgeGroup) %>% mutate(TT = sum(value)) %>%
    ggplot(aes(x = "", y = value/TT,fill = key)) + 
  geom_bar(stat="identity") + 
  facet_wrap(~ AgeGroup) + 
  coord_polar("y", start=0) + 
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = gene_colours,name = NULL) +
  theme(axis.text = element_blank(),axis.ticks = element_blank(),
        axis.line = element_blank(),axis.title = element_blank(),
        legend.key.size = unit(0.2,"cm"),legend.text = element_text(face = "italic")) +
  ggtitle("All clones")

load_data_keep() %>%
    group_by(SardID) %>%
    filter(Age == max(Age)) %>%
    mutate(AgeGroup = cut(Age,c(0,75,80,85,200),c("<75","75-80","80-85",">85"))) %>% 
    group_by(SardID,AgeGroup,Age,Gender) %>%
    summarise(TET2 = sum(ifelse(VAF >= 0.025,grepl("TET2",unique(AAChange.refGene)),F)),
              DNMT3A = sum(ifelse(VAF >= 0.025,grepl("DNMT3A",unique(AAChange.refGene)),F)),
              Splicing = sum(ifelse(VAF >= 0.025,grepl("SF3B1|U2AF1|SRSF2",unique(AAChange.refGene)),F)),
              TotalClones = length(unique(AAChange.refGene)),
              Total = 1) %>%
    group_by(AgeGroup) %>% 
    summarise(TotalTET2 = sum(TET2,na.rm=T),
              TotalDNMT3A = sum(DNMT3A,na.rm=T),
              TotalSplicing = sum(Splicing,na.rm=T),
              TotalClones = sum(TotalClones),
              Total = length(unique(paste(SardID,Age)))) %>%
    gather(key = "key",value = "value",-AgeGroup,-Total,-TotalClones) %>%
    rowwise() %>%
    mutate(p05 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.05),
           p95 = quantile(bootstrap_ci(m = value,n = TotalClones,p = 0.8)/(Total*0.8),0.95)) %>% 
    mutate(key = factor(key, levels = c("TotalDNMT3A","TotalTET2","TotalSplicing"),
                        labels = c("DNMT3A","TET2","Splicing"))) %>% 
    mutate(AgeGroup = factor(AgeGroup,levels = c("<75","75-80","80-85",">85"))) %>%
    mutate(AgeGroup_ = sprintf("%s\n(N=%s)",AgeGroup,Total)) %>% 
    subset(key != 'Splicing') %>%
    group_by(AgeGroup) %>% mutate(TT = sum(value)) %>%
    ggplot(aes(x = "", y = value/TT,fill = key)) + 
  geom_bar(stat="identity") + 
  facet_wrap(~ AgeGroup) + 
  coord_polar("y", start=0) + 
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = gene_colours,name = NULL) +
  theme(axis.text = element_blank(),axis.ticks = element_blank(),
        axis.line = element_blank(),axis.title = element_blank(),
        legend.key.size = unit(0.2,"cm"),legend.text = element_text(face = "italic")) +
  ggtitle("Large clones (VAF > 2.5%)")
```

# Modelling 

We model the counts of mutation $i$ in individual $j$ - clone $c_{ji}$ - at age $t$ as the product of a beta-binomial distribution such that $counts_{c_{ji}} \sim BB(cov_{c_{ji}}(t),\alpha_{c_{ji}}(t),\beta)$, where $cov_{c_{ji}}$ is the coverage for mutation $c_{ji}$ at timepoint $t$, $\beta$ is the technical overdispersion and $\alpha = \frac{\beta p}{1 - p(t)}$, where $p(t)$ is the expected value for the VAF at age $t$. 

Here we work with a relatively simple formulation of clone growth - for a given clone $i$ in individual $j$ we consider there to be a *genetic* growth effect ($b_{genetic_{i}} = b_{gene_i} + b_{site_i}$), an *unknown cause* growth effect ($b_{c_{ji}}$) and a *clone-specific* offset ($u_{ji}$). Only the *growth effects* mentioned here depend on age. These two growth effects and offset are linearly combined and inverse-logit transformed to yield a probability such that $p_{ji}(t) = \frac{ilogit((b_{genetic_{i}} + b_{c_{ji}}) * t + u_{ji})}{2}$. Summarily:

$$
b_{gene_i} \sim N(0,0.1)
\\
b_{site_i} \sim N(0,0.1)
\\
u_{ji} \sim Uniform(-50,0)
\\
b_{c_{ji}} \sim N(0,0.05)
\\
p_i(t) = \frac{ilogit((b_{gene_i} + b_{site_i} + b_{clone_{ji}}) * t + u_{ji})}{2}
\\
\beta \sim N(\mu=\mu_{\beta},\sigma=\sigma_{\beta})
\\
\alpha_{ji}(t) = \frac{\beta p_i}{1 - p_{ji}}
\\
counts \sim BB(coverage,\alpha_{ji}(t),\beta)
$$

## Modelling clone growth as a product of genetic and unknown cause effects

The model specified in `Scripts/bb_gene_site_clone_model.R` was executed using `Scripts/bb_gene_site_clone_model.R`. The model parameters are stored in `models/model_ch.RDS` and will be used for the remainder of this section.

```{r model_load_values}
values_model <- readRDS(model_file_name)
dir.create("figures/",showWarnings = F)
dir.create(sprintf("figures/%s",model_id),showWarnings = F)
c('age_at_clone_foundation','dynamic_coefficients',
  'inferred_trajectories','statistical_analysis') %>%
  sapply( 
    function(x) dir.create(sprintf("figures/%s/%s",model_id,x),showWarnings = F)
  )
```

```{r model_visualization_preparation}
unique_gene <- values_model$sub_data %>%
  select(Gene,gene_numeric) %>% 
  distinct %>%
  arrange(gene_numeric) %>%
  select(Gene) %>%
  unlist
unique_site <- values_model$sub_data %>%
  select(amino_acid_change,site_numeric) %>% 
  distinct %>%
  arrange(site_numeric) %>%
  select(amino_acid_change) %>%
  unlist
unique_site_multiple <- values_model$sub_data %>%
  select(amino_acid_change,site_numeric_multiple) %>% 
  subset(!is.na(site_numeric_multiple)) %>%
  distinct %>%
  arrange(site_numeric_multiple) %>%
  select(amino_acid_change) %>%
  unlist %>% 
  as.character() %>% as.factor

site_samples <- values_model$draws$`11`[,grep('site',colnames(values_model$draws$`11`))]
gene_samples <- values_model$draws$`11`[,grep('gene',colnames(values_model$draws$`11`))]
clone_samples <- values_model$draws$`11`[,grep('clone',colnames(values_model$draws$`11`))]
offset_samples <- values_model$draws$`11`[,grep('u',colnames(values_model$draws$`11`))]
beta_samples <- values_model$draws$`11`[,grep('beta',colnames(values_model$draws$`11`))] %>%
  unlist()

grab_samples <- function(site_numeric,gene_numeric,clone_numeric,sample_size=1000) {
  S <- sample(nrow(site_samples),sample_size,replace = F)
  if (!is.na(site_numeric)) {
    b_site <- site_samples[S,site_numeric]
  } else {
    b_site <- 0
  }
  b_gene <- gene_samples[S,gene_numeric]
  b_clone <- clone_samples[S,clone_numeric]
  u <- offset_samples[S,clone_numeric]
  return(list(site = b_site,clone = b_clone,gene = b_gene,offset = u))
}

X_val <- values_model$sub_data$MUTcount_Xadj

b_gene_values_val <- values_model$b_gene_values[[1]]$values[values_model$b_gene_values[[1]]$labels == 'mean'] %>% 
  matrix(ncol=1) 
b_site_values_val <- values_model$b_site_values[[1]]$values[values_model$b_site_values[[1]]$labels == 'mean'] %>% 
  matrix(ncol=1) 
beta_val <- values_model$beta_values[[1]]$values[values_model$beta_values[[1]]$labels == 'mean'] %>% 
  matrix(ncol=1) 
beta_val_var <- values_model$beta_values[[1]]$values[values_model$beta_values[[1]]$labels == 'var'] %>% 
  matrix(ncol=1)
b_clone_values_val <- values_model$b_clone[[1]]$values[values_model$b_clone[[1]]$labels == 'mean'] %>% 
  matrix(ncol=1)

u_values_val <- values_model$u_values[[1]]$values[values_model$u_values[[1]]$labels == 'mean'] %>%
  matrix(nrow=1) 

ae_gene <- b_gene_values_val[values_model$sub_data$gene_numeric]
ae_site <- b_site_values_val[values_model$sub_data$site_numeric_multiple]
ae_site <- ifelse(is.na(ae_site),0,ae_site)

age_effect_coef <- ae_gene + ae_site
age_effect <- age_effect_coef + b_clone_values_val[values_model$sub_data$clone]
age_effect <- age_effect * (values_model$sub_data$Age - values_model$min_age)

offset_per_individual <- u_values_val[values_model$sub_data$clone]
r <- age_effect + offset_per_individual
mu_val <- inv.logit(r) * 0.5

r_values_full <- data.frame(
  mu_val = mu_val,
  true = values_model$sub_data$MUTcount_Xadj,
  age = values_model$sub_data$Age,
  coverage = values_model$sub_data$TOTALcount,
  individual = values_model$sub_data$SardID,
  site = values_model$sub_data$amino_acid_change,
  site_numeric = values_model$sub_data$site_numeric_multiple,
  gene_numeric = values_model$sub_data$gene_numeric,
  clone_numeric = values_model$sub_data$clone,
  b_gene_005 = c(values_model$b_gene_values[[1]]$values[
    values_model$b_gene_values[[1]]$labels == 'HDPI_low'][values_model$sub_data$gene_numeric]),
  b_gene = c(values_model$b_gene_values[[1]]$values[
    values_model$b_gene_values[[1]]$labels == 'mean'][values_model$sub_data$gene_numeric]),
  b_gene_095 = c(values_model$b_gene_values[[1]]$values[
    values_model$b_gene_values[[1]]$labels == 'HDPI_high'][values_model$sub_data$gene_numeric]),
  b_gene_var = c(values_model$b_gene_values[[1]]$values[
    values_model$b_gene_values[[1]]$labels == 'var'][values_model$sub_data$gene_numeric]),
  b_site_005 = c(values_model$b_site_values[[1]]$values[
    values_model$b_site_values[[1]]$labels == 'HDPI_low'][values_model$sub_data$site_numeric_multiple]),
  b_site = c(values_model$b_site_values[[1]]$values[
    values_model$b_site_values[[1]]$labels == 'mean'][values_model$sub_data$site_numeric_multiple]),
  b_site_095 = c(values_model$b_site_values[[1]]$values[
    values_model$b_site_values[[1]]$labels == 'HDPI_high'][values_model$sub_data$site_numeric_multiple]),
  b_site_var = c(values_model$b_site_values[[1]]$values[
    values_model$b_site_values[[1]]$labels == 'var'][values_model$sub_data$site_numeric_multiple]),
  b_clone_005 = c(values_model$b_clone[[1]]$values[
    values_model$b_clone[[1]]$labels == 'HDPI_low'][values_model$sub_data$clone]),
  b_clone = c(values_model$b_clone[[1]]$values[
    values_model$b_clone[[1]]$labels == 'mean'][values_model$sub_data$clone]),
  b_clone_095 = c(values_model$b_clone[[1]]$values[
    values_model$b_clone[[1]]$labels == 'HDPI_high'][values_model$sub_data$clone]),
  b_clone_var = c(values_model$b_clone[[1]]$values[
    values_model$b_clone[[1]]$labels == 'var'][values_model$sub_data$clone]),
  u_values_005 = c(values_model$u_values[[1]]$values[
    values_model$u_values[[1]]$labels == 'HDPI_low'][values_model$sub_data$clone]),
  u_values = c(values_model$u_values[[1]]$values[
    values_model$u_values[[1]]$labels == 'mean'][values_model$sub_data$clone]),
  u_values_095 = c(values_model$u_values[[1]]$values[
    values_model$u_values[[1]]$labels == 'HDPI_high'][values_model$sub_data$clone])
) %>%
  mutate(
    b_site_005 = ifelse(is.na(b_site_005),0,b_site_005),
    b_site = ifelse(is.na(b_site),0,b_site),
    b_site_095 = ifelse(is.na(b_site_095),0,b_site_095),
    b_site_var = ifelse(is.na(b_site_var),0,b_site_var)
  ) %>%
  mutate(
    coefficient_005 = b_gene_005 + b_site_005,
    coefficient = b_gene + b_site,
    coefficient_095 = b_gene_095 + b_site_095,
    coefficient_var = b_gene_var + b_site_var
  ) %>% 
  as.data.frame() %>% 
  apply(1,FUN = function(x){
    cov <- as.numeric(x[4])
    site_numeric <- as.numeric(x[7])
    gene_numeric <- as.numeric(x[8])
    clone_numeric <- as.numeric(x[9])
    age <- as.numeric(x[3])
    
    S <- grab_samples(site_numeric,gene_numeric,clone_numeric,2500)
    if (!is.na(site_numeric)) {
      b_site <- site_samples[,site_numeric]
    } else {
      b_site <- 0
    }
    b_gene <- gene_samples[,gene_numeric]
    b_clone <- clone_samples[,clone_numeric]
    u <- offset_samples[,clone_numeric]
    b <- b_gene + b_site + b_clone
    mu_val <- inv.logit(unlist((b) * (age - values_model$min_age) + u)) / 2
    mu_val_genetic <- inv.logit(
      unlist((b_gene + b_site) * (age - values_model$min_age) + u)) / 2
    alpha_value <- (beta_samples * mu_val) / (1 - mu_val)
    Q <- extraDistr::rbbinom(n = 2500,
                             size = cov,
                             alpha = alpha_value,
                             beta = beta_samples) %>%
      quantile(c(0.05,0.50,0.95))
    names(Q) <- c("pred_005","pred","pred_095")
    Q_mu_val <- quantile(mu_val,c(0.05,0.50,0.95))
    names(Q_mu_val) <- c("p_005","p_050","p_095")
    Q_mu_val_genetic <- quantile(mu_val_genetic,c(0.05,0.50,0.95))
    names(Q_mu_val_genetic) <- c("p_genetic_005","p_genetic_050","p_genetic_095")
    b_genetic <- b_gene + b_site
    Q_genetic <- quantile(b_genetic,c(0.05,0.5,0.95))
    names(Q_genetic) <- sprintf("b_genetic_%s",c("005","050","095"))
    return(c(x,Q,Q_mu_val,Q_mu_val_genetic,
             b_genetic_mean = mean(b_genetic),
             b_genetic_var = var(b_genetic),
             Q_genetic))
  }) %>% 
  t %>%
  as.data.frame() %>% 
  mutate(genes = sapply(as.character(site),function(x) unlist(strsplit(x,'-'))[[1]])) %>%
  group_by(site) %>%
  mutate(n = length(pred)) %>%
  ungroup() %>%
  mutate_if(function(x) !is.na(sum(as.numeric(as.character(x)))),
            function(x) as.numeric(as.character(x)))

r_values_ <- r_values_full %>% 
  mutate(VAFpred = pred / coverage,
         VAFtrue = true / coverage) %>% 
  subset(coverage > 0) %>%
  group_by(individual,site) %>% 
  mutate(normalizedVAFpred = (VAFpred + 1) / (VAFpred[which(age == min(age))] + 1),
         normalizedVAFtrue = (VAFtrue + 1) / (VAFtrue[which(age == min(age))] + 1)) %>%
  ungroup() %>%
  mutate(individual = as.character(individual)) %>% 
  gather(key = 'key',value = 'value',VAFpred,VAFtrue) %>%
  mutate(
    pred_005 = ifelse(key == 'VAFpred',
                      pred_005 / coverage,
                      NA),
    pred_095 = ifelse(key == 'VAFpred',
                      pred_095 / coverage,
                      NA)
  ) 
```

A simple analysis of the coefficients reveals to us that 68.5% of the variation that our model can explain is due to driver genetics. Looking at specific mutations, we have quite a varied picture in these regards - more than 90% of the variation of CBL and U2AF1 trajectories can be explained by genetics.

```{r variation_explained_plot,fig.width = 1.5,fig.height = 2.5}
explained_variance_total <- r_values_full %>%
  summarise(Genetic = rowSums(cov(data.frame(gene=coefficient*(age-values_model$min_age),
                                             clone=b_clone*(age-values_model$min_age)))) %>% 
              (function(x) x[1] / sum(x))) %>%
  unlist

variation_explained_plot <- r_values_full %>% 
  mutate(genes = str_match(genes,'[A-Z0-9]+')) %>% 
  group_by(genes) %>% 
  summarise(Genetic = rowSums(cov(data.frame(gene=coefficient*(age-values_model$min_age),
                                             clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[1]),
            Unknown = rowSums(cov(data.frame(gene=coefficient*(age-values_model$min_age),
                                             clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[2])
            ) %>% 
  mutate(Genetic = Genetic / (Genetic + Unknown),
         Unknown = Unknown / (Genetic + Unknown)) %>% 
  gather("key" = key,"value",-genes) %>% 
  mutate(model_name = paste0("Driver-associated\ninter-individual\nvariation = ",100*round(explained_variance_total,3),'%')) %>%
  mutate(key = factor(key,levels = rev(unique(key)))) %>% 
  group_by(genes) %>% 
  mutate(G = value[key == "Genetic"]) %>% 
  ggplot(aes(x = reorder(genes,G),y = value,fill = key)) + 
  geom_bar(stat = "identity",color = "black",size = 0) + 
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = c("#79AF97FF","#374E55FF"),name = NULL,breaks = c("Genetic","Unknown")) + 
  #rotate_x_text() + 
  xlab("Driver gene") + 
  ylab("") + 
  facet_grid(~ model_name) + 
  theme(legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.margin=margin(0,0,0,0),
        legend.key.size = unit(0.1,"cm"),
        axis.text.y = element_text(face = 'italic'),
        legend.box.margin=margin(0,0,0,0)) +
  scale_x_discrete(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0,0.01,0)) +
  coord_flip() + 
  guides(fill=guide_legend(nrow=1,byrow=TRUE))

ggsave(useDingbats=FALSE,filename = sprintf("figures/%s/statistical_analysis/variation_explained.pdf",model_id),
       variation_explained_plot,height = 2.5,width = 1.5)
variation_explained_plot
```

```{r,fig.height=4,fig.width=2}
variation_explained_plot <- r_values_full %>% 
  subset(b_site != 0) %>% 
  mutate(genes = str_match(genes,"[A-Z0-9]+"),
         site = str_match(site,"[A-Z0-9]+$")) %>%
  group_by(site,genes) %>% 
  summarise(Gene = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                          site=b_site*(age-values_model$min_age),
                                          clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[1]),
            Site = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                          site=b_site*(age-values_model$min_age),
                                          clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[2]),
            Unknown = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                             site=b_site*(age-values_model$min_age),
                                             clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[3])) %>% 
  # mutate(Gene = Gene / (Gene + Site + Unknown),
  #        Site = Site / (Gene + Site + Unknown),
  #        Unknown = Unknown / (Gene + Site + Unknown)) %>%
  gather("key" = key,"value",-site,-genes) %>% 
  mutate(model_name = paste0("Driver-associated\ninter-individual\nvariation = ",100*round(explained_variance_total,3),'%')) %>%
  mutate(key = factor(key,levels = rev(unique(key)))) %>% 
  group_by(site) %>% 
  mutate(G = sum(value)) %>% 
  group_by(genes) %>% 
  mutate(GG = max(G)) %>% 
  ggplot(aes(x = reorder(site,G),y = value,fill = key)) + 
  geom_bar(stat = "identity",color = "black",size = 0) + 
  geom_hline(yintercept = 0,size = 0.25) +
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = rev(pal_jama()(3)),name = NULL,breaks = c("Gene","Site","Unknown")) + 
  #rotate_x_text() + 
  xlab("Driver site") + 
  ylab("Variance explained") + 
  facet_grid(reorder(genes,-GG) ~ .,scales = "free_y",space = "free") + 
  scale_x_discrete(expand = c(0,0)) +
  coord_flip() + 
  theme(legend.position = "bottom",
        plot.title = element_text(hjust = 0.5),
        legend.margin=margin(0,0,0,0),
        legend.key.size = unit(0.1,"cm"),
        legend.box.margin=margin(0,0,0,0),
        strip.text.y = element_text(angle = 0,face = "italic")) +
  guides(fill=guide_legend(nrow=1,byrow=TRUE))

ggsave(useDingbats=FALSE,filename = sprintf("figures/%s/statistical_analysis/variation_explained_site.pdf",model_id),
       variation_explained_plot,height = 4,width = 2)
variation_explained_plot
```

```{r,fig.height=2,fig.width=2}
var_exp_agg_plot <- r_values_full %>% 
  mutate(genes = str_match(genes,"[A-Z0-9]+"),
         site = str_match(site,"[A-Z0-9]+$")) %>%
  summarise(Gene = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                          site=b_site*(age-values_model$min_age),
                                          clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[1]),
            Site = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                          site=b_site*(age-values_model$min_age),
                                          clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[2]),
            Unknown = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                             site=b_site*(age-values_model$min_age),
                                             clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[3])) %>% 
  # mutate(Gene = Gene / (Gene + Site + Unknown),
  #        Site = Site / (Gene + Site + Unknown),
  #        Unknown = Unknown / (Gene + Site + Unknown)) %>%
  gather("key" = key,"value") %>% 
  mutate(key = factor(key,levels = rev(unique(key)))) %>% 
  mutate(ID = ifelse(key == "Unknown","Unknown","Genetic")) %>% 
  ggplot(aes(x = ID,y = value,fill = key)) + 
  geom_bar(stat = "identity",color = "black",size = 0) + 
  geom_hline(yintercept = 0,size = 0.25) +
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = rev(pal_jama()(3)),name = NULL,breaks = c("Gene","Site","Unknown")) + 
  #rotate_x_text() + 
  xlab("") + 
  ylab("Variance explained") + 
  scale_x_discrete(expand = c(0,0)) +
  coord_flip() + 
  scale_y_continuous(expand = c(0,0)) +
  theme(legend.position = "bottom",
        legend.margin=margin(0,0,0,0),
        legend.key.size = unit(0.1,"cm"),
        legend.box.margin=margin(0,0,0,0)) +
  guides(fill=guide_legend(nrow=1,byrow=TRUE))

var_exp_agg_rec_plot <- r_values_full %>% 
  subset(b_site != 0) %>% 
  mutate(genes = str_match(genes,"[A-Z0-9]+"),
         site = str_match(site,"[A-Z0-9]+$")) %>%
  summarise(Gene = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                          site=b_site*(age-values_model$min_age),
                                          clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[1]),
            Site = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                          site=b_site*(age-values_model$min_age),
                                          clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[2]),
            Unknown = rowSums(cov(data.frame(gene=b_gene*(age-values_model$min_age),
                                             site=b_site*(age-values_model$min_age),
                                             clone=b_clone*(age-values_model$min_age)))) %>% (function(x) x[3])) %>% 
  # mutate(Gene = Gene / (Gene + Site + Unknown),
  #        Site = Site / (Gene + Site + Unknown),
  #        Unknown = Unknown / (Gene + Site + Unknown)) %>%
  gather("key" = key,"value") %>% 
  mutate(key = factor(key,levels = rev(unique(key)))) %>% 
  mutate(ID = ifelse(key == "Unknown","Unknown","Genetic")) %>% 
  ggplot(aes(x = ID,y = value,fill = key)) + 
  geom_bar(stat = "identity",color = "black",size = 0) + 
  geom_hline(yintercept = 0,size = 0.25) +
  theme_gerstung(base_size = 6) + 
  scale_fill_manual(values = rev(pal_jama()(3)),name = NULL,breaks = c("Gene","Site","Unknown")) + 
  #rotate_x_text() + 
  xlab("") + 
  ylab("Variance explained") + 
  scale_x_discrete(expand = c(0,0)) +
  coord_flip() + 
  scale_y_continuous(expand = c(0,0)) +
  theme(legend.position = "bottom",
        legend.margin=margin(0,0,0,0),
        legend.key.size = unit(0.1,"cm"),
        legend.box.margin=margin(0,0,0,0)) +
  guides(fill=guide_legend(nrow=1,byrow=TRUE))

plot_grid(
  var_exp_agg_plot + theme(legend.position = "none") + ggtitle("All sites"),
  var_exp_agg_rec_plot + theme(legend.position = "none") + ggtitle("Recurrent sites"),
  get_legend(var_exp_agg_rec_plot + theme(legend.box.spacing = unit(0,"cm"))),
  rel_heights = c(1,1,0.5),
  ncol = 1
) +
  ggsave(useDingbats=FALSE,filename = sprintf("figures/%s/statistical_analysis/variation_explained_aggregate.pdf",model_id),
         height = 1,width = 2)
```

### Assessing the goodness of fit with the residual effect

The goodness of fit is purely qualitative and is assessed by assessing whether a trajectory has any outliers. If a single outlier is present, we consider this trajectory to be unexplained. Outliers are defined as any point with a tail probability of less than 2.5% according to our model. Following this simple measure of goodness of fit, we estimate that our model explains 92.4% of all trajectories.

```{r model_ch_summary_statistic, fig.height=5,fig.width=6}
beta_value <- values_model$beta_values[[1]][1,1]

statistic_full_tails <- r_values_ %>%
  subset(key == 'VAFpred') %>%
  mutate(
    mu_val_005 = p_005,
    mu_val = p_050,
    mu_val_095 = p_095
  ) %>% 
  mutate(
    tail_prob = pbbinom(
      q = true,size = coverage,
      alpha = (mu_val * beta_value)/(1 - mu_val),beta = beta_value),
    tail_prob_005 = pbbinom(
      q = true,size = coverage,
      alpha = (mu_val_005 * beta_value)/(1 - mu_val_005),beta = beta_value),
    tail_prob_095 = pbbinom(
      q = true,size = coverage,
      alpha = (mu_val_095 * beta_value)/(1 - mu_val_095),beta = beta_value)
    ) %>% 
    mutate(
    tail_prob_genetic = pbbinom(
      q = true,size = coverage,
      alpha = (p_genetic_050 * beta_value)/(1 - p_genetic_050),beta = beta_value)
    )

statistic_full <- statistic_full_tails %>%
  group_by(site,individual,genes) %>% 
  summarise(sum_stat = all(tail_prob > 0.025 & tail_prob < 0.975),
            sum_stat_005 = all(tail_prob_005 > 0.025 & tail_prob_005 < 0.975),
            sum_stat_095 = all(tail_prob_095 > 0.025 & tail_prob_095 < 0.975),
            sum_stat_genetic = all(tail_prob_genetic > 0.025 & tail_prob_genetic < 0.975),
            n_tp = length(unique(age))) %>%
  mutate(split = 'full') %>%
  select(split,site,sum_stat,sum_stat_005,sum_stat_095,gene = genes,individual,sum_stat_genetic,n_tp) %>%
  ungroup() 

total_variants_explained <- statistic_full %>%
  group_by(site) %>%
  summarise(TotalExplainedVariants = sum(sum_stat == T),
            TotalVariants = length(sum_stat),
            TotalExplainedVariants005 = sum(sum_stat_005 == T),
            TotalExplainedVariants095 = sum(sum_stat_095 == T))

statistic_aggregated <- r_values_ %>%
  subset(key == 'VAFpred') %>%
  mutate(
    tail_prob = extraDistr::pbbinom(
    q = true,
    size = coverage,
    alpha = (mu_val * beta_value)/(1 - mu_val),
    beta = beta_value)) %>%
  group_by(site) %>% 
  summarise(sum_stat = pchisq(sum(qnorm(tail_prob)^2),length(tail_prob)),
            gene = genes[1],
            MinimumAge = min(age),
            MaximumAge = max(age),
            MaxValue = max(c(pred_095,value,true/coverage),na.rm = T)) %>%
  mutate(split = 'full') %>%
  ungroup() %>%
  merge(total_variants_explained,by = c("site"))

statistic_table <- statistic_full %>% 
  select(sum_stat) %>%
  mutate(sum_stat = sum_stat == T) %>%
  table

data.frame(
  Explained = statistic_table[2] / sum(statistic_table),
  `NotExplained` = statistic_table[1] / sum(statistic_table)
)

statistics_lists <- list(
  individual_site = statistic_full,
  site = statistic_aggregated
)

r_values <- merge(r_values_full,statistic_full,by = c("site","individual")) %>% 
  subset(sum_stat == T)

r_values %>% 
  select(Site = site,Gene = genes,
         Mean = b_site,Q005 = b_site_005,Q095 = b_site_095,Var = b_site_var) %>%
  distinct %>% 
  write.csv(sprintf("data_output/site_coefficients_%s.csv",model_id),quote = F,row.names = F)
r_values %>%
  select(Site = site,Gene = genes,SardID = individual,
         Mean = b_clone,Q005 = b_clone_005,Q095 = b_clone_095,Var = b_clone_var) %>%
  distinct %>% 
  write.csv(sprintf("data_output/mutation_independent_coefficients_%s.csv",model_id),quote = F,row.names = F)

r_values %>% 
  select(Gene = genes,
         Mean = b_gene,Q005 = b_gene_005,Q095 = b_gene_095,Var = b_gene_var) %>%
  distinct %>% 
  write.csv(sprintf("data_output/gene_coefficients_%s.csv",model_id),quote = F,row.names = F)
```

### Visualising all trajectories

```{r model_ch_dynamic_visualization_all_full,fig.width = 8.3,fig.height = 11.7}
r_values_merged <- merge(statistic_full,
                         r_values_,
                         all_y = F,
                         all.x = F,
                         by = c("site","individual")) %>%
  merge(select(statistic_full_tails,site,individual,tail_prob,tail_prob_005,tail_prob_095,age),
        all = F,
        by = c("site","individual","age")) %>%
  mutate(site_plot = gsub('t|nt','',site)) %>%
  mutate(group = as.factor(paste0(as.character(individual),key,site)))

r_values_merged_prediction_only <- r_values_merged %>%
  subset(key == 'VAFpred') %>% 
  mutate(`Inferred VAF` = value) %>% 
  group_by(individual,group) %>%
  mutate(M = diff(log(c(value[which.min(age)]+1e-5,value[which.max(age)])))) %>%
  group_by(individual) %>%
  mutate(M = max(M,na.rm=T))

M_order <- r_values_merged_prediction_only %>%
  ungroup() %>%
  select(individual,M) %>%
  distinct %>%
  arrange(M) %>%
  na.omit() %>%
  subset(M > 1) %>% 
  select(individual) %>%
  unlist
M_order <- M_order[2:length(M_order)] # keeps it nice and tidy
colour_breaks <- r_values %>% 
  select(gene,site,individual,b_site,b_gene) %>% 
  distinct %>% mutate(gene = str_match(gene,'[A-Z0-9]+')) %>% 
  group_by(gene) %>% summarise(M = mean(b_site+b_gene)) %>% 
  arrange(M) %>%
  ungroup %>%
  select(gene) %>% 
  unlist

r_values_merged_data_only <- r_values_merged %>%
  subset(key == "VAFtrue") %>%
  mutate(`Observed VAF` = value)

ggplot(data = NULL) +
  geom_ribbon(data = r_values_merged_prediction_only,
              aes(ymin = pred_005,
                  ymax = pred_095,
                  x = age,
                  group = group,
                  alpha = 0.1),
              na.rm = T) +
  geom_line(data = r_values_merged_data_only,
            aes(x = age,y = `Observed VAF`,
                colour = str_match(gene,'[A-Z0-9]+'),group = group),
            size = 0.25,linetype = 3) +
  geom_point(data = r_values_merged_data_only,
             aes(x = age,y = `Observed VAF`,
                 colour = str_match(gene,'[A-Z0-9]+')),
             size = 0.25) +
  geom_line(data = r_values_merged_prediction_only,
            aes(x = age,y = `Inferred VAF`,
                colour = str_match(gene,'[A-Z0-9]+'),group = group),
            size = 0.5) +
  facet_wrap(~ factor(individual,levels = M_order),scales="free",ncol = 11) + 
  scale_color_manual(name=NULL,values = gene_colours,breaks = colour_breaks) +
  theme_gerstung(base_size = 6) + 
  theme(strip.text = element_blank(),axis.text = element_blank(),
        legend.key.height = unit(0,"cm")) +
  xlab("Age") + 
  ylab("VAF") + 
  scale_alpha(range = c(0.05,0.05),guide=F) +
  scale_y_continuous(n.breaks = 4,trans = 'log10',
                     labels = function(x) sprintf("%s%%",x*100),
                     limits = c(1e-4,0.5)) +
  scale_x_continuous(breaks = seq(0,120,by = 5)) +
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/inferred_trajectories/plot_all_trajectories_full.pdf",model_id),
         width = 7,height = 8)
```

### Visualising example trajectories

We can see how the previous statistic becomes clearer as soon as we plot these trajectories.

```{r model_ch_dynamic_visualization_all,fig.width = 6.5,fig.height = 3.4}
r_values_merged <- merge(statistic_full,
                         r_values_,
                         all_y = F,
                         all.x = F,
                         by = c("site","individual")) %>%
  merge(select(statistic_full_tails,site,individual,tail_prob,tail_prob_005,tail_prob_095,age),
        all = F,
        by = c("site","individual","age")) %>%
  mutate(site_plot = gsub('t|nt','',site))

r_values_merged_prediction_only <- r_values_merged %>%
  subset(key == 'VAFpred') %>% 
  subset(!is.na(site_numeric))

r_values_merged %>% 
  subset(!is.na(site_numeric)) %>% 
  ggplot(aes(x = age,y = value)) +
  geom_ribbon(alpha = 0.1,
              data = r_values_merged_prediction_only,
              aes(ymin = pred_005,
                  ymax = pred_095,
                  group = as.factor(paste0(as.character(individual),key))),
              na.rm = T) +
  geom_line(aes(colour = key,group = as.factor(paste0(as.character(individual),key))),
            size = 0.25) +
  facet_wrap(~ site_plot,scales = "free",ncol = 8) + 
  scale_color_manual(guide=F,
                     values = c("red4","goldenrod1")) +
  theme_gerstung(base_size = 6) + 
  theme(legend.position = 'bottom') +
  xlab("Age") + 
  ylab("VAF") + 
  scale_x_continuous(n.breaks = 4,
                     breaks = function(x) floor(x[1]) + seq(1,x[2],by = round((x[2] - x[1])/4))) +
  scale_y_continuous(trans = 'log10') +
  ggtitle("Inferred (dark red) and observed (golden) trajectories for recurrent sites") +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/inferred_trajectories/plot_all_trajectories.pdf",model_id),
         width = 6.5,height = 3.4)
```

To better understand this at a higher resolution, by observing the trajectories below we can more finely comprehend what an outlier is in the context of this study, and we can observe how prediction works.

```{r model_dynamic_visualization,warning=F,fig.width = 4,fig.height = 2.5}
full_data %>% subset(amino_acid_change == "SF3B1-K666N") %>% select(SardID) %>% distinct()
site_subset_individuals <- list(
  `DNMT3Ant-R882H` = c(4240,1172,12467),
  `SF3B1-K700E` = c(26841,1905,25209),
  `SF3B1-K666N` = c(2259,3877,4546),
  `JAK2-V617F` = c(11184,3200,2141),
  `SRSF2-P95H` = c(5247,11449,500),
  `U2AF1-Q157R` = c(3877,10094,30271)
)
site_subset <- names(site_subset_individuals)

r_values_merged %>%
  select(-tail_prob_005,-tail_prob_095) %>% 
  pivot_wider(names_from = c("key"),values_from = c("pred_005","pred","pred_095")) %>%
  select(individual,site,sum_stat,true,age,coverage,genes,n,value,site_plot,
         pred_005 = pred_005_VAFpred,
         pred = pred_VAFpred,
         pred_095 = pred_095_VAFpred,
         tail_prob,
         gene) %>%
  na.omit() %>% 
  subset(site %in% site_subset) %>%
  group_by(site) %>%
  filter(individual %in% site_subset_individuals[as.character(site)][[1]]) %>% 
  mutate(site_plot = gsub(x = site_plot,'\n','-')) %>%
  mutate(site_plot = gsub('nt|t','',site_plot)) %>% 
  mutate(site_plot = paste0(
    sprintf('italic("%s")',str_match(genes,'[A-Z0-9]+')),
    str_match(site_plot,'-[a-zA-Z0-9]+'))) %>%
  ggplot(aes(x = age,y = pred/coverage)) +
  geom_ribbon(alpha = 0.1,
              aes(ymin = pred_005,
                  ymax = pred_095,
                  group = as.factor(as.character(individual))),
              na.rm = T) +
  geom_line(aes(colour = gene,
                group = as.factor(as.character(individual))),
            size = 0.5,
            alpha = 0.8) +
  geom_point(aes(y = (true)/coverage,
                 shape = tail_prob < 0.025 | tail_prob > 0.975),
             size = 0.8) +
  facet_wrap(~ site_plot,
             scales = "free",
             ncol = 3,
             labeller = label_parsed) + 
  scale_color_manual(guide = F,values = gene_colours) +
  scale_shape_manual(name = "Outlier",
                     breaks = c(T,F),
                     values = c(4,16),
                     labels = c("Yes","No")) +
  theme_gerstung(base_size = 6) + 
  theme(legend.position = 'bottom',
        legend.key.height = unit(0,"cm"),
        legend.margin = margin(0,0,0,0),
        legend.spacing = unit(0,"cm"),
        legend.background = element_blank()) +
  scale_y_continuous(trans = 'log10') + 
  coord_cartesian(ylim = c(5e-4,0.5)) +
  xlab("Age") + 
  ylab("VAF") +
  scale_alpha(guide = F,range = c(0.4,1.0)) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/inferred_trajectories/example_trajectories.pdf",model_id),
         width = 4,height = 2.5)
```

Observing how trajectories progress assuming that they start one the exact same timepoint allows us to get clear visual impression regarding the variability of these trajectories. A more comprehensive analysis is done further ahead. 

```{r model_dynamic_visualization_genetic_only,fig.width=6,fig.height=3,message=F,warning=F}
pop_size <- 0.5e5
r_values %>% 
  select(coefficient,b_clone,site,genes,u_values) %>% 
  distinct %>% 
  mutate(site = paste0(str_match(genes,'[0-9A-Z]+'),str_match(site,'-[A-Z0-9a-z]+'))) %>%
  mutate(full_effect = coefficient + b_clone) %>%
  apply(1,FUN = function(x) {
    data.frame(
      samples_genetic = trajectory_from_t0(s = as.numeric(x[1]),t0 = 0,x = c(0:20),
                                           N = pop_size,g = 2)/pop_size,
      samples = trajectory_from_t0(s = as.numeric(x[6]),t0 = 0,x = c(0:20),
                                   N = pop_size,g = 2)/pop_size,
      site = x[3],
      genes = x[4],
      x = seq(0,20),
      coefficient = x[1],
      full_effect = x[6],
      row.names = NULL
    ) %>%
      return
  }) %>% 
  do.call(what = rbind) %>%
  mutate(genes = ifelse(genes == "SRSF2",
                        ifelse(site == "SRSF2-P95H",
                               "SRSF2-P95H",
                               "SRSF2-others"),
                        as.character(genes))) %>%
  group_by(x,genes) %>%
  mutate(samples_genetic = mean(samples_genetic)) %>% 
  ungroup() %>% 
  mutate(site = factor(site,levels = sort(as.character(unique(site))))) %>% 
  na.omit() %>% 
  ggplot(aes(x = x,colour = genes)) +
  geom_line(aes(y = samples,group = full_effect),size = 0.25) + 
  geom_line(aes(y = samples_genetic,group = site),linetype = "longdash",size = 0.5,color = 'black') + 
  facet_wrap(~ genes,scales = "free",ncol = 6) +
  theme_gerstung(base_size = 6) +
  scale_y_continuous(trans = 'log10',
                     breaks = c(0.001,0.01,0.1,0.5),
                     labels = scientific(c(0.001,0.01,0.1,0.5)),
                     limits = c(1/pop_size,0.5)) + 
  scale_colour_manual(values = c(gene_colours,
                                 `SRSF2-P95H` = as.character(gene_colours["SRSF2"]),
                                 `SRSF2-others` = as.character(gene_colours["SRSF2"])),guide = F) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_line(size = 0.1)) +
  xlab("Time (years)") + 
  ylab("VAF")
```

### Fitness advantages at different levels of genetic resolution

We can now inspect the coefficients for all genomic resolutions (gene and site) individually and then observe the general picture of their combined effects. We analyse the differences between truncating and non-truncating mutations regarding their fitness advantage, observing that truncating mutations in *CBL* grow faster than non-truncating mutations and that non-truncating mutations in *TP53* grow faster than truncating mutations, with little evidence to allow us to make such statements regarding other genes.

The most relevant aspects are the wide range of dynamics observable for all genes and the seemingly weak influence of site identity in dynamics, with a clear exception being P95H in *SRSF2*, which grows much faster than mutations on the site site for different amino acids (P95A and P95L). Mutations in U2AF1 are unusually fast growers when compared with other genes. Another striking aspect of this analysis is how the prevalence of mutations is quite a poor predictor of their clonal dynamics - the two most frequently mutated genes, *TET2* and *DNMT3A*, do not grow anywhere near as fast as the fastest growing mutations, namely those in *U2AF1*, *PTPN11* and *SRSF2*-P95H.

```{r model_ch_gene_coefficients, fig.width=5.6, fig.height=2.4}
values_b_gene <- values_model$b_gene_values[[1]]
values_b_gene$site_name <- rep(unique_gene,
                               each = length(unique(values_model$b_gene_values[[1]]$labels)))
values_b_gene %>% spread(key = labels,value = values) %>%
  mutate(Truncating = grepl("[A-Z0-9]+t",site_name),
         site_name = str_match(site_name,"[A-Z0-9]+")) %>%
  ggplot(aes(y = mean,x = reorder(site_name,mean),color = site_name,shape = Truncating)) +
  geom_point(size = 2,
             position = position_dodge(width = 1.0)) +
  geom_linerange(aes(ymin = HDPI_low,ymax = HDPI_high),
                 position = position_dodge(width = 1.0)) + 
  geom_hline(yintercept = 0,alpha = 0.2) +
  theme_gerstung(base_size = 6) +
  rotate_x_text() + 
  xlab("Driver gene") + 
  ylab("Effect size") + 
  facet_grid(~ site_name,scales = "free_x",space = "free_x") +
  scale_color_manual(values = gene_colours,
                     guide = F) +
  theme(legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        strip.text = element_blank(),
        axis.text.x = element_text(face = 'italic'),
        legend.key.height = unit(0,"cm")
        ) + 
  scale_shape_discrete(breaks = c(TRUE,FALSE),
                       labels = c("Truncating","Non-truncating"),
                       name = NULL) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/gene_coefficients.pdf",model_id),
         height = 2.4,width = 5.6)
```

```{r truncating_vs_nontruncating_effect,fig.height=2,fig.width=1.5}
tnt_diff <- lapply(c(1:23),function(x) {
  data.frame(
    gene = unique_gene[x],
    samples = grab_samples(site_numeric = 1,gene_numeric = x,clone_numeric = 1,sample_size=2500)[[3]],
    row.names = NULL)}
  ) %>% 
  do.call(what = rbind) %>% 
  mutate(Truncating = grepl('[A-Z0-9]+t',gene),NonTruncating = grepl('[A-Z0-9]+nt',gene)) %>% 
  subset(Truncating | NonTruncating) %>% 
  mutate(gene_root = str_match(gene,'[A-Z0-9]+')) %>% 
  group_by(gene_root) %>% 
  summarise(diff_mean = mean(exp(samples[Truncating]) - exp(samples[NonTruncating])),
            diff_05 = quantile(exp(samples[Truncating]) - exp(samples[NonTruncating]),0.05),
            diff_95 = quantile(exp(samples[Truncating]) - exp(samples[NonTruncating]),0.95),
            diff_50 = median(exp(samples[Truncating]) - exp(samples[NonTruncating])),
            diff_sd = sd(exp(samples[Truncating]) - exp(samples[NonTruncating]))) %>%
  mutate(Sig = sign(diff_05) == sign(diff_95))

tnt_diff %>% 
  ggplot(aes(x = gene_root,y = diff_50,ymin = diff_05,ymax = diff_95,alpha = Sig)) + 
  geom_hline(yintercept = 0,linetype = 2,size = 0.25) + 
  geom_linerange() +
  geom_point() + 
  xlab("") + 
  ylab("Truncating - non-truncating effect") + 
  theme_gerstung(base_size = 6) + 
  scale_alpha_discrete(range = c(0.3,1),
                       guide = F) + 
  theme(axis.text.x = element_text(face = "italic")) + 
  rotate_x_text() + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/tnt_diff.pdf",model_id),
         height = 2,width = 1.5)
```

```{r model_ch_site_coefficients, fig.width=5.6, fig.height=2.4}
values_b_site <- values_model$b_site_values[[1]]
values_b_site$site_name <- rep(unique_site_multiple,
                               each = length(unique(values_model$b_values[[1]]$labels)))
values_b_site$gene <-  sapply(values_b_site$site_name,function(x) (strsplit(as.character(x),'-') %>% unlist)[1])

values_b_site %>% spread(key = labels,value = values) %>%
  ggplot(aes(y = mean,x = substr(site_name,1,30),color = gene)) +
  geom_point() +
  geom_linerange(aes(ymin = HDPI_low,ymax = HDPI_high)) + 
  geom_hline(yintercept = 0,alpha = 0.2) +
  theme_gerstung(base_size = 6) +
  rotate_x_text() + 
  xlab("") + 
  ylab("Effect size") + 
  facet_grid(~ gene,scales = "free_x",space = "free_x") +
  scale_color_manual(values = gene_colours,
                     guide = F) +
    theme(legend.position = "bottom",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        strip.text = element_text(face = "bold")) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/site_coefficients.pdf",model_id),
         height = 2.4,width = 5.6)
```

```{r model_ch_full_effects, fig.width=5.6, fig.height=2.4,warning=F}
min_max_plots <- c(min(c(r_values$coefficient_005,r_values$b_clone_005)),max(c(r_values$coefficient_095,r_values$b_clone_095)))
min_max_plots <- exp(min_max_plots) * 100 - 100
min_max_plots_growth <- sign(min_max_plots) * ceiling(abs(min_max_plots))
min_max_plots_round <- round(min_max_plots,-1)

gene_correspondence_truncation <- values_model$sub_data %>%
  group_by(Gene) %>%
  summarise(Truncating = unique(truncating)) %>%
  select(gene = Gene,Truncating)
  
values_b <- r_values %>%
  select(
    mean = coefficient,
    HDPI_low = coefficient_005,
    HDPI_high = coefficient_095,
    coefficient_var,
    site_name = site,
    gene = genes
  ) %>%
  distinct

gene_values <- values_model$b_gene_values[[1]] %>%
  mutate(
    gene = rep(unique_gene,
               each = length(unique(values_model$b_gene_values[[1]]$labels)))
  ) %>%
  subset(labels == "mean") 
gene_values[gene_values$gene == 'JAK2',] <- c(subset(values_b,gene == 'JAK2')$mean[1],"mean","","JAK2")
gene_values[gene_values$gene == 'IDH2',] <- c(subset(values_b,gene == 'IDH2')$mean[1],"mean","","IDH2")
gene_values <- gene_values %>% 
  mutate(values = as.numeric(values)) %>%
  mutate(gene_pure = str_match(gene,"[A-Z0-9]+")) %>%
  group_by(gene_pure) %>%  
  mutate(gene_average = mean(values)) %>%
  arrange(gene_average,values) %>%
  ungroup() %>% 
  select(gene_average = values,gene,gene_pure) 
gene_order <- gene_values %>% 
  select(gene_pure) %>%
  distinct %>% 
  unlist
sub_gene_order <- gene_order[gene_order %in% str_match(values_model$training_subset$unique_site_multiple,'[A-Z0-9a-z]+')]
gene_values[gene_values$gene == 'JAK2',] <- list(NA,'JAK2','JAK2')
gene_values[gene_values$gene == 'IDH2',] <- list(NA,'IDH2','IDH2')

gene_intervals <- values_model$b_gene_values[[1]] %>%
  mutate(
    gene = rep(unique_gene,
               each = length(unique(values_model$b_gene_values[[1]]$labels)))
  ) %>%
  subset(labels %in% c("HDPI_high","mean","HDPI_low")) %>% 
  spread(key = labels,value = values) %>%
  mutate(Truncating = grepl("[A-Z0-9]+t",gene)) %>%
  mutate(gene = str_match(gene,'[A-Z0-9]+')) %>% 
  select(-variable) %>%
  mutate(Truncating = ifelse(gene %in% c("ASXL1","PPM1D"),T,Truncating))

mutation_aa_order <- values_b %>%
  subset(site_name %in% unique_site_multiple) %>%
  mutate(site_name = as.character(site_name)) %>% 
  mutate(site_name = sapply(site_name,function(x) unlist(strsplit(x,'-'))[2])) %>%
  mutate(site_name = sapply(site_name,function(x) unlist(strsplit(x,';'))[1])) %>%
  mutate(site_name = sapply(site_name,function(x) {
    tmp <- unlist(strsplit(x,':'))
    return(tmp[length(tmp)])})) %>%
  arrange(mean) %>%
  select(site_name) %>%
  unlist

Number_Of_Mutations <- r_values %>%
  select(individual,site) %>%
  distinct() %>%
  group_by(site) %>%
  summarise(n = length(site)) %>%
  ungroup() %>%
  mutate(site_name = site) %>%
  select(-site)

plot_growth_per_year_data <- values_b %>% 
  merge(Number_Of_Mutations,by = 'site_name') %>%
  subset(site_name %in% unique_site_multiple) %>%
  mutate(Truncating = grepl("[A-Z0-9]+t",gene)) %>%
  group_by(gene) %>%
  mutate(mean = exp(mean)*100 - 100,
         HDPI_low = exp(HDPI_low)*100 - 100,
         HDPI_high = exp(HDPI_high)*100 - 100) %>%
  merge(gene_values,by = "gene",all.y=T) %>% 
  group_by(gene) %>%
  mutate(gene_average = as.numeric(gene_average)) %>% 
  mutate(gene_average = ifelse(is.na(gene_average),NA,exp(gene_average)*100 - 100)) %>% 
  ungroup() %>% 
  mutate(site_name = as.character(site_name)) %>% 
  mutate(site_name = sapply(site_name,function(x) unlist(strsplit(x,'-'))[2])) %>%
  mutate(site_name = sapply(site_name,function(x) unlist(strsplit(x,';'))[1])) %>%
  mutate(site_name = sapply(site_name,function(x) {
    tmp <- unlist(strsplit(x,':'))
    return(tmp[length(tmp)])})) %>%
  select(-Truncating) %>% 
  mutate(site_name = ifelse(is.na(site_name),' ',site_name)) %>% 
  merge(gene_correspondence_truncation,by = "gene") %>%
  mutate(gene = factor(as.character(gene_pure),levels = gene_order,
                       ordered = TRUE)) %>%
  mutate(Truncating_facet = ifelse(Truncating,'Trunc.','Non-Trunc.')) %>% 
  mutate(site_name = ifelse(site_name == " ",Truncating_facet,site_name)) %>% 
  group_by(gene,Truncating)

Segment_order <- plot_growth_per_year_data %>% 
  mutate(SiteNumeric = as.numeric(factor(site_name,levels = site_name[order(mean)]))) %>%
  group_by(gene) %>%
  mutate(NT = sum(Truncating==F)) %>%
  group_by(gene,Truncating) %>%
  summarise(SiteNumericMax = max(SiteNumeric),
            SiteNumericMin = min(SiteNumeric),
            NT = NT[1],
            gene_average = gene_average[1]) %>%
  mutate(SiteNumericMax = ifelse(Truncating == T,SiteNumericMax + NT,SiteNumericMax),
         SiteNumericMin = ifelse(Truncating == T,SiteNumericMin + NT,SiteNumericMin)) %>%
  merge(gene_intervals,by = c("gene","Truncating"))

plot_growth_per_year <- plot_growth_per_year_data %>% 
  ggplot(aes(color = gene)) +
  #geom_point(aes(colour = gene,y = gene_average),shape=95,size=9,alpha = 1.0) +
  geom_linerange(aes(ymin = HDPI_low,ymax = HDPI_high,
                     x = reorder(site_name,mean+Truncating*1000))) + 
  geom_point(aes(size = n,
                 y = mean,
                 x = reorder(site_name,mean+Truncating*1000))) +
  geom_hline(yintercept = 0,alpha = 0.2) +
  geom_segment(
    data = Segment_order,
    aes(x = SiteNumericMin-0.5,xend = SiteNumericMax+0.5,
        y = gene_average,yend = gene_average),
    alpha=0.7) +
  geom_rect(
    data = Segment_order,
    inherit.aes = F,
    aes(xmin = SiteNumericMin-0.5,xmax = SiteNumericMax+0.5,
        ymin = exp(HDPI_high)*100-100,ymax = exp(HDPI_low)*100-100,
        fill = gene),colour = NA,
    alpha = 0.2) +
  theme_gerstung(base_size = 6) +
  rotate_x_text() + 
  xlab(" ") + 
  ylab("Driver growth per year") + 
  facet_grid(~ gene,scales = "free_x",space = "free_x") +
  scale_color_manual(values = gene_colours,
                     guide = F) +
  scale_fill_manual(values = gene_colours,
                    guide = F) +
  theme(legend.position = "bottom",
        legend.background = element_rect(colour = NA,fill = "white"),
        panel.grid = element_blank(),
        legend.key.height = unit(0,"cm"),
        axis.text.x = element_text(face = "plain")) +
  #coord_cartesian(ylim = c(-10,60)) + 
  scale_y_continuous(breaks = seq(min_max_plots_round[1],min_max_plots_round[2],by = 20),
                     labels = paste0(seq(min_max_plots_round[1],min_max_plots_round[2],by = 20),'%'),
                     limits = min_max_plots_growth,
                     expand = c(0.05,0,0.0,0)) +
  scale_shape_manual(breaks = c(TRUE,FALSE),
                     values = c(16,17),
                     labels = c("Truncating","Non-truncating"),
                     name = NULL) +
  scale_size(range = c(1,2),
             breaks = c(2,10),
             name = "N") +
  guides(size = guide_legend(nrow = 1,direction = "horizontal")) 

g <- ggplot_gtable(ggplot_build(plot_growth_per_year))
strip_both <- which(grepl('strip-', g$layout$name))
axis_x <- which(grepl('axis-b', g$layout$name))
colors <- unique(gene_colours[gene_values$gene_pure])
gene_labels <- gene_values$gene_pure
k <- 1
m <- 1
for (i in strip_both) {
  j <- which(grepl("text", g$grobs[[i]]$grobs[[1]]$childrenOrder))
  L <- g$grobs[[i]]$grobs[[1]]$children[[j]]$children[[1]]$label
  if (L %in% gene_values$gene_pure) {
    g$grobs[[i]]$grobs[[1]]$children[[j]]$children[[1]]$gp$col <- colors[k]
    g$grobs[[i]]$layout$clip <- "off"
    #nt <- grepl('nt',g$grobs[[i]]$grobs[[1]]$children[[2]]$children[[1]]$label)
    if (k %% 2 == 0) {
      g$grobs[[i]]$heights <- unit.c(unit(6,"pt"),unit(-0.2,"cm"),unit(6,"pt"))
    } else {
      g$grobs[[i]]$heights <- unit.c(unit(6,"pt"),unit(0.2,"cm"),unit(6,"pt"))
    }
    m <- m+1
    k <- k+1
  } else {
    #g$grobs[[i]]$heights <- unit.c(unit(6,"pt"),unit(-,"cm"),unit(6,"pt"))
  }
}
for (i in axis_x) {
  L <- g$grobs[[i]]$children$axis$grobs[[2]]$children[[1]]$label
  gp <- g$grobs[[i]]$children$axis$grobs[[2]]$children[[1]]$gp
  ff <- ifelse(grepl('runc',L),'bold','plain')
  s <- ifelse(grepl('runc',L),gp$fontsize+1,gp$fontsize-1)
  gp <- gpar(
    fontsize = s,
    col = gp$col,
    fontfamily = gp$fontfamily,
    lineheight = gp$lineheight,
    fontface=ff)
  g$grobs[[i]]$children$axis$grobs[[2]]$children[[1]]$gp <- gp
}

grid.newpage()
grid.draw(g) 

ggsave(useDingbats=FALSE,plot = gridExtra::grid.arrange(g),
       sprintf("figures/%s/dynamic_coefficients/driver_coefficients.pdf",model_id),
       width = 5.6,height = 2.4) 
```

```{r model_ch_full_effects_simplified, fig.width=3, fig.height=2,warning=F}
label_p95h <- plot_growth_per_year_data %>%
  subset(gene == 'SRSF2') %>%
  mutate(site_name = ifelse(site_name == 'P95H',site_name,NA))

unknown_cause_growth <- r_values %>%
  select(
    b_clone_005,b_clone,b_clone_095,
    individual,
    site,
    gene
  ) %>%
  distinct %>%
  mutate(gene = str_match(gene,'[A-Z0-9]+'))

plot_growth_per_year <- gene_intervals %>% 
  mutate(gene = factor(gene,levels = gene_order)) %>% 
  ggplot(aes(color = gene)) +
  geom_hline(yintercept = c(0),colour = "black",size = 0.25) +
  geom_hline(yintercept = c(-20,20,40,60,80),size = 0.25,colour = "grey80") +
  geom_tile(
    data = gene_intervals,
    inherit.aes = F,
    aes(x = Truncating,
        y = exp(mean)*100 - 100,
        width = 1,
        color = gene,
        height = 0.01),
    alpha=0.7) +
  geom_tile(
    data = gene_intervals,
    inherit.aes = F,
    aes(x = Truncating,
        y = ((exp(HDPI_high)*100-100) + (exp(HDPI_low)*100-100))/2,
        width = 1,
        height = exp(HDPI_high)*100 - exp(HDPI_low)*100,
        fill = gene),
    colour = NA,
    alpha = 0.2) +
  geom_linerange(
    data = plot_growth_per_year_data,
    aes(ymin = HDPI_low,
        ymax = HDPI_high,
        x = Truncating,
        alpha = as.numeric(site_name == "P95H"),
        group = reorder(site_name,mean)),
    position = position_dodge(width = 0.8),
    size = 0.25) +
  geom_point(
    data = plot_growth_per_year_data,
    aes(size = n,
        y = mean,
        x = Truncating,
        alpha = as.numeric(site_name == "P95H"),
        group = reorder(site_name,mean)),
    position = position_dodge(width = 0.8)) +
  theme_gerstung(base_size = 6) +
  rotate_x_text() + 
  xlab(" ") + 
  ylab("Driver growth per year") + 
  facet_grid(~ factor(gene,levels = gene_order),
             scales = "free_x",
             space = "free_x") + 
  scale_color_manual(values = gene_colours,
                     guide = F) +
  scale_fill_manual(values = gene_colours,
                    guide = F) +
  theme(legend.position = c(0.2,0.8),
        legend.background = element_rect(colour = NA,fill = "white"),
        strip.text = element_text(face = "bold.italic"),
        legend.key.height = unit(0,"cm"),
        axis.text = element_text(face = "plain")) +
  scale_y_continuous(breaks = seq(-20,100,by = 20),
                     labels = paste0(seq(-20,100,by = 20),'%'),
                     limits = min_max_plots_growth,
                     expand = c(0.05,0,0.0,0)) +
  scale_size(range = c(0.5,2),
             breaks = c(2,10,25,50),
             name = "N") +
  scale_x_discrete(breaks = c("FALSE","TRUE","uc"),
                   labels = c("NT","T","UC"),
                   expand = c(0.01,0.01)) +
  scale_alpha(range = c(0.5,1.0),guide = F) +
  guides(size = guide_legend(nrow = 1,direction = "horizontal")) 

library(grid)

g <- ggplot_gtable(ggplot_build(plot_growth_per_year))
strip_both <- which(grepl('strip-', g$layout$name))
axis_x <- which(grepl('axis-b', g$layout$name))
colors <- unique(gene_colours[gene_values$gene_pure])
gene_labels <- gene_values$gene_pure
k <- 1
m <- 1
for (i in strip_both) {
  j <- which(grepl("text", g$grobs[[i]]$grobs[[1]]$childrenOrder))
  L <- g$grobs[[i]]$grobs[[1]]$children[[j]]$children[[1]]$label
  if (L %in% gene_values$gene_pure) {
    g$grobs[[i]]$grobs[[1]]$children[[j]]$children[[1]]$gp$col <- gene_colours[L]
    g$grobs[[i]]$layout$clip <- "off"
    if (k %% 2 == 0) {
      g$grobs[[i]]$heights <- unit.c(unit(6,"pt"),unit(-0.4,"cm"),unit(6,"pt"))
    } else {
      g$grobs[[i]]$heights <- unit.c(unit(6,"pt"),unit(0,"cm"),unit(6,"pt"))
    }
    m <- m+1
    k <- k+1
  } else {
    #g$grobs[[i]]$heights <- unit.c(unit(6,"pt"),unit(-,"cm"),unit(6,"pt"))
  }
}
for (i in axis_x) {
  g$grobs[[i]]$children$axis$grobs$`1`$children[[1]]$gp$font <- ifelse(
    g$grobs[[i]]$children$axis$grobs$`1`$children[[1]]$label == 'UC',
    3,
    2
  ) %>%
    as.integer()
}

grid.newpage()
grid.draw(g) 

ggsave(useDingbats=FALSE,plot = gridExtra::grid.arrange(g),
       sprintf("figures/%s/dynamic_coefficients/driver_coefficients_simplified.pdf",model_id),
       width = 3.2,height = 2) 
```

#### Inspecting driver genetic effects

##### Gender

We find a statistically significant difference between male and female individuals. Since this difference will be dependent on the acquisition of specific mutations, we further look into this. 

```{r differences_sex, fig.width=3,fig.height=2}
totals <- table(distinct(select(full_data,SardID,Gender))$Gender)

N <- full_data %>%
  group_by(Gender) %>%
  summarise(N = length(unique(SardID))) %>%
  arrange(Gender)
  
effect_distribution_sex <- full_data %>%
  select(individual = SardID,Gender) %>% 
  distinct %>%
  merge(r_values %>% select(individual,site,coefficient,b_clone,genes),by = c('individual')) %>%
  group_by(individual) %>%
  mutate(effect = coefficient) %>% 
  distinct() %>% 
  filter(effect == max(effect)) %>%
  mutate(Gender = ifelse(Gender == 'F','Female','Male')) %>%
  mutate(Truncating = !grepl("nt",genes)) %>%
  mutate(Truncating = ifelse(is.na(Truncating),F,T)) %>%
  mutate(genes = str_match(genes,'[A-Z0-9]+')) %>%
  ggplot(aes(x = as.factor(Gender),y = effect)) + 
  geom_boxplot(alpha = 0.8,outlier.size = 0.5,
               size = 0.2) + 
  theme_gerstung(base_size = 6) + 
  geom_signif(comparisons = list(c("Male","Female")),
              textsize = 2,
              size = 0.2,
              map_signif_level = function(x) return(paste("p-value=",format(x,scientifier = T,digits = 2))),
              test = wilcox.test) + 
  ylab("Driver gene effect") + 
  scale_x_discrete(breaks = c("Female","Male"),
                   labels = c(sprintf("Female\n(n=%s)",N$N[1]),
                              sprintf("Male\n(n=%s)",N$N[2]))) +
  xlab("Gender") + 
  theme(axis.text = element_text(size = 6))

gene_mutation_distribution_sex <- load_data_keep() %>%
  select(SardID,Gender,Gene,AAChange.refGene) %>% 
  subset(Gene %in% load_included_genes() | is.na(Gene)) %>% 
  group_by(Gender) %>% 
  mutate(Total = length(unique(SardID))) %>%
  subset(!is.na(Gene)) %>% 
  select(individual = SardID,Gender,Gene,Total) %>% 
  distinct %>%
  mutate(Gender = ifelse(Gender == 'F','Female','Male')) %>%
  mutate(Gene = str_match(Gene,'[A-Z0-9]+')) %>%
  group_by(Gender,Total,Gene) %>% 
  summarise(N = length(unique(individual))) %>%
  mutate(Proportion = N / Total) %>% 
  ggplot(aes(x = Gene,y = Proportion,fill = Gender)) + 
  geom_bar(stat = "identity",position = 'dodge') + 
  theme_gerstung(base_size = 6) +
  scale_fill_lancet(name = NULL) + 
  rotate_x_text() + 
  ylab("Proportion") + 
  xlab("") + 
  theme(axis.text = element_text(size = 6),
        axis.text.x = element_text(face = "italic"),
        legend.position = "top",
        legend.key.height = unit(0.2,"cm"),
        legend.key.width = unit(0.2,"cm")) +
  scale_y_continuous(expand = c(0,0))

effect_distribution_sex %>% 
  ggsave(useDingbats=FALSE,filename = sprintf("figures/%s/dynamic_coefficients/sex_driver_effect.pdf",model_id),
         width = 2,height = 2)

D <- load_data_keep() %>% 
  group_by(SardID,Gene,Gender) %>% 
  summarise(Age = min(Age)) %>% 
  distinct %>% 
  group_by(SardID,Gender,Age) %>% 
  summarise(HasU2AF1 = any(Gene == "U2AF1")) %>%
  mutate(HasU2AF1 = ifelse(is.na(HasU2AF1),F,HasU2AF1))

glm(HasU2AF1 ~ Gender + Age, data = D,
    family = stats::binomial()) %>% summary

plot_grid(effect_distribution_sex,gene_mutation_distribution_sex,
          nrow = 1,rel_widths = c(0.6,1)) + 
  ggsave(useDingbats=FALSE,
         filename = sprintf("figures/%s/dynamic_coefficients/sex_driver_effect_barplot.pdf",model_id),
         width = 3,height = 2)
```

##### Smoking and driver dynamics and CH

We find a statistically significant difference in driver effects between people who smoked in the past and those who have not smoked. We further analyse this to assess whether any relevant difference between mutation acquisition is evident but find no evidence for this (including the previously reported bias for *ASXL1* mutation acquisition in smokers, likely due to the lack of statistical power). We also find no association between past or current smoking experience and mutation acquisition or frequency.  

```{r smoking_gene_composition}
smoking_frequencies_data <- load_smoking_data() %>% 
  mutate(status = ifelse(QuitSmoking == 'N' & CurrentSmoker == 'N',"Never smoked",
                         ifelse(QuitSmoking == 'N' & CurrentSmoker == 'Y',"Smoker",
                                "Previous smoker"))) %>% 
  subset(!is.na(CurrentSmoker)) %>% 
  group_by(SardID) %>%
  summarise(
    status = ifelse(
      sum(status == 'Smoker') > 0,'Smoked during study',
      ifelse(sum(status == 'Previous smoker') > 0, 'Smoked before study',
             'Never smoked')
    )
  ) %>% 
  mutate(status_binary = ifelse(status == 'Never smoked','Never smoked','Smoked')) %>% 
  merge(distinct(select(load_data_keep(),Age,SardID,Gene,Gender,mutation_identifier)),by = "SardID") %>%
  group_by(SardID) %>% 
  mutate(Age = min(Age),
         NMut = length(unique(mutation_identifier))) %>% 
  select(genes=Gene,SardID,status,status_binary,Age,Gender,NMut) %>%
  distinct() %>%
  mutate(genes = str_match(genes,'[A-Z0-9]+')) %>%
  distinct() %>%
  group_by(status_binary) %>% 
  mutate(Total = length(unique(SardID))) %>%
  group_by(genes,status_binary) %>%
  mutate(N = length(unique(SardID)),
         Freq = length(unique(SardID)) / Total) %>%
  distinct

mut_data <- smoking_frequencies_data %>% 
  filter(genes %in% c(NA,names(gene_colours))) %>%
  mutate(Present = 1) %>% 
  spread(key = genes,value = Present,fill = 0) %>%
  select(-`<NA>`) %>% 
  group_by(SardID,status_binary,Gender) %>% 
  mutate(
    Age = min(Age),
    Smoked = ifelse(status_binary == "Smoked",T,F)
  ) %>%
  distinct %>% 
  ungroup() %>% 
  select(-status,-status_binary) %>% 
  group_by(SardID,Smoked,Age,Gender,NMut) %>%
  summarise_all(sum)

mut_data$HasMut <- as.numeric(rowSums(mut_data[,c(9:25)]) > 0)
glm(formula = as.formula(sprintf(
  "Smoked ~ %s + Age + Gender",paste(load_included_genes(),collapse = " + ")
)),
    family = stats::binomial(),
    data = mut_data) %>%
    summary 

glm(formula = HasMut ~ Age + Gender + Smoked,
    family = stats::binomial(),
    data = mut_data) %>%
    summary 

glm(formula = NMut ~ Smoked + Age + Gender,
    family = stats::poisson(),
    data = smoking_frequencies_data %>%
      group_by(SardID,status_binary,Gender) %>%
      summarise(Age = min(Age),
                NMut = NMut[1],
                Smoked = ifelse(status_binary == "Smoked",T,F)) %>%
      distinct) %>%
  summary
```

#### Inspecting unknown cause effects

The advantage of having unknown cause effects is that we have a coefficient which corresponds to effects that affect clone growth at an individual clone level, allowing us to more accurately quantify the genetic growth without the influences of deviations in growth dynamics at an individual clone and quantify deviations from the expected genetic growth. In doing so, we can 1) observe and quantify deviations from genetic effects and inspect the trajectories corresponding to them and 2) see if age has any significant association in these coefficients, which would entail that age has an effect on the growth dynamics of different clones. 

In this section we analyse possible effects influencing the unknown-cause growth effect, both technical - such as the explicit modelling of recurring sites - and biological - such as clonality and within- and between- individual variance.

##### Between- and within-individual variance

We begin by estimating the within and between individual variance. We show here that there is no difference between inter- and intra-individual variance, which would further confirm the independence of clone onset even within individuals. 

```{r within_between_individual_variance,fig.height=1.2,fig.width=1.5}
within_between_individual_variance <- r_values %>%
  select(individual,site,b_clone) %>% 
  distinct %>% 
  mutate(M = mean(b_clone)) %>% 
  group_by(individual) %>% 
  filter(length(unique(site)) > 1) %>% 
  summarise(V = sum((b_clone - mean(b_clone))^2),
            M = sum((mean(b_clone) - M)^2),
            NClones = length(unique(site))) %>%
  ungroup() %>%
  mutate(NIndividuals = length(unique(individual))) %>%
  summarise(WithinIndividualVariance = sum(V) /(sum(NClones)-NIndividuals[1]),
            BetweenIndividualVariance = sum(M)/(NIndividuals[1]-1),
            NIndividuals = NIndividuals[1],
            NClones = sum(NClones)) %>%
  mutate(RatioOfVariances = BetweenIndividualVariance/WithinIndividualVariance) %>%
  mutate(`p-value` = pf(RatioOfVariances,df1 = NIndividuals - 1,
                        df2 = NClones - NIndividuals,lower.tail = F))
ggplot(data = NULL) +
  geom_bar(aes(x = c("Within\nindividuals","Between\nindividuals"),
               y = c(within_between_individual_variance$WithinIndividualVariance,
                     within_between_individual_variance$BetweenIndividualVariance)),
           fill = "grey",
           stat = "identity") + 
  xlab("") + 
  ylab("Variance") + 
  theme_gerstung(base_size = 6) + 
  theme(plot.subtitle = element_text(size = 6)) + 
  scale_y_continuous(breaks = c(0,0.001),
                     expand = c(0,0,0.01,0)) +
  coord_flip() + 
  ggtitle(sprintf("ratio of variances = %.2f\n(p-value = %.2f)",
                  within_between_individual_variance$RatioOfVariances,
                  within_between_individual_variance$`p-value`)) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_var.pdf",model_id),
         height = 1.2,width = 1.5) 
```

##### Truncating and recurring status

It is also important to understand if there are meaningful differences between the groups defined in this study, particularly when it comes to the explicit modelling of a few recurring sites. Here, we also test the possibility of clones with truncating mutations having different UC effects. As shown below, we find no evidence for either of these. 

```{r trucating_recurring,fig.width=2.0,fig.height=1.2}
truncating_recurring_uc <- r_values %>%
  select(individual,site,b_clone) %>%
  distinct %>% 
  mutate(truncating = grepl("[A-Z0-9]+t-",site)) %>% 
  mutate(recurring = site %in% unique_site_multiple) %>%
  mutate(f = paste(
    ifelse(recurring,"Recurring","Non-recurring"),
    ifelse(truncating,"truncating","missense"),
    sep = ' and\n'
  ))
anvar_test <- anova(glm(b_clone ~ f,data = truncating_recurring_uc),
                    test = "F")
truncating_recurring_uc %>%
  ggplot(aes(x = f,y = b_clone)) + 
  geom_boxplot(outlier.size = 0.5,size = 0.25) + 
  xlab("") + 
  ylab("UC effect") +
  theme_gerstung(base_size = 6) + 
  coord_flip() + 
  scale_y_continuous(breaks = c(-0.1,0.0,0.1,0.2),
                     labels = c("-10%","0%","10%","20%")) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_tnt_rec.pdf",model_id),
         height = 1.2,width = 2.0) 
anvar_test
```

##### Co-clonality and unknown cause effect

From our data we cannot observe any association between the unknown cause growth effect of the fastest clone and the number of mutations in each individual while controlling for age. We also conduct this analysis while stratifying by gene. To do so, we compare a multivariate quasipoisson regression which models the number of mutations as dependent on the UC growth effect and age with a null quasipoisson regression which considers the number of mutations in an individual as solely dependent on age. 

We also a stastically significant effect when considering only the UC growth effects of *U2AF1*, particularly that of a -17.2% fold-change.

```{r number_of_mutations, fig.height=1.5,fig.width=1}
n_mut_ind <- load_data() %>%
  group_by(individual = SardID) %>%
  summarise(NMut = length(unique(mutation_identifier)),
            Age = mean(unique(Age)))
n_mut_r_values <- r_values %>%
  merge(n_mut_ind,by = "individual") %>%
  group_by(individual,b_clone,b_genetic = b_genetic_mean,NMut,gene) %>%
  summarise(Age = mean(age)) %>%
  mutate(gene = str_match(gene,"[A-Z0-9]+")) %>% 
  group_by(gene) %>% 
  filter(length(unique(individual)) > 3) %>%
  mutate(b_clone = exp(b_clone)*100 - 100,
         total_effect = exp(b_clone + b_genetic)*100 - 100)

full_model = glm(
  b_clone ~ Age + NMut,
  data = summarise(group_by(n_mut_r_values,individual,NMut,Age),
                   b_clone = mean(b_clone))
)
null_model = glm(
  b_clone ~ Age,
  data = summarise(group_by(n_mut_r_values,individual,NMut,Age),
                   b_clone = mean(b_clone))
)
all_gene_uc_associations <- list(
  `All genes` = list(model = full_model,
                     null_model = null_model,
                     anova = anova(full_model,null_model,test = "Chisq"))
)
for (G in unique(n_mut_r_values$gene)) {
  d <- subset(n_mut_r_values,gene == G) %>%
    group_by(individual,Age,NMut) %>%
    summarise(b_clone = mean(b_clone))
  null_model <- glm(
    b_clone ~ Age,data = d)
  full_model <- glm(
    b_clone ~ Age + NMut,data = d)
  all_gene_uc_associations[[G]] <- list(
    model = full_model,
    null_model = full_model,
    anova = anova(full_model,null_model,test = "Chisq")
  )
}

sig_df <- data.frame(
  gene = names(all_gene_uc_associations),
  p.val = all_gene_uc_associations %>% 
    lapply(function(x) x$anova$`Pr(>Chi)`[2]) %>%
    unlist
)

sig_df$SignificantAfterCorrection <- p.adjust(sig_df$p.val,method = "BH") < 0.05

sig_df
```

##### Age, sex and smoking status

To explain the greater prevalence of CH in older ages, it could be postulated that clone development is accelerated by age, dependent on sex or accelerated by previous smoking experience. Here, we find no evidence for such an effect. As such, it is unlikely that either of these factors plays a role in accelerating clonal dynamics in a senior population.

```{r unknown_cause_sex_smoked_association, fig.width=2,fig.height=1.3}
good_stat_sites <- statistic_full %>%
  subset(sum_stat == T) %>% 
  mutate(X = paste(individual,site)) %>% 
  select(X) %>% 
  unlist
Parameters_full <- r_values_full %>% 
  select(b_clone_mean = b_clone,b_clone_var,b_clone_005,b_clone_095,
         b_gene_mean = b_gene,b_gene_var,b_gene_005,b_gene_095,
         b_site_mean = b_site,b_site_var,b_site_005,b_site_095,
         b_genetic_mean,b_genetic_var,b_genetic_005,b_genetic_095,
         u_values_005,u_values_095,
         individual,age,
         true,coverage,
         site) %>%
  group_by(individual,site) %>%
  mutate(age = min(age[which(round(true/coverage,4) >= 0.005)])) %>%
  mutate(recurring = site %in% unique_site_multiple) %>% 
  select(-true,-coverage) %>%
  distinct %>%
  mutate(gene = str_match(site,'[0-9A-Z]+')) 
Parameters <- Parameters_full %>%
  subset(paste(individual,site) %in% good_stat_sites) 

pearson_correlations <- Parameters %>% 
  group_by(gene) %>% 
  filter(length(unique(paste(individual,site))) > 10) %>%
  summarise(x = list(cor.test(age,b_clone_mean,
                              method = "spearman",exact = F))) %>%
  rowwise() %>%
  mutate(rho = x$estimate,
         p.val = x$p.val) %>%
  select(-x) 

all_gene_pearson_correlation <- cor.test(Parameters$age,Parameters$b_clone_mean)
all_gene_except_tet2_pearson_correlation <- cor.test(
  subset(Parameters,gene != 'TET2')$age,subset(Parameters,gene != 'TET2')$b_clone_mean)
pearson_correlations <- rbind(
  pearson_correlations,
  data.frame(
    gene = "All genes",
    rho = all_gene_pearson_correlation$estimate,
    p.val = all_gene_pearson_correlation$p.value,
    stringsAsFactors = F
  ),
  data.frame(
    gene = "All genes except TET2",
    rho = all_gene_except_tet2_pearson_correlation$estimate,
    p.val = all_gene_except_tet2_pearson_correlation$p.value,
    stringsAsFactors = F
  )
)

smoking_sex_data <- load_smoking_data() %>% 
  group_by(SardID) %>%
  summarise(HasSmoked = any(PackYears>0)) %>%
  merge(distinct(select(full_data,SardID,Gender))) %>%
  select(individual = SardID,HasSmoked,Gender)
min_max_plots <- c(min(c(r_values$coefficient_005,r_values$b_clone)),
                   max(c(r_values$coefficient_095,r_values$b_clone)))
min_max_plots <- exp(min_max_plots) * 100 - 100
min_max_plots_round <- round(min_max_plots,-1)

tmp_2 <- r_values %>%
  mutate(vaf = true/coverage) %>%
  group_by(individual,site,b_clone,b_genetic_mean,genes) %>%
  summarise(vaf_strict = vaf[vaf > 0.005][which.min(age[vaf > 0.005])],
            vaf_median = median(vaf[vaf > 0.005],na.rm=T),
            vaf_last = vaf[which.max(age)],
            vaf = vaf[which.min(age)],
            age = mean(age)) %>% 
  ungroup %>% 
  mutate(genes = str_match(genes,'[A-Z0-9]+')) %>% 
  merge(smoking_sex_data,by = "individual") %>% 
  na.omit

D <- tmp_2 %>%
  na.omit %>% 
  group_by(individual,HasSmoked,Gender,age) %>% 
  summarise(b_clone = mean(b_clone))
anova_age_sex_smoking <- list(
  `All genes` = anova(
    glm(b_clone ~ 1,data = D),
    glm(b_clone ~ age + Gender + HasSmoked,data = D),test = "F")
)

for (G in unique(tmp$genes)) {
  D <- subset(tmp,genes == G) %>%
    na.omit %>% 
    group_by(individual,HasSmoked,Gender,age) %>% 
    summarise(b_clone = mean(b_clone))
  if (nrow(D) > 5) {
    anova_age_sex_smoking[[G]] <- anova(
      glm(b_clone ~ 1,data = D),
      glm(b_clone ~ age + Gender + HasSmoked,data = D),test = "F")
  }
  }

p_values <- anova_age_sex_smoking %>% 
  lapply(function(x) x$`Pr(>F)`[2]) %>% 
  unlist %>%
  as.tibble %>% 
  mutate(Gene = names(anova_age_sex_smoking)) %>%
  mutate(`F-statistic` = sapply(anova_age_sex_smoking,function(x) x$F[2]))

adjusted_pvalues <- p.adjust(c(p_values$value,pearson_correlations$p.val),
                             method = "BH") < 0.05
p_values$Significant <- adjusted_pvalues[1:nrow(p_values)]
pearson_correlations$Significant <- adjusted_pvalues[(nrow(p_values)+1):length(adjusted_pvalues)]

p_values %>%
  select(Gene,`F-statistic`,`p-value` = value,Significant)

uc_age_plot <- tmp_2 %>%
  ggplot(aes(x = age,y = exp(b_clone)*100-100,colour = genes)) + 
  geom_point(alpha = 0.4,size = 0.5) + 
  theme_gerstung(base_size = 6) + 
  xlab("Age at study entry") + 
  ylab("Average UC effect") + 
  scale_color_manual(values = gene_colours,name = "Gene") + 
  guides(colour=guide_legend(ncol=2,byrow=F)) +
  theme(legend.position = 'none') +
  scale_y_continuous(breaks = seq(-10,100,by = 10),
                     labels = paste0(seq(-10,100,by = 10),'%'),
                     expand = c(0.05,0,0.0,0))

uc_vaf_plot <- tmp_2 %>%
  ggplot(aes(x = vaf,y = exp(b_clone)*100-100,colour = genes)) + 
  geom_point(alpha = 0.4,size = 0.5) + 
  theme_gerstung(base_size = 6) + 
  xlab("VAF at study entry") + 
  ylab("UC effect") + 
  scale_color_manual(values = gene_colours,name = "Gene") + 
  guides(colour=guide_legend(ncol=2,byrow=F)) +
  theme(legend.position = 'none') +
  scale_x_continuous(trans = 'log10') +
  scale_y_continuous(breaks = seq(-10,100,by = 10),
                     labels = paste0(seq(-10,100,by = 10),'%'),
                     expand = c(0.05,0,0.0,0))

uc_vaf_last_plot <- tmp_2 %>%
  ggplot(aes(x = vaf_last,y = exp(b_clone)*100-100,colour = genes)) + 
  geom_point(alpha = 0.4,size = 0.5) + 
  theme_gerstung(base_size = 6) + 
  xlab("VAF at study exit") + 
  ylab("UC effect") + 
  scale_color_manual(values = gene_colours,name = "Gene") + 
  guides(colour=guide_legend(ncol=2,byrow=F)) +
  theme(legend.position = 'none') +
  scale_x_continuous(trans = 'log10') +
  scale_y_continuous(breaks = seq(-10,100,by = 10),
                     labels = paste0(seq(-10,100,by = 10),'%'),
                     expand = c(0.05,0,0.0,0))

uc_sex_smoker_plot <- tmp_2 %>% 
  mutate(
    L = paste(ifelse(Gender == "M","Male","Female"),
              ifelse(HasSmoked == T,"\nhas smoked","\nnever smoked"),
              sep = '')
  ) %>%
  ggplot(aes(x = L,y = b_clone, colour = Gender)) + 
  geom_boxplot(size = 0.2,outlier.size = 0.5) + 
  theme_gerstung(base_size = 6) + 
  coord_flip() +
  xlab("") + 
  ylab("Average UC effect") +
  scale_y_continuous(breaks = c(-0.1,0,0.1,0.2),
                     labels = sprintf("%s%%",c(-10,0,10,20))) +
  scale_colour_aaas(guide = F)

uc_sex_smoker_plot + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_sex_smoked.pdf",model_id),
         height = 1.3,width = 2) 
```

```{r unknown_cause_age_association,fig.width=2,fig.height=1.3}
ggsave(plot = uc_age_plot,useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_age.pdf",model_id),
       height = 1.3,width = 2) 

ggsave(plot = uc_vaf_plot,useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_vs_vaf.pdf",model_id),
       height = 1.3,width = 2) 

ggsave(plot = uc_vaf_last_plot,sprintf("figures/%s/dynamic_coefficients/uc_vs_vaf_last.png",model_id,dpi = 300),
       height = 1.3,width = 2) 

uc_age_plot
uc_vaf_plot
uc_vaf_last_plot

tmp_ <-  tmp_2 %>% 
  subset(vaf > 0 & vaf_last > 0) %>%
  mutate(log_vaf = log(vaf),
         log_vaf_last = log(vaf_last),
         log_vaf_median = log(vaf_median),
         log_vaf_strict = log(vaf_strict),
         D = vaf_last / vaf)

cor.test(tmp_$log_vaf,tmp_$b_clone) 

cor.test(tmp_$log_vaf_last,tmp_$b_clone)

glm(b_clone ~ log_vaf + log_vaf : genes,data = tmp_) %>%
  summary

tmp_ %>% 
  ggplot(aes(x = vaf, y = b_clone,colour = genes)) + 
  geom_point(size = 0.25) + 
  theme_gerstung(base_size = 6) + 
  scale_colour_manual(values = gene_colours,guide = F) + 
  geom_smooth(method = 'lm',colour = 'black',formula = y ~ x,
              size = 0.1) + 
  scale_x_continuous(trans = 'log10')
```

```{r uc_age_pearson,fig.width=1.3,fig.height=1}
pearson_correlations %>%
  select(Gene = gene,`Spearman's rho` = rho,`p-value` = p.val,`Significant after correction` = Significant) %>%
  mutate(`p-value` = p.adjust(`p-value`,"BH"))
Parameters %>% 
  subset(gene == "TET2") %>%
  mutate(fraction = exp(b_clone_mean)-1) %>%
  ggplot(aes(x = age, y = fraction,colour = gene)) +
  geom_point(size = 0.5) + 
  geom_point(size = 0.5,alpha = 0.01,colour = "black") + 
  geom_smooth(method = "lm",formula = y ~ x,colour = "black",size = 0.25) +
  ylab("TET2 UC effect") +
  xlab("Age at mutation detection") + 
  scale_colour_manual(values = gene_colours,guide = F) + 
  theme_gerstung(base_size = 6) +
  ggtitle(sprintf("Spearman's rho = %.2f\n(p-value = %.2e)",
                  pearson_correlations$rho[8],pearson_correlations$p.val[8])) + 
  scale_y_continuous(breaks = c(-0.1,0,0.1,0.2),
                     labels = sprintf("%s%%",c(-10,0,10,20))) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_age_pearson.pdf",model_id),
         height = 1.3,width = 2) 
```

##### Driver effect vs. unknown cause effect at the clone level

Here we compare how the predicted driver effect - that which is explained by driver genetics alone - with the observed effect - that which is also explained by the UC effect. We can clearly see that growth is mostly driver by driver genetics.

```{r compare_gene_unknown_cause,fig.height=2.4,fig.width=2.4}
plot_min_max_effects <- c(
  min(Parameters$b_clone_005+Parameters$b_genetic_005,Parameters$b_genetic_005),
  max(Parameters$b_clone_095+Parameters$b_genetic_095,Parameters$b_genetic_095)
)
plot_min_max_effects <- exp(plot_min_max_effects) * 100 - 100
Parameters <- Parameters %>%
  mutate(gene_lighter = paste(gene,'lighter',sep = '_'))

ggplot(Parameters,aes(x = exp(b_genetic_mean + b_clone_mean) * 100 - 100,
                      xmin = exp(b_genetic_005 + b_clone_005) * 100 - 100,
                      xmax = exp(b_genetic_095 + b_clone_095) * 100 - 100,
                      y = exp(b_genetic_mean) * 100 - 100,
                      ymin = exp(b_genetic_005) * 100 - 100,
                      ymax = exp(b_genetic_095) * 100 - 100,
                      colour = gene)) + 
  geom_errorbarh(size = 0.25,
                 aes(colour = gene_lighter)) +
  geom_linerange(size = 0.25,
                 aes(colour = gene_lighter)) +
  geom_abline(slope = 1,linetype = "longdash") +
  geom_point(size = 0.5) +
  theme_gerstung(base_size = 6) + 
  scale_y_continuous(breaks = seq(-20,100,by = 20),
                     labels = paste0(seq(-20,100,by = 20),'%'),
                     expand = c(0.05,0,0.0,0)) +
  scale_x_continuous(breaks = seq(-20,100,by = 20),
                     labels = paste0(seq(-20,100,by = 20),'%'),
                     expand = c(0.05,0,0.0,0)) +
  coord_cartesian(xlim = plot_min_max_effects,
                  ylim = plot_min_max_effects) + 
  xlab("Observed growth per year") +
  ylab("Predicted driver growth per year") +
  scale_colour_manual(values = c(gene_colours,gene_colours_lighter),
                      guide = F) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_vs_driver.pdf",model_id),
         height = 1.7,width = 1.7)
```

```{r mutation_independent_per_gene,fig.width=1.3,fig.height=2,dpi=300}
mutation_independent_coefficients <- r_values %>%
  group_by(individual,site,b_clone,genes) %>%
  summarise(age = min(age)) %>% 
  ungroup %>% 
  mutate(genes = str_match(genes,'[A-Z0-9]+'))
mutation_independent_coefficients %>% 
  ggplot(aes(x = factor(genes,levels = rev(gene_order)),
             y = exp(b_clone)*100-100,colour = genes)) + 
  geom_jitter(alpha = 0.5,width = 0.3,height = 0,size = 0.5) +
  geom_boxplot(outlier.colour = NA,fill = NA,
               colour = "black",size = 0.25) +
  theme_gerstung(base_size = 6) + 
  xlab("Driver gene") + 
  ylab("UC effect") + 
  scale_color_manual(values = gene_colours,name = "Gene") + 
  guides(colour=guide_legend(ncol=2,byrow=F)) +
  theme(legend.position = 'none',
        axis.text.y = element_text(face = 'italic'),
        axis.title.y = element_blank()) +
  scale_y_continuous(breaks = seq(-20,80,by = 10),
                     labels = paste0(seq(-20,80,by = 10),'%'),
                     expand = c(0.05,0,0.0,0)) +
  rotate_x_text() +
  coord_flip() +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_gene.pdf",model_id),
         height = 1.85,width = 1.3)
```

```{r mutation_independent_per_gene_distribution,fig.width=1.9,fig.height=1.6,dpi=300}
mutation_independent_coefficients %>% 
  ggplot(aes(x = exp(b_clone)*100-100,colour = genes)) + 
  geom_density(outlier.colour = NA,fill = NA,size = 0.25) +
  theme_gerstung(base_size = 6) + 
  xlab("Unknown cause effect") + 
  ylab("Density") + 
  scale_color_manual(values = gene_colours,name = "Gene") + 
  guides(colour=guide_legend(ncol=2,byrow=F)) +
  theme(legend.position = 'none') +
  scale_x_continuous(breaks = seq(-20,80,by = 10),
                     labels = paste0(seq(-20,80,by = 10),'%'),
                     expand = c(0.05,0,0.0,0)) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/uc_gene_density.pdf",model_id),
         height = 1,width = 1.5)
```

##### Driver effect vs. unknown cause effect at the gene level

An analysis of the unknown cause variance on a per gene level is valuable, but by accounting for the average total effect we can see there appears to be a linear trend for most genes. The exceptions are, very clearly, *JAK2*, whose UC variance is much higher than expected, and *SRSF2*-P95H (notice that *SRSF2* here is stratified here into P95H (*SRSF2*-P95H) and into other variants(*SRSF2*-other)), *PTPN11* and *U2AF1*, whose variance is much lower than expected. Additionally, we can test whether *TET2* and *DNMT3A*, the most frequently occurring genes, have different variances. Indeed, we can observe that *TET2* has variance which is 2.3 times higher than that of *DNMT3A*.

```{r variance_to_mean,fig.width=3.5,fig.height=3}
Parameters %>% 
  subset(paste(individual,site) %in% good_stat_sites) %>% 
  mutate(gene = ifelse(
    gene == "SRSF2",
    ifelse(site == "SRSF2-P95H","SRSF2-P95H","SRSF2-other"),
    gene)) %>% 
  group_by(gene,gene_lighter) %>% 
  mutate(total_effect = exp(b_clone_mean + b_genetic_mean)*100-100,
         b_clone_mean = exp(b_clone_mean)*100 - 100) %>%
  summarise(std = sd(b_clone_mean),
            m = mean(total_effect),
            N = length(unique(paste(site,individual))),
            m05 = quantile(total_effect,0.1),
            m95 = quantile(total_effect,0.9)) %>% 
  subset(N > 3) %>% 
  mutate(std_10 = std*sqrt(N-1)/sqrt(qchisq(1-0.1,N-1)),
         std_90 = std*sqrt(N-1)/sqrt(qchisq(1-0.9,N-1))) %>%
  arrange(-N) %>% 
  ggplot(aes(x = m,xmin = m05,xmax = m95,
             y = std,ymin = std_10,ymax = std_90,
             colour = gene,label = gene,size = N)) + 
  geom_errorbarh(height = 0,
                 size = 0.4,
                 aes(colour = gene_lighter)) +
  geom_linerange(size = 0.4,
                 aes(colour = gene_lighter)) +
  geom_point() + 
  geom_label_repel(size = 2,
                   fontface = "italic",
                   label.size = unit(0,"cm"),
                   label.padding = unit(0.05,"cm"),
                   label.r = unit(0,"cm"),
                   aes(fill = gene_lighter),
                   colour = 'black',
                   segment.size = 0.25,
                   min.segment.length = 0.2,
                   alpha = 0.8) +
  theme_gerstung(base_size = 6) + 
  scale_x_continuous(limits = c(0,NA),
                     breaks = c(0,10,20,30,40,50),
                     labels = c("0%","10%","20%","30%","40%","50%"),
                     expand = c(0.12,0,0.05,0)) + 
  scale_y_continuous(limits = c(0,NA),
                     breaks = c(0,2,4,6,8,10),
                     labels = c("0%","2%","4%","6%","8%","10%"),
                     expand = c(0.12,0,0.05,0)) + 
  scale_colour_manual(values = c(gene_colours,gene_colours_lighter,
                                 `SRSF2-P95H`=as.character(gene_colours["SRSF2"]),
                                 `SRSF2-other`=as.character(gene_colours["SRSF2"])),
                      guide = F) + 
  scale_fill_manual(values = c(gene_colours,gene_colours_lighter,
                               `SRSF2-P95H`=as.character(gene_colours["SRSF2"]),
                               `SRSF2-other`=as.character(gene_colours["SRSF2"])),
                      guide = F) + 
  scale_size(breaks = c(3,10,100,200),
             labels = c("3 (min)",10,100,200),
             name = "Number of\nindividuals") + 
  theme(legend.position = c(0.5,1),
        legend.justification = c(0.5,1),
        legend.direction = "horizontal",
        legend.key.width = unit(0,"cm")) + 
  xlab("Mean of the total growth effect") + 
  ylab("Standard deviation of the UC effect") +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/var_vs_mean.pdf",model_id),
         height = 3,width = 3.5)

var.test(Parameters$b_clone_mean[Parameters$gene == "TET2"],
         Parameters$b_clone_mean[Parameters$gene == "DNMT3A"])
```

##### Investigating individuals with high or low average unknown cause effect or individuals with a wide range of unknown cause effects

Now that we teased apart most effects which are either technical or biological, we can focus on specific trajectories within each individual. Observing how the unknown cause effects are distributed in each individual with two or more clones also allows us to draw additional conclusions regarding their growth. Immediately apparent is how varied these distributions can be (as evidenced in the first column), with: 

  1. Individuals showing quite a high average unknown cause effect (second column). Here we find mostly individuals with no site-specific effects, making the unknown cause coefficient accountable for the dynamics which are not attributable to the gene (such as the case in individual **11959**), but we also find cases where possibly synergistic effects may be visible, such as in individual **158** (between *DNMT3A*−R635W and *TET2*−R1455fs) or individual **35863** (between *CBL*−c.1228−1G>A and *IDH2*−R140Q competing with *TET2*−I1873T). The very similar UC effects - deviations from the growth explained by driver genetics - in both cases would point towards the co-clonality of these mutations.
  
  2. Individuals showing quite a low average unknown cause effect (third column) - these examples are perhaps the hardest to explain. The most likely explanation for three of these cases - **1904**, **2046** and **3549** we observe the presence of 4 or more mutations, which may point towards an overcroweded clonal landscape. In **1094** there is also evidence for a quickly rising *SRSF2*-P95H clone - indeed, this clone practically suffers no variation in growth dynamics - which may explain the very negative variations observable in other genes.
  
  3. Individuals showing a high range of unknown cause effects (fourth column) - in these cases we observe an overlap with individuals who have a high average unknown cause effect (**2045** and **11959**), which would explain their large range. Other cases can be explained by competition effects: in individual **28831** a very fast clone (*PPM1D*−T483fs) appears to be dominating the clonal landscape and decelerating the *ASXL1*−Q546X clone, in individual **35863** a *IDH2*−R140Q clone is growing quite quickly and leading to the deceleration and decrease of the *TET2*−I1873T clone, and in individual **5201** two very fast growing clones - *TET2*−I274fs and *TET2*−L1332P - are quickly overpowering the clonal landscape and leading to the deceleration of *TET2*−S1481fs, a clone whose magnitude (26% VAF at the last timepoint) was likely due to the absence of competition. 

```{r mutation_independent_effects}
single_occurring_aa <- full_data$amino_acid_change[!full_data$single_occurring] %>%
  unique
mutation_independent_individual_data <- r_values %>%
  group_by(individual) %>% 
  filter(any((site %in% single_occurring_aa))) %>% 
  select(individual,age,coefficient,b_clone) %>%
  group_by(individual,coefficient,b_clone) %>%
  summarise(age = max(age)) %>%
  distinct() %>%
  gather(key = 'key',value = 'value',-individual,-age) %>%
  group_by(key,individual) %>%
  summarise(m=min(value),med=median(value),M=max(value),N = length(unique(value))) %>% 
  mutate(Range = M - m) %>% 
  subset(N > 1) %>% 
  subset(key != 'coefficient') %>%
  mutate(individual = factor(individual,ordered = T)) %>%
  mutate(individual = reorder(individual,med)) %>% 
  mutate(numeric_individual = as.numeric(individual))

mutation_independent_variation_plots <- list(
  top = list(
    data = mutation_independent_individual_data %>%
      arrange(med) %>%
      tail(5) %>%
      ungroup() %>%
      select(individual) %>%
      unlist %>% 
      rev()),
  iqr = list(
    data = mutation_independent_individual_data %>%
      arrange(Range) %>%
      tail(5) %>%
      ungroup() %>%
      select(individual) %>%
      unlist %>% 
      rev()),
  bottom = list(
    data = mutation_independent_individual_data %>%
      arrange(-med) %>%
      tail(5) %>%
      ungroup() %>%
      select(individual) %>%
      unlist %>% 
      rev())
)

plotted_individuals <- mutation_independent_variation_plots %>%
  lapply(function(x) as.character(x$data)) %>% 
  do.call(what = rbind)

mutation_independent_individual_data <- mutation_independent_individual_data %>%
  arrange(individual) %>% 
  mutate(plotted = individual %in% plotted_individuals) %>%
  mutate(plot_individual = ifelse(
    plotted == T | c(T,plotted[seq(1,length(plotted))-1])==T,0.6,0.4) %>% cumsum)
plot_breaks_labels <- mutation_independent_individual_data %>%
  transmute(breaks = plot_individual,labels = ifelse(plotted,as.character(individual),""))

mutation_independent_rectangles <- list(
  top = mutation_independent_individual_data %>%
    filter(individual %in% mutation_independent_variation_plots$top$data) %>%
    arrange(individual) %>%
    summarise(x = plot_individual[3],
              y = (max(M) + min(m))/2,
              height = max(M) - min(m)),
  bottom = mutation_independent_individual_data %>%
    filter(individual %in% mutation_independent_variation_plots$bottom$data) %>%
    arrange(individual) %>%
    summarise(x = plot_individual[3],
              y = (max(M) + min(m))/2,
              height = max(M) - min(m)),
  range = mutation_independent_individual_data %>% 
    filter(individual %in% mutation_independent_variation_plots$iqr$data) %>%
    transmute(x = plot_individual,height = Range,y = (M+m)/2)
)

color_code_plots <- list(top = "hotpink3",
                         bottom = "lightskyblue",
                         iqr = "orange")

mutation_independent_individual_dist <- mutation_independent_individual_data %>% 
  ggplot(aes(x = plot_individual,y = med,ymin=m,ymax=M,
             alpha = plotted,
             color = med > 0)) +
  geom_linerange(size = 0.25) + 
  geom_point(size = 0.5) +
  geom_hline(yintercept = 0,size = 0.25,colour = "white") +
  geom_tile(data = mutation_independent_rectangles$range,
            aes(x = x,y = y,height = height + 0.02, width = 0.3),
            fill = NA,
            colour = color_code_plots$iqr,
            linetype = 1,
            size = 0.3,
            inherit.aes = F) +
  geom_tile(data = mutation_independent_rectangles$top,
            aes(x = x,y = y,height = height + 0.03, width = 2.8),
            fill = NA,
            colour = color_code_plots$top,
            linetype = 1,
            size = 0.3,
            inherit.aes = F) +
  geom_tile(data = mutation_independent_rectangles$bottom,
            aes(x = x,y = y,height = height + 0.02, width = 2.8),
            fill = NA,
            colour = color_code_plots$bottom,
            linetype = 1,
            size = 0.3,
            inherit.aes = F) +
  theme_gerstung(base_size = 6) + 
  rotate_x_text() + 
  scale_color_aaas(guide = F) + 
  scale_x_continuous(breaks = plot_breaks_labels$breaks,
                     labels = plot_breaks_labels$labels,
                     expand = c(0.003,0)) +
  scale_alpha_discrete(guide = F,range = c(0.3,1)) +
  ylab("UC growth effect\ndistribution") + 
  xlab("id")
```

```{r divergent_trajectories, fig.width=6,fig.height=6,warning=F}
plot_titles <- list(top = "high average UC",
                    iqr = "high UC range",
                    bottom = "low average UC")
W <- 1
formatted_data_variation <- r_values %>% 
  mutate(site = gsub("t|nt","",site)) %>% 
  group_by(site) %>% 
  mutate(site = str_split(site,';')[[1]][1]) %>% 
  mutate(site = ifelse(grepl(':',site),
                       paste(str_match(gene,'[A-Z0-9]+'),str_split(site,':')[[1]][3],sep='-'),
                       site)) %>%
  select(site,individual,pred,coverage,pred_005,pred_095,
         age,site_numeric,gene_numeric,clone_numeric,b_clone)

mi <- values_model$min_age
for (element in names(mutation_independent_variation_plots)) {
  mutation_independent_variation_plots[[element]]$plots <- list()
  for (ind in mutation_independent_variation_plots[[element]]$data) {
    tmp <- formatted_data_variation %>% 
      subset(individual == ind) %>%
      mutate(plot = ifelse(age == max(age),b_clone,NA)) %>%
      mutate(pred = ifelse(pred==0,pred + 0.01,pred),
             pred_005 = ifelse(pred_005==0,pred_005+0.01,pred_005),
             pred_095 = ifelse(pred_095==0,pred_095+0.01,pred_095)) # to avoid warning messages
    age_ranges <- seq(min(tmp$age),max(tmp$age),length.out = 20)
    trajectories <- tmp %>%
      select(site,site_numeric,gene_numeric,clone_numeric,b_clone) %>%
      distinct %>%
      apply(1,function(x) {
        S <- grab_samples(as.numeric(x[2]),as.numeric(x[3]),as.numeric(x[4]),2500)
        b <- S$site + S$clone + S$gene
        u <- S$offset
        s <- age_ranges %>%
          lapply(
            function(t) {
              s <- trajectory_from_parameters_unadjusted(b,u,t-mi)
              o <- c(mean(s,na.rm=T),quantile(s,c(0.05,0.95),na.rm=T))
              names(o) <- c("trajectory_mean","trajectory_005","trajectory_095")
              data.frame(
                trajectory_mean = o[1],
                trajectory_005 = o[2],
                trajectory_095 = o[3],
                site = as.character(x[1]),
                b_clone = as.numeric(x[5]),
                age = t
              ) %>%
              return()
            }
          ) %>%
          do.call(what = rbind)
        return(s) 
      }) %>%
      do.call(what = rbind) %>%
      mutate(site = reorder(site,trajectory_mean))
    tmp_2 <- trajectories %>%
      filter(age == max(age)) %>% 
      select(site,b_clone,age,trajectory_mean)
    pos <- c(0,max(tmp$pred_095/tmp$coverage)*0.95)
    mutation_independent_variation_plots[[element]]$plots[[as.character(ind)]] <- tmp %>% 
      ggplot(aes(x = age,
                 group = site)) + 
      geom_ribbon(
        data = trajectories,
        inherit.aes = F,
        aes(fill = site,x = age,y = trajectory_mean,
            ymin = trajectory_005,ymax = trajectory_095),size = 0.25,alpha = 0.5,
        position = position_dodge(width = W)) +
      geom_line(
        data = trajectories,
        inherit.aes = F,
        aes(colour = site,x = age,y = trajectory_mean,group = site),
        size = 0.25,alpha = 0.3,
        position = position_dodge(width = W)) +
      geom_text_repel(data = tmp_2,
                       position = position_dodge(width = W),
                       hjust = -0.8,
                       aes(x = age,
                           y = trajectory_mean,
                           label = sprintf('%s%%',100*round(exp(b_clone)-1,3)),
                           colour = site),
                       direction = "y",
                       segment.size = 0.25,
                       min.segment.length = 0.05,
                       force = 1,
                       show.legend = F,
                       inherit.aes = F,
                       size = 2,
                       segment.colour = "black",
                       alpha=0.9) + 
      scale_shape(guide=F) + 
      ggtitle(sprintf("id=%s (%s)",ind,plot_titles[[element]])) + 
      scale_color_aaas(name = NULL) + 
      scale_fill_aaas(name = NULL) + 
      theme_gerstung(base_size = 6) + 
      xlab("Age") + 
      ylab("VAF") + 
      theme(legend.position = "bottom",
            legend.key.size = unit(0.05,"cm"),
            legend.box.spacing = unit(0.01,"cm"),
            plot.title = element_text(face = 'bold',color = color_code_plots[[element]])) +
      scale_y_continuous(trans = 'log10',
                         labels = function(x) sprintf("%s%%",x*100),
                         expand = c(0,0,0.01,0.01)) +
      coord_cartesian(ylim = c(0.001,pos[2]/0.95)) +
      scale_x_continuous(expand = c(0.05,0,0.4,0),
                         breaks = seq(ceiling(min(tmp$age)),floor(max(tmp$age)),3)) +
      guides(colour=guide_legend(nrow=3))
  }
}

plot_grid(
  mutation_independent_individual_dist + coord_flip(),
  plot_grid(
    plot_grid(plotlist = mutation_independent_variation_plots$top$plots,ncol = 1),
    plot_grid(plotlist = mutation_independent_variation_plots$bottom$plots,ncol = 1),
    plot_grid(plotlist = mutation_independent_variation_plots$iqr$plots,ncol = 1),
    nrow = 1),
  rel_widths = c(0.3,1.0),
  nrow = 1) + 
  ggsave(useDingbats=FALSE,
         filename = sprintf("figures/%s/inferred_trajectories/divergent_trajectories.pdf",model_id),
         width = 6,height = 5)
```

### Inferring ages at onset

The ages at clone foundation are determined using the inferred median values for each parameter for each individual. To do this, we solve $t_0 = \frac{log(\frac{1}{n}) - \alpha + log(\frac{g}{\beta}) - 1}{\beta}$, with $n = 50,000$, $g = 2$, $\beta=growth\ coefficient$ and $alpha = individual\ offset$. We exclude from analysis trajectories with outliers and negative growth coefficients. We use the posterior samples to calculate the age at onset estimates and cap our estimates between 0 and the first point at which the clone was detected. What we observe quite clearly is a clear relationship between the observed growth per year and the age at onset for different clones - indeed, slower clones have a tendency to show up much earlier in life than faster clones. 

```{r age_at_clone_foundation_intervals,fig.height=3.5,fig.width=4.7,warning=F}
c_names <- colnames(values_model$draws$`11`)

u_c_names <- grep('^u',c_names,perl = T)
b_gene_c_names <- grep('^b_gene',c_names,perl = T) 
b_site_c_names <- grep('^b_site',c_names,perl = T)
b_clone_c_names <- grep('^b_clone',c_names,perl = T)

maximum_probability_value <- function(x) {
  d <- density(x)
  return(d$x[which.max(d$y)])
}

sample_ages <- function(sampling_idxs,pop_size,gen_time,sample_size=1000) {
  sampling_idxs %>% 
    apply(1,function(x) {
      site_numeric <- as.numeric(x[3])
      gene_numeric <- as.numeric(x[4])
      clone_numeric <- as.numeric(x[5])
  
      S <- grab_samples(site_numeric,gene_numeric,clone_numeric,sample_size)
      b_site <- S[[1]]
      b_clone <- S[[2]]
      b_gene <- S[[3]]
      u <- S[[4]]
      
      t1 <- as.numeric(x[6])
      age_distribution <- t0_adjusted(u,b_site+b_gene+b_clone,
                                      gen_time,pop_size*2) + values_model$min_age
  
      tmp <- data.frame(
        age_distribution = age_distribution,
        effect_distribution = b_clone + b_site + b_gene,
        offset_distribution = u,
        site = x[2],
        individual = gsub(" ","",x[1]),
        gene = str_match(x[2],'[A-Z0-9]+'),
        t1 = t1
      ) 
      colnames(tmp) <- c("age_distribution","effect_distribution","offset_distribution","site","individual","gene","t1")
      return(tmp)
      }) %>%
    return
}

characterise_age_samples <- function(age_samples) {
  age_samples %>% 
    lapply(function(x) {
      A <- x$age_distribution
      t1 <- x$t1[1]
      A <- ifelse(A < -1, -1, A)
      A <- ifelse(A > t1, t1, A)
      A <- A[!is.na(A)]
      M <- mean(A)
      Qs <- quantile(A,c(0.05,0.10,0.16,0.5,0.84,0.90,0.95))
      if (length(A) == 0) {
        MAP <- NA
      } else {
        D <- density(A)
        MAP <- D$x[which.max(D$y)]
      }
      return(
        c(Mean = M,Qs,MAP = MAP)
      )
    }) %>% 
    do.call(what = rbind) %>% 
    as_tibble() %>%
    return
}

sampling_idxs <- r_values_full %>%
  merge(select(statistic_full,site,sum_stat,individual)) %>% 
  group_by(
    individual,site,
    site_numeric,
    gene_numeric,
    clone_numeric,
    sum_stat
  ) %>% 
  summarise(
    minAge = min(age[which(round(true/coverage,4) >= 0.005)]),
    maxAge = max(age),
    maxVAF = max(true/coverage),
    minVAF = min((true/coverage)[which(round(true/coverage,4) >= 0.005)])
  ) %>%
  arrange(site) %>%
  select(individual,site,site_numeric,gene_numeric,clone_numeric,minAge,maxAge,maxVAF,minVAF,sum_stat)

sampling_idxs %>% 
  apply(1,function(x) {
    o <- do.call(cbind,grab_samples(as.numeric(x[3]),as.numeric(x[4]),as.numeric(x[5]),2500)) %>%
      as.data.frame
    colnames(o) <- c("site_effect_samples","gene_effect_samples",
                     "unknown_effect_samples","offset_samples")
    o$SardID <- as.numeric(x[1])
    o$site <- x[2]
    o$gene <- str_match(o$site,'[A-Z0-9]+[nt]*')
    return(as.data.frame(o))
    }) %>%
  do.call(what = rbind) %>%
  saveRDS(file = "data_output/draws_df.RDS")

list_of_age_sampling <- list()
list_of_age_estimates <- list()
for (gen_time in c(1,2,5,10,13,20)) {
  for (pop_size in c(10e3,5e4,1e5,2e5,600e3)) {
    tmp_samples <- sample_ages(sampling_idxs,pop_size,gen_time,2500)
    tmp_estimates <- characterise_age_samples(tmp_samples)
    colnames(tmp_estimates) <- c("Mean","Q05","Q10","Q16","Q50","Q84","Q90","Q95","MAP")
    list_of_age_sampling[[sprintf('%s_%s',pop_size,gen_time)]] <- tmp_samples
    list_of_age_estimates[[sprintf('%s_%s',pop_size,gen_time)]] <- tmp_estimates
  }
}
all_samples_for_age_estimation <- list_of_age_sampling$`50000_2`
age_estimates <- list_of_age_estimates$`50000_2`

lapply(names(list_of_age_estimates),
       function(x) {
           tmp_str <- str_split(x,'_')[[1]]
           pop_size <- as.numeric(tmp_str[1])
           gen_time <- as.numeric(tmp_str[2])
           X <- list_of_age_estimates[[x]]
           X <- table(floor(X$Q50))
           o <- data.frame(
               X,
               pop_size,
               gen_time
           )
           colnames(o) <- c("Q50","Count","pop_size","gen_time")
           return(o)
       }) %>% 
  do.call(what = rbind) %>% 
  as_tibble() %>% 
  group_by(gen_time,pop_size) %>% 
  summarise(Freq = sum(Count[as.numeric(as.character(Q50)) < 0])/sum(Count)) %>% 
  ggplot(aes(x = as.factor(pop_size),
             y = reorder(sprintf('%.3f (%s)',1/gen_time,gen_time),gen_time),
             fill = Freq)) + 
  geom_tile() + 
  geom_text(aes(label = round(Freq,4))) + 
  theme_gerstung(base_size = 6) + 
  ylab("tau (gen/year)") + 
  xlab("Population size") + 
  scale_fill_material(reverse = F,
                      name = "No. of clones appearing\nbefore birth") + 
  scale_x_discrete(expand = c(0,0)) + 
  scale_y_discrete(expand = c(0,0)) + 
  theme(legend.position = "bottom") + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/age_at_clone_foundation/Q50.pdf",model_id),
         height = 3.5,width = 3.5)

map_age_clones <- list()
N_examples <- length(all_samples_for_age_estimation)
pdf(height = (floor(N_examples/10) + 1)*2,
    width = 10*2,
    file = sprintf("figures/%s/age_at_clone_foundation/samples.pdf",model_id))
par(mfrow = c(floor(N_examples/10) + 1,10),
    mar = c(2,2,3,1))
for (S in all_samples_for_age_estimation) {
  X <- na.omit(S$age_distribution)
  X <- ifelse(X < -1, -1 ,X)
  if (length(X) > 0) {
    D <- density(X)
    map_age_clones[[length(map_age_clones)+1]] <- data.frame(
      MAP = D$x[which.max(D$y)],
      Individual = S$individual[1],
      Site = S$site[1])
    plot(D,
         main = sprintf('SardID = %s\n%s',S$individual[1],S$site[1]),
         xlim = c(-1,max(full_data$Age)))
  }
}
dev.off()

Parameters_Age <- cbind(
  sampling_idxs %>% as_tibble(),age_estimates
) %>% 
  merge(Parameters_full,by = c("individual","site")) %>%
  mutate(
    b_clone_016 = b_clone_mean - sqrt(b_clone_var),
    b_clone_084 = b_clone_mean + sqrt(b_clone_var),
    b_genetic_016 = b_genetic_mean - sqrt(b_genetic_var),
    b_genetic_084 = b_genetic_mean + sqrt(b_genetic_var)
  ) 

M <- Parameters_Age %>%
  subset(b_genetic_mean + b_clone_mean > 0) %>% 
  subset(b_genetic_mean + b_clone_mean == max(b_genetic_mean + b_clone_mean) | 
           (b_genetic_mean + b_clone_mean) == min(b_genetic_mean + b_clone_mean,1)) %>%
  arrange(- b_genetic_mean - b_clone_mean)

all_samples_for_age_estimation_gene <- lapply(
  unique(str_match(levels(values_model$sub_data$Gene),'[0-9A-Z]+')),
  function(x) list())
names(all_samples_for_age_estimation_gene) <- levels(values_model$sub_data$Gene) %>%
  str_match('[0-9A-Z]+') %>% 
  unique
for (S in all_samples_for_age_estimation) {
  G <- as.character(S$gene[1])
  idx <- length(all_samples_for_age_estimation_gene[[G]]) + 1
  all_samples_for_age_estimation_gene[[G]][[idx]] <- S
}
pdf(height = (floor(length(all_samples_for_age_estimation_gene)/5) + 1)*2,
    width = 10*2,
    file = sprintf("figures/%s/age_at_clone_foundation/samples_gene.pdf",model_id))
par(mfrow = c(floor(length(all_samples_for_age_estimation_gene)/5) + 1,5),
    mar = c(2,2,3,1))
for (S in all_samples_for_age_estimation_gene) {
  S <- do.call(rbind,S)
  X <- na.omit(S$age_distribution)
  X <- ifelse(X < -1, -1 ,X)
  if (length(X) > 0) { 
    plot(density(X),
         main = S$gene[1],
         xlim = c(-1,max(full_data$Age)))
  }
}
dev.off()

pdf(height = (floor(length(all_samples_for_age_estimation_gene)/5) + 1)*2,
    width = 5*2,
    file = sprintf("figures/%s/age_at_clone_foundation/samples_gene_above_minus_1.pdf",model_id))
par(mfrow = c(floor(length(all_samples_for_age_estimation_gene)/5) + 1,5),
    mar = c(2,2,3,1))
for (S in all_samples_for_age_estimation_gene) {
  S <- do.call(rbind,S)
  Samp <- S$age_distribution
  Samp <- Samp[Samp > -1]
  if (length(Samp) > 0) {
    plot(density(na.omit(Samp)),
         main = S$gene[1],
         xlim = c(-1,max(full_data$Age)))
  }
}
dev.off()

N <- 0.5e5
detection_threshold <- 0.002
age_at_detection <- 80
fitness <- seq(0,max(Parameters_Age$b_genetic_mean + Parameters_Age$b_clone_mean)+0.1,length.out = 1000)
detection_limits <- data.frame(
  t0 = sapply(fitness,function(x) onset_from_detection(
    b = x,vaf = detection_threshold,
    age_at_detection = age_at_detection,
    g = 2,N = N)),
  fitness = exp(fitness)*100-100)
detection_limits <- rbind(
  detection_limits,
  data.frame(t0 = c(max(detection_limits$t0,na.rm=T),200,200,-1),
             fitness = c(max(detection_limits$fitness,na.rm = T),200,-40,-40))
) 

median_age_at_detection <- median(Parameters_Age$minAge)
age_threshold_accounting <- Parameters_Age %>%
  subset(sum_stat) %>%
  rowwise() %>%
  transmute(
    limit = onset_from_detection(
      b_clone_mean+b_genetic_mean,minVAF,minAge,2,N),
    age_at_onset = Q50,
    coef = b_genetic_mean + b_clone_mean,
    gene = gene
  ) %>%
  subset(coef > 0) %>% 
  mutate(plausible = limit >= 0)

table(age_threshold_accounting$plausible) / nrow(age_threshold_accounting)

Parameters_Age %>%
  subset(b_clone_mean + b_genetic_mean > 0 & sum_stat == T) %>%
  mutate(gene_lighter = paste(gene,'lighter',sep = '_')) %>% 
  ggplot(aes(x = Q50,
             xmin = Q10,xmax = Q90,
             y = exp(b_genetic_mean + b_clone_mean)*100-100,
             ymin = exp(b_genetic_005 + b_clone_005)*100-100,
             ymax = exp(b_genetic_095 + b_clone_095)*100-100,
             colour = gene)) +
  geom_polygon(data = detection_limits,
               inherit.aes = F,
               fill = "grey95",
               aes(x = t0,y = fitness)) +
  geom_linerange(size = 0.25,aes(colour = gene_lighter)) + 
  geom_errorbarh(size = 0.25,aes(colour = gene_lighter)) +
  geom_point(size = 0.5) + 
  geom_line(data = detection_limits[1:(nrow(detection_limits)-4),],
            aes(x = t0,y = fitness),
            size = 0.25,
            inherit.aes = F) + 
  scale_y_continuous(
    breaks = c(-20,0,20,40,60,80,100,120),
    labels = paste0(c(-20,0,20,40,60,80,100,120),'%')
  ) +
  scale_x_continuous(breaks = c(0,20,40,60,80,100),
                     expand = c(0.01,0)) +
  coord_cartesian(xlim = c(0,85),
                  ylim = c(-15,80)) +
  ylab("Observed growth per year") +
  xlab("Age at onset") +
  theme_gerstung(base_size = 6) + 
  scale_color_manual(values = c(gene_colours,gene_colours_lighter),
                     guide = F) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/age_at_clone_foundation/growth_vs_age_at_onset.pdf",model_id),
         width = 3.2,height = 2.3)

Parameters_Age %>% 
  select(
    SardID = individual,
    Site = site,
    Gene = gene,
    CloneID = clone_numeric,
    AgeAtStudyEntry = minAge,
    maxVAF = maxVAF,
    EstimatedAgeAtOnset_Mean = Mean,
    EstimatedAgeAtOnset_Q05 = Q05,
    EstimatedAgeAtOnset_Q10 = Q10,
    EstimatedAgeAtOnset_Q16 = Q16,
    EstimatedAgeAtOnset_Q50 = Q50,
    EstimatedAgeAtOnset_Q84 = Q84,
    EstimatedAgeAtOnset_Q90 = Q90,
    EstimatedAgeAtOnset_Q95 = Q95,
    UnknownCauseEffect_mean = b_clone_mean,
    UnknownCauseEffect_Var = b_clone_var,
    UnknownCauseEffect_Q05 = b_clone_005,
    UnknownCauseEffect_Q95 = b_clone_095,
    GeneEffect_mean = b_gene_mean,
    GeneEffect_Var = b_gene_var,
    GeneEffect_Q05 = b_gene_005,
    GeneEffect_Q95 = b_gene_095,
    SiteEffect_mean = b_site_mean,
    SiteEffect_Var = b_site_var,
    SiteEffect_Q05 = b_site_005,
    SiteEffect_Q95 = b_site_095,
    GeneticEffect_mean = b_genetic_mean,
    GeneticEffect_Var = b_genetic_var,
    GeneticEffect_Q05 = b_genetic_005,
    GeneticEffect_Q95 = b_genetic_095,
    Offset_Q05 = u_values_005,
    Offset_Q95 = u_values_095,
    GoodTrajectory = sum_stat
    ) %>%
  mutate(FullEffectMean = GeneEffect_mean + UnknownCauseEffect_mean + SiteEffect_mean) %>%
  write_tsv("data_output/ParametersAgeAtOnset.tsv",quote = F)
```

Here we have perhaps one of the most striking results - while a few clones show up consistently through life, others - namely those in *U2AF1* and *PTPN11* have a clearly distinct signature regarding the age at onset of clones harbouring mutations in these genes. Indeed, these clones appear only much later in life and while for *PTPN11* this begins at around 20 years of age, for *U2AF1* the first clones appear only at 40. Of particular interest as well is *SRSF2*-P95H, which shows a distinct age signature that is remarkably different not only from other genes but also from other sites in *SRSF2*, with mutations appearing only much later in life.

```{r ages_density_bar,fig.height=2.3,fig.width=2.5}
valid_trajectories <- sampling_idxs %>% 
  ungroup %>% 
  subset(sum_stat == T) %>% 
  mutate(combination = paste(site,individual)) %>% 
  select(combination) %>% 
  unlist

estimates_from_samples <- all_samples_for_age_estimation_gene %>%
  lapply(function(x) {
    X <- lapply(x,function(y) {
      if (paste(y$site[1],y$individual[1]) %in% valid_trajectories) {
        t1 <- y$t1[1]
        X <- y$age_distribution
        X <- X[!is.na(X)]
        X <- ifelse(X < -1,-1,X)
        X <- ifelse(X > t1,t1,X)
        return(X)
      } 
    }) %>% 
      unlist %>%
      as.vector
    gene <- x[[1]]$gene[1]
    data.frame(
      gene = gene,
      values = X,
      row.names = NULL,
      N = length(x)
    ) %>% 
      return
    }) %>% 
  do.call(what = rbind)

stratified_srsf2_samples <- all_samples_for_age_estimation %>%
  lapply(function(x) {
    if (x$gene[1] == "SRSF2" & paste(x$site[1],x$individual[1]) %in% valid_trajectories) {
      x %>%
        mutate(age_distribution = ifelse(age_distribution > t1,t1,age_distribution)) %>%
        return()
    } else return(NULL)
    }) %>%
  do.call(what = rbind) %>%
  as.data.frame %>%
  mutate(strata = ifelse(site == "SRSF2-P95H","SRSF2-P95H","SRSF2-other")) %>%
  mutate(age_distribution = ifelse(age_distribution < -1,-1,age_distribution)) %>% 
  group_by(strata) %>% 
  mutate(N = length(unique(paste(individual,site)))) %>%
  ungroup %>%
  select(gene = strata,values = age_distribution,N)

estimates_from_samples <- estimates_from_samples %>%
  subset(gene != "SRSF2") %>%
  rbind(stratified_srsf2_samples)

estimates_from_samples_parameters <- estimates_from_samples %>%
  group_by(gene,N) %>%
  summarise(Q05 = quantile(values,0.05,na.rm=T),
            Q50 = quantile(values,0.5,na.rm=T),
            Q95 = quantile(values,0.95,na.rm=T),
            LowMass = sum(values <= 0,na.rm = T)/length(values))

gene_order <- estimates_from_samples_parameters %>%
  arrange(-LowMass) %>%
  select(gene) %>%
  unlist

max_growth_fitness_genes <- Parameters_Age %>%
  mutate(gene = ifelse(gene == "SRSF2",
                       ifelse(site == "SRSF2-P95H",
                              "SRSF2-P95H",
                              "SRSF2-other"),
                       gene)) %>% 
  group_by(gene) %>%
  summarise(Growth = exp(mean(b_gene_mean + b_clone_mean + b_site_mean)) - 1) %>%
  rowwise() %>%
  mutate(MaxAgeAtDetection = c(seq(0,110,by = 0.1)[
    which(diff(trajectory_from_t0(age_at_detection,seq(0,110,by = 0.1),
                                  Growth,N=0.5e5,g = 2) > (0.5e5 * detection_threshold)) != 0)])) %>% 
  mutate(AgeAtSampling = age_at_detection) %>%
  merge(estimates_from_samples_parameters,by = "gene") %>%
  mutate(MaxAgeAtDetection_plot = ifelse(MaxAgeAtDetection > Q95,
                                         Q95,
                                         MaxAgeAtDetection))

max_growth_fitness_sites <- Parameters_Age %>%
  group_by(gene,site,individual,Q50) %>%
  summarise(Growth = exp(mean(b_gene_mean + b_clone_mean + b_site_mean)) - 1) %>%
  rowwise() %>%
  apply(1,function(x) {
    A <- seq(0,110,by = 0.1)
    idx <- which(diff(
      trajectory_from_t0(age_at_detection,A,
                         as.numeric(as.character(x[5])),
                         N=0.5e5,g = 2) > (0.5e5 * detection_threshold)) != 0)
    O <- A[idx]
    if (length(O) == 1) {
      return(c(x,MaxAgeAtOnset = O))
      } else {
        return(c(x,MaxAgeAtOnset = NA))
      }
  }) %>% t %>%
  as.data.frame %>%
  mutate(Growth = as.numeric(as.character(Growth)),
         MaxAgeAtOnset = as.numeric(as.character(MaxAgeAtOnset)),
         AgeAtOnset = as.numeric(as.character(Q50))) %>%
  mutate(TimeToDetection = 80 - MaxAgeAtOnset)

max_growth_fitness <- max_growth_fitness_sites %>%
  mutate(gene = ifelse(gene == "SRSF2",
                       ifelse(site == "SRSF2-P95H",
                              "SRSF2-P95H",
                              "SRSF2-other"),
                       as.character(gene))) %>% 
  group_by(gene) %>%
  summarise(Threshold = median(AgeAtOnset+TimeToDetection,na.rm=T))

max_growth_fitness_sites %>%
  group_by(gene) %>%
  summarise(AverageTimeToDetectionByGene = median(TimeToDetection,na.rm=T)) %>%
  arrange(AverageTimeToDetectionByGene)

max_growth_fitness_sites %>%
  summarise(AverageTimeToDetectionAllSites = median(TimeToDetection,na.rm=T))
max_growth_fitness_sites %>%
  subset(gene == "SRSF2") %>%
  mutate(isP95H = ifelse(site == "SRSF2-P95H","P95H","Other mutations")) %>%
  group_by(isP95H) %>%
  summarise(AverageTimeToDetection = median(TimeToDetection,na.rm=T))

for (x in tmp$MaxGrowth) {
  which(diff(trajectory_from_t0(age_at_detection,seq(-10,2000,by = 0.1),tmp$MaxGrowth[1],N=0.5e5,g = 2) > (0.5e5 * detection_threshold)) != 0)
}

labs <- estimates_from_samples_parameters %>%
  mutate(gene = factor(gene,levels = gene_order)) %>%
  arrange(gene) %>% 
  apply(
  1,
  function(x) parse(text = sprintf('italic("%s")~"(n=%d)"',x[1],as.numeric(x[2])))
) %>%
  do.call(what = c)

estimates_from_samples %>% 
  subset(!is.na(values)) %>%
  group_by(gene) %>%
  filter(values >= quantile(values,0.025) & values <= quantile(values,0.975)) %>%
  ggplot(aes(x = factor(gene,levels = gene_order),
             y = values,
             fill = gene)) +
  geom_violin(#scale = "width", 
              size = 0) +
  geom_tile(aes(y = Threshold,
                x = gene,
                width = 0.8,height = 1),
            inherit.aes = F,
            data = max_growth_fitness,
            colour = "grey10",
            size = 0.25,
            fill = "white") + 
  geom_point(colour = "black",
             alpha = 0.5,
             shape = 16,
             size = 0.3,
             position = position_jitter(width = 0.2),
             data = mutate(subset(Parameters_Age,sum_stat),
                           gene = ifelse(
                             gene == "SRSF2",
                             ifelse(site == "SRSF2-P95H","SRSF2-P95H","SRSF2-other"),
                             gene)),
             aes(x = gene,y = Q50)) +
  theme_gerstung(base_size = 6) + 
  coord_flip() + 
  ylab("Projected age at onset") + 
  xlab("Gene") + 
  theme(axis.text.y = element_text(face = "italic")) + 
  scale_fill_manual(values = c(gene_colours,
                               `SRSF2-P95H` = as.character(gene_colours["SRSF2"]),
                               `SRSF2-other` = as.character(gene_colours["SRSF2"])),
                    guide = F) + 
  scale_y_continuous(#limits = c(0,90),
                     breaks = c(0,20,40,60,80),expand = c(0,0,0.01,0)) + 
  scale_x_discrete(labels = labs,breaks = gene_order) + 
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/age_at_onset_violin.pdf",model_id),width = 2.5,height = 2.3)

estimates_from_samples %>% 
  subset(!is.na(values)) %>%
  group_by(gene) %>%
  summarise(Q025 = quantile(values,0.025),
            Q975 = quantile(values,0.975)) %>%
  ggplot(aes(x = factor(gene,levels = gene_order),
             y = (Q975 + Q025)/2,
             width = 0.8,
             height = Q975 - Q025,
             fill = gene)) +
  geom_tile() +
  geom_tile(aes(y = Threshold,
                x = gene,
                width = 0.8,height = 1),
            inherit.aes = F,
            data = max_growth_fitness,
            colour = "grey10",
            size = 0.25,
            fill = "white") + 
  geom_point(colour = "white",
             inherit.aes = F,
             alpha = 0.5,
             shape = 16,
             size = 0.3,
             position = position_jitter(width = 0.2,height = 0),
             data = mutate(subset(Parameters_Age,sum_stat),
                           gene = ifelse(
                             gene == "SRSF2",
                             ifelse(site == "SRSF2-P95H","SRSF2-P95H","SRSF2-other"),
                             gene)),
             aes(x = gene,y = Q50)) +
  theme_gerstung(base_size = 6) + 
  coord_flip() + 
  ylab("Corrected age at onset") + 
  xlab("Gene") + 
  theme(axis.text.y = element_text(face = "italic")) + 
  scale_fill_manual(values = c(gene_colours,
                               `SRSF2-P95H` = as.character(gene_colours["SRSF2"]),
                               `SRSF2-other` = as.character(gene_colours["SRSF2"])),
                    guide = F) + 
  scale_y_continuous(#limits = c(0,90),
                     breaks = c(0,20,40,60,80),expand = c(0,0,0.01,0)) + 
  scale_x_discrete(labels = labs,breaks = gene_order) + 
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/age_at_onset_bars.pdf",model_id),
         width = 2.5,height = 2.3) +
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/age_at_onset_bars_smaller.pdf",model_id),
         width = 1.9,height = 2)
```

##### Clone size vs. growth rate

```{r clone_size_vs_growth_rate, fig.height=1.,fig.width=1.3}
B <- function(n,f) {
  return(log(n) + log(f))
}

expected <- data.frame(
  x = seq(0.005,0.5,by = 0.005))
expected$y_50e3 <- B(100e3,expected$x)

Parameters_Age %>%
  subset(!is.na(Q50) & sum_stat == T) %>% 
  ggplot(aes(x = maxVAF,y = maxAge*(exp(b_genetic_mean + b_clone_mean)-1),
             colour = Q50 <= -1)) + 
  geom_point(size = 0.25) + 
  geom_line(aes(x = x,y = y_50e3),size = 0.5,data = expected,inherit.aes = F) +
  theme_gerstung(base_size = 6) + 
  scale_colour_manual(values = c("lightblue","goldenrod"),labels = c("Feasible","Unfeasible"),guide = F) + 
  theme(legend.position = "bottom",legend.key.size = unit(0,"cm"),
        legend.box.spacing = unit(0,"cm")) + 
  xlab("VAF") + 
  ylab("Age*annual growth") + 
  scale_x_continuous(trans = 'log10') +
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/size_vs_growth_exp.pdf",model_id),width = 1.3,height = 1.0)
```

#### Assessing the technical determinants of the age at onset

Here we inspect how different values for the number of generations per year and population size affects our interval estimates and begin by seeing the fraction of clones which cannot be explained.

```{r,fig.height=2,fig.width=5}
dff <- names(list_of_age_estimates) %>% 
  lapply(function(x) {
    population_gen_time <- strsplit(x,split = "_")
    population <- as.numeric(population_gen_time[[1]][1])
    gen_time <- as.numeric(population_gen_time[[1]][2])
    O <- cbind(list_of_age_estimates[[x]],as.data.frame(sampling_idxs),
               population = population,gen_time = gen_time) 
    return(merge(O,Parameters,by = c("individual","site")))
    }) %>%
  do.call(what = rbind) %>%
  subset(sum_stat == T) 

dff <- dff %>%
  mutate(age_at_onset_unadjusted = t0(alpha = (u_values_095 + u_values_005)/2,
                                      beta = b_clone_mean + b_genetic_mean,
                                      n = population))

dff %>%
  group_by(population,gen_time) %>%
  summarise(UnfeasibleCorrected = sum(Q50 <= -1),
            Unfeasible = sum(age_at_onset_unadjusted <= -1),
            N = sum(!is.na(Q50)),
            N_ = sum(!is.na(age_at_onset_unadjusted))) %>%
  mutate(p05 = qbeta(0.05,UnfeasibleCorrected,N-UnfeasibleCorrected),
         p95 = qbeta(0.95,UnfeasibleCorrected,N-UnfeasibleCorrected)) %>%
  group_by(population) %>%
  mutate(Unfeasible = c(Unfeasible[1],rep(NA,length(Unfeasible)-1))) %>%
  ggplot(aes(x = as.factor(gen_time),y = UnfeasibleCorrected / N,fill = "A")) + 
  geom_bar(stat = "identity") + 
  geom_bar(aes(x = "E",y = Unfeasible / N_,fill = "B"),stat = "identity") +
  geom_linerange(aes(ymin = qbeta(0.05,Unfeasible,N-Unfeasible),
                     ymax = qbeta(0.95,Unfeasible,N-Unfeasible),x="E")) +
  geom_linerange(aes(ymin = p05,ymax = p95)) +
  facet_wrap(~ population,nrow = 1) + 
  theme_gerstung(base_size = 6) + 
  xlab("Generation time") +
  ylab("Fraction of unfeasible clones") + 
  scale_y_continuous(expand = c(0,0,0.1,0)) +
  scale_x_discrete(
    breaks = c("E","1","2","5","10","13","20"),
    limits = c("E","1","2","5","10","13","20")
  ) + 
  scale_fill_manual(values = c("grey40","grey80"),guide = F) +
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/fraction_of_plausible_clones.pdf",model_id),width = 4.5,height = 2)
```

```{r,fig.width=1.3,fig.height=1.45}
dff %>%
  subset(population == 50e3 & gen_time == 2) %>% 
  mutate(gene = ifelse(grepl("DNMT",site),str_match(as.character(site),"[A-Z0-9]+[tn]+"),
                       str_match(as.character(site),"[A-Z0-9]+"))) %>% 
  group_by(population,gen_time,gene) %>%
  summarise(UnfeasibleCorrected = sum(Q50 <= -1),
            N = sum(!is.na(Q50))) %>% 
  mutate(p05 = qbeta(0.05,UnfeasibleCorrected,N-UnfeasibleCorrected),
         p95 = qbeta(0.95,UnfeasibleCorrected,N-UnfeasibleCorrected)) %>%
  ggplot(aes(x = reorder(gene,UnfeasibleCorrected / N),y = UnfeasibleCorrected / N,fill = gene)) + 
  geom_bar(stat = "identity") + 
  geom_linerange(aes(ymin = p05,ymax = p95),size = 0.25) + 
  coord_flip() + 
  scale_fill_discrete(guide = F) +
  scale_fill_manual(guide = F,values = gene_colours) + 
  theme_gerstung(base_size = 6) + 
  theme(axis.text.y = element_text(face = "italic")) + 
  scale_y_continuous(expand = c(0,0),
                     labels = function(x) return(sprintf("%s%%",x * 100))) + 
  ylab("Fraction of clones\nexceeding lifetime") + 
  xlab("") + 
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/fraction_of_plausible_clones_gene.pdf",model_id),
         width = 1.3,height = 1.6)
```

```{r compare_gen_time_pop_size,fig.width=6,fig.height=6}
all_age_estimates <- list_of_age_estimates %>% 
  lapply(function(x) cbind(as.data.frame(sampling_idxs),as.data.frame(x))) %>%
  do.call(what = rbind) 
pop_gen <- rownames(all_age_estimates) %>%
  str_match("[a-z+0-9]+_[0-9]+") %>% 
  str_split("_") %>% 
  do.call(what = rbind)
all_age_estimates$population <- as.numeric(pop_gen[,1])
all_age_estimates$gen_time <- as.numeric(pop_gen[,2])

all_samples_for_age_estimation_gene_compare <- list()

for (parameters in names(list_of_age_sampling)) {
  x <- list()

  for (S in list_of_age_sampling[[parameters]]) {
    gene <- as.character(S$gene[1])
    if (gene == "SRSF2") {
      if (S$site[1] == "SRSF2-P95H") gene <- "SRSF2-P95H"
      else gene <- "SRSF2-other"
    }
    idx <- length(x[[gene]]) + 1
    x[[gene]][[idx]] <- ifelse(S$age_distribution > S$t1,S$t1,S$age_distribution)
  }
  output <- x %>% lapply(function(y) {
    J <- do.call(c,y)
    J <- J[!is.na(J)]
    J <- ifelse(J < 0,0,J)
    J <- quantile(J,c(0.05,0.5,0.95))
    O <- data.frame(
      Q05 = J[1],
      Q50 = J[2],
      Q95 = J[3],
      parameters = parameters
    )
    return(O)}) %>% 
    do.call(what = rbind)
  output$gene <- rownames(output)
  all_samples_for_age_estimation_gene_compare[[parameters]] <- output
}

all_samples_for_age_estimation_gene_compare_df <- all_samples_for_age_estimation_gene_compare %>%
  do.call(what = rbind) %>%
  rowwise() %>%
  mutate(
    pop_size = as.numeric(str_split(parameters,'_')[[1]][1]),
    gen_time = as.numeric(str_split(parameters,'_')[[1]][2])
  ) 

all_samples_for_age_estimation_gene_compare_df %>%
  ggplot(aes(x = as.factor(pop_size),y = Q50,ymin = Q05,ymax = Q95,colour = as.factor(gen_time))) + 
  geom_linerange(position = position_dodge(width = 0.7),
                 size = 0.25) + 
  geom_point(position = position_dodge(width = 0.7),
             size = 0.5) + 
  facet_wrap(~ gene,scales = "free",ncol = 4) + 
  theme_gerstung(base_size = 6) + 
  scale_y_continuous(limits = c(0,100),breaks = c(0,20,40,60,80)) + 
  xlab("Population size") + 
  ylab("Projected age at onset") + 
  scale_color_manual(values = pal_gsea(n = 6)(6),
                     name = "Generations/year") +
  theme(legend.position = "bottom",
        legend.key.size = unit(0.1,"cm")) + 
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/estimate_variability.pdf",model_id),width = 6,height = 6)
```

```{r,fig.height=4,fig.width=7}
ci_cut <- function(x,ci = 0.90) {
  L <- length(x)
  ci_2 <- (1 - ci)/2
  lower_bound <- ci_2
  upper_bound <- 1 - ci_2
  a <- floor(L * lower_bound)
  b <- floor(L * upper_bound)
  o <- rep(F,L)
  o[a:b] <- T
  return(o)
}

all_samples_for_age_estimation_gene_specific <- list()

for (parameters in c("50000_1","50000_2","50000_5","50000_13",
                     "1e+05_1","1e+05_2","1e+05_5","1e+05_13",
                     "2e+05_1","2e+05_2","2e+05_5","2e+05_13")) {
  x <- list()

  for (S in list_of_age_sampling[[parameters]]) {
    gene <- as.character(S$gene[1])
    if (gene == "SRSF2") {
      if (S$site[1] == "SRSF2-P95H") gene <- "SRSF2-P95H"
      else gene <- "SRSF2-other"
    }
    if (gene == "DNMT3A") {
      if (S$site[1] %in% unique_site_multiple) gene <- "DNMT3A-hotspot"
      else gene <- "DNMT3A-other"
    }
    idx <- length(x[[gene]]) + 1
    x[[gene]][[idx]] <- ifelse(S$age_distribution > S$t1,S$t1,S$age_distribution)
  }
  output <- names(x) %>% lapply(function(n) {
    y <- x[[n]]
    J <- do.call(c,y)
    J <- J[!is.na(J)]
    J <- ifelse(J < -1,-1,J)
    O <- data.frame(
      age_dist = J,
      gene = n
    )
    return(O)}) %>% 
    do.call(what = rbind)
  
  output <- output %>%
    group_by(gene) %>%
    arrange(age_dist) %>% 
    filter(ci_cut(age_dist,0.95)) %>%
    summarise(
      x = density(age_dist,bw = 1)$x,y = density(age_dist,bw = 1)$y
    ) 
  output$pop_size <- as.numeric(str_split(parameters,'_')[[1]][1])
  output$gen_time <- as.numeric(str_split(parameters,'_')[[1]][2])
  all_samples_for_age_estimation_gene_specific[[parameters]] <- output
}

all_samples_for_age_estimation_gene_specific_df <- all_samples_for_age_estimation_gene_specific %>%
  do.call(what = rbind)

all_samples_for_age_estimation_gene_specific_df %>%
  mutate(shift = as.numeric(as.factor(pop_size))) %>% 
  group_by(gene) %>%
  mutate(y = y / max(y)) %>%
  ggplot(aes(x = x,y = shift + y,group = paste(pop_size,gen_time),colour = factor(gen_time,levels = c(1,2,5,13)))) +
  geom_line(size = 0.25,alpha = 0.7) + 
  theme_gerstung(base_size = 6) +
  facet_wrap(~ gene,scales = "free_y") + 
  theme(legend.position = "bottom",
        legend.key.height = unit(0.2,"cm"),
        legend.key.width = unit(0.3,"cm")) + 
  xlab("Projected age at onset") + 
  ylab("Density") + 
  scale_y_continuous(breaks = c(1,2,3),labels = c("N=50,000","N=100,000","N=200,000")) +
  scale_colour_aaas(name = "Generations/year") + 
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/age_at_clone_foundation/variability_density.pdf",model_id),width = 7,height = 4)
```

#### Smoking and age at onset

The distribution of the ages at onset between those with and without a past or current smoking experience appears to have no distinct differences.

```{r differences_smoking_coefficient, fig.height=2,fig.width=1.125}
smoking_freq_data <- load_smoking_data() %>% 
  merge(distinct(select(full_data,SardID,phase = Phase,Age,Gender)),by = c("SardID","phase")) %>% 
  group_by(SardID) %>% 
  filter(phase == max(phase)) %>% 
  select(SardID,CurrentSmoker,QuitYears,QuitSmoking,PackYears,Age,Gender) %>% 
  distinct %>%
  mutate(status = ifelse(QuitSmoking == 'N' & CurrentSmoker == 'N',"Never smoked",
                         ifelse(QuitSmoking == 'N' & CurrentSmoker == 'Y',"Smoker",
                                "Previous smoker"))) %>% 
  subset(!is.na(CurrentSmoker)) %>% 
  group_by(SardID,Gender) %>%
  summarise(
    status = ifelse(
      sum(status == 'Smoker') > 0,'During',
      ifelse(sum(status == 'Previous smoker') > 0, 'Before',
             'Never')
    ),
    pack_years = mean(PackYears)
  ) %>% 
  mutate(status = factor(status,levels = c("Never","Before","During"))) %>%
  merge(r_values,by.x = "SardID",by.y = "individual") %>%
  mutate(individual = SardID) %>%
  subset(sum_stat == T) %>% 
  mutate(genes = str_match(genes,'[A-Z0-9]+')) %>%
  group_by(SardID) %>%
  mutate(effect = coefficient) %>%
  filter(effect == max(effect)) %>% 
  select(SardID,effect,status,age,pack_years,Gender) %>%
  mutate(status = ifelse(status == 'Never','No','Yes')) %>%
  distinct() 

n_smokers <- smoking_freq_data %>%
  select(SardID,status) %>%
  distinct %>% 
  group_by(status) %>%
  summarise(N = length(unique(SardID))) %>% 
  arrange(status)

smoking_plot <- smoking_freq_data %>% 
  select(-age) %>%
  distinct %>% 
  ggplot(aes(x = status,y = effect)) + 
  geom_jitter(width = 0.2,size = 0.5) +
  geom_boxplot(alpha = 0.5,
               size = 0.2,outlier.alpha = 0) +
  theme_gerstung(base_size = 6) +
  geom_signif(comparisons = list(c("Yes","No")),
              textsize = 2,
              size = 0.2,
              map_signif_level = function(x) return(paste("p-value=",format(x,scientifier = T,digits = 2))),
              test = wilcox.test) + 
  theme(legend.position = 'bottom',
        axis.text = element_text(size = 6)) + 
  scale_color_manual(values = gene_colours,name = NULL) + 
  xlab("Smoked") + 
  ylab("Driver gene effect") + 
  coord_cartesian(ylim = c(NA,0.52)) +
  scale_x_discrete(breaks = c("No","Yes"),
                   labels = c(sprintf("No\n(n=%s)",n_smokers$N[1]),
                              sprintf("Yes\n(n=%s)",n_smokers$N[2])))
```

```{r smoking_clone_age, fig.width=3, fig.height=2}
INTERVAL <- 5
smoking_clone_ages <- load_smoking_data() %>% 
  mutate(status = ifelse(QuitSmoking == 'N' & CurrentSmoker == 'N',"Never smoked",
                         ifelse(QuitSmoking == 'N' & CurrentSmoker == 'Y',"Smoker",
                                "Previous smoker"))) %>% 
  subset(!is.na(CurrentSmoker)) %>% 
  group_by(SardID) %>%
  filter(phase == max(phase)) %>%
  ungroup() %>%
  merge(Parameters_Age,by.x = "SardID",by.y = "individual",all.x = T) %>%
  select(-age) %>%
  merge(distinct(select(subset(load_data_keep(),Gene %in% load_included_genes() | is.na(Gene)),Gender,SardID,age=Age)),
        by = "SardID",all.y = T) %>%
  mutate(individual = SardID) %>%
  mutate(genes = str_match(gene,'[A-Z0-9]+')) %>%
  group_by(SardID) %>%
  mutate(effect = Q50) %>%
  filter(age == max(age)) %>%
  select(SardID,effect,site,CurrentSmoker, QuitYears,status,age,genes,Gender,b_clone_mean) %>% 
  #subset(effect > 18 & effect < 100) %>% 
  mutate(effect = age - effect) %>% 
  mutate(QuitYears = ifelse(CurrentSmoker == 'Y',0,QuitYears)) %>% 
  mutate(CurrentlySmoking = effect > QuitYears) %>%
  mutate(CurrentlySmoking = ifelse(status == "Never smoked",FALSE,CurrentlySmoking)) %>%
  mutate(AgeGroups = floor((age - effect) / INTERVAL) * INTERVAL)

clone_ages_plot <- smoking_clone_ages %>% 
  ungroup() %>%
  subset(AgeGroups <= 60 | is.na(AgeGroups)) %>% 
  subset(AgeGroups >= 15) %>% 
  mutate(CurrentlySmoking = ifelse(CurrentlySmoking == T,"Yes","No")) %>% 
  ggplot(aes(x = effect,colour = CurrentlySmoking,fill = CurrentlySmoking)) +
  geom_density(alpha = 0.2,size = 0.25) + 
  theme_gerstung(base_size = 6) + 
  scale_x_continuous(expand = c(0,0)) + 
  scale_y_continuous(expand = c(0,0,0.05,0)) + 
  xlab("Age") + 
  ylab("Density of clones appearing") +
  theme(legend.position = "top",
        legend.key.height = unit(0.2,"cm"),
        legend.key.width = unit(0.2,"cm"),
        axis.text = element_text(size = 6)) + 
  scale_colour_simpsons(name = NULL,breaks = c("No","Yes"),
                        labels = c("Not smoking","Smoking")) + 
  scale_fill_simpsons(name = NULL,breaks = c("No","Yes"),
                      labels = c("Not smoking","Smoking"))

plot_grid(smoking_plot,clone_ages_plot,
          rel_widths = c(0.6,1),nrow = 1) +
  ggsave(useDingbats=FALSE,
         filename = sprintf("figures/%s/dynamic_coefficients/smoking.pdf",
                            model_id),
         width = 3,height = 1.5)
```

### Assessing the predictive power of this model

Upon colony sequencing, 13 individuals were sequenced at an additional timepoint. We use this data to understand how well our model performs when predicting for individuals which have been studied for some time. Here we can observe that $MAE(test)=3.5\%$, which is 7 times larger than our clone detection limit ($0.5\%$). As such, it can be risky to make predictions at low detection values, but clones show a fairly stable behaviour as they progress and as implied by the fact that 92.3% of all clones are explained by the simple model we propose here.

```{r predictive_power_6th_tp,fig.width=4,fig.height=2,warning=F,message=F}
beta_value <- values_model$beta_values[[1]][1,1]
r_values_for_prediction <- r_values_full %>%
  select(site,individual,
         coefficient,coefficient_005,coefficient_095,
         b_clone,b_clone_005,b_clone_095,
         u_values,u_values_005,u_values_095) %>%
  distinct %>%
  mutate(site = paste(str_match(site,'[A-Z0-9]+'),str_match(site,'-.*'),sep = ''))

predict_timepoint <- function(Individual,Site,Age,coverage=1000) {
  Site <- as.character(Site)
  Individual <- as.numeric(Individual)
  Age <- as.numeric(Age) - values_model$min_age
  tmp <- r_values_for_prediction %>% 
    subset((individual == Individual) & (site == Site)) %>% 
    mutate(
      p_005 = inv.logit((coefficient_005+b_clone_005) * Age + u_values_005) * 0.5,
      p = inv.logit((coefficient+b_clone) * Age + u_values) * 0.5,
      p_095 = inv.logit((coefficient_095+b_clone_095) * Age + u_values_095) * 0.5) %>% 
    transmute(
      coefficient = coefficient,
      b_clone = b_clone,
      u_values = u_values,
      p = p
    ) 
  S <- rbbinom(1000,coverage,(tmp$p*beta_value)/(1-tmp$p),beta_value)/coverage
  S <- ifelse(S > 0.5, 0.5,S)
  tmp$pred_005 <- quantile(S,0.05,names=F,na.rm = T)
  tmp$pred = mean(S)
  tmp$pred_095 = quantile(S,0.95,names=F,na.rm = T)
  return(tmp)
}
data_6th_tp <- load_data_6th_tp()

predictions_6th_tp <- data_6th_tp %>%
  apply(1,FUN=function(x) {
    c(x[1],x[12],x[4],
      predict_timepoint(x[1],x[12],x[4],1000),
      VAF=as.numeric(x[10]))
    }) %>%
  lapply(unlist) %>% 
  do.call(what=rbind) %>%
  as.data.frame() %>%
  mutate(coefficient = as.numeric(as.character(coefficient)),
         b_clone = as.numeric(as.character(b_clone)),
         u_values = as.numeric(as.character(u_values)),
         p = as.numeric(as.character(p)),
         pred_005 = as.numeric(as.character(pred_005)),
         pred = as.numeric(as.character(pred)),
         pred_095 = as.numeric(as.character(pred_095)),
         VAF = as.numeric(as.character(VAF))) %>%
  na.omit

last_tp <- r_values_full %>%
  mutate(site = paste(str_match(site,'[A-Z0-9]+'),str_match(site,'-.*'),sep = '')) %>% 
  subset(individual %in% data_6th_tp$SardID & site %in% data_6th_tp$amino_acid_change) %>% 
  transmute(SardID=individual,amino_acid_change=site,Age=age,coefficient,b_clone,u_values,
            p=mu_val,pred_005=pred_005/coverage,pred=pred/coverage,pred_095=pred_095/coverage,VAF=true/coverage) %>%
  group_by(SardID,amino_acid_change) %>%
  filter(Age == max(Age)) %>% 
  ungroup() %>%
  rbind(predictions_6th_tp) %>%
  mutate(gene = sapply(amino_acid_change,function(x) unlist(strsplit(x,'-'))[1])) %>%
  mutate(amino_acid_change = sapply(amino_acid_change,function(x) unlist(strsplit(x,'-'))[2])) %>%
  mutate(amino_acid_change = sapply(amino_acid_change,function(x) unlist(strsplit(x,';'))[1])) %>%
  mutate(amino_acid_change = sapply(amino_acid_change,function(x) {
    tmp <- unlist(strsplit(x,':'))
    return(tmp[length(tmp)])})) %>% 
  mutate(amino_acid_change = paste(gene,amino_acid_change,sep = '-')) %>% 
  na.omit() %>% 
  mutate(
    SardID = gsub(" ","",as.character(SardID)),
    amino_acid_change = as.character(amino_acid_change)) %>% 
  group_by(SardID,amino_acid_change) %>%
  mutate(last_timepoint = Age == max(Age)) %>%
  ungroup() %>%
  mutate(Age = as.numeric(Age)) %>% 
  filter(amino_acid_change %in% c("DNMT3A-R882H","GNB1-K57E",
                                  "JAK2-V617F","SF3B1-K666N",
                                  "SF3B1-K700E","SRSF2-P95H",
                                  "SRSF2-P95L","TET2-Q1542X",
                                  "U2AF1-Q157P","U2AF1-Q157R"))
last_tp_only <- last_tp %>%
  group_by(SardID, amino_acid_change) %>% 
  filter(Age == max(Age))

prediction_data <- r_values_full %>%
  mutate(site = paste(str_match(site,'[A-Z0-9]+'),str_match(site,'-.*'),sep = '')) %>% 
  subset(individual %in% data_6th_tp$SardID & site %in% data_6th_tp$amino_acid_change) %>% 
  transmute(SardID=individual,amino_acid_change=site,Age=age,coefficient,b_clone,u_values,
            p=mu_val,pred_005=pred_005/coverage,pred=pred/coverage,pred_095=pred_095/coverage,VAF=true/coverage) %>%
  mutate(gene = sapply(amino_acid_change,function(x) unlist(strsplit(x,'-'))[1])) %>%
  mutate(amino_acid_change = sapply(amino_acid_change,function(x) unlist(strsplit(x,'-'))[2])) %>%
  mutate(amino_acid_change = sapply(amino_acid_change,function(x) unlist(strsplit(x,';'))[1])) %>%
  mutate(amino_acid_change = sapply(amino_acid_change,function(x) {
    tmp <- unlist(strsplit(x,':'))
    return(tmp[length(tmp)])})) %>% 
  mutate(amino_acid_change = paste(gene,amino_acid_change,sep = '-')) %>% 
  na.omit() %>% 
  mutate(
    SardID = gsub(" ","",as.character(SardID)),
    amino_acid_change = as.character(amino_acid_change)) %>% 
  filter(amino_acid_change %in% c("DNMT3A-R882H","GNB1-K57E",
                                  "JAK2-V617F","SF3B1-K666N",
                                  "SF3B1-K700E","SRSF2-P95H",
                                  "SRSF2-P95L","TET2-Q1542X",
                                  "U2AF1-Q157P","U2AF1-Q157R"))

prediction_data %>%
  #filter(any(last_timepoint == T)) %>% 
  ggplot(aes(x = Age,
             y = pred,
             ymin = pred_005, 
             ymax = pred_095,
             group = SardID,
             fill = SardID,
             color = SardID)) + 
  geom_line(aes(y = VAF),linetype = 'dotted',color = 'black',alpha = 0.2,
            size = 0.25) + 
  geom_ribbon(alpha = 0.2,color = NA) + 
  geom_point(aes(y = VAF),alpha = 0.2,color = 'black',
             size = 0.8) + 
  geom_ribbon(data = last_tp,alpha = 0.2,color = NA) + 
  geom_point(data = last_tp_only,aes(y = VAF),shape = 1,
             size = 0.5) + 
  geom_point(data = last_tp_only,aes(y = VAF),color = 'black',shape = 3,
             size = 0.5) + 
  geom_line(data = last_tp,aes(y = VAF),linetype = 'dotted',color = 'black',alpha = 0.2,
            size = 0.25) + 
  geom_line(alpha = 0.6,
            size = 0.25) + 
  geom_line(data = last_tp,linetype = 'longdash',alpha = 0.6,
            size = 0.25) +
  facet_wrap(~ paste(amino_acid_change,sep='-'),nrow = 2,scales = 'free') +
  #scale_y_continuous(trans = 'log10') + 
  theme_gerstung(base_size = 6) + 
  scale_color_discrete(guide=F) + 
  scale_fill_discrete(guide=F) + 
  scale_shape_discrete(guide=F) + 
  scale_y_continuous(trans = 'log10',limits = c(5e-4,0.5),
                     labels = function(x) sprintf("%s%%",x*100)) +
  scale_x_continuous(breaks = seq(0,100,by = 10),limits = c(58,92),expand=c(0,0)) + 
  ylab("VAF") + 
  ggtitle("Predictions for the last timepoint conditioned on the previous timepoints") +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/inferred_trajectories/predicted_trajectories_6th_tp.pdf",model_id),
         width = 5,height = 2)
```

```{r compare_truth_and_predictions,fig.height=2,fig.width=4}
pred_true_change <- last_tp %>%
  group_by(SardID,amino_acid_change,gene) %>%
  summarise(pred_increase = pred[which.max(Age)]/pred[which.max(Age)-1],
            vaf_increase = VAF[which.max(Age)]/VAF[which.max(Age)-1],
            pred_005 = pred_005[which.max(Age)],
            pred = pred[which.max(Age)],
            pred_095 = pred_095[which.max(Age)],
            vaf = VAF[which.max(Age)]) 

plot_min_max <- c(
  min(pred_true_change$pred_increase,pred_true_change$vaf_increase),
  max(pred_true_change$pred_increase,pred_true_change$vaf_increase)
)

performance_scatter_plot <- pred_true_change %>%
  ggplot(aes(x = vaf_increase,
             y = pred_increase,
             color = gene)) + 
  geom_abline(slope = 1,linetype = 3) + 
  geom_point(size=0.8) + 
  xlab("Observed fold change in VAF") + 
  ylab("Predicted fold change in VAF")  + 
  theme_gerstung(base_size = 6) +
  scale_x_continuous(limits = plot_min_max,trans = 'log10') + 
  scale_y_continuous(limits = plot_min_max,trans = 'log10') + 
  scale_color_manual(values=gene_colours,name=NULL) + 
  scale_fill_manual(values=gene_colours,name=NULL) + 
  theme(legend.text = element_text(face = "italic"),
        legend.key.size = unit(0,"cm"))

plot_min_max_ <- c(
  0.005,
  max(pred_true_change$pred_095,pred_true_change$vaf)
)

pred_true_change %>% 
  ggplot(aes(x = vaf,y = pred,ymin = pred_005,ymax = pred_095,colour = gene)) + 
  geom_abline(intercept = 0,slope = 1,linetype = 3) + 
  geom_point(size = 0.5) + 
  geom_linerange(size = 0.25) + 
  coord_cartesian(xlim = plot_min_max_,ylim = plot_min_max_) +
  xlab("Observed VAF") + 
  ylab("Predicted VAF")  + 
  theme_gerstung(base_size = 6) +
  scale_x_continuous(trans = 'log10') + 
  scale_y_continuous(trans = 'log10') + 
  theme(legend.text = element_text(face = "italic"),
        legend.key.size = unit(0,"cm")) +
  scale_colour_manual(values = gene_colours,name = NULL) +
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/inferred_trajectories/performance_6tp.pdf",model_id),
         height = 1.5,width = 2.3)


interesting_cases_6th_tp_plot <- prediction_data %>%
  subset(amino_acid_change %in% c("SF3B1-K666N","SRSF2-P95H")) %>% 
  ggplot(aes(x = Age,
             y = pred,
             ymin = pred_005, 
             ymax = pred_095,
             group = SardID,
             fill = SardID,
             color = SardID)) + 
  geom_ribbon(alpha = 0.2,color = NA) + 
  geom_point(aes(y = VAF),alpha = 0.2,color = 'black',
             size = 0.5) + 
  geom_ribbon(data = subset(last_tp,amino_acid_change %in% c("SF3B1-K666N","SRSF2-P95H")),alpha = 0.2,color = NA) + 
  geom_point(data = subset(last_tp_only,amino_acid_change %in% c("SF3B1-K666N","SRSF2-P95H")),
             aes(y = VAF),shape = 1,size = 0.5) + 
  geom_point(data = subset(last_tp_only,amino_acid_change %in% c("SF3B1-K666N","SRSF2-P95H")),
             aes(y = VAF),color = 'black',shape = 3,size = 0.5) + 
  geom_line(alpha = 0.6,size = 0.25) + 
  geom_line(data = subset(last_tp,amino_acid_change %in% c("SF3B1-K666N","SRSF2-P95H")),
            linetype = 'longdash',alpha = 0.6,size = 0.25) +
  facet_wrap(~ paste(amino_acid_change,sep='-'),nrow = 2,scales = 'free_x') +
  theme_gerstung(base_size = 6) + 
  scale_color_discrete(guide=F) + 
  scale_fill_discrete(guide=F) + 
  scale_shape_discrete(guide=F) + 
  scale_y_continuous(trans = 'log10') +
  scale_x_continuous(breaks = seq(0,100,by = 10),limits = c(58,92),expand=c(0,0)) + 
  ylab("VAF") 

plot_grid(
  plot_grid(performance_scatter_plot + theme(legend.position = "none",
                                             panel.grid = element_blank()),
            interesting_cases_6th_tp_plot,nrow=1,rel_widths = c(4.3/7,2.7/7)),
  plot_grid(get_legend(performance_scatter_plot + theme(legend.position = "bottom")),
            ggplot() + theme_void(),
            rel_widths = c(4.3/7,2.7/7)),
  rel_heights = c(0.9,0.1),
  ncol = 1
) + 
  ggsave(useDingbats=FALSE,
         sprintf("figures/%s/inferred_trajectories/performance_and_predicted_6thpoint_trajectories.pdf",model_id),
         height = 2,width = 4)
```

## Global picture

Using what we have so far we can replicate the general picture derived by [McKerrel et al.](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4542313/).

```{r global_picture,fig.width=5,fig.height=4}
average_age <- Parameters_Age %>%
  group_by(site) %>%
  summarise(MedianAge = mean(Q50),
            Gene = gene[1]) %>% 
  transmute(site = site,
            Gene = Gene,
            MedianAge = MedianAge)

interesting_sites_to_see <- r_values %>% 
  subset(!is.na(site_numeric)) %>% 
  mutate(Gene = str_match(genes,'[A-Z0-9]+')) %>%
  group_by(Gene) %>% 
  filter(site %in% c("JAK2-V617F","DNMT3Ant-R882H","SF3B1-K666N","SRSF2-P95H","U2AF1-Q157R")) %>%
  select(Gene,site,coefficient) %>%
  distinct() 

global_picture_parameters <- merge(average_age,interesting_sites_to_see,
                                   by = c("Gene","site"),
                                   all = F) %>%
  arrange(Gene) %>%
  mutate(site = paste(Gene,str_match(site,'[0-9A-Za-z]+$'),sep = '-'))

TOT_READS <- 1e3
global_picture_out <- list()
for (S in unique(global_picture_parameters$site)) {
  subset_df <- global_picture_parameters[global_picture_parameters$site == S,]
  global_picture_out[[S]] <- data.frame(
    time = seq(floor(subset_df$MedianAge),80),
    draws = trajectory_from_t0(
      x = c(floor(subset_df$MedianAge):80),
      t0 = subset_df$MedianAge,
      s = subset_df$coefficient,N = 0.5e5,g = 2)/0.5e5,
    site = S,
    Gene = subset_df$Gene
  )
}

text_global_picture_out <- global_picture_out %>% 
  do.call(what = rbind) %>%
  group_by(site) %>% 
  summarise(x = max(time),
            y = max(draws),
            label = site[1])

global_picture_out %>% 
  do.call(what = rbind) %>%
  ggplot(aes(x = time,y = draws,colour = Gene,group = site)) + 
  geom_line(size = 1,
            alpha = 0.5) + 
  geom_label_repel(inherit.aes = F,
            size = 2.0,
            alpha = 0.5,
            data = text_global_picture_out,
            hjust = 0,
            vjust = 1,
            fontface = "bold",
            min.segment.length = 0,
            force = 0.1,
            direction = "y",
            label.r = unit(0,"cm"),
            label.size = unit(0,"cm"),
            label.padding = unit(0.05,"cm"),
            aes(x = x,y = y,label = label,fill = str_match(label,'[0-9A-Z]+'))) +
  scale_y_continuous(trans = 'log10',breaks = c(0.001,0.010,0.10,0.20,0.3,0.5),
                     minor_breaks = NULL,
                     expand = c(0,0,0.05,0)) + 
  scale_x_continuous(limits = c(0,100),
                     breaks = c(0,20,40,60,80),
                     minor_breaks = NULL,
                     expand = c(0,0.05)) +
  scale_color_manual(guide=F,values = gene_colours) + 
  scale_fill_manual(guide=F,values = gene_colours) +
  theme_gerstung(base_size = 6) +
  theme(legend.position = 'bottom',
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        axis.line = element_line()) +
  xlab("Age") + 
  ylab("logVAF") +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/inferred_trajectories/fiugre_like_mckerrel2015.pdf",model_id),
         height = 4,width = 5) 
```

### DNMT3A and TET2

Here we present a systematic comparison of *DNMT3A* and *TET2* clones to better assess whether there is any real difference between them.

```{r,fig.height=2,fig.width=3}
Parameters_Age_DNMT3A_TET2 <- subset(Parameters_Age,gene %in% c("DNMT3A","TET2"))

Parameters_Age %>%
  mutate(AgeGroup = cut(maxAge,c(0,60,70,80,90,200),c("<60","60-70","70-80","80-90",">90"))) %>% 
  group_by(individual) %>%
  summarise(TET2 = "TET2" %in% gene,
            DNMT3A = "DNMT3A" %in% gene) %>%
  summarise(TotalTET2 = sum(TET2),
            TotalDNMT3A = sum(DNMT3A),
            TotalBoth = sum(TET2 == T & DNMT3A == T),
            TotalNone = sum(TET2 == F & DNMT3A == F)) %>%
  mutate(TotalTET2 = TotalTET2 - TotalBoth,
         TotalDNMT3A = TotalDNMT3A - TotalBoth) %>% 
  gather(key = "key",value = "value") %>% 
  mutate(key = factor(key,levels = c("TotalTET2","TotalDNMT3A","TotalBoth","TotalNone"))) %>% 
  ggplot(aes(x = "",y = value,fill = key)) + 
  geom_bar(stat = "identity") + 
  coord_polar("y",start = 0) + 
  theme_gerstung(base_size = 6) +
  theme(axis.line = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank()) +
  scale_fill_manual(values = c(TotalTET2 = "goldenrod",TotalDNMT3A = "red4",TotalBoth = "orange",TotalNone = "white"),
                    labels = c("Only TET2","Only DNMT3A","TET2 & DNMT3A",""),
                    name = NULL) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/prevalence/dnmt3a_vs_tet2_prevalence.pdf",model_id),
         height = 2,width = 3) 
```

```{r,fig.height=6,fig.width=3}
plot_grid(
  load_data_keep() %>%
    subset(VAF > 0.005) %>%
    group_by(SardID) %>% 
    mutate(Age = max(Age)) %>% 
    mutate(AgeGroup = cut(Age,c(0,70,80,90,200),c("<70","70-80","80-90",">90"))) %>% 
    group_by(SardID,AgeGroup) %>%
    summarise(TET2 = sum(grepl("TET2",unique(AAChange.refGene))),
              DNMT3A = sum(grepl("DNMT3A",unique(AAChange.refGene))),
              Total = 1) %>%
    group_by(AgeGroup) %>% 
    summarise(TotalTET2 = sum(TET2),
              TotalDNMT3A = sum(DNMT3A),
              Total = sum(Total)) %>%
    gather(key = "key",value = "value",-AgeGroup,-Total) %>%
    group_by(AgeGroup) %>% 
    mutate(p05 = qbeta(0.05,value,Total - value),
           p95 = qbeta(0.95,value,Total - value)) %>% 
    mutate(key = factor(key, levels = c("TotalDNMT3A","TotalTET2","TotalBoth","TotalNone"),
                        labels = c("DNMT3A","TET2","DNMT3A+TET2","None"))) %>% 
    ggplot(aes(x = key,y = value/Total,fill = AgeGroup,group = AgeGroup)) + 
    geom_bar(stat = "identity",position = position_dodge(width = 0.9)) + 
    geom_linerange(aes(ymin = p05,ymax = p95),position = position_dodge(width = 0.9)) +
    theme_gerstung(base_size = 6) + 
    scale_fill_brewer(name = NULL,palette = 3) + 
    theme(axis.title.x = element_blank()) + 
    scale_y_continuous(labels = function(x) return(sprintf("%s%%",x*100)),
                       expand = c(0,0)) + 
    ylab("Prevalence") + 
    ggtitle("Considering all individuals only once\nand all clones"),
  load_data_keep() %>%
    subset(VAF > 0.005) %>% 
    mutate(AgeGroup = cut(Age,c(0,70,80,90,200),c("<70","70-80","80-90",">90"))) %>% 
    group_by(SardID,AgeGroup) %>%
    summarise(TET2 = sum(grepl("TET2",unique(AAChange.refGene))),
              DNMT3A = sum(grepl("DNMT3A",unique(AAChange.refGene))),
              Total = 1) %>%
    group_by(AgeGroup) %>% 
    summarise(TotalTET2 = sum(TET2),
              TotalDNMT3A = sum(DNMT3A),
              Total = sum(Total)) %>%
    gather(key = "key",value = "value",-AgeGroup,-Total) %>%
    group_by(AgeGroup) %>% 
    mutate(p05 = qbeta(0.05,value,Total - value),
           p95 = qbeta(0.95,value,Total - value)) %>% 
    mutate(key = factor(key, levels = c("TotalDNMT3A","TotalTET2","TotalBoth","TotalNone"),
                        labels = c("DNMT3A","TET2","DNMT3A+TET2","None"))) %>% 
    ggplot(aes(x = key,y = value/Total,fill = AgeGroup,group = AgeGroup)) + 
    geom_bar(stat = "identity",position = position_dodge(width = 0.9)) + 
    geom_linerange(aes(ymin = p05,ymax = p95),position = position_dodge(width = 0.9)) +
    theme_gerstung(base_size = 6) + 
    scale_fill_brewer(name = NULL,palette = 3) + 
    theme(axis.title.x = element_blank()) + 
    scale_y_continuous(labels = function(x) return(sprintf("%s%%",x*100)),
                       expand = c(0,0)) + 
    ylab("Prevalence") + 
    ggtitle("Considering all individuals across all timepoints\nand all clones"),
  load_data_keep() %>%
    subset(VAF > 0.005) %>% 
    group_by(SardID) %>% 
    mutate(Age = max(Age)) %>% 
    mutate(AgeGroup = cut(Age,c(0,70,80,90,200),c("<70","70-80","80-90",">90"))) %>% 
    group_by(SardID,AgeGroup) %>%
    summarise(TET2 = any(grepl("TET2",unique(AAChange.refGene))),
              DNMT3A = any(grepl("DNMT3A",unique(AAChange.refGene)))) %>%
    group_by(AgeGroup) %>% 
    summarise(TotalTET2 = sum(TET2),
              TotalDNMT3A = sum(DNMT3A),
              TotalBoth = sum(TET2 == T & DNMT3A == T),
              TotalNone = sum(TET2 == F & DNMT3A == F)) %>%
    mutate(TotalTET2 = TotalTET2 - TotalBoth,
           TotalDNMT3A = TotalDNMT3A - TotalBoth) %>%
    gather(key = "key",value = "value",-AgeGroup) %>%
    group_by(AgeGroup) %>% 
    mutate(Total = sum(value)) %>% 
    mutate(p05 = qbeta(0.05,value,Total - value),
           p95 = qbeta(0.95,value,Total - value)) %>% 
    mutate(key = factor(key, levels = c("TotalDNMT3A","TotalTET2","TotalBoth","TotalNone"),
                        labels = c("DNMT3A","TET2","DNMT3A+TET2","None"))) %>% 
    ggplot(aes(x = key,y = value/Total,fill = AgeGroup,group = AgeGroup)) + 
    geom_bar(stat = "identity",position = position_dodge(width = 0.9)) + 
    geom_linerange(aes(ymin = p05,ymax = p95),position = position_dodge(width = 0.9)) +
    theme_gerstung(base_size = 6) + 
    scale_fill_brewer(name = NULL,palette = 3) + 
    theme(axis.title.x = element_blank()) + 
    scale_y_continuous(labels = function(x) return(sprintf("%s%%",x*100)),
                       expand = c(0,0)) + 
    ylab("Prevalence") + 
    ggtitle("Considering all individuals only once\nand without considering more than one clone per individual"),
  load_data_keep() %>%
    subset(VAF > 0.005) %>% 
    mutate(AgeGroup = cut(Age,c(0,70,80,90,200),c("<70","70-80","80-90",">90"))) %>% 
    group_by(SardID,AgeGroup) %>%
    summarise(TET2 = any(grepl("TET2",unique(AAChange.refGene))),
              DNMT3A = any(grepl("DNMT3A",unique(AAChange.refGene)))) %>%
    group_by(AgeGroup) %>% 
    summarise(TotalTET2 = sum(TET2),
              TotalDNMT3A = sum(DNMT3A),
              TotalBoth = sum(TET2 == T & DNMT3A == T),
              TotalNone = sum(TET2 == F & DNMT3A == F)) %>%
    mutate(TotalTET2 = TotalTET2 - TotalBoth,
           TotalDNMT3A = TotalDNMT3A - TotalBoth) %>%
    gather(key = "key",value = "value",-AgeGroup) %>%
    group_by(AgeGroup) %>% 
    mutate(Total = sum(value)) %>% 
    mutate(p05 = qbeta(0.05,value,Total - value),
           p95 = qbeta(0.95,value,Total - value)) %>% 
    mutate(key = factor(key, levels = c("TotalDNMT3A","TotalTET2","TotalBoth","TotalNone"),
                        labels = c("DNMT3A","TET2","DNMT3A+TET2","None"))) %>% 
    ggplot(aes(x = key,y = value/Total,fill = AgeGroup,group = AgeGroup)) + 
    geom_bar(stat = "identity",position = position_dodge(width = 0.9)) + 
    geom_linerange(aes(ymin = p05,ymax = p95),position = position_dodge(width = 0.9)) +
    theme_gerstung(base_size = 6) + 
    scale_fill_brewer(name = NULL,palette = 3) + 
    theme(axis.title.x = element_blank()) + 
    scale_y_continuous(labels = function(x) return(sprintf("%s%%",x*100)),
                       expand = c(0,0)) + 
    ylab("Prevalence") + 
    ggtitle("Considering all individuals across all timepoints\nand without considering more than one clone per individual"),
  ncol = 1) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/prevalence/dnmt3a_vs_tet2_prevalence.pdf",model_id),
         height = 6,width = 3)
```

```{r,fig.height=3,fig.width=4}
plot_grid(
  full_data %>%
    mutate(AgeGroup = cut(Age,c(0,60,70,80,90,200),c("<60","60-70","70-80","80-90",">90"))) %>% 
    group_by(SardID,AgeGroup) %>%
    summarise(TET2 = sum(grepl("TET2",unique(amino_acid_change))),
              DNMT3A = sum(grepl("DNMT3A",unique(amino_acid_change)))) %>%
    mutate(TET2 = cut(TET2,c(-1,0,1,3,5,20),c("None","1","2-3","4-5","6+"))) %>%
    mutate(DNMT3A = cut(DNMT3A,c(-1,0,1,3,5,20),c("None","1","2-3","4-5","6+"))) %>%
    gather(key = "Gene",value = "value",TET2,DNMT3A) %>%
    group_by(AgeGroup,Gene,value) %>% 
    summarise(N = length(unique(SardID))) %>%
    group_by(AgeGroup,Gene) %>%
    mutate(Total = sum(N)) %>% 
    mutate(p05 = qbeta(0.05,N,Total - N),
           p95 = qbeta(0.95,N,Total - N)) %>% 
    na.omit() %>% 
    mutate(value = factor(value,levels = c("None","1","2-3","4-5","6+"))) %>% 
    ggplot(aes(x = AgeGroup,y = N/Total,fill = value,group = value)) + 
    geom_bar(stat = "identity",position = position_dodge(width = 0.9)) + 
    geom_linerange(aes(ymin = p05,ymax = p95),position = position_dodge(width = 0.9)) +
    theme_gerstung(base_size = 6) + 
    scale_fill_brewer(palette = 5,name = "Number of mutations") + 
    theme(axis.title.x = element_blank()) + 
    scale_y_continuous(labels = function(x) return(sprintf("%s%%",x*100)),
                       expand = c(0,0)) + 
    ylab("Prevalence") + 
    facet_wrap(~ Gene) +
    ggtitle("Considering all individuals across all timepoints"),
  full_data %>%
    group_by(SardID) %>% 
    mutate(Age = max(Age)) %>% 
    mutate(AgeGroup = cut(Age,c(0,60,70,80,90,200),c("<60","60-70","70-80","80-90",">90"))) %>% 
    group_by(SardID,AgeGroup) %>%
    summarise(TET2 = sum(grepl("TET2",unique(amino_acid_change))),
              DNMT3A = sum(grepl("DNMT3A",unique(amino_acid_change)))) %>%
    mutate(TET2 = cut(TET2,c(-1,0,1,3,5,20),c("None","1","2-3","4-5","6+"))) %>%
    mutate(DNMT3A = cut(DNMT3A,c(-1,0,1,3,5,20),c("None","1","2-3","4-5","6+"))) %>%
    gather(key = "Gene",value = "value",TET2,DNMT3A) %>%
    group_by(AgeGroup,Gene,value) %>% 
    summarise(N = length(unique(SardID))) %>%
    group_by(AgeGroup,Gene) %>%
    mutate(Total = sum(N)) %>% 
    mutate(p05 = qbeta(0.05,N,Total - N),
           p95 = qbeta(0.95,N,Total - N)) %>% 
    na.omit() %>% 
    mutate(value = factor(value,levels = c("None","1","2-3","4-5","6+"))) %>% 
    ggplot(aes(x = AgeGroup,y = N/Total,fill = value,group = value)) + 
    geom_bar(stat = "identity",position = position_dodge(width = 0.9)) + 
    geom_linerange(aes(ymin = p05,ymax = p95),position = position_dodge(width = 0.9)) +
    theme_gerstung(base_size = 6) + 
    scale_fill_brewer(palette = 5,name = "Number of mutations") + 
    theme(axis.title.x = element_blank()) + 
    scale_y_continuous(labels = function(x) return(sprintf("%s%%",x*100)),
                       expand = c(0,0)) + 
    ylab("Prevalence") + 
    facet_wrap(~ Gene) +
    ggtitle("Considering all individuals only once"),
  ncol = 1) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/prevalence/dnmt3a_vs_tet2_prevalence_n_mutations.pdf",model_id),
         height = 3,width = 4) 
```

```{r,fig.height=1.2,fig.width=3}
Parameters_Age %>%
  group_by(individual) %>%
  filter(length(unique(site)) > 1) %>% 
  summarise(TET2IsFastest = any((b_clone_mean + b_genetic_mean)[gene == "TET2"]==max(b_clone_mean + b_genetic_mean)),
            DNMT3AIsFastest = any((b_clone_mean + b_genetic_mean)[gene == "DNMT3A"]==max(b_clone_mean + b_genetic_mean)),
            TET2 = "TET2" %in% gene,
            DNMT3A = "DNMT3A" %in% gene) %>%
  #filter(TET2 == T | DNMT3A == T) %>% 
  group_by(TET2,DNMT3A) %>%
  summarise(TET2Fastest = sum(TET2IsFastest),
            DNMT3AFastest = sum(DNMT3AIsFastest),
            N = length(TET2IsFastest)) %>%
  mutate(Other = N - TET2Fastest - DNMT3AFastest) %>% 
  rowwise() %>%
  mutate(label = paste(c("TET2","DNMT3A")[c(TET2,DNMT3A)],collapse = " & ")) %>%
  mutate(label = ifelse(label == "","None",label)) %>%
  mutate(label = factor(label,rev(c("TET2","DNMT3A","TET2 & DNMT3A","None")))) %>% 
  gather(key = "key",value = "value",TET2Fastest,DNMT3AFastest,Other) %>% 
  mutate(key = factor(key,levels = c("TET2Fastest","DNMT3AFastest","Other"),
                      labels = c("TET2","DNMT3A","Other"))) %>% 
  ggplot(aes(x = label,y = value,fill = key)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  theme_gerstung(base_size = 6) +
  theme(legend.position = "bottom",legend.key.size = unit(0.2,"cm")) + 
  scale_fill_aaas(name = "Fastest clone") + 
  scale_y_continuous(expand = c(0,0)) + 
  ylab("Number of individuals\nwith two or more mutations") + 
  xlab("") +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/dnmt3a_vs_tet2_fastest.pdf",model_id),
         height = 1.2,width = 3) 
```

```{r,fig.height=3,fig.width=6}
plot_grid(
  Parameters_Age_DNMT3A_TET2 %>%
    ggplot(aes(x = maxAge,y = maxVAF,colour = gene)) +
    geom_point(size = 0.5) + 
    scale_y_continuous(trans = 'log10') + 
    scale_color_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    xlab("Age") + 
    ylab("Maximum VAF"),
  Parameters_Age_DNMT3A_TET2 %>%
    ggplot(aes(x = maxVAF,y = (b_genetic_mean + b_clone_mean)*maxAge,colour = gene)) +
    geom_point(size = 0.5) + 
    scale_x_continuous(trans = 'log10') +
    scale_color_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    xlab("Maximum VAF") + 
    ylab("Age * Growth rate"),
  Parameters_Age_DNMT3A_TET2 %>%
    ggplot(aes(x = maxVAF,y = b_genetic_mean + b_clone_mean,colour = gene)) +
    geom_point(size = 0.5) + 
    geom_smooth(size = 0.25,method = "lm") +
    scale_color_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    scale_x_continuous(trans = 'log10') + 
    xlab("Maximum VAF") + 
    ylab("Growth rate"),
  Parameters_Age_DNMT3A_TET2 %>%
    ggplot(aes(x = maxVAF,y = b_clone_mean,colour = gene)) +
    geom_point(size = 0.5) + 
    geom_smooth(size = 0.25,method = "lm") +
    scale_color_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    scale_x_continuous(trans = 'log10') + 
    xlab("Maximum VAF") + 
    ylab("Unknown cause effect"),
  Parameters_Age_DNMT3A_TET2 %>%
    ggplot(aes(x = maxVAF,y = b_genetic_mean,colour = gene)) +
    geom_point(size = 0.5) + 
    geom_smooth(size = 0.25,method = "lm") +
    scale_color_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    scale_x_continuous(trans = 'log10') + 
    xlab("Maximum VAF") + 
    ylab("Genetic effect"),
  plot_grid(
    Parameters_Age_DNMT3A_TET2 %>%
      ggplot(aes(x = gene,y = maxVAF,fill = gene)) +
      geom_boxplot(size = 0.5) + 
      scale_y_continuous(trans = 'log10') + 
      scale_fill_manual(values = gene_colours) + 
      theme_gerstung(base_size = 6) + 
      theme(legend.position = "none") + 
      xlab("Gene") + 
      ylab("Maximum VAF") + 
      theme(axis.text.x = element_text(face = "italic")),
    Parameters_Age_DNMT3A_TET2 %>%
      ggplot(aes(x = gene,y = b_genetic_mean + b_clone_mean,fill = gene)) +
      geom_boxplot(size = 0.5) + 
      scale_fill_manual(values = gene_colours) + 
      theme_gerstung(base_size = 6) + 
      theme(legend.position = "none") + 
      xlab("Gene") + 
      ylab("Growth rate") + 
      theme(axis.text.x = element_text(face = "italic")),
    nrow = 1),
  nrow = 2
) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/dnmt3a_vs_tet2_size_growth.pdf",model_id),
         height = 3,width = 6) 

glm(log(maxVAF) ~ b_genetic_mean + b_clone_mean + gene + maxAge,data = Parameters_Age_DNMT3A_TET2) %>%
  summary()

cor.test(Parameters_Age_DNMT3A_TET2$maxVAF,Parameters_Age_DNMT3A_TET2$b_clone_mean,method = "kendall")
```

```{r,fig.height=1.5,fig.width=4}
plot_grid(
  full_data %>%
    subset(grepl("DNMT3A|TET2",Gene)) %>% 
    mutate(Gene = str_match(Gene,"[A-Z0-9]+")) %>%
    group_by(Gene) %>% 
    summarise(N = length(unique(SardID))) %>%
    ggplot(aes(x = Gene,y = N,fill = Gene)) +
    geom_bar(stat = "identity") + 
    scale_fill_manual(values = gene_colours) + 
    scale_y_continuous(expand = c(0,0)) +
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    xlab("Gene") + 
    ylab("Age") + 
    theme(axis.text.x = element_text(face = "italic")),
  full_data %>%
    subset(grepl("DNMT3A|TET2",Gene)) %>% 
    mutate(Gene = str_match(Gene,"[A-Z0-9]+")) %>%
    group_by(SardID,Gene) %>% 
    summarise(Age = min(Age[VAF > 0.005])) %>%
    ggplot(aes(x = Gene,y = Age,fill = Gene)) +
    geom_boxplot(size = 0.5,outlier.size = 0.25) + 
    scale_fill_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    xlab("Gene") + 
    ylab("Age") + 
    theme(axis.text.x = element_text(face = "italic")),
  Parameters_Age_DNMT3A_TET2 %>%
    ggplot(aes(x = gene,y = maxVAF,fill = gene)) +
    geom_boxplot(size = 0.5,outlier.size = 0.25) + 
    scale_y_continuous(trans = 'log10') + 
    scale_fill_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    xlab("Gene") + 
    ylab("Maximum VAF") + 
    theme(axis.text.x = element_text(face = "italic")),
  Parameters_Age_DNMT3A_TET2 %>%
    ggplot(aes(x = gene,y = b_genetic_mean + b_clone_mean,fill = gene)) +
    geom_boxplot(size = 0.5,outlier.size = 0.25) + 
    scale_fill_manual(values = gene_colours) + 
    theme_gerstung(base_size = 6) + 
    theme(legend.position = "none") + 
    xlab("Gene") + 
    ylab("Growth rate") + 
    theme(axis.text.x = element_text(face = "italic")),
  nrow = 1) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/dynamic_coefficients/dnmt3a_vs_tet2_summary.pdf",model_id),
         height = 1.5,width = 4) 
```

## Statistical analysis 

Here we present a more detailed statistical analysis, particularly regarding how much each gene is explained by driver genetics or by the full model we proposed so far. Strikingly, we observe that, while genes such as *U2AF1* are explained almost exclusively by driver genetics, others such as *SRSF2* or *JAK2* only become explainable once we include a term that accounts for variability from the expected driver trajectory.

```{r statistics_model_d,fig.height=2.2,fig.width=1.5}
statistics_data <- statistics_lists$individual_site %>% 
  select(site,sum_stat,gene,individual,sum_stat_genetic,n_tp) %>% 
  mutate(model = "") %>%
  mutate(recurrent = ifelse(site %in% unique_site_multiple,TRUE,FALSE)) %>%
  group_by(model) %>%
  mutate(gene = str_match(gene,"[A-Z0-9]+"),
         truncating = str_match(gene,"[nt]+")) %>%
  mutate(model_label = sprintf("%s\np(s-value < 0.7) = %s",model,round(sum(sum_stat == T)/length(sum_stat),3))) %>%
  mutate(gene_numeric = as.numeric(factor(gene,levels = rev(unique(gene))))) 

gene_numeric_correspondence <- statistics_data %>%
  ungroup() %>%
  select(gene,gene_numeric) %>%
  distinct()

statistics_data_bars <- statistics_data %>%
  ungroup() %>%
  mutate(model_label = sprintf("Trajectories with\nno outliers = %s%%",
                               100*round(sum(sum_stat == T )/length(sum_stat),3))) %>% 
  mutate(good = sum_stat == T,
         not_good = sum_stat == F) %>% 
  gather(key = "key",value = "value",good,not_good) %>%
  mutate(key = factor(key,levels = c("not_good","good"))) %>%
  filter(value == T) %>%
  select(site,sum_stat,gene,individual,gene_numeric,key,value,model_label) %>%
  group_by(gene,model_label) %>%
  mutate(total = length(key)) %>%
  group_by(gene,key,model_label) %>%
  summarise(N = length(site),
            Total = total[1]) %>% 
  mutate(Proportion = N / Total,
         Proportion005 = ifelse(
           key == "good",
           qbeta(p = 0.05,N+1,Total+1-N),
           NA),
          Proportion095 = ifelse(
           key == "good",
           qbeta(p = 0.95,N+1,Total+1-N),
           NA)) %>% 
  mutate(
    Proportion095 = ifelse(Proportion095 < Proportion,
                           Proportion,
                           Proportion095),
    Proportion005 = ifelse(Proportion005 > Proportion,
                           Proportion,
                           Proportion005)
  ) 

statistics_data_bars_simplified <- statistics_data %>%
  ungroup() %>%
  mutate(model_label = sprintf("Trajectories with\nno outliers = %s%%",
                               100*round(sum(sum_stat == T )/length(sum_stat),3))) %>% 
  mutate(good_genetic = sum_stat == T & sum_stat_genetic == T,
         good_clone = sum_stat == T & sum_stat_genetic == F,
         not_good = sum_stat == F) %>% 
  gather(key = "key",value = "value",good_genetic,not_good,"good_clone") %>%
  mutate(key = factor(key,levels = c("not_good","good_clone","good_genetic"))) %>%
  filter(value == T) %>%
  select(site,sum_stat,gene,individual,gene_numeric,key,value,model_label) %>%
  group_by(gene,model_label) %>%
  mutate(total = length(unique(paste(site,individual)))) %>%
  group_by(gene,key,model_label) %>%
  summarise(N = length(site),
            Total = total[1]) %>% 
  mutate(Proportion = N / Total)

gene_order_explained <- statistics_data_bars$gene[statistics_data_bars$key == 'good'][
  order(statistics_data_bars$Proportion[statistics_data_bars$key == 'good'])]

gene_order_explained_simp <- statistics_data_bars_simplified$gene[statistics_data_bars_simplified$key == 'good_genetic'][
  order(statistics_data_bars_simplified$Proportion[statistics_data_bars_simplified$key == 'good_genetic'])]

statistics_data_bars <- statistics_data_bars %>%
  mutate(
    gene = factor(gene,rev(gene_order_explained),ordered = T)
  )

statistics_data_bars_simplified <- statistics_data_bars_simplified %>%
  mutate(
    gene = factor(gene,rev(gene_order_explained_simp),ordered = T)
  )

statistics_data_bars_good <- statistics_data_bars %>%
  subset(key == "good") %>%
  mutate(gene = factor(gene,levels = gene_order_explained))
statistics_data_bars %>% 
  mutate(gene = factor(gene,levels = gene_order_explained)) %>%
  ggplot(aes(x = gene,y = Proportion)) + 
  geom_bar(data = statistics_data_bars_good,
           aes(y = 1),stat = "identity",fill = "grey95",width = 0.8) +
  geom_bar(data = statistics_data_bars_good,
           aes(fill = gene),
           stat="identity",width = 0.8) +
  geom_linerange(aes(ymin = Proportion005-0.0005,ymax = Proportion095+0.0005),
                 colour = "black",size = 0.75) +
  geom_linerange(aes(ymin = Proportion,ymax = Proportion095,colour = gene),
                 size = 0.5) + 
  geom_linerange(aes(ymin = Proportion,ymax = Proportion005),
                 colour = "grey95",
                 size = 0.5) + 
  theme_gerstung(base_size = 6) +
  facet_grid(~ model_label) + 
  #rotate_x_text() + 
  scale_x_discrete(
    expand = c(0.02,0.02)) +
  scale_y_continuous(
    breaks = c(0.00,0.5,1.00),
    expand = c(0.003,0.003)
  ) +
  scale_fill_manual(name = NULL,
                    values = gene_colours,
                    guide = F) +
  scale_colour_manual(name = NULL,
                      values = gene_colours,
                      guide = F) +
  theme(legend.position = "bottom",
        axis.text.y = element_text(face = 'italic'),
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(0,0,0,0),
        legend.key.size = unit(0.25,"cm")) + 
  xlab("Driver gene") + 
  ylab("Fraction of clones growing\nat constant rate") + 
  coord_flip() +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/statistical_analysis/statistics_barplot.pdf",model_id),
         height = 2.4,width = 1.9) 

statistics_data_bars_simplified %>% 
  mutate(
    gene_colour = ifelse(
      key == "not_good",
      "1NONE",
      ifelse(
        key == "good_clone",
        sprintf("%s_lighter",gene),
        as.character(gene)
      )
    )
  ) %>%
  mutate(
    gene_colour = factor(gene_colour,levels = c("1NONE",
                                                names(gene_colours_lighter),
                                                names(gene_colours)))
  ) %>%
  ggplot(aes(x = gene,y = Proportion,fill = gene_colour)) + 
  geom_bar(stat = 'identity',color='black',size = 0.1,
           width = 0.8) +
  theme_gerstung(base_size = 6) +
  facet_grid(~ model_label) + 
  #rotate_x_text() + 
  scale_x_discrete(
    expand = c(0.02,0.02)) +
  scale_y_continuous(
    breaks = c(0.00,0.5,1.00),
    expand = c(0,0,0.003,0)
  ) +
  scale_fill_manual(guide = F,values = c(`1NONE` = "white",gene_colours,gene_colours_lighter)) +
  scale_alpha() +
  theme(legend.position = "bottom",
        axis.text.y = element_text(face = 'italic'),
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(0,0,0,0),
        legend.key.size = unit(0.25,"cm")) + 
  xlab("Driver gene") + 
  ylab("Fraction of trajectories") + 
  coord_flip() +
  #guides(fill = guide_legend(direction = "horizontal")) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/statistical_analysis/statistics_barplot_discriminated.pdf",model_id),
         height = 2.2,width = 1.9) 

statistics_data_bars_simplified %>% 
  select(-Proportion,-model_label) %>%
  spread(key = key,value = N) %>%
  mutate(not_good = ifelse(is.na(not_good),0,not_good)) %>%
  ungroup %>% 
  select(-gene) %>%
  summarise_all(sum) %>%
  transmute(
    `Unexplained trajectories` = not_good / Total,
    `Trajectories explained by genetics alone` = good_genetic / Total,
    `Trajectories explained by the full model` = (good_genetic + good_clone) / Total
  )

stats_nmut <- statistics_data %>%
  merge(n_mut_ind,by = "individual") %>%
  group_by(individual) %>%
  select(site,sum_stat,individual,NMut,Age) %>%
  distinct() %>%
  mutate(gene = str_match(site,"[A-Z0-9]+"))

glm(NMut ~ sum_stat + Age,stats_nmut,family = stats::poisson()) %>%
  summary

stats_nmut %>%
  group_by(NMut) %>% 
  summarise(N=sum(sum_stat),Total=length(sum_stat),
            PropGood = sum(sum_stat)/length(sum_stat),
            PropGood_05 = qbeta(0.05,sum(sum_stat)+1,sum(1-sum_stat)+1),
            PropGood_95 = qbeta(0.95,sum(sum_stat)+1,sum(1-sum_stat)+1)) %>% 
  ggplot(aes(x = NMut,y = PropGood)) +
  geom_bar(stat="identity",fill = "grey80",colour = NA) + 
  geom_errorbar(aes(ymin = PropGood_05,ymax = PropGood_95),
                width = 0.4) + 
  theme_gerstung(base_size = 6) + 
  ylab("Fraction of explained trajectories") + 
  xlab("Number of mutations") + 
  scale_x_continuous(breaks = seq(1,12)) + 
  scale_colour_manual(values = gene_colours,
                      name = NULL,
                      guide = F) + 
  theme(legend.position = "bottom") + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/statistical_analysis/statistics_vs_mut_no.pdf",model_id),
         height = 1.5,width = 2) 
```

```{r explained_n_tp,fig.height=1.5,fig.width=1.4}
statistics_data %>% 
  group_by(n_tp) %>% 
  summarise(N = sum(sum_stat),Total = length(sum_stat),
            q05 = qbeta(0.05,N,Total-N),
            q95 = qbeta(0.95,N,Total-N)) %>% 
  ggplot(aes(x = n_tp,y = N/Total)) + 
  geom_bar(stat = "identity",fill = "grey80") + 
  geom_linerange(aes(ymin = q05,ymax = q95)) + 
  theme_gerstung(base_size = 6) + 
  ylab("Fraction of explained trajectories") + 
  xlab("Number of available\ntimepoints") + 
  scale_y_continuous(expand = c(0,0)) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/statistical_analysis/statistics_vs_n_tp.pdf",model_id),
         height = 1.5,width = 1.4)

stats_nmut %>%
  group_by(NMut) %>% 
  summarise(PropGood = sum(sum_stat)/length(sum_stat),
            PropGood_05 = qbeta(0.05,sum(sum_stat)+1,sum(1-sum_stat)+1),
            PropGood_95 = qbeta(0.95,sum(sum_stat)+1,sum(1-sum_stat)+1)) %>% 
  ggplot(aes(x = NMut,y = PropGood)) +
  geom_bar(stat="identity",fill = "grey80",colour = NA) + 
  geom_errorbar(aes(ymin = PropGood_05,ymax = PropGood_95),
                width = 0.4) + 
  theme_gerstung(base_size = 6) + 
  ylab("Fraction of explained trajectories") + 
  xlab("Number of\nmutations") + 
  scale_x_continuous(breaks = seq(1,12)) + 
  scale_y_continuous(expand = c(0,0)) +
  scale_colour_manual(values = gene_colours,
                      name = NULL,
                      guide = F) + 
  theme(legend.position = "bottom") + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/statistical_analysis/statistics_vs_mut_no.pdf",model_id),
         height = 1.5,width = 1.4) 
```

```{r,fig.height=1.9,fig.width=4}
statistics_data %>% 
  subset(gene %in% c("DNMT3A","TET2","ASXL1","PPM1D","TP53","SF3B1")) %>% 
  group_by(n_tp,gene) %>% 
  summarise(N = sum(sum_stat),Total = length(sum_stat),
            q05 = qbeta(0.05,N,Total-N),
            q95 = qbeta(0.95,N,Total-N)) %>% 
  mutate(gene = factor(gene,levels = c("TET2","DNMT3A","SF3B1","ASXL1","TP53","PPM1D"))) %>% 
  subset(n_tp > 2) %>%
  ggplot(aes(x = factor(n_tp,levels=c(5,4,3,2)))) + 
  geom_bar(aes(y = Total/Total),stat = "identity",fill = "red4",
           colour = "black",size = 0.25) +
  geom_bar(aes(y = N/Total),stat = "identity",fill = "green3",
           colour = "black",size = 0.25) + 
  geom_text(aes(y = 1,label = sprintf("%.2f%% (%s/%s)",N/Total*100,N,Total)),
            size=2,hjust=-0.02) +
  theme_gerstung(base_size = 6) + 
  ylab("Fraction of explained trajectories") + 
  xlab("Number of available timepoints") + 
  scale_y_continuous(expand = c(0,0),breaks = c(0,25,50,75,100)/100) + 
  facet_wrap(~ gene,ncol=2) +
  coord_flip(ylim = c(NA,1.6)) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/statistical_analysis/statistics_vs_n_tp_by_gene.pdf",model_id),
         height = 1.9,width = 4) 
```

```{r save_r_value,fig.width=18,fig.height=10,echo=F}
r_values %>% 
  write.csv(sprintf("data_output/r_values_%s.csv",model_id))
```

These statistical measures also allow us to pinpoint trajectories that do not follow the simple dynamics proposed in this work. We use these to forward explanations regarding the poor performance of the model in these instances and categorise, finding relevant cases for each. Here we present three distinct possibilities for the unexplainability of a trajectory, namely competition (the fitness of a clone suddenly decreases as a fitter clone appears), sampling artefacts or drift (clones are detected at a single timepoint, which may be attributable to drift or to a sequencing artefact) and size-dependent deceleration (clones start growing more slowly as their size increases).

```{r bad_statistics_examples,fig.width=4,fig.height=1.5}
data_bad_statistics <- r_values_full %>%
  merge(select(statistics_data,-gene_numeric),by = c("site","individual")) %>% 
  group_by(individual) %>%
  filter(any(sum_stat == F)) %>% 
  mutate(bad = sum_stat == F) %>%
  mutate(site_ = site) %>% 
  mutate(site = as.character(site)) %>% 
  mutate(site = sapply(site,function(x){
    tmp <- unlist(strsplit(x,'-'))
    paste(tmp[2:length(tmp)],collapse = '-')
  })) %>%
  mutate(site = sapply(site,function(x) unlist(strsplit(x,';'))[1])) %>%
  mutate(site = sapply(site,function(x) {
    tmp <- unlist(strsplit(x,':'))
    return(tmp[length(tmp)])})) %>%
  mutate(site = paste(gene,site,sep='-'))

competition_example <- data_bad_statistics %>% 
  subset(individual == 22109) %>%
  mutate(Label = "Competition")

sequencing_artefact_example <- data_bad_statistics %>%
  subset(individual == 12824) %>%
  mutate(Label = "Sequencing artefact\nor drift")

size_dependent_deceleration_example <- data_bad_statistics %>%
  subset(individual == 11959 & sum_stat == F) %>%
  mutate(Label = "Size-dependent\ndeceleration")

plot_df <- list(
  competition_example,
  sequencing_artefact_example,
  size_dependent_deceleration_example
)
mi <- values_model$min_age
p_list <- lapply(plot_df,function(x) {
  m <- floor(min(x$age))
  age_range <- seq(min(x$age),max(x$age),length.out = 20)
  trajectory_df <- x %>%
    ungroup %>%
    select(site_numeric,gene_numeric,clone_numeric,site,sum_stat) %>%
    distinct %>%
    apply(1,function(x) {
      S <- grab_samples(as.numeric(x[1]),as.numeric(x[2]),
                        as.numeric(x[3]),2500)
      b <- S$site + S$clone + S$gene
      u <- S$offset
      traj <- lapply(age_range-mi,function(t) trajectory_from_parameters_unadjusted(b,u,t)) %>%
        do.call(what = cbind) %>%
        apply(2,function(s) {return(c(mean(s),quantile(s,c(0.05,0.95))))}) %>%
        t
      data.frame(
        trajectory = traj[,1],
        trajectory_005 = traj[,2],
        trajectory_095 = traj[,3],
        time = age_range,
        site = x[4],
        sum_stat = as.logical(x[5])
      ) %>%
        return
    }) %>%
    do.call(what = rbind) 

  ggplot(x) + 
    geom_ribbon(data = trajectory_df,
                aes(x = time,y = trajectory,fill = site,
                    ymin = trajectory_005,ymax = trajectory_095),
                inherit.aes = F,
                position = position_dodge(width = 0.5),
                alpha = 0.3) + 
    geom_line(data = trajectory_df,
              size = 0.25,
              aes(x = time,y = trajectory,colour = site,linetype = sum_stat),
              position = position_dodge(width = 0.5),
              inherit.aes = F) + 
    geom_linerange(
      aes(x = age,
          ymin = qbeta(p=0.05,true+1,coverage+1-true),
          ymax = qbeta(p=0.95,true+1,coverage+1-true),
          colour = site),
      size = 0.25,position = position_dodge(width = 0.5)) + 
    geom_point(aes(x = age,y = true/coverage,colour = site),
               position = position_dodge(width = 0.5),
               size=0.5) + 
    scale_y_continuous(trans = 'log10',
                       breaks = c(0.005,0.01,0.03,0.1,0.3,0.5),
                       labels = sprintf("%s%%",c(0.005,0.01,0.03,0.1,0.3,0.5)*100)) +
    scale_x_continuous(breaks = seq(m,m+50,by=3)) +
    theme_gerstung(base_size = 6) + 
    scale_colour_lancet(name = NULL) + 
    scale_fill_lancet(name = NULL) +
    xlab("Age") + 
    ylab("VAF") + 
    facet_wrap(~ Label) +
    theme(legend.position = "bottom",
          legend.key.size = unit(0,"cm"),
          legend.margin = margin(),legend.spacing = unit(0,"cm")) + 
    guides(colour = guide_legend(nrow = 2,byrow = T)) + 
    coord_cartesian(ylim = c(0.002,NA))
}) 

plot_grid(plotlist = p_list,nrow=1,align = T,axis = "hv") + 
  ggsave(useDingbats=F,sprintf("figures/%s/inferred_trajectories/bad_trajectories_subset.pdf",model_id),
         height = 1.5,
         width = 4)
```

```{r saving_everything_exporting}
save(r_values,statistics_data,Parameters_Age,r_values_full,
     file = sprintf('data_output/vaf_modelling_coefficients_%s.Rdata',model_id))
```

# Post-hoc analysis

## Is the U2AF1 and SRSF2 prevalence we observe in our cohort normal or explained by the old age?

```{r}
n_sites <- Parameters_Age %>% 
  subset(site == "SRSF2-P95H" | gene == "U2AF1") %>%
  subset(sum_stat == T) %>%
  mutate(gene = ifelse(gene == "SRSF2","SRSF2-P95H",gene)) %>% 
  group_by(gene) %>%
  summarise(N = length(unique(individual)))

n_muts_uniform <- estimates_from_samples %>% 
  subset(gene %in% c("SRSF2-P95H","U2AF1")) %>% 
  subset(!is.na(values)) %>%
  group_by(gene) %>%
  summarise(Q025 = quantile(values,0.025),
            Q975 = quantile(values,0.975)) %>%
  merge(n_sites,by = "gene") %>%
  mutate(rate = N / (Q975 - Q025)) %>%
  mutate(over_70_years = rate * 70)
```

### *U2AF1*

To assess this we need to get the prevalence of U2AF1 of myeloid origin, which we obtain by marginalising over AML, CMML, MDS and CH.

The incidence of U2AF1 mutations in our cohort is quite higher than what would be expected in a sample of individuals over 70. As such, we can, to some extent, validate our findings regarding the late appearance of most U2AF1 clones.

```{r checking_U2AF1_prevalence,fig.height=1,fig.width=3}
u2af1_CH_list <- list(
  "Jaiswal 2014" = c(n_mutations = 0,ch_individuals = 805,
                     total_individuals = 17182),
  "Genovese 2014" = c(n_mutations = 0,ch_individuals = 2113,
                      total_individuals = 11845),
  "McKerrell 2015" = c(n_mutations = NA,ch_individuals = 105,
                       total_individuals = 4219),
  "Zink 2017" = c(n_mutations = 0, ch_individuals = 1403,
                  total_individuals = 11262),
  "Acuna-Hidalgo 2017" = c(n_mutations = 1, ch_individuals = 192,
                           total_individuals = 2007),
  # "Coombs 2017" = c(n_mutations = 1, ch_individuals = 192,
  #                   total_individuals = 6898), # therapy-related
  "Young 2016" = c(n_mutations = 0,ch_individuals = 19,
                   total_individuals = 20),
  "Young 2019" = c(n_mutations = 1,ch_individuals = 63,
                   total_individuals = 69), # 3 more U2AF1 mutations in AML
  # "Desai 2018" = c(n_mutations = 0,ch_individuals = 56,
  #                  total_individuals = 181), # 6 U2AF1 mutations in MAL
  "Midic 2020" = c(n_mutations = 1,ch_individuals = 13,
                   total_individuals = 50)
)

p_u2af1_c_ch <- sum(sapply(u2af1_CH_list,function(x) x[1]),na.rm = T) / sum(sapply(u2af1_CH_list,function(x) x[2]))
p_u2af1_c_mds <- 21.7/100 # Wang 2020 (Sci. rep.)
p_u2af1_c_cml <- 13.1/100 # Abu Kar 2013 (Haematologica)
p_u2af1_c_aml <- 9/100 # Bamopoulos 2020 (Leukaemia)
p_ch <- sum(sapply(u2af1_CH_list,function(x) x[2]),na.rm = T) / sum(sapply(u2af1_CH_list,function(x) x[3]))
p_mds <- 0.0002
p_cmml <- 0.34 / 1e5
p_aml <- 19.8 / 1e5

# marginalise P(U2AF1)
o <- p_u2af1_c_ch * p_ch + p_u2af1_c_mds * p_mds + p_u2af1_c_aml * p_aml + p_u2af1_c_cml * p_cmml

S_u2af1 <- load_data() %>% 
  select(Gene,SardID) %>% 
  distinct %>% 
  summarise(`U2AF1 in our cohort` = length(unique(SardID[Gene == 'U2AF1'])) / length(unique(SardID)),
            `U2AF1 mutation in CH (other studies)` = p_u2af1_c_ch,
            `U2AF1 mutation in MDS` = p_u2af1_c_mds,
            `U2AF1 mutation in CMML` = p_u2af1_c_cml,
            `U2AF1 mutation in AML` = p_u2af1_c_aml,
            `U2AF1 mutation of myeloid origin` = o,
            `U2AF1 in our cohort at constant mutation rate` = n_muts_uniform$over_70_years[2] / length(unique(SardID))) %>%
  gather(key = "Condition",value = "Expected prevalence") %>%
  mutate(Condition = factor(
    Condition,
    levels = c("U2AF1 in our cohort at constant mutation rate",
               "U2AF1 in our cohort","U2AF1 mutation of myeloid origin",
               "U2AF1 mutation in CH (other studies)","U2AF1 mutation in MDS",
               "U2AF1 mutation in CMML","U2AF1 mutation in AML")))

S_u2af1

S_u2af1 %>%
  ggplot(aes(x = Condition,y = `Expected prevalence`)) +
  geom_bar(stat = "identity",fill = gene_colours[["U2AF1"]]) + 
  coord_flip() + 
  theme_gerstung(base_size = 6) + 
  geom_text(aes(label = sprintf("%s%%",round(`Expected prevalence`*100,2))),
            size = 2.6,hjust = -0.1) +
  xlab("") +
  scale_y_continuous(expand = c(0,0,0.4,0),
                     labels = function(x) sprintf("%s%%",x*100)) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/prevalence/u2af1.pdf",model_id),
         height = 1,width = 3)
```

### *SRSF2*-P95H

To assess this we need to get the prevalence of SRSF2-P95H in blood cells, which we obtain by marginalising over AML, CMML, MDS and CH.

In Genovese *et al.*, 3 individuals with SRSF2-P95H CH contracted some form of myeloid malignancy upon follow-up (ages: 73, 71, 85), one of which died of leukaemia ~2 years later (age: 85).

In Desai *et al.*, 13 individuals with SRSF2-P95H CH contracted AML upon follow-up. No individuals which did not develop AML had SRSF2-P95H mutations.

The incidence of SRSF2-P95H mutations in our cohort, much like that of U2AF1 mutations, is quite higher than what would be expected in a sample of individuals over 70. As such, we can, to some extent, validate our findings regarding the late appearance of most SRSF2-P95H clones.

```{r checking_SRSF2_P95H_prevalence,fig.height=1,fig.width=3}
srsf2_p95h_CH_list <- list(
  "Jaiswal 2014" = c(n_mutations = 4,ch_individuals = 805,
                     total_individuals = 17182),
  "Genovese 2014" = c(n_mutations = 2,ch_individuals = 2113,
                      total_individuals = 11845), 
  "McKerrell 2015" = c(n_mutations = 5,ch_individuals = 105,
                       total_individuals = 4219),
  "Zink 2017" = c(n_mutations = 2, ch_individuals = 1403,
                  total_individuals = 11262),
  "Acuna-Hidalgo 2017" = c(n_mutations = 1, ch_individuals = 192,
                           total_individuals = 2007),
  "Young 2016" = c(n_mutations = NA,ch_individuals = 19,
                   total_individuals = 20),
  "Young 2019" = c(n_mutations = 1,ch_individuals = 63,
                   total_individuals = 69), 
  # "Desai 2018" = c(n_mutations = 0,ch_individuals = 56,
  #                  total_individuals = 181),
  "Midic 2020" = c(n_mutations = 0,ch_individuals = 13,
                   total_individuals = 50))

srsf2_p95h_CH_ages <- rbind(
  data.frame(ages = rep(NA,4),
             study = "Jaiswal 2014",
             average_age = 57),
  data.frame(ages = rep(NA,2),
             study = "Genovese 2014",
             average_age = 55),
  data.frame(ages = c(73,77,80,82,88),
             study = "McKerrell 2015",
             average_age = NA),
  data.frame(ages = c(68,70),
             study = "Zink 2017",
             average_age = NA),
  data.frame(ages = c(65),
             study = "Acuna-Hidalgo 2017",
             average_age = NA),
  data.frame(ages = c(NA),
             study = "Young 2016",
             average_age = NA),
  data.frame(ages = c(NA),
             study = "Young 2019",
             average_age = 75),
  data.frame(ages = c(NA),
             study = "Desai 2018",
             average_age = NA),
  data.frame(ages = c(NA),
             study = "Midic 2020",
             average_age = NA)
)

p_srsf2_p95h_c_ch <- sum(sapply(srsf2_p95h_CH_list,function(x) x[1]),na.rm = T) / sum(sapply(srsf2_p95h_CH_list,function(x) x[2]))
p_srsf2_p95h_c_mds <- 11.5/100 # Jafari 2018 (for SRSF2 mutations)
p_srsf2_p95h_c_cml <- 39.8/100 # Jafari 2018 (for SRSF2 mutations)
p_srsf2_p95h_c_aml <- 69/1119 # Bamopoulos 2020
p_srsf2_p95h_c_aml <- 10/100 # Bamopoulos 2020 (SRSF2 mutations)
p_ch <- sum(sapply(srsf2_p95h_CH_list,function(x) x[2]),na.rm = T) / sum(sapply(srsf2_p95h_CH_list,function(x) x[3]))
p_mds <- 0.0002
p_cmml <- 0.34 / 1e5
p_aml <- 19.8 / 1e5

# marginalise P(SRSF2_p95h)
o <- p_srsf2_p95h_c_ch * p_ch + p_srsf2_p95h_c_mds * p_mds + p_srsf2_p95h_c_aml * p_aml + p_srsf2_p95h_c_cml * p_cmml

S <- load_data() %>% 
  select(Gene,SardID,amino_acid_change) %>% 
  distinct %>% 
  summarise(`SRSF2-P95H in our cohort` = length(unique(SardID[amino_acid_change == 'SRSF2-P95H'])) / length(unique(SardID)),
            `SRSF2-P95H mutation in CH (other studies)` = p_srsf2_p95h_c_ch,
            `SRSF2 mutation in MDS` = p_srsf2_p95h_c_mds,
            `SRSF2 mutation in CMML` = p_srsf2_p95h_c_cml,
            `SRSF2 mutation in AML` = p_srsf2_p95h_c_aml,
            `SRSF2 mutation of myeloid origin` = o,
            `SRSF2-P95H in our cohort at constant mutation rate` = n_muts_uniform$over_70_years[1] / length(unique(SardID))) %>%
  gather(key = "Condition",value = "Expected prevalence") %>%
  mutate(Condition = factor(
    Condition,
    levels = c("SRSF2-P95H in our cohort at constant mutation rate",
               "SRSF2-P95H in our cohort","SRSF2 mutation of myeloid origin",
               "SRSF2-P95H mutation in CH (other studies)","SRSF2 mutation in MDS",
               "SRSF2 mutation in CMML","SRSF2 mutation in AML")))

S
```

### Figure for *U2AF1* and *SRSF2*-P95H

```{r,fig.height=2,fig.width=4}
plot_grid(
  S %>% 
    ggplot(aes(x = Condition,y = `Expected prevalence`)) +
    geom_bar(stat = "identity",fill = gene_colours[["SRSF2"]]) + 
    coord_flip() + 
    theme_gerstung(base_size = 6) + 
    geom_text(aes(label = sprintf("%s%%",round(`Expected prevalence`*100,2))),
              size = 2.6,hjust = -0.1) +
    xlab("") +
    scale_y_continuous(expand = c(0,0,0.25,0),
                       labels = function(x) sprintf("%s%%",x*100)),
  S_u2af1 %>% 
    ggplot(aes(x = Condition,y = `Expected prevalence`)) +
    geom_bar(stat = "identity",fill = gene_colours[["U2AF1"]]) + 
    coord_flip() + 
    theme_gerstung(base_size = 6) + 
    geom_text(aes(label = sprintf("%s%%",round(`Expected prevalence`*100,2))),
              size = 2.6,hjust = -0.1) +
    xlab("") +
    scale_y_continuous(expand = c(0,0,0.25,0),
                       labels = function(x) sprintf("%s%%",x*100)),
  align = "v",
  ncol = 1) + 
  ggsave(useDingbats=FALSE,sprintf("figures/%s/prevalence/u2af1_srsf2.pdf",model_id),
         height = 2,width = 4)
```

### Considering age

The problem with the previous analysis is that we do not account for the dependence that prevalence has on age. As such, in this next code block we simply calculate what the prevalence of leukaemia and AML would be after 60 years of age and compare it with the prevalence of *U2AF1* and *SRSF2*-P95H - this gives us a maximum (yet unlikely) prevalence possible for *U2AF1* and *SRSF2-P95H* mutations in blood cells.

This analysis shows that even if all the leukaemias were characterised by mutations in either *SRSF2*-P95H or *U2AF1* we would never see a prevalence this high for either driver, implying that the main cause for this high prevalence is age rather than myeloid disease and that the reason why we observe it not connected with other cases of mutations in these drivers being flagged as leukaemia.

```{r}
aml_incidence <- read_csv("data/aml_incidence.csv",progress = F)[1:12,] %>%
  summarise(prevalence_after_60 = sum(female_incidence + male_incidence)/100000)

leukaemia_incidence <- read_csv("data/leukaemia_incidence.csv",progress = F)[1:12,] %>%
  summarise(prevalence_after_60 = sum(female_incidence + male_incidence)/100000)

data.frame(`Prevalence` = c(round(100*S_u2af1[c(1,7),2],2),
                            round(100*S[c(1,7),2],2),
                            100*aml_incidence$prevalence_after_60,100*leukaemia_incidence$prevalence_after_60) %>%
             sprintf(fmt = "%s%%"),
           id = c("U2AF1 in our cohort",
                  "U2AF1 in our cohort at constant mutation rate",
                  "SRSF2-P95H in our cohort",
                  "SRSF2-P95H in our cohort at constant mutation rate",
                  "AML incidence (at 60yo)",
                  "All leukaemia incidence (at 60yo)"))
```

## Is growth rate associated with dNdS in AML and MDS?

```{r,fig.height=1.5,fig.width=2.5}
all_samples <- sampling_idxs %>%
  group_by(site) %>%
  mutate(re = length(unique(individual)) > 1) %>%
  subset(re == T) %>% 
  subset(sum_stat == T) %>%
  apply(1,function(x) {
    site_numeric <- as.numeric(x[3])
    gene_numeric <- as.numeric(x[4])
    clone_numeric <- as.numeric(x[5])
    site <- x[2]
    S <- grab_samples(site_numeric,gene_numeric,clone_numeric,2500)
    return(data.frame(samples = S$site + S$clone + S$gene,site = site))
  }) %>%
  do.call(what = rbind) %>%
  mutate(gene = str_match(site,"[A-Z0-9]+"),
         aachange = gsub("-","",str_match(site,"-.+"))) %>%
  mutate(site = paste(gene,aachange,sep = "-")) %>% 
  group_by(site) %>%
  summarise(mean = mean(samples),
            s = sd(samples,0.025))

dnds <- read.csv("data/dnds_aml_mds.csv") %>%
  mutate(site = paste(gene,gsub("\\*","X",aachange),sep = "-"))

dd <- all_samples %>% 
  merge(dnds,by = "site",all = F) %>% 
  as.data.frame

dd %>% 
  ggplot(aes(y = dnds,x = exp(mean)-1,colour = gene,shape = condition)) + 
  geom_point(size = 1) + 
  geom_errorbarh(aes(xmin = exp(mean-s*1.96)-1,xmax = exp(mean+s*1.96)-1),size = 0.1) + 
  geom_linerange(aes(ymin = CI_low,ymax = CI_high),size = 0.1) + 
  geom_smooth(aes(group = NA),method = "lm",colour = "black",alpha = 0.2,size = 0.5) + 
  scale_colour_manual(values = gene_colours,guide = F) + 
  ylab("dN/dS") +
  xlab("Annual driver growth rate in CH") +
  scale_x_continuous(labels = function(x) sprintf("%s%%",x*100)) + 
  scale_y_continuous(trans = 'log10') + 
  theme_gerstung(base_size = 6) +
  scale_shape(guide = F) +
  ggsave(useDingbats=FALSE,
         filename = sprintf("figures/%s/dynamic_coefficients/dnds_association.pdf",model_id),
         width = 2.5,height = 1.5)

lm(formula = dnds ~ mean,data = dd) %>%
  summary

lm(formula = dnds ~ mean + condition,data = dd) %>%
  summary
```

## Is growth rate associated with AML risk?

```{r,fig.height=1.5,fig.width=2.5}
all_samples <- sampling_idxs %>%
  apply(1,function(x) {
    site_numeric <- as.numeric(x[3])
    gene_numeric <- as.numeric(x[4])
    clone_numeric <- as.numeric(x[5])
    site <- x[2]
    S <- grab_samples(site_numeric,gene_numeric,clone_numeric,2500)
    return(data.frame(samples = S$site + S$clone + S$gene,site = site))
  }) %>%
  do.call(what = rbind) %>%
  mutate(gene = str_match(site,"[A-Z0-9]+")) %>%
  group_by(gene) %>%
  summarise(mean = mean(samples),
            s = sd(samples,0.025))

all_gene_vaf <- full_data %>%
  group_by(SardID,amino_acid_change) %>%
  mutate(medianVAF = median(VAF[VAF >= 0.005]),
         randomVAF = sample(VAF[VAF >= 0.005],1)) %>%
  group_by(SardID) %>% 
  filter(Age %in% range(Age)) %>%
  mutate(Age = ifelse(Age == min(Age),"Entry","Exit")) %>%
  select(-relative_timepoint,-Phase,-MUTcount_Xadj,-WTcount,-TOTALcount) %>%
  spread(key = Age,value = VAF) %>%
  mutate(Gene = str_match(Gene,'[A-Z0-9]+')) %>%
  group_by(Gene) %>%
  summarise(meanlogVAF = log(median(Entry,na.rm=T)),
            medianMedianVAF = log(median(medianVAF)),
            medianRandomVAF = log(median(randomVAF)),
            meanlogVAFStrict = log(median(Entry[Entry == 0 | Entry >= 0.005],na.rm=T)),
            meanlogVAFLate = log(median(Exit,na.rm=T)),
            IncipientRatio = sum(Entry == 0,na.rm=T)/length(Entry)) %>%
  select(gene = Gene,meanlogVAF,meanlogVAFLate,IncipientRatio,medianMedianVAF,medianRandomVAF)

dd <- all_samples %>% 
  merge(read.csv("data/aml-risk-coef.csv"),by = "gene",all = F) %>% 
  merge(all_gene_vaf,by = "gene",all = F) %>%
  as.data.frame

dd %>% 
  ggplot(aes(y = coef,x = exp(mean)-1,colour = gene)) + 
  geom_point(size = 1) + 
  geom_errorbarh(aes(xmin = exp(mean-s*1.96)-1,xmax = exp(mean+s*1.96)-1),size = 0.1) + 
  geom_linerange(aes(ymin = coef - sd*1.96,ymax = coef + sd*1.96),size = 0.1) + 
  geom_smooth(method = "lm",colour = "black",alpha = 0.2,size = 0.5) + 
  scale_colour_manual(values = gene_colours,guide = F) + 
  xlab("Annual driver growth rate in CH") +
  ylab("AML risk") +
  scale_x_continuous(labels = function(x) sprintf("%s%%",x*100)) + 
  theme_gerstung(base_size = 6) +
  coord_cartesian(xlim = c(NA,0.5)) +
  ggsave(useDingbats=FALSE,
         filename = sprintf("figures/%s/dynamic_coefficients/aml_risk.pdf",model_id),
         width = 2.5,height = 1.5)

j <- cor.test(dd$mean,dd$coef)

sprintf("R^2 = %s (p-value = %s)",round(j$estimate^2,4),round(j$p.value,4))
```

### Can VAF during study be a possible confounder?

In Abelson & Collord *et al.* (2018), clone size has a sizeable impact on AML risk, so it is important to assess this.

```{r,fig.height=2.5,fig.width=6}
n <- nrow(dd)
correlations <- cor(dd[,c(2,5,11,12)]) %>%
  as.data.frame
colnames(correlations) <- c("Average growth\nperyear","AML risk",
                            "VAF at\nstudy entry","VAF at\nstudy exit")
correlations$ids <- c("Average growth\nperyear","AML risk",
                      "VAF at\nstudy entry","VAF at\nstudy exit")
correlations_long <- correlations %>%
  gather(key = "key",value = "R",-ids)
correlations_long <- correlations_long %>%
  mutate(t.val = abs(R)*sqrt(n-2)/(sqrt(1-(R^2)))) %>%
  mutate(p.val = (1-pt(t.val,df = n-2))*2) %>% 
  mutate(sig = cut(p.val,c(0,0.001,0.005,0.05,1),c("***","**","*","")))
  
plot_grid(
  plot_grid(
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),y = exp(meanlogVAFLate),colour = gene)) +
      geom_smooth(size = 0.25,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black",
                  linetype = 3) +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("Median VAF at study exit") + 
      scale_y_continuous(trans = 'log10') + 
      scale_x_continuous(trans = 'log10') +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),y = mean,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("Average growth per year") + 
      scale_x_continuous(trans = 'log10') +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAFLate),mean,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study exit") + 
      ylab("Average growth per year") + 
      scale_x_continuous(trans = 'log10') +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("AML risk") + 
      scale_x_continuous(trans = 'log10') +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAFLate),coef,colour = gene)) +
      geom_smooth(size = 0.25,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black",
                  linetype = 3) +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study exit") + 
      ylab("AML risk") + 
      scale_x_continuous(trans = 'log10') +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = mean,coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Average growth per year") + 
      ylab("AML risk") + 
      scale_colour_manual(values = gene_colours)),
  plot_grid(
    NULL,
    correlations_long %>%
      ggplot(aes(x = ids,y = key,fill = R)) + 
      geom_tile(colour = "black",size = 0.25) + 
      geom_text(aes(label = sig),vjust = 0.75) + 
      theme_gerstung(base_size = 6) + 
      scale_fill_gradient2(low = "blue",high = "red4",mid = "white",midpoint = 0,
                           limits = c(-1,1)) + 
      theme(axis.title = element_blank(),
            legend.key.width = unit(0.1,"cm"),
            legend.key.height = unit(0.4,"cm")) + 
      scale_x_discrete(expand = c(0,0)) + 
      scale_y_discrete(expand = c(0,0)) + 
      rotate_x_text(angle = 90),
    NULL,ncol = 1,rel_heights = c(0.3,1,0.3)),
  nrow = 1,rel_widths = c(2,1)) + 
  ggsave(useDingbats=FALSE,
         filename = sprintf("figures/%s/dynamic_coefficients/aml_risk_full.pdf",model_id),
         width = 6,height = 2.5)
```

```{r,fig.height=1.3,fig.width=6}
n <- nrow(dd)
correlations <- cor(dd[,c(2,5,11)]) %>%
  as.data.frame
colnames(correlations) <- c("Average growth\nperyear","AML risk",
                            "VAF at\nstudy entry")
correlations$ids <- c("Average growth\nperyear","AML risk",
                      "VAF at\nstudy entry")
correlations_long <- correlations %>%
  gather(key = "key",value = "R",-ids)
correlations_long <- correlations_long %>%
  mutate(t.val = abs(R)*sqrt(n-2)/(sqrt(1-(R^2)))) %>%
  mutate(p.val = (1-pt(t.val,df = n-2))*2) %>% 
  mutate(sig = cut(p.val,c(0,0.001,0.005,0.05,1),c("***","**","*","")))

plot_grid(
  plot_grid(
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),y = mean,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("Average growth per year") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("A") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("AML risk") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("B") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = mean,coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Average growth per year") + 
      ylab("AML risk") + 
      ggtitle("C") +
      scale_colour_manual(values = gene_colours),
    nrow = 1),
  plot_grid(
    correlations_long %>%
      ggplot(aes(x = ids,y = key,fill = R)) + 
      geom_tile(colour = "black",size = 0.25) + 
      geom_text(aes(label = sig),vjust = 0.75) + 
      theme_gerstung(base_size = 6) + 
      scale_fill_gradient2(low = "blue",high = "red4",mid = "white",midpoint = 0,
                           limits = c(-1,1)) + 
      theme(axis.title = element_blank(),
            legend.key.width = unit(0.1,"cm"),
            legend.key.height = unit(0.25,"cm")) + 
      scale_x_discrete(expand = c(0,0)) + 
      scale_y_discrete(expand = c(0,0)) + 
      ggtitle("D") +
      rotate_x_text(angle = 90),ncol = 1),
  nrow = 1,rel_widths = c(2,1)) + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/aml_risk_early_vaf.png",model_id),
         width = 6,height = 1.3)
```

```{r,fig.height=2.6,fig.width=6}
n <- nrow(dd)
correlations <- cor(dd[,c(2,5,11,14)]) %>%
  as.data.frame
colnames(correlations) <- c("Average growth\nperyear","AML risk",
                            "VAF at\nstudy entry","Median\naverage VAF")
correlations$ids <- c("Average growth\nperyear","AML risk",
                      "VAF at\nstudy entry","Median\naverage VAF")
correlations_long <- correlations %>%
  gather(key = "key",value = "R",-ids)
correlations_long <- correlations_long %>%
  mutate(t.val = abs(R)*sqrt(n-2)/(sqrt(1-(R^2)))) %>%
  mutate(p.val = (1-pt(t.val,df = n-2))*2) %>% 
  mutate(sig = cut(p.val,c(0,0.001,0.005,0.05,1),c("***","**","*","")))

plot_grid(
  plot_grid(
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),y = exp(medianMedianVAF),colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("Median average VAF") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("A") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),y = mean,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("Average growth per year") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("B") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("AML risk") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("C") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(medianMedianVAF),y = mean,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median average VAF") + 
      ylab("Average growth per year") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("D") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(medianMedianVAF),y = coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median average VAF") + 
      ylab("AML risk") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("E") +
      scale_colour_manual(values = gene_colours),
    NULL,
    nrow = 2),
  plot_grid(
    ggplot() + theme_nothing(),
    correlations_long %>%
      ggplot(aes(x = ids,y = key,fill = R)) + 
      geom_tile(colour = "black",size = 0.25) + 
      geom_text(aes(label = sig),vjust = 0.75) + 
      theme_gerstung(base_size = 6) + 
      scale_fill_gradient2(low = "blue",high = "red4",mid = "white",midpoint = 0,
                           limits = c(-1,1)) + 
      theme(axis.title = element_blank(),
            legend.key.width = unit(0.1,"cm"),
            legend.key.height = unit(0.25,"cm")) + 
      scale_x_discrete(expand = c(0,0)) + 
      scale_y_discrete(expand = c(0,0)) + 
      ggtitle("F") +
      rotate_x_text(angle = 90),ncol = 1,rel_heights = c(0.5,1)),
  nrow = 1,rel_widths = c(2,1)) + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/aml_risk_early_vaf_median_vaf.png",model_id),
         width = 6,height = 2.6)
```

```{r,fig.height=2.6,fig.width=6}
n <- nrow(dd)
correlations <- cor(dd[,c(2,5,11,15)]) %>%
  as.data.frame
colnames(correlations) <- c("Average growth\nperyear","AML risk",
                            "VAF at\nstudy entry","Median\nrandom VAF")
correlations$ids <- c("Average growth\nperyear","AML risk",
                      "VAF at\nstudy entry","Median\nrandom VAF")
correlations_long <- correlations %>%
  gather(key = "key",value = "R",-ids)
correlations_long <- correlations_long %>%
  mutate(t.val = abs(R)*sqrt(n-2)/(sqrt(1-(R^2)))) %>%
  mutate(p.val = (1-pt(t.val,df = n-2))*2) %>% 
  mutate(sig = cut(p.val,c(0,0.001,0.005,0.05,1),c("***","**","*","")))

plot_grid(
  plot_grid(
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),y = exp(medianRandomVAF),colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("Median random VAF") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("A") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),y = mean,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("Average growth per year") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("B") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(meanlogVAF),coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median VAF at study entry") + 
      ylab("AML risk") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("C") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(medianRandomVAF),y = mean,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median random VAF") + 
      ylab("Average growth per year") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("D") +
      scale_colour_manual(values = gene_colours),
    dd %>%
      ggplot(aes(x = exp(medianRandomVAF),y = coef,colour = gene)) +
      geom_smooth(size = 0.5,alpha = 0.5,method = "lm",formula = y ~ x,colour = "black") +
      geom_point(size = 0.5) + 
      theme_gerstung(base_size = 6) +
      theme(legend.position = "none") + 
      xlab("Median random VAF") + 
      ylab("AML risk") + 
      scale_x_continuous(trans = 'log10') +
      ggtitle("E") +
      scale_colour_manual(values = gene_colours),
    NULL,
    nrow = 2),
  plot_grid(
    ggplot() + theme_nothing(),
    correlations_long %>%
      ggplot(aes(x = ids,y = key,fill = R)) + 
      geom_tile(colour = "black",size = 0.25) + 
      geom_text(aes(label = sig),vjust = 0.75) + 
      theme_gerstung(base_size = 6) + 
      scale_fill_gradient2(low = "blue",high = "red4",mid = "white",midpoint = 0,
                           limits = c(-1,1)) + 
      theme(axis.title = element_blank(),
            legend.key.width = unit(0.1,"cm"),
            legend.key.height = unit(0.25,"cm")) + 
      scale_x_discrete(expand = c(0,0)) + 
      scale_y_discrete(expand = c(0,0)) + 
      ggtitle("F") +
      rotate_x_text(angle = 90),ncol = 1,rel_heights = c(0.5,1)),
  nrow = 1,rel_widths = c(2,1)) + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/aml_risk_early_vaf_random_vaf.png",model_id),
         width = 6,height = 2.6)
```

# Session details

```{r session_details}
cat(system("git log -n1", intern=TRUE), sep="\n")
```

```{r session_objects}
l <- ls()
data.frame(variable=l, Reduce("rbind",lapply(l, function(x) data.frame(classes=paste(class(get(x)),collapse = ','), size=format(object.size(get(x)), units="auto")))))
```

```{r session_packages}
sessionInfo()
```