---
title: "SBIR project - revisions"
author: "Jeremy Li"
date: "9/23/2020"
output: html_document
editor_options: 
  chunk_output_type: inline
header-includes:
  - \usepackage{eqnarray}
---

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

# Requested revisions: results using different imputation algorithms 

We imputed chromosome 21 using GLIMPSE for sequence data and IMPUTE5 for array data.
We recompute concordance on chromosome 21 from the original data and the new results. 

```{r message = FALSE}
# load in requisite packages
library(tidyverse)
library(ggthemes)
library(patchwork)
library(rms)
library(glue)
library(knitr)
library(kableExtra)

# populations
pops = read_delim("../data_revisions/all_samples_with_pops.list", delim = " ")

# load in precomputed data
concordance = read_delim("../data_revisions/all-gc-counts.txt",
                          col_names = c("cell_line", "sample", "stat_type", "set_id", "af_midpoint", "rr_hom_matches", "ra_het_matches", "aa_hom_matches",
                   "rr_hom_mismatches", "ra_het_mismatches", "aa_hom_mismatches", "dosage_rsquared", "n_genotypes", "site_type", "tech", "imputation_program"), 
                   delim = " ") %>% 
   mutate(
    # characterize into ilmn 0.5x, 1x, BGI, array
    Experiment = case_when(
      grepl("-0-2$", sample) ~ "Experiment D",
      grepl("-0-[0-9]-[01]-0$", sample) ~ "Experiment A",
      grepl("-1-[0-9]-[01]-0$", sample) ~ "Experiment B",
      grepl("_[1-3]$", sample) ~ "Experiment E (array)",
      grepl("-1$", sample) ~ "ilmn_1x_repl",
    )
  ) %>% 
  inner_join(pops, by = "cell_line") %>% 
  mutate(
    nrc_numerator = (ra_het_matches + aa_hom_matches),
    nrc_denominator =  (ra_het_matches + aa_hom_matches + rr_hom_mismatches + ra_het_mismatches + aa_hom_mismatches),
    nrc = nrc_numerator / nrc_denominator
  )

```

We can summarize results and plot, like so:

```{r, fig.width = 6, fig.height = 4}
concordance_summary = concordance %>%  
  group_by(af_midpoint, Experiment, stat_type, site_type, tech, imputation_program, super_pop) %>%
  summarise(mean_nrc = mean(nrc)) %>% 
  rename(
    # rename to human-readable
    `Super population` = super_pop,
  )

concordance_summary %>% 
  filter(stat_type == "GCsAF", site_type == "allsites") %>% 
  ggplot(aes(x = af_midpoint, y = mean_nrc, color = Experiment)) + 
  facet_grid(~`Super population`) + geom_point() + 
  geom_line(aes(linetype = factor(imputation_program))) + scale_color_tableau() + theme_few() + 
  labs(
    title = "Average NRC by non-reference allele frequency, unfiltered on chr21",
    x = "Non-Reference Allele Frequency",
    y = "Non-Reference Concordance"
  ) + scale_x_log10()

ggsave("../paper/src/figs/revisions-chr21-nrc.pdf", width = 10, height = 7)
```

from which we see that the qualitative trend of the sequence imputation being more accurate than the array imputation holding, with a fixed effect of the imputation program being used given an assay type/Experiment. 

Interestingly, if one looks at filtered SNPs, the gap between concordance for the two array imputation algorithms (IMPUTE5 vs minimac4) widens significantly:

```{r, fig.width = 6, fig.height = 4}
concordance_summary %>% 
  filter(stat_type == "GCsAF", site_type == "passing") %>% 
  ggplot(aes(x = af_midpoint, y = mean_nrc, color = Experiment)) + 
  facet_grid(~`Super population`) + geom_point() + 
  geom_line(aes(linetype = factor(imputation_program))) + scale_color_tableau() + theme_few() + 
  labs(
    title = "Average NRC by non-reference allele frequency, filtered for high confidence on chr21",
    x = "Non-Reference Allele Frequency",
    y = "Non-Reference Concordance"
  ) + scale_x_log10()

```

particularly at higher allele frequencies. 

# IQRs for imputation r2s

It is informative to include IQRs for the imputation r2s in addition to the means. 

In particular we want to find the IQRs for the common variants by cohort. 
We must compute this from the site-level r2 data:
```{r}

# site-level r2 data
unbinned_r2s =  read_delim(
  "../data/dos_r2_biallelic_snps.gz",
  delim = " ",
  col_names = c("cohort", "dos_r2", "af")) %>% 
  mutate(
    maf = if_else(af > 0.5, 1 - af, af)
  )

unbinned_r2s
```

from which we can compute 

```{r}
r2_means = unbinned_r2s %>% 
  filter(maf > 0.05) %>% 
  group_by(cohort) %>% 
  summarise(
    mean_r2 = mean(dos_r2), 
    r2_quantile_25 = quantile(dos_r2, 1/4),
    r2_quantile_75 = quantile(dos_r2, 3/4),
  )  %>% 
  mutate(
    # display 
    mean_r2_display =  paste(round(mean_r2, 4), " (", round(r2_quantile_25, 4), "-", round(r2_quantile_75, 4), ")", sep = ""),
    `Super population` = if_else(grepl("afr", cohort), "AFR", "EUR"),
    Experiment = case_when(
      cohort == 'afr-0-ilmn' ~ "A",
      cohort == "afr-1-bgi" ~ "D",
      cohort == "afr-1-ilmn" ~ "B",
      cohort == "afr-array-ilmn2" ~ "E",
      cohort == "eur-0-ilmn" ~ "A",
      cohort == "eur-1-bgi" ~ "D",
      cohort == "eur-1-ilmn" ~ "B",
      cohort == "eur-array-ilmn2" ~ "E"
    )
    )

r2_means
```

