The longitudinal dynamics and natural history of clonal haematopoiesis

Simulations and modelling

1 Validating inference through simulations

Here we present the modelling work that validates this approach. We simulate several Fischer-Wright dynamics with passenger and driver mutations across a grid of different driver mutation rates (\(0.1*10^{-9}\) to \(40.0*10^{-9}\)) and fitness values (\(0.001\) to \(0.040\)) per generation for a fixed population of 200000 HSC. We assume that HSCs divide 13 times per year, yielding for 100 years, 1300 generations. An example of some of these runs is presented below. We also assume that 10 driver and 200 passenger sites exist for each of these mutation rates and fitness values.

N_DRIVERS <- 50
all_files <- list.files(path = "/hps/nobackup/research/gerstung/josegcpa/projects/05VAF_DYNAMICS/hsc_output",
                        full.names = T) %>%
  lapply(function(x) list.files(x,pattern = 'csv',full.names = T)) %>%
  unlist
pop_size <- 2e5

We wish to determine if a model such as the one posited by us and in a similar scenario can accurately infer the fitness effects. As such, we need to simulate the same data collection circumstances, namely:

• Sample one individual between 3 and 5 times with the earliest timepoint being sampled from a Gamma distribution fitted to the first timepoints from the collected data and the following times separated by approximately 3 years (40 generations)
• Coverage is also simulated according to what was observed in our data, simmulating this using a Gamma distribution fitted to the coverage for all timepoints
• Sampling at each timepoint for each individual is done by sampling from a beta binomial distribution with the sampled coverage as the number of trials, the probabilities as the allele frequencies of the driver mutations, and a fixed technical dispersion (with \(\alpha = \frac{p\beta}{1-p}\) and \(\beta=dispersion\)). We only account for drivers because we are working with targetted sequencing, which focuses on mutations known to be associated with CH

1.1 Inspecting and fitting distributions for age and coverage

plot(main="Coverage density",density(full_data$TOTALcount))

full_data %>%
  select(SardID,Age) %>%
  distinct %>% 
  group_by(SardID) %>%
  filter(Age == min(Age)) %>%
  ungroup %>%
  select(Age) %>%
  unlist %>% 
  density %>%
  plot(main = "Age at first timepoint density")

min_age_data <- full_data %>%
  select(SardID,Age) %>%
  distinct %>% 
  group_by(SardID) %>%
  filter(Age == min(Age)) %>%
  ungroup %>%
  select(Age) %>%
  unlist

coverage_distr <- fitdistrplus::fitdist(full_data$TOTALcount,
                                        method = "mme",
                                        #probs = c(0.1,0.9),
                                        distr = "gamma")
min_age_distr <- fitdistrplus::fitdist(min_age_data,
                                       method = "mge",
                                       distr = "gamma")

## Warning in fitdistrplus::fitdist(min_age_data, method = "mge", distr = "gamma"):
## maximum GOF estimation has a default 'gof' argument set to 'CvM'

plot(min_age_distr)

plot(coverage_distr)

sample_age <- function(n) {
  s <- rgamma(n,min_age_distr$estimate[1],min_age_distr$estimate[2])
  s <- ifelse(s < min(min_age_data),min(min_age_data),s)
  s <- ifelse(s > max(min_age_data),max(min_age_data),s)
  return(s)
}

sample_coverage <- function(n) {
  s <- rgamma(n,coverage_distr$estimate[1],coverage_distr$estimate[2])
  s <- ifelse(s < min(full_data$TOTALcount),
              min(full_data$TOTALcount),s)
  s <- ifelse(s > max(full_data$TOTALcount),
              max(full_data$TOTALcount),s)
  return(round(s))
}

While these fits are not perfect they are good enough to help us simulate the conditions under the experimental data was obtained.

1.2 Loading the sampling

all_data <- list(
  N200 = mget(load("models/simulation_model_200k_13.RData")),
  N100 = mget(load("models/simulation_model_100k_5.RData")),
  N50 = mget(load("models/simulation_model_50k_1.RData"))
  )

predicted_variables <- variable_summaries(do.call(rbind,all_data$N200$draws)) %>% 
  as_tibble() %>%
  subset(grepl('genetic',variable)) %>% 
  mutate(variable_no = as.numeric(str_match(variable,'[0-9]+(?=\\])'))) %>% 
  mutate(true_fitness = unique(all_data$N200$simulated_samples$fitness)[variable_no]) %>% 
  mutate(true_fitness_per_year = true_fitness * 13 / age_mult) %>%
  spread(key = labels,value = values)
plot_min_max <- c(min(predicted_variables$HDPI_low,predicted_variables$true_fitness_per_year),
                  max(predicted_variables$HDPI_high,predicted_variables$true_fitness_per_year))
predicted_coefficients_plot <- predicted_variables %>%
  ggplot(aes(y = `0.50`,ymin = `HDPI_low`,ymax = `HDPI_high`,x = true_fitness_per_year)) + 
  geom_abline(slope=1,intercept = 0,linetype=2) +
  #geom_smooth(method = 'lm',formula = y ~ x,color = 'grey4',linetype=2,size=1,se=F) +
  geom_point(size = 0.5) + 
  geom_linerange(size = 0.2) +
  theme_gerstung(base_size = 6) + 
  theme(axis.title = element_text(size = 6)) + 
  xlab("Simulated genetic\nselection coefficient") + 
  ylab("Inferred genetic\nselection coefficient")

predicted_coefficients_plot

cor_test_inference <- cor.test(predicted_variables$true_fitness_per_year,predicted_variables$`0.50`)
cor_test_inference

## 
##  Pearson's product-moment correlation
## 
## data:  predicted_variables$true_fitness_per_year and predicted_variables$`0.50`
## t = 15.743, df = 23, p-value = 8.273e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9026530 0.9809485
## sample estimates:
##       cor 
## 0.9565979

cat(paste("Rsquared =",cor_test_inference$estimate^2))

## Rsquared = 0.915079532301844

all_data$N200$Data %>%
  ggplot(aes(x = Gen,y = VAF,group = paste(Individual,Mutation))) +
  geom_line(alpha = 0.5) +
  theme_gerstung() +
  theme(legend.position = "bottom") + 
  facet_wrap(~fitness,nrow = 3,scales = "free") + 
  ggtitle("Samples from our simulated cohort for each fitness value (N=200,000; g=13)") + 
  scale_y_continuous(limits = c(0,0.5)) + 
  scale_x_continuous(c(min(Data$Gen),max(Data$Gen)))

all_data$N100$Data %>%
  ggplot(aes(x = Gen,y = VAF,group = paste(Individual,Mutation))) +
  geom_line(alpha = 0.5) +
  theme_gerstung() + 
  theme(legend.position = "bottom") + 
  facet_wrap(~fitness,nrow = 3,scales = "free") + 
  ggtitle("Samples from our simulated cohort for each fitness value  (N=100,000; g=13)") + 
  scale_y_continuous(limits = c(0,0.5)) + 
  scale_x_continuous(c(min(Data$Gen),max(Data$Gen)))

all_data$N50$Data %>%
  ggplot(aes(x = Gen,y = VAF,group = paste(Individual,Mutation))) +
  geom_line(alpha = 0.5) +
  theme_gerstung() + 
  theme(legend.position = "bottom") + 
  facet_wrap(~fitness,nrow = 3,scales = "free") + 
  ggtitle("Samples from our simulated cohort for each fitness value  (N=50,000; g=1)") + 
  scale_y_continuous(limits = c(0,0.5)) + 
  scale_x_continuous(c(min(Data$Gen),max(Data$Gen)))

We now need to set a model similar to the one described below, which is similar to the one used on the Sardinian cohort. The main difference here is that the model here presented does not have three distinct levels of a genetic effect (gene, domain and site) and has a single effect whose prior is the addition of the three previous priors (\(N(0,0.1) + N(0,0.1) + N(0,0.1) = N(0,0.14)\)). We also choose to not have mutations with \(VAF > 0.45\) for all timepoints because this is an unrealistic scenario which we do not observe in our data.

1.3 Assessing the quality of the fit

1.3.1 MCMC convergence diagnostic

1.3.1.1 N=200,000; g=13

mcmc_trace(all_data$N200$draws,regex_pars = "genetic")

rg_stat <- coda::gelman.diag(all_data$N200$draws)
ess <- coda::effectiveSize(all_data$N200$draws)
plot(rg_stat$psrf[,1],main = "Distribution of Rubin-Gelman Statistic",ylab = "Rubin Gelman Statistic",xlab = "Variable index")

plot(ess,main = "Effective sample sizes for all variables",ylab = "Effective sample size",xlab = "Variable index",ylim = c(0,max(ess)),axes=F)
axis(side = 2, at=seq(0,max(ess),by = 500),las=2)
axis(side = 1)

1.3.1.2 N=100,000; g=5

mcmc_trace(all_data$N100$draws,regex_pars = "genetic")

rg_stat <- coda::gelman.diag(all_data$N100$draws)
ess <- coda::effectiveSize(all_data$N100$draws)
plot(rg_stat$psrf[,1],main = "Distribution of Rubin-Gelman Statistic",ylab = "Rubin Gelman Statistic",xlab = "Variable index")

plot(ess,main = "Effective sample sizes for all variables",ylab = "Effective sample size",xlab = "Variable index",ylim = c(0,max(ess)),axes=F)
axis(side = 2, at=seq(0,max(ess),by = 500),las=2)
axis(side = 1)

