prpf8.Rmd

---
title: "Prpf8"
output: 
  bookdown::word_document2:
    toc: true
    number_sections: false
    reference_docx: "output_template.docx"
---

```{r setup, echo=FALSE, message=FALSE, warning = FALSE, results = "hide"}
knitr::opts_chunk$set(echo = FALSE, fig.width = 8, fig.height = 5, cache = TRUE)
library(tidyverse)
library(bayesplot)
library(tidybayes)
library(ggdist)
library(brms)
library(patchwork)
options(mc.cores = 4, brms.backend = "cmdstanr")
theme_set(cowplot::theme_cowplot())
devtools::load_all()

options(width = 160)

document_output <- isTRUE(getOption('knitr.in.progress'))

if(document_output) {
  table_format <- function(x, caption = NULL, align = NULL) { 
    knitr::kable(x, caption = caption, align = align) 
  }
} else {
  table_format <- function(x, caption = NULL, align = NULL) { x }
}

housekeeping <- c("Ubb", "Gapdh", "Actb") 
pp_checks <- FALSE # Uncomment this to show all the pp checks
```

For each of the genotype-timepoint combinations we will follow the same structure: show exploratory plots, then summary of the fitted model and finally look at comparisons. The details will be explained for the d17 - 4w case and the explanation not repeated for the additional cases.

# d17 4w

```{r}
f <- here::here("private_data/qPCR/qPCR Prpf8 d17 4w final.xlsx")
```

```{r}
a1 <- read_qpcr(f, "qPCR set a1", range = "A5:AR31")
a2 <- read_qpcr(f, "qPCR set a2", range = "A5:Z22")
a3 <- read_qpcr(f, "qPCR set a3", range = "A5:AU22")
a4 <- read_qpcr(f, "qPCR set a4", range = "A5:AC56")

d17_4w <- rbind(a1, a2, a3, a4)
```


```{r d17-4w-box-housekeeping, fig.cap = "Cq values of the housekeeping genes stratified by genotype and sex. Each dot is a single technical replicate."}
genes_boxplot(d17_4w, housekeeping)
```


```{r d17-rw-parcoord-housekeeping, fig.cap = "Average Cq number across technical replicates for the housekeeping genes in each biological replicate. Each line connects measurements from the same biological replicate."}
genes_parcoord(d17_4w, housekeeping)
```


The set of plots below shows Cq numbers for all the measured primers stratified by sex and genotype. Each dot is a single technical replicate.

```{r, fig.width=10, fig.height=6}
to_compare <- setdiff(unique(d17_4w$primer_short), housekeeping)
genes_boxplot(d17_4w, to_compare, step_size = 8)
```

## Fitting the model

```{r}
#d17_4w_model <- d17_4w 
 d17_4w_model <- data_for_model(d17_4w, housekeeping)


priors = c(prior(normal(0,2), class = "sd"),
           prior(normal(0,8), class = "b"))

fit_d17_4w <- brm(
    bf(Cq ~ primer_short * genotype + sex*primer_short + (1 || run + animal_no + animal_no:primer_short + litter:primer_short),
       sigma ~ (1 || primer_short + run)
    ),
    data = d17_4w_model,
    control = list(max_treedepth = 11),
    prior = priors, refresh = 500,
    file = here::here("local_temp_data/fit_d17_4w"), file_refit = "on_change")


fit_d17_4w

nuts_params(fit_d17_4w) %>% filter(Parameter == "treedepth__") %>%
  group_by(Value) %>% summarise(n())
```


Posterior predictive checks to asses model fit.

```{r}
if(pp_checks) {
  show_pp_checks(d17_4w_model, fit_d17_4w)
}
```

## Comparisons - alternative forms

The basic quantity of interest are the ratios of two products, we treat one product as numerator and one as denominator. When detecting circRNAs we take the numerator to be the primer that should be specific to the circRNA and denominator as the primer that should correspond to the linear form, i.e. the larger ratio for circRNA, the higher proportion of the products are circRNA. For detecting alternative exons, we take the primer including the alternative exons as the numerator, i.e. the the larger ratio for alternative exon, the more products contain the alternative exon. 

