Notebook_HistoricalGrowth_PoorFits.Rmd

---
title: "The longitudinal dynamics and natural history of clonal haematopoiesis"
subtitle: "Historical growth and its determinants; poor fits and their determinants"
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)
```

# Historical growth effect

We here define historical growth effect as the growth effect calculated only from posterior samples which yield age at onset estimates which are larger than 0.

```{r setup_data}
pop_size <- 0.5e5
source("Scripts/vaf_dynamics_functions.R")
source("Scripts/prepare_data.R")
model_file_name <- "models/model_ch.RDS"
dir.create("figures/model_ch/other_associations")
set.seed(42)
model_id <- "model_ch"
load(file = sprintf('data_output/vaf_modelling_coefficients_%s.Rdata',model_id))
values_model <- readRDS(file = model_file_name)

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()

sampling_idxs <- r_values_full %>%
  subset(b_site_095 + b_gene_095 + b_clone_095 > 0) %>%
  group_by(
    individual,site,
    site_numeric,
    gene_numeric,
    clone_numeric
  ) %>% 
  summarise(
    minAge = min(age[which(round(true/coverage,4) >= 0.005)]),
    maxAge = max(age),
    maxVAF = max(true/coverage)
  ) %>%
  arrange(site)

c_names <- colnames(values_model$draws$`11`)

good_sites_individual <- statistics_data %>%
  ungroup %>% 
  select(site,individual,sum_stat) %>%
  filter(sum_stat) %>%
  select(-sum_stat)
```

```{r calculate_historical_estimates}
historical_estimates <- sampling_idxs %>%
  apply(1,function(x) {
    individual <- x[1]
    site_numeric <- as.numeric(x[3])
    gene_numeric <- as.numeric(x[4])
    clone_numeric <- as.numeric(x[5])
    minAge <- as.numeric(x[6])
    
    if (is.na(site_numeric)) {
      b_site <- 0} else {
        b_site <- site_samples[,site_numeric]
      }
    b_gene <- gene_samples[,gene_numeric]
    b_clone <- clone_samples[,clone_numeric]
    b <- b_gene + b_clone + b_site
    u_values <- offset_samples[,clone_numeric]
    ages_at_onset <- t0_adjusted(u_values,b,g=2,n=pop_size) + values_model$min_age
    ages_at_onset_corrected <- ifelse(ages_at_onset < -1,-1,ages_at_onset)
    ages_at_onset_corrected <- ifelse(ages_at_onset > minAge,minAge,ages_at_onset)
    biologically_feasible <- ages_at_onset > -1 & ages_at_onset < minAge
    historical_b <- b[biologically_feasible]
    Q <- quantile(historical_b,c(0.05,0.5,0.95),na.rm = T)
    Q_observed <- quantile(b,c(0.05,0.5,0.95),na.rm = T)
    names(Q) <- c("Q05","Q50","Q95")
    return(data.frame(
                      ages_Q50 = quantile(ages_at_onset_corrected,0.5,na.rm = T)[1],
                      observed_growth_q05=Q_observed[1],
                      observed_growth_q50=Q_observed[2],
                      observed_growth_q95=Q_observed[3],
                      observed_growth_mean=mean(b,na.rm=T),
                      historical_growth_q05=Q[1],
                      historical_growth_q50=Q[2],
                      historical_growth_q95=Q[3],
                      historical_growth_mean=mean(historical_b,na.rm=T)))
  }) %>%
  do.call(what=rbind) 
historical_estimates <- cbind(
  as.data.frame(sampling_idxs),as.data.frame(historical_estimates))
historical_estimates$gene <- str_match(historical_estimates$site,"[A-Z0-9]+")
historical_estimates <- historical_estimates %>%
  filter(
    paste(individual,as.character(site)) %in% paste(
      good_sites_individual$individual,good_sites_individual$site))
historical_estimates %>%
  select(SardID = individual,
         Site = site,
         Gene = gene,
         ObservedGrowthQ05 = observed_growth_q05,
         ObservedGrowthQ50 = observed_growth_q50,
         ObservedGrowthQ95 = observed_growth_q95,
         ObservedGrowthMean = observed_growth_mean,
         HistoricalGrowthQ05 = historical_growth_q05,
         HistoricalGrowthQ50 = historical_growth_q50,
         HistoricalGrowthQ95 = historical_growth_q95,
         HistoricalGrowthMean = historical_growth_mean) %>% 
  write.csv("data_output/historical_growth.csv")
```

```{r,fig.height=0.65,fig.width=1.4}
idxs <- sampling_idxs %>% 
  subset(individual == 4088 & site == "DNMT3At-E578X")

tmp_data <- full_data %>%
  subset(SardID == idxs$individual & amino_acid_change == idxs$site)

#idxs <- sampling_idxs[1,]

S <- data.frame(
  u = offset_samples[,idxs$clone_numeric],
  b_gene = gene_samples[,idxs$gene_numeric],
  b_site = site_samples[,idxs$site_numeric],
  b_clone = clone_samples[,idxs$clone_numeric],
  idx = seq(1,nrow(clone_samples)),
  maxAge = idxs$maxAge) %>%
  apply(1,function(x) {
    u <- as.numeric(x[1])
    b_gene <- as.numeric(x[2])
    b_site <- as.numeric(x[3])
    b_site <- ifelse(is.na(b_site),0,b_site)
    b_clone <- as.numeric(x[4])
    max_age <- as.numeric(x[6]) + 25
    b <- b_gene+b_site+b_clone
    traj <- suppressWarnings(trajectory_from_parameters_2(
      b = b,u = u,g = 2,N = pop_size,x = seq(-20,max_age,by = 0.5) - values_model$min_age))
    age_at_onset <- t0_adjusted(u,b,g = 2,n = pop_size)
    return(data.frame(
      age = traj$x + values_model$min_age,
      VAF = traj$y,
      t0 = age_at_onset + values_model$min_age,
      ID = as.numeric(x[5])
    ))
  }) %>%
  do.call(what = rbind) %>%
  mutate(plausible = t0 >= 0) %>%
  mutate(VAF = ifelse(is.na(VAF),0,VAF))