1.3.1.3 N=50,000; g=1

mcmc_trace(all_data$N50$draws,regex_pars = "genetic")

rg_stat <- coda::gelman.diag(all_data$N50$draws)
ess <- coda::effectiveSize(all_data$N50$draws)
plot(rg_stat$psrf[,1],main = "Distribution of Rubin-Gelman Statistic",ylab = "Rubin Gelman Statistic",xlab = "Variable index")

plot(ess,main = "Effective sample sizes for all variables",ylab = "Effective sample size",xlab = "Variable index",ylim = c(0,max(ess)),axes=F)
axis(side = 2, at=seq(0,max(ess),by = 500),las=2)
axis(side = 1)

We can see clearly from analysing the draws for the genetic coefficients that the MCMC has converged. Additionally, the above plots shows Rubin-Gelman statistic values are very close to 1 for all variables and effective sample sizes are adequate. These 3 separate analysis combined together indicate that convergence has been achieved and that our estimates are acceptable.

1.3.2 Goodness of fit (residual effect)

To assess the goodness of fit we assume that, for each trajectory, the probability density of our true counts under the inferred model - \(p\) - should follow a uniform distribution if the fit is perfect. As such, the normal normal quantiles of \(p\) should follow a normal distribution and the sum of their squares should follow a Chi-squared distribution. This allows us to have a "goodness of fit" for each trajectory varying between 0 (not good) and 1 (perfect), which can be stated as \(s=CDF_{\chi^2}(\sum_{i=1}^{n}Q_{N(0,1)}(CDF_{BBinom(coverage,\alpha_{inferred},\beta_{inferred})}(counts))^2)\) (here \(p=CDF_{BBinom(coverage,\alpha_{inferred},\beta_{inferred})}(counts)\)). Using this, we now define a residual effect (RE) as \(RE = 1 - s\) and 0.7 as a threshold for declaring a trajectory as good (\(RE\lt0.7\)) or bad (\(RE\geq0.7\)). Unless noted otherwise, we here use only the simulations with N = 200,000 and g = 13.

get_r_values <- function(draws,Data,Age,gen_time) {
  all_predicted_variables <- variable_summaries(do.call(rbind,draws)) %>% 
    as_tibble() %>%
    mutate(variable_no = as.numeric(str_match(variable,'[0-9]+(?=\\])'))) %>% 
    mutate(true_fitness = unique(simulated_samples$fitness)[variable_no]) %>%
    mutate(variable_no = as.numeric(as.character(variable_no))) %>%
    mutate(true_fitness_per_year = gen_time * true_fitness) %>%
    spread(key = labels,value = values)
  genetic_effect_inferred <- all_predicted_variables %>%
    subset(grepl('genetic',variable)) %>%
    arrange(variable_no)
  clone_effect_inferred <- all_predicted_variables %>%
    subset(grepl('clone',variable)) %>%
    arrange(variable_no)
  individual_effect_inferred <- all_predicted_variables %>%
    subset(grepl('individual',variable)) %>%
    arrange(variable_no)
  
  beta_val <- all_predicted_variables %>%
    subset(variable == 'beta_coefficient')
  
  min_age <- min(Data$Gen / gen_time * age_mult)
  
  r_values <- data.frame(
    fitness = Data$fitness,
    fitness_factor = Data$fitness_factor,
    individual_factor = Data$individual_factor,
    genetic_effect = genetic_effect_inferred$mean[Data$fitness_factor],
    genetic_effect_005 = genetic_effect_inferred$`0.05`[Data$fitness_factor],
    genetic_effect_095 = genetic_effect_inferred$`0.95`[Data$fitness_factor],
    clone_effect = clone_effect_inferred$mean[Data$individual_factor],
    clone_effect_005 = clone_effect_inferred$`0.05`[Data$individual_factor],
    clone_effect_095 = clone_effect_inferred$`0.95`[Data$individual_factor],
    individual_offset = individual_effect_inferred$mean[Data$individual_factor],
    individual_offset_005 = individual_effect_inferred$`0.05`[Data$individual_factor],
    individual_offset_095 = individual_effect_inferred$`0.95`[Data$individual_factor],
    true_clone_age = Data$GenAtCloneFoundation/gen_time,
    true_count = Data$sample,
    coverage = Data$coverage,
    age = Age,
    Mutation = Data$Mutation,
    Individual = Data$Individual,
    min_age = min_age) %>%
    as_tibble() %>% 
    mutate(mu_val = inv.logit((genetic_effect + clone_effect)*age + individual_offset)/2) %>%
    mutate(
      tail_prob = extraDistr::pbbinom(
      q = true_count,
      size = coverage,
      alpha = (mu_val * beta_val$mean)/(1 - mu_val),
      beta = beta_val$mean)) %>% 
    group_by(individual_offset) %>%
    mutate(sum_stat = all(tail_prob > 0.025 & tail_prob < 0.975))
  return(r_values)
}

r_values_N200 <- get_r_values(
  draws = all_data$N200$draws,
  Data = all_data$N200$Data,
  Age = all_data$N200$Age,
  gen_time = 13)

r_values_N100 <- get_r_values(
  draws = all_data$N100$draws,
  Data = all_data$N100$Data,
  Age = all_data$N100$Age,
  gen_time = 5)

r_values_N50 <- get_r_values(
  draws = all_data$N50$draws,
  Data = all_data$N50$Data,
  Age = all_data$N50$Age,
  gen_time = 1)

gof_plot_values <- r_values_N200 %>%
  na.omit() %>% 
  select(Individual,fitness,sum_stat) %>%
  distinct() %>%
  mutate(bad = sum_stat == F) %>%
  mutate(bad = factor(bad,levels = c(T,F),labels = c("One or more\noutliers","No outliers")))

## Adding missing grouping variables: `individual_offset`

explained_variants_plot <- gof_plot_values %>%
  mutate(fitness = cut(fitness,c(0.000,0.005,0.010,0.015,0.020,0.026),
                       c("0.001-0.005","0.005-0.010","0.010-0.015","0.015-0.020","0.020-0.025"))) %>% 
  group_by(fitness) %>%
  mutate(Total = length(unique(individual_offset))) %>% 
  group_by(fitness,bad) %>%
  summarise(N = length(unique(individual_offset)),
            Total = Total[1]) %>%
  ggplot(aes(x = fitness,y = N/Total,fill = bad)) + 
  geom_bar(alpha=1,
           stat = "identity") +
  theme_gerstung(base_size = 8) + 
  scale_x_discrete(expand = c(0,0)) + 
  scale_y_continuous(expand = c(0,0)) + 
  scale_fill_lancet(name=NULL) + 
  theme(legend.position = "bottom") + 
  ylab("Fraction") + 
  xlab("Fitness") 

explained_variants_plot_bar <- gof_plot_values %>% 
  group_by(bad) %>%
  summarise(N = length(bad)) %>%
  ungroup() %>%
  mutate(Freq = N / sum(N)) %>% 
  ggplot(aes(x = 1,
             y = Freq,
             fill = bad)) + 
  geom_bar(stat="identity") + 
  theme_gerstung() + 
  scale_x_discrete(expand = c(0,0)) + 
  scale_y_continuous(expand = c(0,0)) + 
  scale_fill_manual(values = rev(ggsci::pal_lancet()(2)),
                  breaks = c(F,T),
                  name = NULL,
                  labels = c("Good fit","Poor fit")) +
  theme(legend.position = "bottom") + 
  xlab("") + 
  ylab("Proportion of trajectories")

explained_variants_plot

The graphic shown above shows that a fairly good amount of trajectories (\(\sim87\%\)) are explained by our model. We exclude the non-explained trajectories from the calculation of age at clone foundation.

1.3.3 Genetic coefficients

predicted_variables <- variable_summaries(do.call(rbind,all_data$N200$draws)) %>% 
  as_tibble() %>%
  subset(grepl('genetic',variable)) %>% 
  mutate(variable_no = as.numeric(str_match(variable,'[0-9]+(?=\\])'))) %>% 
  mutate(true_fitness = unique(all_data$N200$simulated_samples$fitness)[variable_no]) %>% 
  mutate(true_fitness_per_year = true_fitness * 13 / age_mult) %>%
  spread(key = labels,value = values)
plot_min_max <- c(min(predicted_variables$HDPI_low,predicted_variables$true_fitness_per_year),
                  max(predicted_variables$HDPI_high,predicted_variables$true_fitness_per_year))
predicted_coefficients_plot <- predicted_variables %>%
  ggplot(aes(y = `0.50`,ymin = `HDPI_low`,ymax = `HDPI_high`,x = true_fitness_per_year)) + 
  geom_abline(slope=1,intercept = 0,linetype=2) +
  #geom_smooth(method = 'lm',formula = y ~ x,color = 'grey4',linetype=2,size=1,se=F) +
  geom_point(size = 0.5) + 
  geom_linerange(size = 0.2) +
  theme_gerstung(base_size = 6) + 
  theme(axis.title = element_text(size = 6)) + 
  xlab("Simulated genetic\nselection coefficient") + 
  ylab("Inferred genetic\nselection coefficient")

predicted_coefficients_plot

cor_test_inference <- cor.test(predicted_variables$true_fitness_per_year,predicted_variables$`0.50`)
cor_test_inference

## 
##  Pearson's product-moment correlation
## 
## data:  predicted_variables$true_fitness_per_year and predicted_variables$`0.50`
## t = 15.743, df = 23, p-value = 8.273e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9026530 0.9809485
## sample estimates:
##       cor 
## 0.9565979