Those are the ratios products we use:

```{r}
comparisons_raw <-  readxl::read_excel(here::here("private_data/qPCR/qPCR selected amplicons.xlsx"), sheet = "qPCR amplicons to be compared", range = "T2:W40") %>%
  filter(!is.na(label))

comparisons_d17_4w <- comparisons_raw %>% 
  process_comparisons(data_model = d17_4w_model)
  
eff_variants <- data.frame(eff_numerator = c(1.5, 1.8, 2), eff_denominator = c(2, 1.8, 1.5), eff_label = factor(1:3, labels = c("Favor denom", "Neutral", "Favor num")))

comparisons_d17_4w_eff <- comparisons_d17_4w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())

comparisons_d17_4w %>% select(label, variant, numerator, denominator) %>% table_format()

```


```{r}
comparisons_pred_d17_4w <- predict_comparisons(fit_d17_4w, d17_4w, comparisons_d17_4w_eff, c("wt", "wt"), c("d17", "het"))
ratios_observed_d17_4w <- compute_ratios_observed(d17_4w_model, comparisons_d17_4w)
```


Below is a set of somewhat complex plots showing both A) the observed ratios of products within a biological sample for all of the comparisons modelled (left) and B) model-derived "ratios of ratios" between genotypes (right). 

Assuming equal PCR efficiency, the logarithm of ratios of two products can be approximated directly from the observed data by computing the difference in Cq numbers between the numerator and denominator products. In the left panel, each dot is a biological replicate and represents the difference between the mean Cq numbers of corresponding technical replicates. Since larger Cq implies lower expression, the larger the difference Cq (denominator) - Cq (numerator) the higher proportion of the numerator product was observed.

Using the statistical model, we can further compute the ratios of those ratios between genotypes in a hypothetical "noise free" experiment, after accounting for various sources of biological and technical variation. In the right panel we show the posterior credible intervals (95% - thin, 50% - thick) and means (points) for the ratio of ratios. Here ratio of ratios > 1 means that the d17 mice have larger numerator/denominator ratio than the wt mice. (note the log scale on the vertical axis)

```{r}
plot_detailed_comparisons(comparisons_d17_4w, comparisons_pred_d17_4w, ratios_observed_d17_4w)
```


```{r, fig.cap = "Ratio of product ratios between genotypes derived from the linear mixed model. Ratio of each product (horizontal axis) to a corresponding canonical product is computed for each genotype and the ratio of those ratios is treated as the estimand of interest (vertical axis, log scale). We show posterior credible intervals (95% - thin, 50% - thick) and means (points). Here ratio of ratios > 1 means that the mutant mice have larger proportion of circular/alternative products than the wt mice. The model assumes equal transcription efficiency of both product and accounts for various sources of biological (animal, litter, sex) and technical (replicate, run) variation."}
label_variant_ord_fun <- function(x) {
  unique(-10 * grepl("^(Chd7|Rbfox1|Zfp521)", x) + grepl("^Gapdh", x))
}

plot_all_d17_4w <- comparisons_pred_d17_4w %>%
  filter(variant == 1 | label == "Rims2 circ 595",
         !grepl("^Cnksr3", label)) %>%
  compute_label_variant() %>%
  mutate(label_variant = fct_reorder(label_variant, .x = label_variant, 
                                     .fun = label_variant_ord_fun)) %>%
  plot_all_comparisons()
  
plot_all_d17_4w
ggsave(here::here("local_temp_data/comparisons_d17_4w.svg"), plot = plot_all_d17_4w)

```


Computing the ratio of ratios requires us to guess a value for PCR efficiency. In the above plots we assumed both products have the same efficiency of 1.8. To check robustness, we take two additional cases - either the numerator or the denominator product have substantially larger efficiency - 1.5 vs 2.0. Many of the comparisons stay almost unaffected, but some do change the sign of the difference assuming different efficiency.

```{r, fig.width=8, fig.height = 8}
plot_comparisons_sensitivity(comparisons_pred_d17_4w)
```