and we can write this out 

```{r}
r2_means %>%
  select(`Super population`, Experiment, mean_r2_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_r2_display) %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  add_header_above(c(" ", "Experiment" = 4)) %>% 
  cat(file = "../paper/src/tabs/mean_r2s_common.tex")
```



And if we want to bin, 

```{r}
# binning 
seq_breaks = c(0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.15, 0.2, 0.25, 0.5)
seq_break_labels = c(0.005,0.015,0.025,0.035,0.045,0.055,0.065,0.075,0.085,0.095,0.125,0.175,0.225,0.375)

bin_r2s = function(df, breaks, r2_column, af_column, grouping_col){
  #' bin_r2s
  #' 
  #' Take a DF with an allele frequency column and an r2 column and bin according to the breaks
  #' 
  #' @param df dataframe w/ the requisite data & columns
  #' @param breaks list of breakpoints
  #' @param r2_column name of the column w/ the r2 data to break
  #' @param af_column name of the allele frequency on which to break
  #' @param grouping_col name of col on which to group on in addition to the allele frequency bin
  
  # add column denoting the bin before aggregating
  binned_df = df %>% 
    mutate(
      af_bin = cut(
        x = !!af_column,
        breaks = breaks,
        include.lowest = TRUE,
        right = TRUE, # interval closed on right, open on left
        # labels = TRUE
      )
    ) %>% 
    group_by(
      !!grouping_col, af_bin
    ) %>% 
    summarise(
      mean_r2 = mean(!!r2_column, na.rm = TRUE),
      n_snps = n(),
      
      # display 
     r2_quantile_25 = quantile(!!r2_column, 1/4),
     r2_quantile_75 = quantile(!!r2_column, 3/4),
     
     mean_r2_display =  paste(round(mean_r2, 4), " (", round(r2_quantile_25, 4), "-", round(r2_quantile_75, 4), ")", sep = "")
    ) %>% 
    ungroup() %>% 
    group_by(!!grouping_col) %>% 
  mutate(
    af_bin_index = row_number(),
    af_midpoint = seq_break_labels[af_bin_index]
  ) %>% 
    ungroup()
  
  return(binned_df)
}


binned_r2s =  bin_r2s(unbinned_r2s, 
                     breaks = seq_breaks,
                     r2_column = quo(dos_r2),
                     af_column = quo(maf),
                     grouping_col = quo(cohort)
                     )

binned_r2s = binned_r2s %>%
  mutate(
    super_pop = if_else(grepl("afr", cohort), "AFR", "EUR"),
    assay_type = case_when(
      cohort == 'afr-0-ilmn' ~ "ilmn_0.5x",
      cohort == "afr-1-bgi" ~ "bgi_1x",
      cohort == "afr-1-ilmn" ~ "ilmn_1x",
      cohort == "afr-array-ilmn2" ~ "gsa",
      cohort == "eur-0-ilmn" ~ "ilmn_0.5x",
      cohort == "eur-1-bgi" ~ "bgi_1x",
      cohort == "eur-1-ilmn" ~ "ilmn_1x",
      cohort == "eur-array-ilmn2" ~ "gsa"
    ), 
    
    # for display 
    mean_r2_display = paste(round(mean_r2, 4), " (", round(r2_quantile_25, 4), "-", round(r2_quantile_75, 4), ")", sep = "")
  )

binned_r2s
```

# Tests for significance 

We can test for difference in means 

```{r}
unbinned_r2s = unbinned_r2s %>% 
  mutate(
    super_pop = if_else(grepl("afr", cohort), "AFR", "EUR"),
    assay_type = case_when(
      cohort == 'afr-0-ilmn' ~ "ilmn_0.5x",
      cohort == "afr-1-bgi" ~ "bgi_1x",
      cohort == "afr-1-ilmn" ~ "ilmn_1x",
      cohort == "afr-array-ilmn2" ~ "gsa",
      cohort == "eur-0-ilmn" ~ "ilmn_0.5x",
      cohort == "eur-1-bgi" ~ "bgi_1x",
      cohort == "eur-1-ilmn" ~ "ilmn_1x",
      cohort == "eur-array-ilmn2" ~ "gsa"
    )
    )

# results for arrays 
afr_array_r2s = unbinned_r2s %>% filter(maf > 0.05) %>% filter(super_pop == "AFR", assay_type == "gsa") %>% pull(dos_r2)
eur_array_r2s = unbinned_r2s %>% filter(maf > 0.05) %>% filter(super_pop == "EUR", assay_type == "gsa") %>% pull(dos_r2)

print("Two-sample paired t-tests for difference in means between sequence and array r2s")
for (i in c("ilmn_0.5x", "ilmn_1x", "bgi_1x")) {
    eurr2s = unbinned_r2s %>% filter(super_pop == "EUR", assay_type == i) %>% pull(dos_r2)
    print(paste("p value for eur", i, "vs array", t.test(eurr2s, eur_array_r2s)$p.value))

    afrr2s = unbinned_r2s %>% filter(super_pop == "AFR", assay_type == i) %>% pull(dos_r2)
    print(paste("p value for afr", i, "vs array", t.test(afrr2s, afr_array_r2s)$p.value))

}
```

these have absurdly large sample sizes (millions), so of course the p-values will be tiny. 