cat(paste("Rsquared =",cor_test_inference$estimate^2))

## Rsquared = 0.915079532301844

The coefficients closely agree with the data, with high correlation (\(R^2>0.87\)) between the inferred and the true coefficients.

1.3.4 Age at clone foundation

The advantage of these simulations is that we have a very precise value for when the mutation first appeared. This allows us to assess how closely we can estimate the age at clone foundation.

calculate_age_estimates <- function(draws,r_values,min_age) {
  c_names <- colnames(draws$`11`)
  
  offset_names <- grep('^individual_offset',c_names,perl = T)
  genetic_c_names <- grep('^genetic_coef',c_names,perl = T)
  clone_c_names <- grep('^clone_specific_coef',c_names,perl = T)
  
  genetic_samples <- draws$`11`[,genetic_c_names]
  clone_samples <- draws$`11`[,clone_c_names]
  offset_samples <- draws$`11`[,offset_names]
  beta_samples <- draws$`11`[,length(c_names)] %>%
    unlist()
  
  sampling_idxs <- r_values %>%
    filter(sum_stat == T) %>% 
    group_by(
      fitness,fitness_factor,
      Individual,individual_factor,
      true_clone_age,
      Mutation
    ) %>% 
    summarise(
      minAge = min(age[which(round(true_count/coverage,4) > 0.005)]),
      maxAge = max(age),
      maxVAF = max(true_count/coverage),
      minVAF = min((true_count/coverage)[which(round(true_count/coverage,4) >= 0.005)])
    ) %>% 
    filter(!is.infinite(minAge)) %>%
    ungroup %>%
    mutate(true_clone_age_intervals = floor(true_clone_age/10) * 10)
  
  age_data_list <- list()
  i <- 0
  for (gen_time in c(0.5,1,2,5,10,13,20)) {
    for (population in c(10e3,5e4,1e5,2e5,600e3)) {
        i <- i + 1
    age_estimates <- sampling_idxs %>% 
      apply(1,FUN = function(x){
        genetic_numeric <- as.numeric(as.character(x[2]))
        clone_numeric <- as.numeric(as.character(x[4]))
    
        b_genetic <- genetic_samples[,genetic_numeric]
        b_clone <- clone_samples[,clone_numeric]
        u <- offset_samples[,clone_numeric]
        t1 <- as.numeric(x[7]) + min_age
        t1 <- t1 / age_mult
        b <- b_genetic + b_clone
        age_distribution <- t0_adjusted(u,b_genetic+b_clone,gen_time / age_mult,population) + min_age
        age_distribution <- age_distribution / age_mult
        historical_b <- b[age_distribution > -1]
        age_distribution <- ifelse(age_distribution < -1,-1,age_distribution)
        age_distribution <- ifelse(age_distribution > t1, t1, age_distribution)
        age_distribution <- age_distribution[!is.na(age_distribution)]
        D <- density(age_distribution)
        MAP <- D$x[which.max(D$y)]
        return(c(MAP,
                 mean(age_distribution),
                 quantile(age_distribution,c(0.05,0.10,0.16,0.5,0.84,0.90,0.95)),
                 mean(historical_b,na.rm = T),
                 quantile(historical_b,c(0.05,0.5,0.95),na.rm=T),
                 mean(b,na.rm = T),
                 quantile(b,c(0.05,0.5,0.95),na.rm=T)))
      }) %>% 
      t %>% 
      as_tibble() %>%
      mutate(
        gen_time = gen_time,
        population = population
      )
    
    colnames(age_estimates) <- c(
      "MAP","Mean","Q05","Q10","Q16","Q50","Q84","Q90","Q95",
      "Mean_HG","Q05_HG","Q50_HG","Q95_HG",
      "Mean_OG","Q05_OG","Q50_OG","Q95_OG",
      "gen_time","population")
    age_data_list[[i]] <- age_estimates
    }
  }
  return(list(age_data_list = age_data_list,sampling_idxs = sampling_idxs))
}

age_data_list_sampling_idxs <- calculate_age_estimates(
  draws = all_data$N200$draws,r_values = r_values_N200,min_age = min(all_data$N200$Data$Gen)/13)

## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

age_data_list_sampling_idxs_N100 <- calculate_age_estimates(
  draws = all_data$N100$draws,r_values = r_values_N100,min_age = min(all_data$N100$Data$Gen)/5)

## Warning in min(age[which(round(true_count/coverage, 4) > 0.005)]): no non-
## missing arguments to min; returning Inf

age_data_list_sampling_idxs_N100$age_data_full <- age_data_list_sampling_idxs_N100$age_data_list %>% 
  lapply(function(x) cbind(x,age_data_list_sampling_idxs_N100$sampling_idxs)) %>%
  do.call(what = rbind) %>%
  as_tibble()
age_data_list_sampling_idxs_N50 <- calculate_age_estimates(
  draws = all_data$N50$draws,r_values = r_values_N50,min_age = min(all_data$N50$Data$Gen)/1)
age_data_list_sampling_idxs_N50$age_data_full <- age_data_list_sampling_idxs_N50$age_data_list %>% 
  lapply(function(x) cbind(x,age_data_list_sampling_idxs_N50$sampling_idxs)) %>%
  do.call(what = rbind) %>%
  as_tibble()

age_data_list <- age_data_list_sampling_idxs$age_data_list
sampling_idxs <- age_data_list_sampling_idxs$sampling_idxs

age_data_full <- age_data_list %>% 
  lapply(function(x) cbind(x,sampling_idxs)) %>%
  do.call(what = rbind) %>%
  as_tibble()
age_data <- age_data_list[[29]] %>%
  cbind(sampling_idxs)

plot_min_max <- c(
  min(age_data$Q05,
      age_data$true_clone_age),
  max(age_data$Q05,
      age_data$true_clone_age)
)

estimated_ages_plot <- age_data %>% 
  ggplot(aes(x = true_clone_age,y = Q50)) + 
  geom_abline(slope=1,intercept=0,linetype=2,color = 'grey4') + 
  geom_point(aes(x = true_clone_age),
             colour = "grey",
             size = 0.5,alpha = 0.8) + 
  geom_linerange(aes(x = true_clone_age,
                     ymin = Q05,
                     ymax = Q95),
                 colour = "grey",
                 size = 0.2,
                 alpha = 0.8) + 
  geom_boxplot(
    aes(group = true_clone_age_intervals),
    alpha = 0.8,
    outlier.alpha = 0,
    size = 0.2
  ) + 
  scale_color_gradient(low = "pink",high = "red4",name = "Fitness") + 
  theme_gerstung(base_size = 6) + 
  ylab("Inferred age\nat clone foundation") +
  xlab("Simulated age\nat clone foundation") + 
  theme(legend.position = 'right') +
  coord_cartesian(ylim = plot_min_max) +
  scale_shape_discrete(name = "Generation time, population size") 

parameters_df <- r_values_N200 %>%
  select(genetic_effect,genetic_effect_005,genetic_effect_095,
         clone_effect,clone_effect_005,clone_effect_095,
         individual_offset,Mutation,Individual) %>%
  distinct

parameters_age_df <- merge(age_data,parameters_df,by = c("Mutation","Individual")) %>%
  mutate(minAge = minAge + min(all_data$N200$Data$Gen)/13,
         maxAge = maxAge + min(all_data$N200$Data$Gen)/13)
detection_threshold <- 0.002
N <- 2e5
age_threshold_accounting <- parameters_age_df %>%
  rowwise() %>%
  transmute(
    limit = onset_from_detection(clone_effect + genetic_effect,minVAF,mean(minAge,maxAge),g = 13,N = N),
    limit_005 = onset_from_detection(clone_effect_005 + genetic_effect_005,minVAF,mean(minAge,maxAge),13,N),
    limit_095 = onset_from_detection(clone_effect_095 + genetic_effect_095,minVAF,mean(minAge,maxAge),13,N),
    age_at_onset = Q50,
    coef = clone_effect + genetic_effect
  ) %>%
  subset(coef > 0) %>%
  mutate(plausible = limit >= 0,
         plausible_005 = limit_005 >= 0,
         plausible_095 = limit_095 >= 0)

estimated_ages_plot

table(age_threshold_accounting$plausible) / nrow(age_threshold_accounting)

## 
##     FALSE      TRUE 
## 0.1033058 0.8966942

cor_test_ages <- cor.test(age_data$true_clone_age,
                          age_data$Q50)

cat(paste("Mean absolute error =",mean(abs(age_data$Q50 - age_data$true_clone_age)),'\n'))

## Mean absolute error = 8.72497246458656

cat(paste("Rsquared =",cor_test_ages$estimate^2,'\n'))

## Rsquared = 0.658406026289533

cat(paste('p-value =',cor_test_ages$p.value))

## p-value = 8.90015725717923e-60

age_data_post <- age_data
age_data_post$residuals_005 <- age_data_post$true_clone_age - age_data_post$Q05
age_data_post$residuals <- age_data_post$true_clone_age - age_data_post$Q50
age_data_post$residuals_095 <- age_data_post$true_clone_age - age_data_post$Q95
age_data_post <- age_data_post %>%
  mutate(residuals_min = ifelse(residuals_095 > residuals_005,residuals_005,residuals_095),
         residuals_max = ifelse(residuals_095 < residuals_005,residuals_005,residuals_095)) 