The figure is summarized also in the following table - `low` and `high` are the bounds of the 95% posterior credible intervals, `sign` is the main direction of the effect and `robust` indicates whether the sign would change assuming different efficiency.

```{r}
summ_d17_4w <- summarise_comparisons(comparisons_pred_d17_4w)
summ_d17_4w %>% table_format()
```

To highlight - the qualitative conclusions for the following comparisons are not robust to big differences in amplification efficiency.

```{r}
summ_d17_4w %>% 
  filter(!robust) %>%
  group_by(genotypes) %>%
  summarise(non_robust_comparisons = paste0(label_variant, collapse = ", ")) %>% table_format()
```



## Comparisons - linear forms

```{r}

lin_label_variant_ord_fun <- function(x) {
  res <- grepl("^Ntm", as.character(unique(x)))
  as.numeric(res)
}



comparisons_lin_raw <- tibble::tribble(
  ~label,           ~denominator, ~numerator_1, ~numerator_2,
  "Chd7/Gapdh",    "Gapdh",   "MR523/524", "MR526/527",
  "Rbfox1/Gapdh",  "Gapdh",   "MR513/514", "MR515/516", 
  "Zfp521/Gapdh",  "Gapdh",   "MR518/519", "MR520/521",
  "Ntm circ./Gapdh",     "Gapdh",   "MR494/495", NA,
  "Ntm lin./Gapdh","Gapdh",   "MR492/493", NA,
) 
  
  

comparisons_lin_d17_4w <- process_comparisons(comparisons_lin_raw, d17_4w_model)

comparisons_lin_d17_4w_eff <-  comparisons_lin_d17_4w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())

comparisons_lin_d17_4w %>% select(label, variant, numerator, denominator) %>% table_format()

```
```{r}
comparisons_lin_pred_d17_4w <- predict_comparisons(fit_d17_4w, d17_4w, comparisons_lin_d17_4w_eff, c("wt", "wt"), c( "d17", "het")) 

comparisons_lin_pred_d17_4w$label_variant <- fct_reorder(.f = comparisons_lin_pred_d17_4w$label_variant, .x = comparisons_lin_pred_d17_4w$label_variant, 
                                     .fun = lin_label_variant_ord_fun)

ratios_observed_d17_4w <- compute_ratios_observed(d17_4w_model, comparisons_lin_d17_4w)
```


```{r}
plot_detailed_comparisons(comparisons_lin_d17_4w, comparisons_lin_pred_d17_4w, ratios_observed_d17_4w, n_to_show = 4)
```


```{r, fig.cap = "Ratio of product ratios between genotypes derived from the linear mixed model. Ratio of each product to Gapdh (horizontal axis) is computed for each genotype and the ratio of those ratios is treated as the estimand of interest (vertical axis, log scale). We show posterior credible intervals (95% - thin, 50% - thick) and means (points). Here ratio of ratios > 1 means that the mutant mice have larger proportion of the product than the wt mice. The model assumes equal transcription efficiency of both products and accounts for various sources of biological (animal, litter, sex) and technical (replicate, run) variation."}
plot_all_comparisons(comparisons_lin_pred_d17_4w)
```


```{r}
plot_comparisons_sensitivity(comparisons_lin_pred_d17_4w, ncol = 2)
```

```{r}
summ_lin_d17_4w <- summarise_comparisons(comparisons_lin_pred_d17_4w)
summ_lin_d17_4w

to_include <- summ_lin_d17_4w$label_variant != "Ntm circ./Gapdh"

min_change <- exp(-min(summ_lin_d17_4w$low[to_include]))
max_change <- exp(max(summ_lin_d17_4w$high[to_include]))
```
For all of the linear comparisons (except Ntm circ.), we can rule out more then `r min_change` fold decrease
and more than `r max_change` fold increase in expression relative to Gapdh.

# d17 8w


```{r}
f_d17_8w <- here::here("private_data/qPCR/qPCR Prpf8 d17 8w final.xlsx")
d17_8w <- rbind(
  read_qpcr(f_d17_8w, "qPCR set 1", range = "A5:BA22"),
  read_qpcr(f_d17_8w, "qPCR set 2", range = "A5:BA22")
)
```