x <- rbind(
  S %>% 
    subset(plausible == T) %>%
    group_by(age,plausible) %>%
    summarise(q05 = quantile(VAF,0.05,na.rm = T),
              q50 = quantile(VAF,0.5,na.rm = T),
              q95 = quantile(VAF,0.95,na.rm = T)) %>%
    as.data.frame,
  S %>% 
    group_by(age) %>%
    summarise(plausible = as.logical(0),
              q05 = quantile(VAF,0.05,na.rm = T),
              q50 = quantile(VAF,0.5,na.rm = T),
              q95 = quantile(VAF,0.95,na.rm = T)) %>%
    as.data.frame) %>%
  subset(q95 > 0)

plot_good_bad <- x %>% 
  group_by(age) %>% 
  subset(q95 > 0) %>% 
  mutate(q05 = ifelse(
    all(c(T,F) %in% plausible),
    ifelse(plausible == F,min(q95),q05),
    q05)) %>% 
  ggplot(aes(x = age,ymin = q05,ymax = q95,fill = plausible)) + 
  geom_ribbon(alpha = 0.3) + 
  geom_line(aes(y = q50,colour = plausible),size = 0.25) +
  geom_vline(xintercept = 0,linetype = 3,
             size = 0.5) +
  theme_gerstung(base_size = 6) + 
  xlab("Age") + 
  ylab("VAF") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.x = element_blank()) + 
  scale_fill_manual(values = colorspace::lighten(c("red4","goldenrod"),amount = 0.2),
                    guide = F) + 
  scale_colour_manual(values = colorspace::darken(c("red4","goldenrod"),amount = 0.2),
                      guide = F) + 
  scale_y_continuous(trans = 'log10') + 
  scale_x_continuous(breaks = 0)

plot_plausible <- x %>% 
  ggplot(aes(x = age,y = q50,group = plausible,colour = plausible)) + 
  geom_line(size = 0.25) + 
  geom_point(data = tmp_data,aes(x = Age,y = VAF),inherit.aes = F,size = 0.25) +
  geom_vline(xintercept = 0,linetype = 3,
             size = 0.5) +
  theme_gerstung(base_size = 6) + 
  xlab("Age") + 
  ylab("VAF") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.x = element_blank()) +
  scale_colour_manual(values = c("red4","goldenrod"),guide = F) + 
  scale_y_continuous(trans = 'log10') + 
  scale_x_continuous(breaks = 0,limits = c(min(x$age),max(x$age)))

plot_grid(plot_good_bad,plot_plausible) + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/good_bad_traj.pdf",model_id),
         width = 1.4,height = 0.55,
         useDingbats = F)
```

```{r,fig.width=0.9,fig.height=0.8}
plot_good_bad + 
    geom_point(data = tmp_data,aes(x = Age,y = VAF),inherit.aes = F,size = 0.25) + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/good_bad_traj_combo.pdf",model_id),
         width = 0.8,height = 0.8,
         useDingbats = F)
```

```{r,fig.height=0.5,fig.width=0.6}
idxs <- sampling_idxs %>% 
  subset((individual == 4088 & site == "DNMT3At-E578X") | (individual == 3877 & site == "U2AF1-Q157R"))

tmp_data <- full_data %>%
  subset(SardID %in% idxs$individual & amino_acid_change %in% idxs$site)

#idxs <- sampling_idxs[1,]

S <- data.frame(
  u = c(offset_samples[,idxs$clone_numeric]),
  b_gene = c(gene_samples[,idxs$gene_numeric]),
  b_site = c(site_samples[,idxs$site_numeric]),
  b_clone = c(clone_samples[,idxs$clone_numeric]),
  idx = c(replicate(2,seq(1,nrow(clone_samples)))),
  maxAge = c(sapply(idxs$maxAge, function(x) rep(x,nrow(offset_samples)))),
  site = c(sapply(idxs$site, function(x) rep(x,nrow(offset_samples)))),
  id = c(sapply(idxs$individual, function(x) rep(x,nrow(offset_samples))))) %>%
  apply(1,function(x) {
    u <- as.numeric(x[1])
    b_gene <- as.numeric(x[2])
    b_site <- as.numeric(x[3])
    b_site <- ifelse(is.na(b_site),0,b_site)
    b_clone <- as.numeric(x[4])
    max_age <- as.numeric(x[6]) + 10
    b <- b_gene+b_site+b_clone
    traj <- suppressWarnings(trajectory_from_parameters_2(
      b = b,u = u,g = 2,N = pop_size,x = seq(-20,max_age,by = 0.5) - values_model$min_age))
    age_at_onset <- t0_adjusted(u,b,g = 2,n = pop_size)
    return(suppressWarnings(data.frame(
      age = traj$x + values_model$min_age,
      VAF = traj$y,
      t0 = age_at_onset + values_model$min_age,
      ID = as.numeric(x[5]),
      site = x[7],
      individual = x[8]
    )))
  }) %>%
  do.call(what = rbind) %>%
  mutate(plausible = t0 >= 0) %>%
  mutate(VAF = ifelse(is.na(VAF),0,VAF))