ggplot(age_data_post,aes(x = fitness,y = residuals)) +
  geom_jitter(width = 0.0002,height = 0,size = 0.5) + 
  geom_boxplot(aes(group = cut(fitness,c(0,0.010,Inf))),
               alpha = 0.8,outlier.alpha = 0) + 
  theme_gerstung(base_size = 6)

The estimate for the ages, while not as good as the estimate for the coefficients, appears to be quite close to the true value. Additionally, on average, there is a difference of 6 years between our estimates and the true values. Residuals do not appear to be normally distributed, but a wider spread is observable for lower ages at clone foundation as per the figure above. Below we will observe that this is mostly caused by an unexpectedly high amount of time spent during the stochastic growth regime.

Next, we verify how much changing population sizes and time per generation affects estimates. To assess this we have calculated the values for different population sizes - 50,000 and 200,000 HSC - and for different generation times - 2 and 13 generations per year. For this, we observe that the median difference between all 4 possible combinations, is \(\sim19\) for estimates for the upper bound, \(\sim17\) for the median value, and \(\sim15\) for the lower bound.

age_data_full %>%
  group_by(Individual,true_clone_age) %>%
  summarise(m_005 = min(Q05),
            M_005 = max(Q05),
            m = min(Q50),
            M = max(Q50),
            m_095 = min(Q95),
            M_095 = max(Q95)) %>%
  na.omit() %>%
  mutate(Range_005 = M_005 - m_005,
         Range = M - m,
         Range_095 = M_095 - m_095) %>%
  ungroup() %>%
  summarise(Range_005Mean = mean(Range_005),
            RangeMean = mean(Range),
            Range_095Mean = mean(Range_095))

residual_plot <- age_data_full %>%
  group_by(Individual) %>%
  filter(length(unique(true_clone_age)) == 1) %>%
  ggplot(aes(x = reorder(Individual,true_clone_age),y = true_clone_age - Q50,
             color = paste(population))) + 
  geom_point(size = 0.5) +
  xlab("Index") +
  ylab("Residuals") +
  scale_colour_lancet(name = "Population size") +
  theme_gerstung(base_size = 6) +
  theme(axis.text.x = element_blank(),
        legend.position = 'bottom',
        legend.key.size = unit(0,"cm")) + 
  facet_wrap(~ gen_time,ncol = 1)

true_clone_age_plot <- age_data_full %>%
  group_by(Individual) %>%
  filter(length(unique(true_clone_age)) == 1) %>%
  ggplot(aes(x = reorder(Individual,true_clone_age),y = true_clone_age)) + 
  geom_point(size = 0.4) +
  theme_gerstung(base_size = 6) +
  theme(axis.text.x = element_blank()) + 
  xlab("") + 
  ylab("Age at onset")

age_data_full$label <- factor(
    sprintf("N = %i, g = %f",age_data_full$population,age_data_full$gen_time,sep=' ')) 
age_data_full %>% 
  ggplot(aes(x = true_clone_age,y = Q50,
             color = paste(population))) + 
  geom_abline(slope = 1) +
  geom_point(position = position_jitter(height = 0,width = 0),size = 0.4) +
  theme_gerstung(base_size = 6) +
  theme() + 
  xlab("Simulated age at onset") + 
  scale_colour_lancet(guide = F) +
  ylab("Inferred age at onset") + 
  theme(legend.position = "bottom",
        legend.key.size = unit(0,"cm")) + 
  geom_boxplot(aes(group = paste(true_clone_age,paste(gen_time,population,sep='\n')),
                   colour = paste(population,sep='\n')),
               size = 0.25,outlier.size = 0.5) + 
  facet_wrap(~ label,ncol = 5,dir = "v")

age_data_full %>%
  group_by(population,gen_time) %>%
  summarise(MAE = sum(abs(true_clone_age - Q50))/length(true_clone_age)) %>%
  mutate(Y = population == 2e5 & gen_time == 13) %>% 
  ggplot(aes(x = as.factor(population),y = as.factor(gen_time),fill = MAE,label = round(MAE,2))) + 
  geom_tile(aes(colour = Y),size = 0.5) + 
  geom_label(size = 2.6,label.size = unit(0,"cm"),label.r = unit(0,"cm"),fill = "white") + 
  theme_gerstung(base_size = 6) + 
  scale_fill_material(palette = "teal",name = "MAE") + 
  scale_colour_manual(values = c(NA,"black"),guide = F) + 
  xlab("Population size") + 
  ylab("Generations per year")

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

plot_grid(true_clone_age_plot,residual_plot,
          ncol = 1,align = "v",axis = "tlrb",
          rel_heights = c(0.15,1))

1.3.5 Rank analysis

We here begin by assessing that we can clearly observe mutations appearing consistently over life for all simulation conditions.

dir.create(path = "figures/rank_analysis",showWarnings = F)
plot_grid(
  all_data$N50$Data %>% 
    ggplot(aes(x = GenAtCloneFoundation)) +
    geom_density(width = 2) + 
    xlab("Age at onset") +
    theme_gerstung(base_size = 6) + 
    ggtitle("N = 50,000; Tg = 1"),
  all_data$N100$Data %>% 
    ggplot(aes(x = GenAtCloneFoundation / 5)) +
    geom_density(width = 2) + 
    xlab("Age at onset") +
    theme_gerstung(base_size = 6) + 
    ggtitle("N = 100,000; Tg = 5"),
  all_data$N200$Data %>% 
    ggplot(aes(x = GenAtCloneFoundation / 13)) +
    geom_density(width = 2) + 
    xlab("Age at onset") +
    theme_gerstung(base_size = 6) + 
    ggtitle("N = 200,000; Tg = 13"),
  nrow = 1
) + 
  ggsave("figures/rank_analysis/distributions_real.pdf",height = 1,width = 4)

## Warning: Ignoring unknown parameters: width

## Warning: Ignoring unknown parameters: width

## Warning: Ignoring unknown parameters: width

1.3.5.1 Rank distributions

ran <- function(x) {
  m <- 0#min(x,na.rm=T)
  M <- max(x,na.rm=T)
  d <- M-m
  if (d == 0) {
    d <- 1
  }
  return((x - m) / d)
}
extended_rank <- function(x,y) {
  return(order(order(x,y)))
}
distances_age_at_onset <- rbind(
  mutate(age_data_list_sampling_idxs_N50$age_data_full,true_pop_size = 50e3,true_gen_time = 1),
  mutate(age_data_list_sampling_idxs_N100$age_data_full,true_pop_size = 100e3,true_gen_time = 5),
  mutate(select(age_data_full,-label),true_pop_size = 200e3,true_gen_time = 13)
  ) %>%
  group_by(population,gen_time,true_pop_size,true_gen_time) %>%
  mutate(Rank = extended_rank(Q50,Q95), # helps us solve ties
         RankTrue = rank(true_clone_age,ties.method = "min")) %>%
  mutate(DistToFirst = Q50 - min(Q50)) %>% 
  mutate(DistToFirstTrue = true_clone_age - min(true_clone_age)) %>%
  na.omit() %>%
  mutate(Rank_ = ran(Rank),
         RankTrue_ = ran(RankTrue),
         DistToFirst_ = ran(DistToFirst)) %>%
  arrange(Q50) %>% 
  mutate(DistToPrevious = c(NA,diff(Q50))) %>%
  arrange(true_clone_age) %>%
  mutate(DistToPreviousTrue = c(NA,diff(true_clone_age))) %>%
  mutate(ExpectedDistToPrevious = rep(mean(DistToPrevious,na.rm=T),length(DistToFirst))) %>%
  mutate(ExpectedDistToFirst = cumsum(ExpectedDistToPrevious)) %>%
  mutate(RankExpected = rank(ExpectedDistToFirst,ties.method = "min")) %>%
  mutate(RankExpected_ = ran(RankExpected))

For all conditions we observe that using the correct estimates for population size and generation time renders a relatively uniform distribution of mutations through life for the distance to the earliest clone. The same is observable when analysing the rank of ages at onset. Importantly, there is a range of more or less likely solutions provided different combinations of population size and number of generations per year. This is unsurprising since we calculte the age at onset as \(age\ at\ onset = \frac{log(\frac{gn}{\beta}) - \alpha - 1}{\beta}\), where \(\alpha\) and \(\beta\) are the offset and growth coefficient, respectively, \(g\) is the number of generations per year and \(n\) is the population size. Thus, the same solution can be if \(g\) is multiplied by \(x\) while \(n\) is divided by \(x\) and vice-versa.

distance_distribution_N50 <- distances_age_at_onset %>%
  subset(true_pop_size == 50e3) %>%
  gather(key = "key",value = "value",DistToFirst,DistToFirstTrue,ExpectedDistToFirst) %>% 
  mutate(key = c(DistToFirst = "Inferred",DistToFirstTrue = "Simulated",ExpectedDistToFirst = "Uniform")[key]) %>% 
  ggplot(aes(x = value,colour = key)) +
  geom_density(size = 0.25) +
  facet_wrap(gen_time ~ population,
             ncol = 5) + 
  theme_gerstung(base_size = 6) + 
  xlab("Distance to earliest clone") + 
  ylab("Density") + 
  scale_colour_lancet(name = NULL) +
  scale_y_continuous(n.breaks = 3,expand = c(0,0)) +
  scale_x_continuous(breaks = c(0,30,60,90)) +
  theme(legend.position = "bottom",
        strip.text = element_text(margin = margin()),
        panel.spacing.x = unit(0.1,"cm"))