```{r}
genes_boxplot(d17_8w, housekeeping)
```

```{r}
genes_boxplot(d17_8w, sort(setdiff(unique(d17_8w$primer_short), housekeeping)))
```

## Fitting the model

Simpler model as we have very few runs includ . Including sex also makes the model diverge, as there are not enough data to inform all the parameters.

```{r}
d17_8w_model <- data_for_model(d17_8w, housekeeping) 

fit_d17_8w <-   brm(
    bf(Cq ~ primer_short * genotype + (1 || animal_no + animal_no:primer_short + litter:primer_short),
       sigma ~ (1 || primer_short)
    ),
    data = d17_8w_model,
    adapt_delta = 0.95,
    prior = priors, refresh = 500,
    file = here::here("local_temp_data/fit_d17_8w"), file_refit = "on_change")
  
fit_d17_8w
```

However, looking at the PP checks for sex we see that the variability between sexes is already accounted for by the model as it includes a per-animal and primer varying intercept.

```{r}
if(pp_checks) {
  show_pp_checks(d17_8w_model, fit_d17_8w)
}
```


## Comparisons - alternative forms


```{r}
comparisons_d17_8w <- comparisons_raw %>% 
  process_comparisons(data_model = d17_8w_model)
  

comparisons_d17_8w_eff <- comparisons_d17_8w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())
```


```{r}
comparisons_pred_d17_8w <- predict_comparisons(fit_d17_8w, d17_8w, comparisons_d17_8w_eff, c( "wt"), c("d17"))
ratios_observed_d17_8w <- compute_ratios_observed(d17_8w_model, comparisons_d17_8w)
```




```{r}
plot_detailed_comparisons(comparisons_d17_8w, comparisons_pred_d17_8w, ratios_observed_d17_8w)
```


```{r}
plot_all_d17_8w <-  comparisons_pred_d17_8w %>%     
  filter(variant == 1, !grepl("^Cnksr3", label)) %>%
  compute_label_variant() %>%
  mutate(label_variant = fct_reorder(label_variant, .x = label_variant, 
                                     .fun = label_variant_ord_fun)) %>%
  plot_all_comparisons()

plot_all_d17_8w
ggsave(here::here("local_temp_data/comparisons_d17_8w.svg"), plot = plot_all_d17_8w)

```

```{r, fig.width=8, fig.height = 8}
plot_comparisons_sensitivity(comparisons_pred_d17_8w)
```

```{r}
summ_d17_8w <- summarise_comparisons(comparisons_pred_d17_8w)
summ_d17_8w %>% table_format()
```

To highlight - the qualitative conclusions for the following comparisons are not robust to big differences in amplification efficiency.

```{r}
summ_d17_8w %>% 
  filter(!robust) %>%
  group_by(genotypes) %>%
  summarise(non_robust_comparisons = paste0(label_variant, collapse = ", ")) %>% table_format()
```



## Comparisons - linear forms

```{r}
comparisons_lin_d17_8w <- process_comparisons(comparisons_lin_raw, d17_8w_model)

comparisons_lin_d17_8w_eff <-  comparisons_lin_d17_8w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())
```


```{r}
comparisons_lin_pred_d17_8w <- predict_comparisons(fit_d17_8w, d17_8w, comparisons_lin_d17_8w_eff, c("wt"), c("d17"))

comparisons_lin_pred_d17_8w$label_variant <- fct_reorder(.f = comparisons_lin_pred_d17_8w$label_variant, .x = comparisons_lin_pred_d17_8w$label_variant, 
                                     .fun = lin_label_variant_ord_fun)


ratios_observed_d17_8w <- compute_ratios_observed(d17_8w_model, comparisons_lin_d17_8w)
```


```{r}
plot_detailed_comparisons(comparisons_lin_d17_8w, comparisons_lin_pred_d17_8w, ratios_observed_d17_8w, n_to_show = 4)
```