x <- S %>% 
  group_by(age,site,individual) %>%
  summarise(plausible = as.logical(0),
            q05 = quantile(VAF,0.05,na.rm = T),
            q50 = quantile(VAF,0.5,na.rm = T),
            q95 = quantile(VAF,0.95,na.rm = T)) 

plot_slow_fast <- x %>% 
  ggplot(aes(x = age,y = q50,colour = site)) + 
  geom_ribbon(aes(ymin = q05,ymax = q95,fill = site),
              colour = NA,alpha = 0.5) +
  geom_line(size = 0.25) + 
  #geom_point(data = tmp_data,aes(x = Age,y = VAF),inherit.aes = F,size = 0.25) +
  geom_vline(xintercept = 0,linetype = 3,
             size = 0.5) +
  theme_gerstung(base_size = 6) + 
  xlab("Age") + 
  ylab("VAF") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  scale_colour_manual(values = c("goldenrod","lightblue"),guide = F) + 
  scale_fill_manual(values = c("goldenrod","lightblue"),guide = F) + 
  scale_y_continuous(trans = 'log10') + 
  scale_x_continuous(breaks = 0,limits = c(min(x$age),max(x$age))) + 
  theme(axis.title.x = element_blank())

plot_slow_fast + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/slow_fast_traj.pdf",model_id),
         width = 0.6,height = 0.5,
         useDingbats = F)
```

```{r,fig.height=0.65,fig.width=1.4}
idxs <- sampling_idxs %>% 
  subset(grepl("TP53",site))# %>% subset(individual == 4088 & site == "DNMT3At-E578X")
idxs <- idxs[1,]
tmp_data <- full_data %>%
  subset(SardID == idxs$individual & amino_acid_change == idxs$site)

S <- data.frame(
  u = offset_samples[,idxs$clone_numeric],
  b_gene = gene_samples[,idxs$gene_numeric],
  b_site = site_samples[,idxs$site_numeric],
  b_clone = clone_samples[,idxs$clone_numeric],
  idx = seq(1,nrow(clone_samples)),
  maxAge = idxs$maxAge) %>%
  apply(1,function(x) {
    u <- as.numeric(x[1])
    b_gene <- as.numeric(x[2])
    b_site <- as.numeric(x[3])
    b_site <- ifelse(is.na(b_site),0,b_site)
    b_clone <- as.numeric(x[4])
    max_age <- as.numeric(x[6]) + 25
    b <- b_gene+b_site+b_clone
    traj <- suppressWarnings(trajectory_from_parameters_2(
      b = b,u = u,g = 2,N = pop_size,x = seq(-20,max_age,by = 0.5) - values_model$min_age))
    age_at_onset <- t0_adjusted(u,b,g = 2,n = pop_size)
    return(data.frame(
      age = traj$x + values_model$min_age,
      VAF = traj$y,
      t0 = age_at_onset + values_model$min_age,
      ID = as.numeric(x[5])
    ))
  }) %>%
  do.call(what = rbind) %>%
  mutate(plausible = t0 >= 0) %>%
  mutate(VAF = ifelse(is.na(VAF),0,VAF))

x <- rbind(
  S %>% 
    subset(plausible == T) %>%
    group_by(age,plausible) %>%
    summarise(q05 = quantile(VAF,0.05,na.rm = T),
              q50 = quantile(VAF,0.5,na.rm = T),
              q95 = quantile(VAF,0.95,na.rm = T)) %>%
    as.data.frame,
  S %>% 
    group_by(age) %>%
    summarise(plausible = as.logical(0),
              q05 = quantile(VAF,0.05,na.rm = T),
              q50 = quantile(VAF,0.5,na.rm = T),
              q95 = quantile(VAF,0.95,na.rm = T)) %>%
    as.data.frame) %>%
  subset(q95 > 0)

plot_good_bad <- x %>% 
  group_by(age) %>% 
  subset(q95 > 0) %>% 
  mutate(q05 = ifelse(
    all(c(T,F) %in% plausible),
    ifelse(plausible == F,min(q95),q05),
    q05)) %>% 
  ggplot(aes(x = age,ymin = q05,ymax = q95,fill = plausible)) + 
  geom_ribbon() + 
  geom_vline(xintercept = 0,linetype = 3,
             size = 0.5) +
  theme_gerstung(base_size = 6) + 
  xlab("Age") + 
  ylab("VAF") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) + 
  scale_fill_manual(values = colorspace::lighten(c("red4","goldenrod"),amount = 0.01),
                    guide = F) + 
  scale_y_continuous(trans = 'log10') + 
  scale_x_continuous(breaks = 0)

plot_plausible <- x %>% 
  ggplot(aes(x = age,y = q50,group = plausible,colour = plausible)) + 
  geom_line(size = 0.25) + 
  geom_point(data = tmp_data,aes(x = Age,y = VAF),inherit.aes = F,size = 0.25) +
  geom_vline(xintercept = 0,linetype = 3,
             size = 0.5) +
  theme_gerstung(base_size = 6) + 
  xlab("Age") + 
  ylab("VAF") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  scale_colour_manual(values = c("red4","goldenrod"),guide = F) + 
  scale_y_continuous(trans = 'log10') + 
  scale_x_continuous(breaks = 0,limits = c(min(x$age),max(x$age)))

plot_grid(plot_good_bad,plot_plausible) + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/good_bad_traj_2.pdf",model_id),
         width = 1.4,height = 0.55,
         useDingbats = F)