distance_distribution_N100 <- distances_age_at_onset %>%
  subset(true_pop_size == 100e3) %>%
  gather(key = "key",value = "value",DistToFirst,DistToFirstTrue,ExpectedDistToFirst) %>% 
  mutate(key = c(DistToFirst = "Inferred",DistToFirstTrue = "Simulated",ExpectedDistToFirst = "Uniform")[key]) %>% 
  ggplot(aes(x = value,colour = key)) +
  geom_density(size = 0.25) +
  facet_wrap(gen_time ~ population,
             ncol = 5) + 
  theme_gerstung(base_size = 6) + 
  xlab("Distance to earliest clone") + 
  ylab("Density") + 
  scale_colour_lancet(name = NULL) +
  scale_y_continuous(n.breaks = 3,expand = c(0,0)) +
  scale_x_continuous(breaks = c(0,30,60,90)) +
  theme(legend.position = "bottom",
        strip.text = element_text(margin = margin()),
        panel.spacing.x = unit(0.1,"cm"))

distance_distribution_N200 <- distances_age_at_onset %>%
  subset(true_pop_size == 200e3) %>%
  gather(key = "key",value = "value",DistToFirst,DistToFirstTrue,ExpectedDistToFirst) %>% 
  mutate(key = c(DistToFirst = "Inferred",DistToFirstTrue = "Simulated",ExpectedDistToFirst = "Uniform")[key]) %>% 
  ggplot(aes(x = value,colour = key)) +
  geom_density(size = 0.25) +
  facet_wrap(gen_time ~ population,
             ncol = 5) + 
  theme_gerstung(base_size = 6) + 
  xlab("Distance to earliest clone") + 
  ylab("Density") + 
  scale_colour_lancet(name = NULL) +
  scale_y_continuous(n.breaks = 3,expand = c(0,0)) +
  scale_x_continuous(breaks = c(0,30,60,90)) +
  theme(legend.position = "bottom",
        strip.text = element_text(margin = margin()),
        panel.spacing.x = unit(0.1,"cm"))

plot_grid(
  distance_distribution_N50 + ggtitle("N = 50,000; Tg = 1") + theme(legend.position = "none"),
  distance_distribution_N100 + ggtitle("N = 100,000; Tg = 5") + theme(legend.position = "none"),
  distance_distribution_N200 + ggtitle("N = 200,000; Tg = 13") + theme(legend.position = "none"),
  ggplot() + theme_nothing(),
  get_legend(distance_distribution_N200 + theme(legend.key.height = unit(0.1,"cm"))),
  ggplot() + theme_nothing(),
  rel_heights = c(1,0.05),
  ncol = 3) + 
  ggsave("figures/rank_analysis/distance_distributions.pdf",width = 7,height = 4)

rank_distribution_N50 <- distances_age_at_onset %>%
  subset(true_pop_size == 50e3) %>%
  ggplot(aes(alpha = population == 50e3 & gen_time == 1)) +
  geom_line(aes(x = Rank_,y = DistToFirst,colour = "Inferred")) +
  geom_line(aes(x = RankTrue_,y = DistToFirstTrue,colour = "Simulated")) +
  geom_line(aes(x = RankExpected_,y = ExpectedDistToFirst,colour = "Uniform"),
            size = 0.25) +
  facet_wrap(gen_time ~ population,
             ncol = 5) + 
  theme_gerstung(base_size = 6) + 
  xlab("Relative rank") + 
  ylab("Distance to first point") + 
  scale_colour_lancet(name = NULL) +
  scale_x_continuous(breaks = c(0,1)) +
  theme(legend.position = "bottom",
        strip.text = element_text(margin = margin()),
        panel.spacing.x = unit(0.1,"cm")) + 
  scale_y_continuous(breaks = c(0,25,50,75),expand = c(0,0)) +
  scale_alpha_discrete(range = c(0.2,1),guide = F) 

## Warning: Using alpha for a discrete variable is not advised.

rank_distribution_N100 <- distances_age_at_onset %>%
  subset(true_pop_size == 100e3) %>%
  ggplot(aes(alpha = population == 100e3 & gen_time == 5)) +
  geom_line(aes(x = Rank_,y = DistToFirst,colour = "Inferred")) +
  geom_line(aes(x = RankTrue_,y = DistToFirstTrue,colour = "Simulated")) +
  geom_line(aes(x = RankExpected_,y = ExpectedDistToFirst,colour = "Uniform"),
            size = 0.25) +
  facet_wrap(gen_time ~ population,
             ncol = 5) + 
  theme_gerstung(base_size = 6) + 
  xlab("Relative rank") + 
  ylab("Distance to first point") + 
  scale_colour_lancet(name = NULL) +
  scale_x_continuous(breaks = c(0,1)) +
  theme(legend.position = "bottom",
        strip.text = element_text(margin = margin()),
        panel.spacing.x = unit(0.1,"cm")) + 
  scale_y_continuous(breaks = c(0,25,50,75),expand = c(0,0)) +
  scale_alpha_discrete(range = c(0.2,1),guide = F) 

## Warning: Using alpha for a discrete variable is not advised.

rank_distribution_N200 <- distances_age_at_onset %>%
  subset(true_pop_size == 200e3) %>%
  ggplot(aes(alpha = population == 200e3 & gen_time == 13)) +
  geom_line(aes(x = Rank_,y = DistToFirst,colour = "Inferred")) +
  geom_line(aes(x = RankTrue_,y = DistToFirstTrue,colour = "Simulated")) +
  geom_line(aes(x = RankExpected_,y = ExpectedDistToFirst,colour = "Uniform"),
            size = 0.25) +
  facet_wrap(gen_time ~ population,
             ncol = 5) + 
  theme_gerstung(base_size = 6) + 
  xlab("Relative rank") + 
  ylab("Distance to first point") + 
  scale_colour_lancet(name = NULL) +
  scale_x_continuous(breaks = c(0,1)) +
  theme(legend.position = "bottom",
        strip.text = element_text(margin = margin()),
        panel.spacing.x = unit(0.1,"cm")) + 
  scale_y_continuous(breaks = c(0,25,50,75),expand = c(0,0)) +
  scale_alpha_discrete(range = c(0.2,1),guide = F) 

## Warning: Using alpha for a discrete variable is not advised.

plot_grid(
  rank_distribution_N50 + ggtitle("N = 50,000; Tg = 1") + theme(legend.position = "none"),
  rank_distribution_N100 + ggtitle("N = 100,000; Tg = 5") + theme(legend.position = "none"),
  rank_distribution_N200 + ggtitle("N = 200,000; Tg = 13") + theme(legend.position = "none"),
  ggplot() + theme_nothing(),
  get_legend(rank_distribution_N200),
  ggplot() + theme_nothing(),
  rel_heights = c(1,0.05),
  ncol = 3) + 
  ggsave("figures/rank_analysis/rank_distributions.pdf",width = 7,height = 4)

1.3.5.2 Calculating statistical distance between inferred and simulated ages at onset

dist_subset <- distances_age_at_onset %>%
  subset(population == true_pop_size & gen_time == true_gen_time)
d <- c(ran(dist_subset$DistToFirstTrue[dist_subset$true_pop_size == 200e3]),
       ran(dist_subset$DistToFirstTrue[dist_subset$true_pop_size == 100e3]),
       ran(dist_subset$DistToFirstTrue[dist_subset$true_pop_size == 50e3]))
d <- d[!is.na(d)]
age_at_onset_distribution <- fitdistrplus::fitdist(d,"beta",method = "mse")

hist_entropy <- function(x,M=100) {
  x <- x[!is.na(x)]
  D <- density(x,bw = M/10)
  x <- data.frame(dx = D$x,dy = D$y) %>%
    mutate(groups = cut(dx,seq(-1,M,length.out = 10))) %>%
    group_by(groups) %>%
    summarise(S = sum(dy))
  x <- x$S
  x <- x / sum(x)
  x <- x[x > 0]
  return(-sum(x * log(x)))
}

chisq_statistic <- distances_age_at_onset %>%
  #subset((fitness * true_gen_time) <= 0.05) %>% 
  group_by(population,gen_time,true_pop_size,true_gen_time) %>%
  na.omit() %>% 
  summarise(
    ChiSq = sum(qnorm(ran(DistToFirst)*(1-1e-6) + 1e-8)^2,na.rm=T),
    ChiSqTrue = sum(qnorm(ran(DistToFirst)*(1-1e-6) + 1e-8)^2,na.rm=T),
    Tau = cor(Rank,RankTrue,method = "kendall"),
    N = sum(!is.na(DistToPrevious)),
    ks_unif = suppressWarnings(ks.test(ran(DistToFirst),"punif",0,1)$statistic),
    ks_beta = suppressWarnings(ks.test(ran(DistToFirst),"pbeta",age_at_onset_distribution$estimate[1],age_at_onset_distribution$estimate[2])$statistic),
    HE = hist_entropy(ran(DistToFirst),M=1)) 

plot_stat_data <- chisq_statistic %>% 
  group_by(true_pop_size,true_gen_time) %>% 
  mutate(R = rank(ks_unif,ties.method = "min")) %>%
  ungroup %>% 
  mutate(Y = (gen_time == true_gen_time) & (population == true_pop_size)) 