```{r}


plot_all_lin_d17 <- rbind(comparisons_lin_pred_d17_4w %>% mutate(time = "4 weeks"),
      comparisons_lin_pred_d17_8w %>% mutate(time = "8 weeks")) %>% 
  plot_all_comparisons(expand = c(0.1,2), breaks = c(0.1,0.25,0.5,1,2)) + facet_wrap(~time)

plot_all_lin_d17

ggsave(here::here("local_temp_data/comparisons_lin_d17.svg"), plot = plot_all_lin_d17)
```


```{r}
plot_comparisons_sensitivity(comparisons_lin_pred_d17_8w, ncol = 2)
```

```{r}
summ_lin_d17_8w <- summarise_comparisons(comparisons_lin_pred_d17_8w)
summ_lin_d17_8w

min_decrease <- exp(-max(summ_lin_d17_8w$high[summ_lin_d17_8w$label_variant != "Ntm lin./Gapdh"]))
```

All of the linear comparisons except for Ntm linear point downards and for those, we can rule out less then `r min_decrease` fold decrease in expression relative to Gapdh.

# N 4w


```{r}
f_N_4w <- here::here("private_data/qPCR/qPCR Prpf8 N 4w final.xlsx")
N_4w <- rbind(
  read_qpcr(f_N_4w, "qPCR set a1", range = "A5:AR31", mutation = "N"),
  read_qpcr(f_N_4w, "qPCR set a2", range = "A5:Z22", mutation = "N"),
  read_qpcr(f_N_4w, "qPCR set a3", range = "A5:AU22", mutation = "N"),
  read_qpcr(f_N_4w, "qPCR set a4", range = "A6:T31", mutation = "N")
)
```


```{r}
genes_boxplot(N_4w, housekeeping)
```
```{r}
genes_boxplot(N_4w, sort(setdiff(unique(N_4w$primer_short), housekeeping)))
```

## Fitting the model

```{r}
N_4w_model <- data_for_model(N_4w, housekeeping)

fit_N_4w <-   brm(
    bf(Cq ~ primer_short * genotype + (1 || animal_no + animal_no:primer_short + litter:primer_short),
       sigma ~ (1 || primer_short)
    ),
    data = N_4w_model,
    adapt_delta = 0.95,
    max_treedepth = 11,
    prior = priors, refresh = 500,
    file = here::here("local_temp_data/fit_N_4w"), file_refit = "on_change")
  
fit_N_4w
```

```{r}
if(pp_checks) {
  show_pp_checks(N_4w_model, fit_N_4w)
}
```

## Comparisons - alternative forms


```{r}
comparisons_N_4w <- comparisons_raw %>% 
  process_comparisons(data_model = N_4w_model)
  

comparisons_N_4w_eff <- comparisons_N_4w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())
```


```{r}
comparisons_pred_N_4w <- predict_comparisons(fit_N_4w, N_4w, comparisons_N_4w_eff, c( "wt", "wt"), c("N", "het"))
ratios_observed_N_4w <- compute_ratios_observed(N_4w_model, comparisons_N_4w)
```




```{r}
plot_detailed_comparisons(comparisons_N_4w, comparisons_pred_N_4w, ratios_observed_N_4w)
```


```{r}
plot_all_N_4w <-  comparisons_pred_N_4w   %>%
  filter(variant == 1, !grepl("^Cnksr3", label)) %>%
  compute_label_variant() %>%
  mutate(label_variant = fct_reorder(label_variant, .x = label_variant, 
                                     .fun = label_variant_ord_fun)) %>%
  plot_all_comparisons()

plot_all_N_4w
ggsave(here::here("local_temp_data/comparisons_N_4w.svg"), plot = plot_all_N_4w)

```

```{r, fig.width=8, fig.height = 8}
plot_comparisons_sensitivity(comparisons_pred_N_4w)
```

```{r}
summ_N_4w <- summarise_comparisons(comparisons_pred_N_4w)
summ_N_4w %>% table_format()
```

To highlight - the qualitative conclusions for the following comparisons are not robust to big differences in amplification efficiency.

```{r}
summ_N_4w %>% 
  filter(!robust) %>%
  group_by(genotypes) %>%
  summarise(non_robust_comparisons = paste0(label_variant, collapse = ", ")) %>% table_format()
```



## Comparisons - linear forms