```

## Comparing historical with observed growth

We can see clearly a bias in terms of genes where the difference between historical and observed growth is concerned - particularly, the phase of clonal dynamics we are able to observe for *TP53*, *DNMT3A* and *BRCC3* is far too slow to explain the observed VAF.

```{r historical_growth_plots,fig.height=2,fig.width=4}
aging_order <- c(
  "DNMT3A","BRCC3","TP53","GNB1","ASXL1","TET2","CBL",
  "CTCF","PPM1D","JAK2","KRAS","SRSF2-other","SF3B1","IDH1",
  "IDH2","PTPN11","U2AF1","SRSF2-P95H")
min_max <- c(
  exp(min(historical_estimates$observed_growth_mean,historical_estimates$historical_growth_mean,na.rm=T))-1,
  exp(max(historical_estimates$observed_growth_mean,historical_estimates$historical_growth_mean,na.rm=T))-1
)

growth_per_year_plot <- historical_estimates %>% 
  ggplot(aes(x = exp(observed_growth_mean)-1,
             y = exp(historical_growth_mean)-1,
             colour = gene)) + 
  geom_abline(slope = 1,
              size = 0.2,
              linetype=2) +
  geom_point(size = 0.2) +
  scale_y_continuous(limits = min_max,
                     breaks = seq(-0.2,1,by = 0.2),
                     expand = c(0.02,0),
                     labels = sprintf("%s%%",seq(-0.2,1,by = 0.2)*100)) + 
  scale_x_continuous(limits = min_max,
                     breaks = seq(-0.2,1,by = 0.2),
                     expand = c(0.02,0),
                     labels = sprintf("%s%%",seq(-0.2,1,by = 0.2)*100)) + 
  theme_gerstung(base_size = 6) + 
  scale_colour_manual(values = gene_colours,guide = F) +
  xlab("Observed growth per year") + 
  ylab("Historical growth per year")

growth_per_year_gene_plot_data <- historical_estimates %>%
  mutate(gene = ifelse(gene == "SRSF2",
                       ifelse(site == "SRSF2-P95H",
                              "SRSF2-P95H",
                              "SRSF2-other"),
                       gene)) %>%
  group_by(gene) %>%
  summarise(average_historical_growth = mean(exp(historical_growth_mean)-1,na.rm=T),
            average_observed_growth = mean(exp(observed_growth_mean)-1,na.rm=T))

growth_per_year_gene_plot <- growth_per_year_gene_plot_data %>% 
  ggplot(aes(x = factor(gene,levels = aging_order),
             y = (average_observed_growth + average_historical_growth)/2,
             colour = str_match(gene,"[A-Z0-9]+"))) +  
  geom_linerange(aes(ymin = average_observed_growth,ymax = average_historical_growth),
                 colour = "black") + 
  geom_text(aes(y = average_historical_growth,
                label = sprintf("%s%%",100*round(average_observed_growth/average_historical_growth,3))),
            colour = "black",
            size = 2.7,
            hjust = -0.2) + 
  geom_point(aes(y = average_observed_growth),size = 1) + 
  geom_point(aes(y = average_historical_growth),size = 1) + 
  geom_point(aes(y = average_historical_growth),size = 0.5,colour = 'white') + 
  theme_gerstung(base_size = 6) + 
  xlab("") + 
  ylab("Growth per year") +
  scale_y_continuous(breaks = seq(-0.2,1,by = 0.1),
                     labels = sprintf("%s%%",seq(-0.2,1,by = 0.1)*100),
                     expand = c(0,0,0.05,0),limits = c(0,NA)) + 
  scale_colour_manual(values = gene_colours,guide = F) + 
  rotate_x_text(face = "italic")

plot_grid(growth_per_year_plot,growth_per_year_gene_plot,nrow=1,align = "h") +
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/historical_growth_figure.pdf",model_id),
         height = 2,
         width = 3.5,
         useDingbats = F)

growth_per_year_gene_plot + 
  theme(axis.text.x = element_text(face = "plain",angle = 0,hjust = 0.5),
        axis.text.y = element_text(face = "italic")) + 
  coord_flip() + 
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/historical_growth_figure_larger.pdf",model_id),
         width = 1.7,height = 2,useDingbats = F)
```

### Comparison between the expected and observable growth for phylogenetic trees and longitudinal data

We can estimate the annual growth that one would expect and the growth one would be able to observe in the later years of a clone from phylogenetic trees alone. This allows us to draw a comparison between the EPS trajectory defined by a tree and the ratio between observed and plausible - or historical - growth. In this figure we present the combined evidence from simulations and real data (both longitudinal and phylogenetic) that clonal deceleration is observable across most clones, with a few genes being characterised by their larger deceleration - we estimate that *DNMT3A*, one of the most prevalent genes in CH, loses in average at least half of its annual growth rate.

```{r early_late_ratio,fig.height=1.4,fig.width=2.5}
mitchell_coefficients <- read.csv("data_output/mitchell_bnpr_coefficients.csv") %>%
  subset(w_sum < 10) %>%
  group_by(id,clade) %>% 
  select(id,clade,
         early = cp_1,late = cp_2,expected = nonlinear_vaf) %>%
  mutate(plotting_id = "Real data") %>% 
  mutate(source = "Mitchell et al. 2021") %>%
  mutate(colour = ifelse(grepl("DNMT3A",clade),"DNMT3A",
                         ifelse(grepl("Clade",clade),"No known driver","Other driver"))) 

tree_coefficients <- read.csv("data_output/tree_bnpr_coefficients.csv") %>%
  transmute(id = individual,clade = plot_label,
            early = b_early,late = b_early + b_late,expected = expected_late) %>%
  mutate(plotting_id = "Real data") %>% 
  mutate(source = "Fabre et al. 2021") %>%
  mutate(colour = str_match(clade,"[A-Z0-9]+")) %>%
  mutate(colour = ifelse(grepl("UD",colour),"No known driver",colour))