plot_grid(
  plot_stat_data %>% 
    subset(true_pop_size == 50e3) %>%
    ggplot(aes(x = as.factor(population),y = as.factor(gen_time), fill = ks_unif)) +
    geom_tile(colour = NA) +
    geom_tile(aes(colour = Y),size = 0.25,fill = NA) +
    geom_text(aes(label = R),size = 2.3) +
    xlab("Population size") +
    ylab("Generations/year") +
    theme_gerstung(base_size = 6) +
    scale_fill_material(palette = "red",name = "D",n.breaks = 4) +
    scale_x_discrete(expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0)) +
    theme(legend.position = "bottom",
          legend.key.height = unit(0.2,"cm"),
          legend.key.width = unit(0.4,"cm"),
          strip.text = element_text(margin = margin()),
          panel.spacing = unit(0.1,"cm")) + 
    scale_color_manual(breaks = c(F,T),values = c(NA,"black"),guide = F) +
    ggtitle("N = 50,000; Tg = 1"),
  plot_stat_data %>% 
    subset(true_pop_size == 100e3) %>%
    ggplot(aes(x = as.factor(population),y = as.factor(gen_time), fill = ks_unif)) +
    geom_tile(colour = NA) +
    geom_tile(aes(colour = Y),size = 0.25,fill = NA) +
    geom_text(aes(label = R),size = 2.3) +
    xlab("Population size") +
    ylab("Generations/year") +
    theme_gerstung(base_size = 6) +
    scale_fill_material(palette = "red",name = "D",n.breaks = 4) +
    scale_x_discrete(expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0)) +
    theme(legend.position = "bottom",
          legend.key.height = unit(0.2,"cm"),
          legend.key.width = unit(0.4,"cm"),
          strip.text = element_text(margin = margin()),
          panel.spacing = unit(0.1,"cm")) + 
    scale_color_manual(breaks = c(F,T),values = c(NA,"black"),guide = F) +
    ggtitle("N = 100,000; Tg = 5"),
  plot_stat_data %>% 
    subset(true_pop_size == 200e3) %>%
    ggplot(aes(x = as.factor(population),y = as.factor(gen_time), fill = ks_unif)) +
    geom_tile(colour = NA) +
    geom_tile(aes(colour = Y),size = 0.25,fill = NA) +
    geom_text(aes(label = R),size = 2.3) +
    xlab("Population size") +
    ylab("Generations/year") +
    theme_gerstung(base_size = 6) +
    scale_fill_material(palette = "red",name = "D",n.breaks = 4) +
    scale_x_discrete(expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0)) +
    theme(legend.position = "bottom",
          legend.key.height = unit(0.2,"cm"),
          legend.key.width = unit(0.4,"cm"),
          strip.text = element_text(margin = margin()),
          panel.spacing = unit(0.1,"cm")) + 
    scale_color_manual(breaks = c(F,T),values = c(NA,"black"),guide = F) +
    ggtitle("N = 200,000; Tg = 13"),
  ncol = 3
) + 
  ggsave("figures/rank_analysis/statistical_distance_d.pdf",height = 1.5,width = 5)

1.3.6 Inspecting poor age at onset predictions

Below we inspect the clones for which the difference between predicted and true age at clone foundation differed more than by 20 years.

root <- "/hps/nobackup/research/gerstung/josegcpa/projects/05VAF_DYNAMICS/hsc_output_200k_13/hsc_%s_%s/r00%s.csv"

individuals_mutations <- age_data_post %>% 
  mutate(difference = abs(true_clone_age - Mean)) %>% 
  subset(difference > 20) %>% 
  arrange(difference)
file_ids <- str_split(individuals_mutations$Individual,pattern = '_') %>%
  do.call(what = rbind) %>%
  as.data.frame()
colnames(file_ids) <- c("fitness","mutation_rate","replicate")
file_ids <- file_ids %>% 
  mutate(mutation_rate = as.numeric(str_match(file_ids$mutation_rate,'[0-9.]+')),
         div = as.numeric(str_match(file_ids$mutation_rate,'[0-9.]+$')) - 9
  ) %>%
  mutate(mutation_rate = mutation_rate / (10^div)) %>% 
  select(-div)
file_ids$Mutation <- individuals_mutations$Mutation
file_ids$bad <- T

bad_predictions <- file_ids %>%
  select(-Mutation) %>%
  apply(1,function(x) {
    x <- as.numeric(x)
    if (x[2] < 1) {
      mr <- as.character(format(x[2],nsmall=3))
    } else if (x[2] >= 5) {
      mr <- as.character(format(x[2]))
    } else {
      mr <- as.character(paste0(0,format(x[2],nsmall=3)))
    }
    while (nchar(x[1]) < 6) {
      x[1] <- paste0(x[1],"0")
      }
    file_name <- sprintf(root,x[1],mr,x[3])
    tmp <- read_tsv(file_name,
                    col_names = c("Gen","Count","Clone","Mutation"),
                    col_types = list(col_number(),col_number(),
                                     col_number(),col_number()),
                    progress = F) 
    tmp <- tmp %>%
      mutate(fitness = x[1],
             mutation_rate = x[2],
             replicate = x[3])
    return(tmp)
  }) %>%
  do.call(what = rbind) %>% 
  mutate_all(as.numeric) %>% 
  filter(Mutation > 0) %>%
  group_by(Mutation,Gen,fitness,mutation_rate,replicate) %>%
  summarise(Count = max(Count)) %>%
  group_by(Mutation) %>%
  filter(length(Count) > 10) %>% 
  mutate(Driver = Mutation <= (N_DRIVERS)) %>% 
  mutate(Alpha = ifelse(Driver,0.9,0.05)) %>% 
  ungroup()

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

## Warning in FUN(newX[, i], ...): NAs introduced by coercion

bad_predictions <- merge(
  bad_predictions,file_ids,by = c("fitness","mutation_rate","Mutation","replicate"),all.x=T) %>%
  mutate(bad = ifelse(is.na(bad),F,bad))

bad_predictions <- bad_predictions %>% 
  arrange(Mutation,Gen) %>% 
  group_by(Mutation) %>%
  filter(length(Gen) > 1) %>%
  mutate(Break = Gen - c(0,Gen[1:max(0,(length(Gen)-1))]) > 20) %>% 
  mutate(Component = cumsum(Break)) %>%
  ungroup() %>%
  mutate(Mutation_C = paste(Mutation,Component,sep = '_')) 

bad_predictions %>% 
  group_by(mutation_rate,fitness,replicate) %>%
  filter(Mutation %in% c(1:50,sample(unique(Mutation),10))) %>%
  ggplot(aes(x = Gen,y = Count,color = Driver,alpha=Alpha,
             group = paste(Driver,Mutation_C,replicate),
             size = as.numeric(Driver))) + 
  geom_hline(aes(yintercept = (1 - (1+fitness)^-1) / (1 - (1+fitness)^-pop_size) * pop_size),
             size = 0.25) +
  geom_line(aes(linetype = bad)) + 
  scale_y_continuous(trans = 'log',
                     breaks = c(10^seq(0,5)),
                     limits = c(10,2e5)) + 
  scale_x_continuous(limits = c(0,1300)) + 
  scale_alpha(guide = F) + 
  theme_gerstung(base_size = 6) + 
  scale_color_d3(name = NULL,labels = c("Passenger","Driver")) + 
  theme(legend.position = "bottom",strip.text = element_text(margin = margin())) + 
  scale_size(labels = c(FALSE,TRUE),
             breaks = c(FALSE,TRUE),
             limits = c(FALSE,TRUE),
             range = c(0.25,0.5),
             guide = F) + 
  ylab("Count") + 
  facet_wrap(~ fitness + mutation_rate + replicate,nrow = 3,scales = 'free') + 
  ggtitle("Trajectories with large errors for the inferred age at clone foundation",
          "Mean absolute error > 20 years")

## Warning: Removed 1817 row(s) containing missing values (geom_path).

1.4 Historical growth

L <- c(min(age_data_full$Q50_HG,age_data_full$Q50_OG),
       max(age_data_full$Q50_HG,age_data_full$Q50_OG))
age_data_full %>% 
  ggplot(aes(x = Q50_OG,y = Q50_HG, colour = fitness)) + 
  geom_abline(slope = 1,size = 0.25) + 
  geom_point(size = 0.5) + 
  facet_grid(gen_time ~ population,scales = "free") + 
  theme_gerstung(base_size = 6) + 
  ylab("Historical growth") + 
  xlab("Observed growth") + 
  scale_x_continuous(limits = L) + 
  scale_y_continuous(limits = L) + 
  scale_color_material(palette = "deep-purple") +
  theme(legend.position = "bottom",
        legend.key.height = unit(0.2,"cm"),
        legend.title = element_text(vjust = 1))

## Warning: Removed 34 rows containing missing values (geom_point).

age_data_full %>%
  group_by(fitness,population,gen_time) %>%
  summarise(ExpectedHistoricalGrowth = mean(exp(Q50_HG) - 1,na.rm=T),
            ExpectedObservedGrowth = mean(exp(Q50_OG) - 1,na.rm=T)) %>%
  ggplot(aes(x = fitness * 13,
             y = (ExpectedHistoricalGrowth+ExpectedObservedGrowth)/2,
             ymin = ExpectedObservedGrowth,
             ymax = ExpectedHistoricalGrowth,
             colour = as.factor(population),
             width = 0.2)) + 
  geom_abline(slope = 1,size = 0.25) +
  geom_point(size = 0.3,shape = 3,position = position_dodge(width = 0.01)) + 
  geom_linerange(size = 0.2,position = position_dodge(width = 0.01)) +
  facet_grid(. ~ gen_time) + 
  theme_gerstung(base_size = 6) +
  xlab("Fitness") + 
  ylab("Growth effect") + 
  coord_flip() +
  scale_colour_manual(values = c("red4","red","goldenrod","forestgreen","blue"),
                      name = "N") +
  theme(legend.position = "bottom",
        legend.key.height = unit(0.2,"cm"),
        legend.title = element_text(vjust = 1))