```{r}
comparisons_lin_N_4w <- process_comparisons(comparisons_lin_raw, N_4w_model)

comparisons_lin_N_4w_eff <-  comparisons_lin_N_4w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())
```


```{r}
comparisons_lin_pred_N_4w <- predict_comparisons(fit_N_4w, N_4w, comparisons_lin_N_4w_eff, c("wt","wt"), c("N","het"))

comparisons_lin_pred_N_4w$label_variant <- fct_reorder(.f = comparisons_lin_pred_N_4w$label_variant, .x = comparisons_lin_pred_N_4w$label_variant, 
                                     .fun = lin_label_variant_ord_fun)


ratios_observed_N_4w <- compute_ratios_observed(N_4w_model, comparisons_lin_N_4w)
```


```{r}
plot_detailed_comparisons(comparisons_lin_N_4w, comparisons_lin_pred_N_4w, ratios_observed_N_4w, n_to_show = 4)
```


```{r}
plot_all_comparisons(comparisons_lin_pred_N_4w)
```


```{r}
plot_comparisons_sensitivity(comparisons_lin_pred_N_4w, ncol = 2)
```

```{r}
summ_lin_N_4w <- summarise_comparisons(comparisons_lin_pred_N_4w)
summ_lin_N_4w

to_include <- summ_lin_N_4w$label_variant != "Ntm circ./Gapdh"


min_change <- exp(-min(summ_lin_N_4w$low[to_include]))
max_change <- exp(max(summ_lin_N_4w$high[to_include]))
```

For all of the linear comparisons (not including Ntm circ), we can rule out more then `r min_change` fold decrease
and more than `r max_change` fold increase in expression relative to Gapdh.


# N 12w


```{r}
f_N_12w <- here::here("private_data/qPCR/qPCR Prpf8 N 12w final.xlsx")
N_12w <- rbind(
  read_qpcr(f_N_12w, "qPCR set 1", range = "A5:BA31", mutation = "N"),
  read_qpcr(f_N_12w, "qPCR set 2", range = "A5:BA31", mutation = "N")
)
```


```{r}
genes_boxplot(N_12w, housekeeping)
```
```{r}
genes_parcoord(N_12w, housekeeping)
```


```{r}
genes_boxplot(N_12w, sort(setdiff(unique(N_12w$primer_short), housekeeping)))
```

## Fitting the model

```{r}
N_12w_model <- data_for_model(N_12w, housekeeping)

fit_N_12w <-   brm(
    bf(Cq ~ primer_short * genotype + (1 || animal_no + animal_no:primer_short),
       sigma ~ (1 || primer_short)
       ),
    data = N_12w_model,
    control = list(max_treedepth = 11),
    prior = priors, refresh = 500,
    file = here::here("local_temp_data/fit_N_12w"), file_refit = "on_change")
  
fit_N_12w
``` 

```{r}
if(pp_checks) {
  show_pp_checks(N_12w_model, fit_N_12w)
}
```

## Comparisons - alternative forms


```{r}
comparisons_N_12w <- comparisons_raw %>% 
  process_comparisons(data_model = N_12w_model)
  

comparisons_N_12w_eff <- comparisons_N_12w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())
```


```{r}
comparisons_pred_N_12w <- predict_comparisons(fit_N_12w, N_12w, comparisons_N_12w_eff, c( "wt", "wt"), c("N", "het"))
ratios_observed_N_12w <- compute_ratios_observed(N_12w_model, comparisons_N_12w)
```




```{r}
plot_detailed_comparisons(comparisons_N_12w, comparisons_pred_N_12w, ratios_observed_N_12w)
```


```{r}
plot_all_N_12w <-  comparisons_pred_N_12w %>% 
  filter(variant == 1, !grepl("^Cnksr3", label)) %>%
  compute_label_variant() %>%
  mutate(label_variant = fct_reorder(label_variant, .x = label_variant, 
                                     .fun = label_variant_ord_fun)) %>%
  plot_all_comparisons()

plot_all_N_12w
ggsave(here::here("local_temp_data/comparisons_N_12w.svg"), plot = plot_all_N_12w)

```