simulation_coefficients_full <- read.csv("data_output/simulated_bnpr_coefficients.csv") %>%
  subset(w_sum < 5)
  
simulation_coefficients <-  simulation_coefficients_full %>% 
  subset(col == "Low variance") %>% 
  mutate(plotting_id = ifelse(late_growth/growth_rate < 0.8,
                              "Saturating",
                              "Near\nconstant")) %>% 
  select(id = name,clade = driver,early = cp_1,late = cp_late,expected = expected_nl,plotting_id) %>%
  mutate(source = "Simulated") %>%
  mutate(colour = "No known driver")

all_bnpr_coefficients <- rbind(as.data.frame(mitchell_coefficients),
                               as.data.frame(tree_coefficients),
                               as.data.frame(simulation_coefficients))

x_order <- c("Near\nconstant","Saturating",
             "Real data","Obs./plaus.\ngrowth ratio") %>%
  rev

rbind(
  all_bnpr_coefficients %>% 
    group_by(source) %>% 
    summarise(N = sum((late / expected) < 0.44,na.rm = T),
              Total = length(id)) %>% 
    mutate(Ratio = N / Total),
  historical_estimates %>% 
    mutate(ratio = (exp(observed_growth_mean)-1) / (exp(historical_growth_mean)-1)) %>% 
    summarise(N = sum(ratio < 0.44,na.rm=T),Total = length(ratio)) %>% 
    mutate(Ratio = N / Total) %>%
    mutate(source = "Fabre et al. 2021 (longitudinal)") %>%
    select(source,N,Total,Ratio)) %>% 
  select(Source = source,`Number of saturating trajectories` = N,
         `Total number of trajectories` = Total,Proportion = Ratio) %>%
  as.data.frame

all_bnpr_coefficients %>%
  mutate(late_expected = late / expected) %>%
  mutate(late_expected = ifelse(late_expected < -1,-1,late_expected)) %>% 
  ggplot(aes(x = plotting_id,y = late_expected)) +
  geom_point(
    data = growth_per_year_gene_plot_data,inherit.aes = F,
    mapping = aes(x = "Obs./plaus.\ngrowth ratio",
                  y = average_observed_growth/average_historical_growth,
                  group = reorder(gene,average_observed_growth/average_historical_growth),
                  colour = str_match(gene,"[A-Z0-9]+"),
                  shape = "Fabre et al. 2021"),
    size = 0.3,
    position = position_dodge(width = 0.75)) +
  geom_point(aes(shape = source,colour = colour),size = 0.25,
             position = position_jitter(width = 0.2,height = 0),
             alpha = 0.9) + 
  geom_boxplot(size = 0.25,fill = NA,outlier.color = NA) +
  geom_hline(yintercept = 1.0,size = 0.25) + 
  theme_gerstung(base_size = 6) + 
  scale_shape_manual(values = c("circle","triangle","cross"),name = NULL) +
  scale_x_discrete(breaks = x_order,limits = x_order) +
  xlab("") + 
  ylab("Growth ratio") + 
  scale_colour_manual(values = c(gene_colours,
                                 `Other driver` = "#18A558",
                                 `No known driver` = "grey50"),guide = F) + 
  coord_flip(ylim = c(NA,NA)) + 
  scale_y_continuous(breaks = c(-1,0,0.5,seq(1,10)),
                     labels = c("<-1","0","0.5",seq(1,10))) +
  theme(legend.position = "top",
        legend.key.size = unit(0,"cm"),legend.box.spacing = unit(0.1,"cm")) +
  guides(shape = guide_legend(nrow = 2)) +
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/bnpr_historical_growth.pdf",model_id),
         useDingbats = F,
         width = 1.9,height = 1.2)
```

```{r,fig.height=0.8,fig.width=1.9}
all_bnpr_coefficients %>%
  mutate(late_expected = late / expected) %>%
  mutate(late_expected = ifelse(late_expected < -1,-1,late_expected)) %>% 
  subset(source != "Simulated") %>%
  ggplot(aes(x = plotting_id,y = late_expected)) +
  geom_point(
    data = growth_per_year_gene_plot_data,inherit.aes = F,
    mapping = aes(x = "Obs./plaus.\ngrowth ratio",
                  y = average_observed_growth/average_historical_growth,
                  group = reorder(gene,average_observed_growth/average_historical_growth),
                  colour = str_match(gene,"[A-Z0-9]+"),
                  shape = "Fabre et al. 2021"),
    size = 0.3,
    position = position_dodge(width = 0.75)) +
  geom_point(size = 0.25,
             position = position_jitter(width = 0.2,height = 0),
             alpha = 0.9,
             colour = "grey") + 
  geom_boxplot(size = 0.25,fill = NA,outlier.color = NA) +
  geom_hline(yintercept = 1.0,size = 0.25) + 
  theme_gerstung(base_size = 6) + 
  scale_shape_manual(values = c("circle","triangle","cross"),guide = F) +
  scale_x_discrete(breaks = x_order[1:2],limits = x_order[1:2]) +
  xlab("") + 
  ylab("Growth ratio") + 
  scale_colour_manual(values = c(gene_colours,
                                 `Other driver` = "#18A558",
                                 `No known driver` = "grey50"),guide = F) + 
  scale_y_continuous(breaks = c(-1,0,0.5,1,2),
                     labels = c("<-1","0","0.5","1","2")) +
  theme(legend.position = "none",
        legend.key.size = unit(0,"cm"),legend.box.spacing = unit(0.1,"cm")) +
  guides(shape = guide_legend(nrow = 2)) +
  coord_flip() +
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/bnpr_historical_growth_simplified.pdf",model_id),
         useDingbats = F,
         width = 1.9,height = 0.8)