age_data_full %>% 
  group_by(fitness,population,gen_time) %>%
  summarise(ExpectedHistoricalGrowth = mean(exp(Q50_HG) - 1,na.rm=T),
            ExpectedObservedGrowth = mean(exp(Q50_OG) - 1,na.rm=T)) %>%
  mutate(D = ExpectedHistoricalGrowth - ExpectedObservedGrowth) %>% 
  mutate(fitness = factor(fitness,
                          levels = unique(fitness),
                          labels = paste("fitness =",unique(fitness)))) %>% 
  ggplot(aes(x = as.factor(population),y = as.factor(gen_time),fill = D)) + 
  geom_tile(colour = "white") + 
  theme_gerstung(base_size = 6) +
  facet_wrap(~ fitness,ncol = 5) +
  scale_fill_gradientn(colours = c("goldenrod","red4"),
                       name = "Hist. - obs. growth") +
  xlab("N") +
  ylab("Generations/year") + 
  theme(legend.position = "bottom",
        legend.key.height = unit(0.2,"cm"),
        legend.title = element_text(vjust = 1)) + 
  scale_x_discrete(expand = c(0,0)) +
  scale_y_discrete(expand = c(0,0))

1.5 Generating a final summary figure

fixation_probility <- function(s,N) {
  return(
    1 / (N * s)
  )
}

counts <- all_data$N200$Data %>%
  group_by(fitness,mutation_rate,Replicate) %>%
  summarise(L = length(unique(Mutation_C)),
            M = max(VAF),
            minGen = abs(GenAtCloneFoundation - 600)[1]) %>%
  filter(L < 2) %>%
  filter(fitness >= 0.005) 
chosen_one_full <- counts %>% 
  arrange(minGen,M,L)

good_examples <- chosen_one_full#[c(129),]

plot_list <- list()
for (i in 1:nrow(good_examples)) {
  chosen_one <- good_examples[i,]
  params <- c(f = chosen_one$fitness,
              m = chosen_one$mutation_rate/(10^-9),
              r = chosen_one$Replicate)

  if (params[2] < 1) {
        params[2] <- as.character(format(params[2],nsmall=3))
      } else if (params[2] >= 4.5) {
        params[2] <- format(params[2])
      } else {
        params[2] <- as.character(paste0(0,format(params[2],nsmall=3)))
      }
  
  while (nchar(params[1]) < 6) {
    params[1] <- paste0(params[1],"0")
  }
  
  file <- sprintf(root,params[1],params[2],as.numeric(params[3]))

  df <- read_tsv(file,
                 col_names = c("Gen","Count","Clone","Mutation"),
                 col_types = list(col_number(),col_number(),
                                  col_number(),col_number()),
                 progress = F)

  fit_adv <- as.numeric(str_match_all(file,"[0-9.]+")[[1]][4])
  population_threshold <- fixation_probility(fit_adv,pop_size)
  time_stochastic <- 1 / fit_adv
  
  example_trajectory <- df %>% 
    filter(Mutation > 0) %>%
    group_by(Mutation,Gen) %>%
    summarise(Count = sum(Count)) %>%
    group_by(Mutation) %>%
    filter(length(Count) > 1) %>% 
    mutate(Driver = Mutation <= (N_DRIVERS)) %>% 
    mutate(Alpha = ifelse(Driver,0.9,0.02)) %>% 
    mutate(Break = Gen - c(0,Gen[1:max(0,(length(Gen)-1))]) > 5) %>% 
    mutate(Component = cumsum(Break)) %>%
    ungroup() %>%
    mutate(Mutation = paste(Mutation,Component,sep = '_')) %>%
    group_by(Mutation) %>%
    mutate(Count = Count/2e5) %>%
    mutate(Deterministic = ifelse(Count > population_threshold & Driver == T,
                                  Count,
                                  NA)) 
  
  tmp <- example_trajectory %>%
    subset(Driver == T) %>% 
    group_by(Mutation) %>% 
    summarise(minGen = min(Gen),
              C = Count[which.min(Gen)]) %>%
    filter(C > 5e-6)
  
  example_trajectory <- rbind(
    example_trajectory,
    tmp %>% 
      apply(1,function(x) example_trajectory %>% 
              subset(Gen == x[2] & Mutation == x[1]) %>% 
              mutate(Gen = Gen - 5,Count = 5e-6)) %>% 
      do.call(what = rbind))
  
  dead_clones <- example_trajectory %>% 
    group_by(Mutation) %>%
    mutate(GenDiff = max(Gen) - min(Gen)) %>%
    ungroup %>% 
    filter(min(Count) == 5e-6,
           max(Count) < population_threshold,
           max(Gen) - min(Gen) > 5) %>%
    filter(Driver == T) %>% 
    filter(GenDiff == max(GenDiff)) %>%
    mutate(Count = ifelse(Gen < Gen[min(which(Count == 5e-6))],5e-6,Count)) 
  
  big_clone <- example_trajectory$Mutation[which.max(example_trajectory$Count)]
  big_clone_fixation <- example_trajectory$Gen[
    which.min(which(example_trajectory$Count > population_threshold))]
  
  if (T %in% example_trajectory$Driver) {
    simulated_gen_at_foundation <- example_trajectory %>%
      subset(Driver == T) %>%
      summarise(m = min(Gen),
                M = max(Count)) %>%
      ungroup() %>% 
      subset(M == max(M))
    inferred_gen_at_foundation <- example_trajectory %>%
      subset(Driver == T & Mutation == simulated_gen_at_foundation$Mutation) %>% 
      summarise(m = min(Gen) - 5,
                M = max(Count)) %>%
      ungroup() %>% 
      subset(M == max(M)) %>% 
      mutate(InferredAgeAtCloneFoundation = m - time_stochastic) 
  
    age_at_clone_foundation_info <- data.frame(
      x = c(inferred_gen_at_foundation$InferredAgeAtCloneFoundation,
            simulated_gen_at_foundation$m),
      label = c("Inferred age at clone foundation",
                "Simulated age at clone foundation"),
      y_label = c(1e-5,1e-4),
      y_points = c(6.5e-6)
    )
    example_trajectory_neutral <- example_trajectory %>% 
      subset(Driver == F)
    example_trajectory_driver <- example_trajectory %>%
      subset(Driver == T)
    
    tmp <- example_trajectory_driver %>% 
      subset(!is.na(Deterministic) & Mutation == inferred_gen_at_foundation$Mutation) 
    if (nrow(tmp) > 0) {
      
      driver_fit_coefficient <- lsfit(tmp$Gen / 13,logit(tmp$Count*0.999999))$coefficients
      points_meeting_threshold <- c(
        x=(
          logit(
            population_threshold) - driver_fit_coefficient[1])/driver_fit_coefficient[2],
        y=population_threshold
        )
  
      points_fit <- data.frame(
        x = c(
          points_meeting_threshold[1]-time_stochastic/13,
          points_meeting_threshold[1],1300/13),
        y = c(
          5e-6,points_meeting_threshold[2],
          inv.logit(100 * driver_fit_coefficient[2] + driver_fit_coefficient[1])),
        colour = "Fit"
      )
      points_fit_plot <- rbind(
        data.frame(
          x = seq(points_fit$x[1],points_fit$x[2],length.out = 100),
          y = seq(points_fit$y[1],points_fit$y[2],length.out = 100),
          colour = "Fit"
        ),
        data.frame(
          x = seq(points_fit$x[2],100,by = 0.1),
          y = inv.logit(
            seq(points_fit$x[2],
                100,by = 0.1) * driver_fit_coefficient[2] + driver_fit_coefficient[1]),
          colour = "Fit"
        )
      )
    
      example_trajectory_plot <- example_trajectory_neutral %>% 
        ggplot(aes(x = Gen / 13,
                   y = Count,
                   alpha=Alpha,
                   color = Driver,
                   group = Mutation)) + 
        geom_line(aes(size = as.numeric(Driver)),alpha = 1,colour = "#ffe08c") + 
        geom_hline(yintercept = population_threshold,size = 0.5) +
        geom_hline(yintercept = 5e-6,linetype="dotted",size=0.5) +
        # geom_line(data = dead_clones,
        #           size = 1.5,
        #           colour = 'orange',
        #           alpha = 1) +
        geom_line(
          data = points_fit_plot,
          inherit.aes = F,
          aes(x = x,y = y,colour = colour),
          #color = '#d68b82',
          size = 1
        ) +
        geom_point(
          data = example_trajectory_driver,
          alpha = 1
        ) + 
        geom_line(
          data = example_trajectory_driver,
          size = 0.4,
          alpha = 1,
          linetype = 5
        ) + 
        annotate(geom = 'label',
                 x = 20 / 13,
                 y = 10^((log10(population_threshold)+log10(1/2e5))/2),
                 label = "Stochastic\ngrowth",
                 hjust = 0,
                 size = 1.9,
                 label.size = unit(0,"cm"),
                 label.r = unit(0,"cm"),
                 alpha = 0.8) +
        annotate(geom = 'label',
                 x = 100 / 13,
                 size = 1.9,
                 y = sqrt(population_threshold),
                 label = "Deterministic\ngrowth",
                 hjust = 0,
                 label.size = unit(0,"cm"),
                 label.r = unit(0,"cm"),
                 alpha = 0.8) +
        annotate(geom = 'label',
                 x = 1300 / 13,
                 y = population_threshold,
                 label = "drift threshold = 1/(Nf)",
                 size = 1.9,
                 hjust = 1.1,
                 vjust = -0.2,
                 label.size = unit(0,"cm"),
                 label.r = unit(0,"cm"),
                 label.padding = unit(0,"cm"),
                 alpha = 0.8) +
        annotate(geom = 'label',
               x = 1300 / 13,
               size=1.9,
               y = population_threshold/100,
               label = "clones that do not\ncross the drift\nthreshold tend to\ndisappear",
               hjust = 1.1,
               vjust = -0.2,
               label.size = unit(0,"cm"),
               label.r = unit(0,"cm"),
               label.padding = unit(0,"cm"),
               alpha = 0.8) +
        annotate(
          geom = 'label',
          x = 600 / 13,
          size = 1.9,
          y = 0.03,
          label=(paste(expression("counts ~ e"^"ft"*""))),
          parse = TRUE,
          hjust = 1,
          vjust = -0.2,
          label.size = unit(0,"cm"),
          label.r = unit(0,"cm"),
          fontface="italic",
          alpha = 0.8) +
        geom_errorbar(
          inherit.aes = F,
          data = data.frame(
            ymin = c(min(example_trajectory$Count),population_threshold),
            ymax = c(population_threshold,1.0),
            x = c(20 / 13,100 / 13)
          ),
          aes(ymin = ymin, 
              ymax = ymax,
              x = x,
              group = x),
          color = 'black',
          size = 0.3) + 
        scale_y_continuous(trans = 'log10',
                           breaks = c(5e-6,c(10^c(-5:0))),
                           limits = c(5e-6,1),
                           labels = c("Single cell",10^c(-5:0)),
                           expand = c(0.01,0.0,0.0,0.0)) + 
        scale_x_continuous(expand = c(0,1)) +
        scale_alpha(guide = F) + 
        theme_gerstung(base_size = 6) + 
        scale_color_manual(name = NULL,
                           breaks = c(T,F,"Fit"),
                           drop = F,
                           values = c("#ab2503","#00468BFF","#d68b82"),
                           labels = c("Driver","Neutral","Fit")) + 
        theme(legend.position = "bottom",
              axis.title = element_text(size = 6)) + 
        scale_size(breaks = c(FALSE,T),
                   limits = c(FALSE,T),
                   labels = c("Neutral","Driver"),
                   range = c(0.4,0.45),
                   name = NULL) +
        ylab("VAF") +
        xlab("t (years)")
      LL <- length(unique(example_trajectory_driver$Mutation))
        if (LL <= 4 & LL > 1) {
          plot_list[[file]] <- example_trajectory_plot
        }
    }
  }
}