```{r, fig.width=8, fig.height = 8}
plot_comparisons_sensitivity(comparisons_pred_N_12w)
```

```{r}
summ_N_12w <- summarise_comparisons(comparisons_pred_N_12w)
summ_N_12w %>% table_format()
```

To highlight - the qualitative conclusions for the following comparisons are not robust to big differences in amplification efficiency.

```{r}
summ_N_12w %>% 
  filter(!robust) %>%
  group_by(genotypes) %>%
  summarise(non_robust_comparisons = paste0(label_variant, collapse = ", ")) %>% table_format()
```



## Comparisons - linear forms

```{r}
comparisons_lin_N_12w <- process_comparisons(comparisons_lin_raw, N_12w_model)

comparisons_lin_N_12w_eff <-  comparisons_lin_N_12w %>%
  crossing(eff_variants) %>%
  mutate(comparison_id = 1:n())
```


```{r}
comparisons_lin_pred_N_12w <- predict_comparisons(fit_N_12w, N_12w, comparisons_lin_N_12w_eff, c("wt", "wt"), c("N", "het"))

comparisons_lin_pred_N_12w$label_variant <- fct_reorder(.f = comparisons_lin_pred_N_12w$label_variant, .x = comparisons_lin_pred_N_12w$label_variant, 
                                     .fun = lin_label_variant_ord_fun)


ratios_observed_N_12w <- compute_ratios_observed(N_12w_model, comparisons_lin_N_12w)
```


```{r}
plot_detailed_comparisons(comparisons_lin_N_12w, comparisons_lin_pred_N_12w, ratios_observed_N_12w, n_to_show = 4)
```


```{r}

plot_all_lin_N <- rbind(comparisons_lin_pred_N_4w %>% mutate(time = "4 weeks"),
      comparisons_lin_pred_N_12w %>% mutate(time = "12 weeks")) %>% 
  mutate(time = factor(time, levels = c("4 weeks", "12 weeks"))) %>%
  plot_all_comparisons(expand = c(0.1,2), breaks = c(0.1,0.25,0.5,1,2)) + facet_wrap(~time)

plot_all_lin_N

ggsave(here::here("local_temp_data/comparisons_lin_N.svg"), plot = plot_all_lin_N)
```


```{r}
plot_comparisons_sensitivity(comparisons_lin_pred_N_12w, ncol = 2)
```

```{r}
summ_lin_N_12w <- summarise_comparisons(comparisons_lin_pred_N_12w)
summ_lin_N_12w
```

# Mathematical model


Steibel et al.: Mixed models, normal error https://www.sciencedirect.com/science/article/pii/S0888754309000986
The same, but Poisson-LogNormal https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0071448



Below, assuming we have "noise free" $Cq$ values.

The ratio of the products $a,b$ in sample $i$ is:

$$
v_i = \frac{E_a^{S_a - Cq_{a,i}}}{E_b^{S_b - Cq_{b,i}}}
$$

Where $E_a, E_b$ are the PCR efficiencies and $S_a, S_b$ - the cycle that would be reached when amplifying a single molecule (similarly an abstract "noise free" value).

We are interested in 

$$
r = \frac{v_1}{v_2} \\
\log{r} = \log v_1 - \log v_2 = \left(  (\log E_a)(S_a - Cq_{a,1}) - (\log E_b)(S_b - Cq_{b,1}) \right) - 
\left(  (\log E_a)(S_a - Cq_{a,2}) - (\log E_b)(S_b - Cq_{b,2}) \right) = \\
(\log E_a)(Cq_{a,2} - Cq_{a,1}) - (\log E_b)(Cq_{b,2}-Cq_{b,1})
$$


$$
\log{r} = \log v_1 - \log v_2 = \left(  (\log E_a)(S_a - Cq_{a,1}) - (\log E_b)(S_b - Cq_{b,1}) \right) - \left(  (\log E_a)(S_a - Cq_{a,2}) - (\log E_b)(S_b - Cq_{b,2}) \right) = 
(\log E_a)(Cq_{a,2} - Cq_{a,1}) - (\log E_b)(Cq_{b,2}-Cq_{b,1})
$$