```

### Comparison between expected and observed VAF 

```{r,fig.height=1,fig.width=1}
mitchell_vaf <- read.csv("data_output/mitchell_bnpr_coefficients.csv") %>%
  subset(w_sum < 10) %>%
  group_by(id,clade) %>% 
  select(id,clade,
         observed_vaf=vaf_from_clade,expected_vaf) %>%
  mutate(plotting_id = "Real data") %>% 
  mutate(source = "Mitchell et al. 2021") %>%
  mutate(colour = ifelse(grepl("DNMT3A",clade),"DNMT3A",
                         ifelse(grepl("Clade",clade),"No known driver","Other driver"))) 

tree_vaf <- read.csv("data_output/tree_bnpr_coefficients.csv") %>%
  subset(source == "Colonies") %>% 
  select(id = individual,clade = plot_label,
         observed_vaf,expected_vaf) %>%
  mutate(plotting_id = "Real data") %>% 
  mutate(source = "Fabre et al. 2021") %>%
  mutate(colour = str_match(clade,"[A-Z0-9]+")) %>%
  mutate(colour = ifelse(grepl("UD",colour),"No known driver",colour)) %>%
  mutate(id = factor(id))

rbind(as.data.frame(mitchell_vaf),as.data.frame(tree_vaf)) %>%
  mutate(expected_vaf = expected_vaf) %>% 
  mutate(expected_vaf = ifelse(expected_vaf > 1,1,expected_vaf)) %>% 
  ggplot(aes(x = observed_vaf,y = expected_vaf)) + 
  geom_point(size = 0.5) + 
  geom_abline(slope = 1,size = 0.25) + 
  theme_gerstung(base_size = 6) + 
  scale_y_continuous(breaks = c(0,0.5,1.0)) +
  scale_x_continuous(breaks = c(0,0.5,1.0)) + 
  xlab("Observed clone freq.") + 
  ylab("Extrapolated\nclone frequency") 
```

```{r,fig.height=0.7,fig.width=1}
rbind(as.data.frame(mitchell_vaf),as.data.frame(tree_vaf)) %>%
  mutate(expected_vaf = expected_vaf) %>% 
  mutate(expected_vaf = ifelse(expected_vaf > 1,1,expected_vaf)) %>% 
  ggplot(aes(y = expected_vaf/observed_vaf,x = "")) + 
  geom_hline(yintercept = 1,size = 0.25) +
  geom_point(position = position_jitter(height = 0.2),size = 0.25,colour = "grey") + 
  geom_boxplot(outlier.color = NA,alpha = 0.5,size = 0.25) +
  theme_gerstung(base_size = 6) + 
  theme(axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank()) +
  ylab("Expected/observed\nclone size fold change") +
  scale_y_continuous(trans = 'log10',
                     breaks = c(0.1,1e-2,1,10,100),labels = c("0.1","0.01","1","10","100")) + 
  coord_flip() +
  ggsave(filename = sprintf("figures/%s/dynamic_coefficients/expected_vs_observed_vaf.pdf",model_id),
         useDingbats = F,
         width = 1.1,height = 0.75)
```

## Investigating the determinants of the difference between historical and observed growth effect

```{r historical_ratios}
unique_genes <- historical_estimates %>%
  group_by(gene) %>%
  filter(length(unique(paste(individual,site))) > 10) %>%
  select(gene) %>%
  distinct %>%
  unlist
historical_estimates <- historical_estimates %>%
  mutate(
    ratio_q05 = observed_growth_q05/historical_growth_q05,
    ratio_q50 = observed_growth_q50/historical_growth_q50,
    ratio_q95 = observed_growth_q95/historical_growth_q95,
    ratio_mean = observed_growth_mean/historical_growth_mean
  )
```

### Distribution by gene

```{r historical_distribution_gene,fig.height=2,fig.width=2.5}
historical_estimates %>%
  group_by(gene) %>%
  mutate(M = unlist(quantile(ratio_q50,c(0.5),na.rm = T))) %>%
  mutate(ratio_q50 = ifelse(ratio_q50 < -0.5,-0.5,ratio_q50)) %>% 
  ggplot(aes(x = reorder(gene,M),y = ratio_q50)) + 
  geom_hline(yintercept = 1,linetype = 2,size = 0.25) + 
  geom_point(size = 0.5,position = position_jitter(width = 0.1)) +
  geom_boxplot(outlier.size = 0,alpha = 0.7,size = 0.25) + 
  coord_flip() + 
  theme_gerstung(base_size = 6) +
  scale_y_continuous(breaks = c(-0.5,0,1),labels = c("<-0.5","0","1")) + 
  theme(axis.text.y = element_text(face = "italic")) + 
  ylab("Ratio of historical and observed annual fold change") + 
  xlab("Gene")
```

### Association with recurrence and truncation status 

```{r historical_recurrence_truncation,fig.width=2.5,fig.height=1}
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

tmp <- historical_estimates %>%
  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","non-truncating"),
    sep = ' and\n'
  ))

anvar_test <- anova(glm(ratio_q50 ~ f,data = tmp),
                    test = "F")
tmp %>%
  mutate(ratio_q50 = ifelse(ratio_q50 < -1,-1,ratio_q50)) %>%
  ggplot(aes(x = f,y = ratio_q50)) + 
  geom_boxplot(outlier.size = 0.5,size = 0.25) + 
  ylab("Historical and observed annual fold change") + 
  xlab("") +
  theme_gerstung(base_size = 6) + 
  coord_flip(ylim = c(-1,NA)) + 
  scale_y_continuous(breaks = c(-1,0,0.5,1),labels = c("<-1",0,0.5,1)) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/other_associations/uc_tnt_rec.pdf",model_id),
         height = 1.2,width = 2.0) 