## Warning in min(Gen): no non-missing arguments to min; returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in min(Gen): no non-missing arguments to min; returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

## Warning in max(GenDiff): no non-missing arguments to max; returning -Inf

## Warning: Problem with `mutate()` input `Count`.
## ℹ no non-missing arguments to min; returning Inf
## ℹ Input `Count` is `ifelse(Gen < Gen[min(which(Count == 5e-06))], 5e-06, Count)`.

## Warning in min(which(Count == 5e-06)): no non-missing arguments to min;
## returning Inf

main_results <- plot_grid(
  example_trajectory_plot + theme(legend.position = c(0.5,0.9),legend.direction = "horizontal",
                                  legend.key.height = unit(0.2,"cm"),
                                  axis.text = element_text(size = 6)),
  plot_grid(
    predicted_coefficients_plot + 
      xlab("Simulated driver effect") + 
      ylab("Inferred driver effect") + 
      theme(axis.text = element_text(size = 6)),
    estimated_ages_plot + 
      xlab("Simulated age at onset") + 
      ylab("Inferred age at onset") + 
      theme(axis.text = element_text(size = 6)),
    ncol = 1
  ),
  nrow = 1,
  rel_widths = c(7,3.5)
)

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

main_results

ggsave(main_results,filename = "figures/simulations/simulation_plot.pdf",height = 3,width = 6,
       useDingbats=FALSE)
ggsave(example_trajectory_plot + theme(legend.position = "top"),
       filename = "figures/simulations/example_trajectory_plot.pdf",height = 3,width = 3)

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

## Warning in numnotnull("lwd"): NAs introduced by coercion

ggsave(predicted_coefficients_plot,
       filename = "figures/simulations/predicted_coefficients_plot.pdf",height = 2,width = 2)
ggsave(explained_variants_plot + theme(legend.position = "top"),
       filename = "figures/simulations/explained_variants_plot.pdf",height = 2,width = 2)
ggsave(estimated_ages_plot + theme(legend.position = "right"),
       filename = "figures/simulations/estimated_ages_plot.pdf",height = 2,width = 2)

2 Session details

cat(system("git log -n1", intern=TRUE), sep="\n")

## commit 5bcb05932603c8c4a9756c112d35c89f3c41c92b
## Author: josegcpa <josegcpa@ebi.ac.uk>
## Date:   Wed Aug 11 19:08:29 2021 +0100
## 
##     removed redundant parts of the analysis (regarding age at onset)

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

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux
## 
## Matrix products: default
## BLAS:   /hps/nobackup/research/gerstung/josegcpa/software/R/3.6.3/lib64/R/lib/libRblas.so
## LAPACK: /hps/nobackup/research/gerstung/josegcpa/software/R/3.6.3/lib64/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
##  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
##  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] survminer_0.4.8   survival_3.2-7    reghelper_1.0.1   scatterpie_0.1.5 
##  [5] ggtree_2.0.4      ape_5.4-1         dendextend_1.14.0 default_1.0.0    
##  [9] extraDistr_1.9.1  ggrepel_0.8.2     ggsci_2.9         cowplot_1.1.0    
## [13] gtools_3.8.2      openxlsx_4.2.2    bayesplot_1.7.2   forcats_0.5.0    
## [17] stringr_1.4.0     dplyr_1.0.2       purrr_0.3.4       readr_1.4.0      
## [21] tidyr_1.1.2       tibble_3.0.4      tidyverse_1.3.0   ggpubr_0.4.0     
## [25] ggplot2_3.3.2     greta_0.3.1       reticulate_1.16  
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_1.4-1    ggsignif_0.6.0      ellipsis_0.3.1     
##   [4] rio_0.5.16          ggridges_0.5.2      base64enc_0.1-3    
##   [7] fs_1.5.0            dichromat_2.0-0     rstudioapi_0.11    
##  [10] farver_2.0.3        listenv_0.8.0       fansi_0.4.1        
##  [13] lubridate_1.7.9     xml2_1.3.2          splines_3.6.3      
##  [16] codetools_0.2-16    knitr_1.30          polyclip_1.10-0    
##  [19] jsonlite_1.7.1      km.ci_0.5-2         broom_0.7.1        
##  [22] dbplyr_1.4.4        tfruns_1.4          ggforce_0.3.2      
##  [25] mapproj_1.2.7       BiocManager_1.30.10 compiler_3.6.3     
##  [28] httr_1.4.2          rvcheck_0.1.8       backports_1.1.10   
##  [31] assertthat_0.2.1    Matrix_1.2-18       lazyeval_0.2.2     
##  [34] cli_2.1.0           tweenr_1.0.1        htmltools_0.5.0    
##  [37] prettyunits_1.1.1   tools_3.6.3         coda_0.19-4        
##  [40] gtable_0.3.0        glue_1.4.2          reshape2_1.4.4     
##  [43] maps_3.3.0          Rcpp_1.0.5          carData_3.0-4      
##  [46] cellranger_1.1.0    vctrs_0.3.4         nlme_3.1-144       
##  [49] xfun_0.18           globals_0.13.0      ps_1.4.0           
##  [52] rvest_0.3.6         lifecycle_0.2.0     rstatix_0.6.0      
##  [55] future_1.19.1       zoo_1.8-8           MASS_7.3-51.5      
##  [58] scales_1.1.1        hms_0.5.3           parallel_3.6.3     
##  [61] yaml_2.2.1          curl_4.3            gridExtra_2.3      
##  [64] KMsurv_0.1-5        stringi_1.5.3       tensorflow_2.2.0   
##  [67] tidytree_0.3.3      zip_2.1.1           pals_1.6           
##  [70] rlang_0.4.8         pkgconfig_2.0.3     evaluate_0.14      
##  [73] lattice_0.20-38     labeling_0.4.2      treeio_1.10.0      
##  [76] tidyselect_1.1.0    plyr_1.8.6          magrittr_1.5       
##  [79] R6_2.4.1            generics_0.0.2      DBI_1.1.0          
##  [82] pillar_1.4.6        haven_2.3.1         whisker_0.4        
##  [85] foreign_0.8-75      withr_2.3.0         fitdistrplus_1.1-1 
##  [88] abind_1.4-5         modelr_0.1.8        crayon_1.3.4       
##  [91] car_3.0-10          survMisc_0.5.5      rmarkdown_2.4      
##  [94] viridis_0.5.1       progress_1.2.2      readxl_1.3.1       
##  [97] data.table_1.13.0   blob_1.2.1          reprex_0.3.0       
## [100] digest_0.6.26       xtable_1.8-4        munsell_0.5.0      
## [103] viridisLite_0.3.0