anvar_test
```

```{r historical_recurrence_truncation_dnmt3a_hs,fig.width=2.5,fig.height=1}
tmp_dnmt3a <- tmp %>%
  subset(gene == "DNMT3A") %>%
  mutate(f = ifelse(grepl("Recurring",f),"Hot-spot","Normal mutation"))

wilcox.test(ratio_q50 ~ f, data = tmp_dnmt3a)

tmp_dnmt3a %>%
  ggplot(aes(x = f,y = ratio_q50)) + 
  geom_boxplot(outlier.size = 0.5,size = 0.25) + 
  ylab("Observed vs. historical annual fold change") + 
  xlab("") +
  theme_gerstung(base_size = 6) + 
  coord_flip(ylim = c(-1,NA)) + 
  scale_y_continuous() +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/other_associations/uc_tnt_rec_dnmt3a.pdf",model_id),
         height = 1.2,width = 2.5) 
```

### Association with gender and smoking

We tested whether age, gender and smoking could explain a greater amount of variance than age alone and found no evidence for this.

```{r historical_gender_smoking,fig.width=2,fig.height=1.3}
smoking_gender_data <- load_smoking_data() %>% 
  group_by(SardID) %>%
  summarise(HasSmoked = any(PackYears>0)) %>%
  merge(distinct(select(full_data,SardID,Gender))) %>%
  select(individual = SardID,HasSmoked,Gender)

historical_estimates_gender_smoking <- merge(
  historical_estimates,smoking_gender_data,
  by = c("individual")) 

tmp <- historical_estimates_gender_smoking %>%
  mutate(min_age = minAge) %>% 
  filter(!is.na(ratio_q50) & !is.na(HasSmoked) & !is.na(Gender) & !is.na(min_age))
all_models <- list(
  `All genes` = list(
    m = glm(ratio_q50 ~ HasSmoked + Gender + min_age,data = tmp),
    null_m = glm(ratio_q50 ~ min_age,data = tmp)))

for (g in unique_genes) {
  x <- tmp %>% 
    subset(gene == g)
  all_models[[g]] <- list(
    m = glm(ratio_q50 ~ HasSmoked + Gender + min_age,data = x),
    null_m = glm(ratio_q50 ~ min_age,data = x))
}

significance_df <- all_models %>% 
  lapply(
    function(x) anova(x$null_m,x$m,test="F")[2,c(5,6)]
  ) %>%
  do.call(what = rbind)
significance_df$gene <- names(all_models)
significance_df$SignificantAfterCorrection <- p.adjust(significance_df$`Pr(>F)`,method = "BH") < 0.05

significance_df

tmp %>% 
  mutate(
    L = paste(ifelse(Gender == "M","Male","Female"),
              ifelse(HasSmoked == T,"\nhas smoked","\nnever smoked"),
              sep = '')
  ) %>%
  mutate(ratio_q50 = ifelse(ratio_q50 < -1,-1,ratio_q50)) %>%
  ggplot(aes(x = L,y = ratio_q50, colour = Gender)) + 
  geom_boxplot(size = 0.2,outlier.size = 0.5) + 
  theme_gerstung(base_size = 6) + 
  coord_flip(ylim = c(-1,NA)) + 
  xlab("") + 
  ylab("Ratio of historical and observed annual fold change") +
  scale_y_continuous() +
  scale_colour_aaas(guide = F)
```

### dNdS in AML and MDS

```{r dnds_aml_mds,fig.height=2,fig.width=3.5}
all_samples <- sampling_idxs %>%
  group_by(site) %>%
  mutate(re = length(unique(individual)) > 1) %>%
  subset(re == T) %>% 
  apply(1,function(x) {
    individual <- x[1]
    site_numeric <- as.numeric(x[3])
    gene_numeric <- as.numeric(x[4])
    clone_numeric <- as.numeric(x[5])
    minAge <- as.numeric(x[6])
    
    if (is.na(site_numeric)) {
      b_site <- 0} else {
        b_site <- site_samples[,site_numeric]
      }
    b_gene <- gene_samples[,gene_numeric]
    b_clone <- clone_samples[,clone_numeric]
    b <- b_gene + b_clone + b_site
    u_values <- offset_samples[,clone_numeric]
    ages_at_onset <- t0_adjusted(u_values,b,g=2,n=pop_size) + values_model$min_age
    ages_at_onset_corrected <- ifelse(ages_at_onset < -1,-1,ages_at_onset)
    ages_at_onset_corrected <- ifelse(ages_at_onset > minAge,minAge,ages_at_onset)
    biologically_feasible <- ages_at_onset > -1 & ages_at_onset < minAge
    historical_b <- b[biologically_feasible] %>%
      unlist
    return(tibble(
      historical_samples = historical_b,
      aachange = gsub("-","",str_match(x[2],"-.+")),
      gene = str_match(x[2],"[A-Z0-9]+")
    ))
  }) %>%
  do.call(what = rbind) %>%
  group_by(gene,aachange) %>%
  summarise(mean = mean(historical_samples,na.rm = T),
            s = sd(historical_samples,na.rm = T)) %>%
  mutate(site = sprintf("%s-%s",gene,aachange))

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

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

dd %>% 
  #subset(gene == "DNMT3A") %>% 
  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,height = 0) + 
  geom_linerange(aes(ymin = CI_low,ymax = CI_high),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("dN/dS") +
  scale_x_continuous(labels = function(x) sprintf("%s%%",x*100)) + 
  scale_y_continuous(trans = 'log10') +
  theme_gerstung(base_size = 6) +
  #coord_cartesian(xlim = c(NA,0.5)) +
  ggsave(useDingbats=FALSE,
         filename = sprintf("figures/%s/dynamic_coefficients/dnds_association_hist.pdf",model_id),
         width = 3.5,height = 2)

lm(formula = dnds ~ mean,data = dd,subset = dd$condition == "AML") %>%
  summary

lm(formula = dnds ~ mean,data = dd,subset = dd$condition == "MDS") %>%
  summary

lm(formula = dnds ~ mean,data = dd,subset = dd$condition == "AML" & dd$gene == "DNMT3A") %>%
  summary

lm(formula = dnds ~ mean,data = dd,subset = dd$condition == "MDS" & dd$gene == "DNMT3A") %>%
  summary

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

### AML risk

```{r aml_risk,fig.height=2,fig.width=3.5}
all_samples <- sampling_idxs %>%
  apply(1,function(x) {
    individual <- x[1]
    site_numeric <- as.numeric(x[3])
    gene_numeric <- as.numeric(x[4])
    clone_numeric <- as.numeric(x[5])
    minAge <- as.numeric(x[6])
    
    if (is.na(site_numeric)) {
      b_site <- 0} else {
        b_site <- site_samples[,site_numeric]
      }
    b_gene <- gene_samples[,gene_numeric]
    b_clone <- clone_samples[,clone_numeric]
    b <- b_gene + b_clone + b_site
    u_values <- offset_samples[,clone_numeric]
    ages_at_onset <- t0_adjusted(u_values,b,g=2,n=pop_size) + values_model$min_age
    ages_at_onset_corrected <- ifelse(ages_at_onset < -1,-1,ages_at_onset)
    ages_at_onset_corrected <- ifelse(ages_at_onset > minAge,minAge,ages_at_onset)
    biologically_feasible <- ages_at_onset > -1 & ages_at_onset < minAge
    historical_b <- b[biologically_feasible] %>%
      unlist
    return(tibble(
      historical_samples = historical_b,
      gene = str_match(x[2],"[A-Z0-9]+")
    ))
  }) %>%
  do.call(what = rbind) %>%
  group_by(gene) %>%
  summarise(mean = mean(historical_samples,na.rm = T),
            s = sd(historical_samples,na.rm = T))

dd <- all_samples %>% 
  merge(read.csv("data/aml-risk-coef.csv"),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,height = 0) + 
  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_hist.pdf",model_id),
         width = 3.5,height = 2)

cor.test(dd$mean,dd$coef)

cor.test(dd$mean,dd$coef,method = "kendall")
```

# Poor fits

## Investigating the determinants of poor fits

### Association with recurrence and truncation status

```{r poor_fits_recurrence_truncation,fig.width=2.5,fig.height=1}
tmp <- statistics_data %>%
  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(sum_stat ~ f,data = tmp,family = stats::binomial()),
                    test = "Chisq")
tmp %>%
  group_by(f) %>% 
  summarise(Proportion = sum(sum_stat)/length(sum_stat)) %>%
  ggplot(aes(x = f,y = Proportion)) + 
  geom_bar(outlier.size = 0.5,size = 0.25,stat = "identity") + 
  ylab("Proportion of fits with no outliers") + 
  xlab("") +
  theme_gerstung(base_size = 6) + 
  coord_flip() + 
  scale_y_continuous(expand = c(0,0)) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/other_associations/poor_tnt_rec.pdf",model_id),
         height = 1.2,width = 2.0) 
anvar_test
```

### Associations with smoking status and gender

```{r poor_fits_gender_smoking,fig.width=2,fig.height=1.3}
good_stat_smoking_gender <- merge(
  statistics_data,smoking_gender_data,
  by = c("individual")) 

tmp <- good_stat_smoking_gender %>%
  filter(!is.na(sum_stat) & !is.na(HasSmoked) & !is.na(Gender))
all_models <- list(
  `All genes` = list(
    m = glm(sum_stat ~ HasSmoked + Gender,data = tmp,
            family = stats::binomial()),
    null_m = glm(sum_stat ~ 1,data = tmp,
                 family = stats::binomial())))

for (g in unique_genes) {
  x <- tmp %>% 
    subset(gene == g)
  all_models[[g]] <- list(
    m = glm(sum_stat ~ HasSmoked + Gender,data = x,
            family = stats::binomial()),
    null_m = glm(sum_stat ~ 1,data = x,
                 family = stats::binomial()))
}

significance_df <- all_models %>% 
  lapply(
    function(x) anova(x$null_m,x$m,test="Chisq")[2,c(4,5)]
  ) %>%
  do.call(what = rbind)
significance_df$gene <- names(all_models)
significance_df$SignificantAfterCorrection <- p.adjust(significance_df$`Pr(>Chi)`,method = "BH") < 0.05

significance_df

tmp %>% 
  mutate(
    L = paste(ifelse(Gender == "M","Male","Female"),
              ifelse(HasSmoked == T,"\nhas smoked","\nnever smoked"),
              sep = '')
  ) %>%
  group_by(L,Gender) %>%
  summarise(proportion = sum(sum_stat) / length(sum_stat)) %>%
  ggplot(aes(x = L,y = proportion, fill = Gender)) + 
  geom_bar(stat = "identity") + 
  theme_gerstung(base_size = 6) + 
  coord_flip() +
  xlab("") + 
  ylab("Proportion of fits with no outliers") +
  scale_y_continuous(expand = c(0,0)) +
  scale_fill_aaas(guide = F) +
  ggsave(useDingbats=FALSE,sprintf("figures/%s/other_associations/poor_smoke_sex.pdf",model_id),
         height = 1.2,width = 2.0) 
```

# 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()
```