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

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

## Results

Some results from the SBIR project. We first load in data which we'll proceed to analyze and visualize. 

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

# load in precomputed data
scores = read_csv("../data/scores.csv")
coverage = read_csv("../data/coverage.csv")
concordance = read_csv("../data/concordance.csv")
binned_r2s = read_csv("../data/binned_r2s.csv")
random_samples = read_csv("../data/all.random.samples", col_names = c("sample"))
bgi_cl = read_csv("../data/bgi-cell-lines.csv", col_names = c("cell_line"))
```

Recall that we sequenced 120 individuals (60 AFR and 60 EUR) samples in triplicate at 0.5x and 1x on Illumina machines and a subset of 30 each from those populations on BGI machines at 1x (with no replicates). We also did NA12878 30 times on Illumina machines at a target coverage of 1x. 

There's a couple views we can take of the data --- first, the most obvious one is just look at every single sample that we sequenced. Another view that we want to take a look at is to, for each cell line, take a random sample of one individual from the replicates for each "assay type" (of which there are 4: Illumina 1x, Illumina 0.5x, Illumina GSA, BGI 1x) and then compare the same stats for the same cell lines across assay types. 

We treated the 1000 Genomes Phase 3 (1KGP3) release on GRCh37 as the gold standard "truth" against which we computed accuracy metrics --- for imputation for all assay types we performed imputation in a leave-one-out manner. 

We will proceed by giving a quick outline of the data that we have, before proceeding to an outline of "effective coverage" and presenting cohort-level results, after which we will discard replicates and look at a single iteration of a cell line for each assay type and examine apples-to-apples results. 

## Data 
We have coverage data on our samples:
```{r}
head(coverage)
```

We have concordance information at the allele frequency bin level broken down by a number of things including

* imputation quality (whether a genotype call was considered `PASS`ing) 
* variant type (SNPs or indels)
* a given non-reference allele frequency bin, where the cuts are at breakpoints 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, 1 according to the allele frequency in 1KGP3. 
* types of concordance (e.g. details on whether it's RR or RA sites in the gold standard that were called wrong, etc)
* population
* assay type
* and more... 

```{r}
glimpse(concordance)
```

We also have polygenic risk scores for breast cancer (brca) and coronary artery disease (CAD) which were calculated off the imputed dosages for each sample: 

```{r}
head(scores)
```

We have a list of a randomly selected samples, one for each assay type so that we can compare head-to-head stats across assay types:

```{r}
head(random_samples)
```

as well as the imputation $r^2$s as calculated for each assay type between those individuals and the gold standard data, binned by minor allele frequency:

```{r}
head(binned_r2s)
```

We also have a list of the cell lines that BGI ran (recall 58 out of 60 of them succeeded, so it's useful to have a note of which didn't).

```{r}
head(bgi_cl)
```


Let's start out from the bottom; _i.e._, details of what we got out of the sequencer. 

## Samples

As mentioned previously, we have 60 AFR and 60 EUR individuals in this dataset: specifically, 

```{r}
concordance %>% distinct(cell_line, super_pop, pop)
```

we can group these by population and superpopulation and print them out into a latex table for the manuscript
```{r}
concordance %>%
  distinct(cell_line, super_pop, pop) %>% 
  group_by(pop, super_pop) %>% 
  tally() %>% 
  arrange(super_pop, -n) %>% 
  rename( #rename for human readable
    Population = pop,
    `Super Population` = super_pop,
  ) %>% 
  kable("latex", booktabs = TRUE, linesep = "") %>% 
  cat(file = "../paper/src/tabs/sample-pops.tex")
```

We also write out the number of samples for each cell line for each experiment which passed QC. 

```{r}
# count how many samples for each experiment passed QC 
samples_passing_qc = concordance %>% 
  distinct(sample, cell_line, pop, super_pop, assay_type) %>% 
  group_by(cell_line, pop, super_pop, assay_type) %>% 
  tally() %>% 
  pivot_wider(names_from = assay_type, values_from = n) %>% 
  rename(
    # rename to human-readable
    `Cell line` = cell_line,
    Population = pop,
    `Super population` = super_pop,
    `Exp. A` = ilmn_0.5x,
    `Exp. B` = ilmn_1x,
    `Exp. C` = ilmn_1x_repl,
    `Exp. D` = bgi_1x,
    `Exp. E` = ilmn_GSA
  ) %>% 
  select(`Cell line`, Population, `Super population`, sort(current_vars())) 
  
# replace NA with more human readable. depending on the experiment either 0 it out (for those which ran but did not succeed) or replace NA with a dot, which indicates its not part of the experiment
samples_passing_qc_print = samples_passing_qc %>% 
  mutate(`Exp. D` = as.character(`Exp. D`)) %>% 
  mutate(`Exp. D` = if_else(`Cell line` %in% bgi_cl$cell_line, `Exp. D`, ".")) %>% 
  replace_na(
    list(
      `Exp. A` = 0,
      `Exp. B` = 0,
      `Exp. C` = ".",
      `Exp. D` = 0,
      `Exp. E` = 0
    )
  ) %>% 
  arrange(
    `Super population`
  )

# write out 
samples_passing_qc_print %>% 
  kable(
    "latex", 
    booktabs = TRUE, 
    longtable = TRUE, 
    caption = "This table shows the number of samples for each cell line in each experiment which passed QC. A dot indicates that the particular cell line was not part of the experiment.",
    label = "supptab:samples-passing-qc"        
  ) %>% 
  kable_styling(
    latex_options = c("repeat_header", "HOLD_position")
    ) %>% 
  cat(file = "../paper/src/tabs/samples-passing-qc.tex")
```



## Effective coverage

Even though we have target coverages to which we sequence samples, there's always variation in what you actually get back. For example, we can look at the expected number of bases sequenced for each sample and the actual number of bases sequenced for each sample. We've precomputed the overall coverage (bases sequenced / size of genome) and the same for deduped coverage, deduped mapped coverage, and more: this is stored in `coverage`:

```{r}
head(coverage)
```
which we recode to include Experiment notation:

```{r}
coverage = coverage %>% 
  mutate(
    # replace assay type with Experiment X
    Experiment = case_when(
      assay_type == "ilmn_0.5x" ~ "A",
      assay_type == "ilmn_1x" ~ "B",
      assay_type == "ilmn_1x_repl" ~ "C",
      assay_type == "bgi_1x" ~ "D"
    ),
    
    # clean up names
    `Mapped coverage` = coverage_deduped_mapped,
    `Effective coverage` = eff_cov
    )
```



and we can plot out the target vs realized coverages:

```{r fig.width=7, fig.height=10}
realized_plot = coverage %>%
  ggplot(aes(x = coverage_all, fill = assay_type)) + 
  geom_density(alpha = 0.7) + theme_few() + scale_fill_tableau() + 
  labs(title = "Overall realized coverage") 

deduped_plot = coverage %>% 
  ggplot(aes(x = coverage_deduped, fill = assay_type)) + 
  geom_density(alpha = 0.7) + theme_few() + scale_fill_tableau() + 
  labs(title = "Overall deduped coverage") 

mapped_deduped_plot = coverage %>% 
  ggplot(aes(x = coverage_deduped_mapped, fill = assay_type)) + 
  geom_density(alpha = 0.7) + theme_few() + scale_fill_tableau() + 
  labs(title = "Overall mapped deduped coverage")

realized_plot / deduped_plot / mapped_deduped_plot
```

However, the raw number of bases is not the only thing that matters; ideally, we performing shotgun sequencing where reads are randomly distributed across the whole genome. We can illustrate how closely this hews to reality by counting the number of sites which are imputed which are actually covered by _at least one read_ and seeing how close that number is to the expected value of the number of these sites covered by at least one read. 

We can obtain this by treating the number of reads $k$ on a site as a random variable with a Poisson distribution (whose PMF $f(k; \lambda) = \lambda^k e^{-\lambda} / k!$) where the target coverage is $\lambda$. Then 

$$
\begin{aligned}
P(k >0; \lambda) &= 1 -  \sum_{k=1}^{\infty}f(k; \lambda) \\
 &= 1 - f(k = 0; \lambda) \\
 &= 1 - e^{-\lambda}\
 \end{aligned}.
$$

so for example a given site on a genome sequenced at a target coverage of $\lambda = 1$ the probability that it is has at least one read on it is $P(k > 0 ; 1) = 1 - e ^ {-1} = 0.632120559$. Then assuming the $n$ sites we impute are IID we would expect to see around $n(1 - e^{-\lambda})$ sites covered w/ at least one read. From this the idea of an "effective coverage" follows naturally for a sequenced genome: _i.e._, the value of $\lambda_{\mathrm{eff}}$ such that $f_{\mathrm{covered}} = 1 - e^{-\lambda_{\mathrm{eff}}}$ holds, where $f_{\mathrm{covered}}$ is the fraction of sites actually covered with at least one read. 

More formally, for a single site, we can model whether or not it was covered by at least one read as a Bernoulli process, with its being covered w/ at least one read to be considered a "success" --- as previously derived, this would follow a Bernoulli distribution with probability of success being $p = 1 - e^{-\lambda}$. Then the number of pileup sites covered by at least one read over the whole genome (denoted $X$) can be modeled by a binomial distribution with a PMF 

$$
\begin{align}
  f(k, n, p) = \mathrm{Pr}(k; n, p) = Pr(X = k) = {n \choose k}p^k (1 - p)^{n-k} = {n \choose k}(1 - e^{-\lambda}) ( e^{-\lambda})^{n - k} = {n \choose k}(1 - e^{-\lambda})e^{\lambda(k -n)}
\end{align}
$$
where we plug in $p = 1 - e ^ {-\lambda}$. Then $\mathrm{E}[X] = np = n(1 - e^{-\lambda})$ and $\mathrm{Var}(X) = np(1 - p) = n(1-e^{-\lambda})e^{-\lambda}$. Then defining $f_{\mathrm{covered}} =  X / n$ we have $\mathrm{E}[f_{\mathrm{covered}}] = (1 - e^{-\lambda})$. 

We therefore obtained $f_{\mathrm{covered}}$ for each sample and calculated the $\lambda_{\mathrm{eff}} = -\ln(1 - f_{\mathrm{covered}})$. Plotting, we find that the effective coverage is consistently lower than the target coverage, even though our overall realized coverage is higher than the target coverage, indicating that sequencing oftentimes deviates from truly spatially uniform sampling: 

```{r echo = FALSE}
# density plots of pileup sites covered 
pileup_plot1 = coverage %>% 
  ggplot(aes(x = pileup_covered, fill = assay_type)) + 
  scale_fill_tableau() + geom_density(alpha = 0.7) + theme_few() + 
  labs(title = "Pileup sites covered by at least one read by assay type")

pileup_plot2 = coverage %>% 
  ggplot(aes(x = eff_cov, fill = assay_type)) + scale_fill_tableau() + 
  geom_density(alpha = 0.7) + theme_few() + 
  labs(title = "Effective coverage by assay type")

pileup_plot1 / pileup_plot2
```

Let's write out these histograms to plots as well

```{r, fig.width = 6, fig.height = 4}
# nominal mapped coverage 
coverage %>% 
  ggplot(aes(x = `Mapped coverage`, fill = Experiment)) + 
  geom_density(alpha = 0.7) + 
  scale_fill_tableau() + theme_few() + 
  labs(
    title = "Distribution of nominal mapped coverage by experiment"
  )

ggsave("../paper/src/figs/density-mapped-coverage.pdf")
```

```{r, fig.width = 6, fig.height = 4}
# effective  coverage 
coverage %>% 
  ggplot(aes(x = `Effective coverage`, fill = Experiment)) + 
  geom_density(alpha = 0.7) + 
  scale_fill_tableau() + theme_few() + 
  labs(
    title = "Distribution of effective coverage by experiment"
  )

ggsave("../paper/src/figs/density-effective-coverage.pdf")
```


Indeed, we can also plot the realized overall mapped bases against the effective coverage:
```{r echo = FALSE, fig.width = 5, fig.height = 4}
coverage %>% 
  ggplot(aes(x = `Mapped coverage`, y = `Effective coverage`, color = Experiment)) + 
  geom_point(alpha= 0.7) + theme_few() + scale_color_tableau() + 
  labs(title = "Mapped vs effective coverage") + geom_abline()

ggsave("../paper/src/figs/mapped_vs_effcov.pdf")
```
and see that that there's nontrivial variation in effective coverage; however, given a sequencing technology type (_e.g._, Illumina sequencing), effective coverage tracks overall coverage quite well: 

```{r}
# calculate pearsons correlation between the points plotted above^
coverage %>% filter(assay_type != "bgi_1x") %>% summarise(`r between overall coverage and effective coverage` = cor(coverage_all, eff_cov)) 
```

We can summarize the effective and nominal coverages by experiment: 

```{r}
# group by assay type, then get mean and std
coverage_summary = coverage %>% 
  group_by(assay_type) %>% 
  summarise(
    mean_mapped_coverage = mean(coverage_deduped_mapped),
    std_mapped_coverage = sqrt(var(coverage_deduped_mapped)),
    
    # eff cov
    mean_eff_cov = mean(eff_cov),
    std_eff_cov = sqrt(var(eff_cov)),
    
    n = n()
  )

coverage_summary
```

which summarizes the coverage trends by assay type. 
## Cohort trends

We can first examine cohort-level trends by assay and effective coverage --- we want to what extent this construct effective coverage actually trends with performance. But first, a word about non-reference concordance:

#### Non-reference Concordance (NRC)
Since by and large variants are found at low minor allele frequencies in 1KGP3, overall concordance can be a little less than informative since low MAF variants can be imputed to hom ref with high confidence as the prior is overwhelmingly against there being variation at that site. 

We therefore consider non-reference concordance, and we define it as follows: taking the variant calls at the same site between two normed VCFs, we consider all possible pairwise combinations of genotype calls (order doesn't matter):

| |0|1|2|
|-|-|-|-|
|0|a|b|c|
|1|d|e|f|
|2|g|h|i|

then we write non-reference concordance $\mathrm{NRC} = (e + i) / (b + c + d + e + f + g + h + i)$; in other words, it's the concordance at all sites where the pair of genotypes is not 0 and 0. 

We can compute this from our concordance data for each sample, taking care to group by whether the site is imputed confidently and also whether the variant is an indel or a SNP: 

```{r}
# sequence non-reference concordance by sample
seq_conc_bysample = concordance %>% 
  filter(tech == "sequence") %>% 
  group_by(
    cell_line, sample, pop, super_pop, assay_type, eff_cov, coverage_deduped_mapped, site_type, stat_type
    ) %>% 
  summarise(
    # non-reference and overall concordance
    nrc = 1 - sum(nrd_numerator) / sum(nrd_denominator),
    conc = sum(total_matches) / sum(total_sites)
    )
head(seq_conc_bysample)
```

we can then plot the broken-down sample-level NRCs against effective coverage for SNPs, stratified by population and whether to filter to passing sites. Let's first make a copy of the dataframe and rename it with more human-readable labels so we can plot it without having to deal with manually relabeling labels on plots

```{r}
# make a copy of by-sample concordance, replacing the labels to be plotted with human-friendly values
seq_conc_bysample_plot = seq_conc_bysample %>% 
  mutate(
    # replace assay type with Experiment X
    Experiment = case_when(
      assay_type == "ilmn_0.5x" ~ "A",
      assay_type == "ilmn_1x" ~ "B",
      assay_type == "ilmn_1x_repl" ~ "C",
      assay_type == "bgi_1x" ~ "D"
    ),
    # just capitalize it
    NRC = nrc,
    Concordance = conc,
    # unfiltered status 
    `Site Type` = case_when(
      site_type == "allsites" ~ "Unfiltered",
      site_type == "passing" ~ "Filtered"
    ),
    # coverage columns
    `Nominal Mapped Coverage` = coverage_deduped_mapped,
    `Effective Coverage` = eff_cov
  )
seq_conc_bysample_plot
```



```{r fig.width=7, fig.height=6}
# plot NRC by effective coverage for SNPs
seq_conc_bysample_plot %>% filter(stat_type == "GCsAF") %>% 
  ggplot(aes(x = `Effective Coverage`, y = NRC, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  theme_few() + scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) + 
  labs(title = "Non-reference concordance by effective coverage for SNPs")

ggsave("../paper/src/figs/effcov_nrc_snps.pdf")
```

and the same for indels:

```{r fig.width=7, fig.height=6}
# filtered NRC by indels 
seq_conc_bysample_plot %>% filter(stat_type == "GCiAF") %>%
  ggplot(aes(x = `Effective Coverage`, y = nrc, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  theme_few() + scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) + 
  labs(title = "Non-reference concordance by effective coverage for indels")

ggsave("../paper/src/figs/effcov_nrc_indels.pdf")
```

in both of which we observe firstly that SNPs are consistently imputed with greater accuracy with respect to NRC, imputation quality is relatively well-calibrated (_i.e._, the NRC at sites confidently imputed, "passing", is consistently higher than at unfiltered sites), AFR samples tend to be imputed with lower accuracy with respect to NRC, and finally that the NRC depends quite a bit on effective coverage. 

The same pattern holds for overall concordance as well -- for SNPs, we see
`
```{r fig.width=7, fig.height=6}
# plot overall concordance for snps
seq_conc_bysample_plot %>% 
  filter(stat_type == "GCsAF") %>% 
  ggplot(aes(x = `Effective Coverage`, y = conc, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  theme_few() + 
  scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) + 
  labs(title = "Overall concordance by effective coverage for SNPs")

ggsave("../paper/src/figs/effcov_oconc_snps.pdf")
```

and the same for indels

```{r fig.width=7, fig.height=6}
# plot overall concordance for indels
seq_conc_bysample_plot %>% 
  filter(stat_type == "GCiAF") %>% 
  ggplot(aes(x = `Effective Coverage`, y = conc, color = Experiment)) + 
  geom_point(alpha = 0.7) + theme_few() + scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) +
  labs(title = "Overall concordance by effective coverage for indels")

ggsave("../paper/src/figs/effcov_oconc_snps.pdf")
```


For *nominal coverage* as measured by deduped mapped bases, we can plot the same for SNPs for NRC:

```{r fig.width=7, fig.height=6}
# plot NRC for snps for unfiltered snps
seq_conc_bysample_plot %>% 
  filter(stat_type == "GCsAF") %>% 
  ggplot(aes(x = `Nominal Mapped Coverage`, y = nrc, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  theme_few() + 
  scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) + 
  labs(title = "NRC by nominal mapped coverage for SNPs")

ggsave("../paper/src/figs/nomcov_nrc_snps.pdf")
```

and indels

```{r fig.width=7, fig.height=6}
# plot NRC for indels for nominal  coverage
seq_conc_bysample_plot %>% 
  filter(stat_type == "GCiAF") %>% 
  ggplot(aes(x = `Nominal Mapped Coverage`, y = nrc, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  theme_few() + 
  scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) + 
  labs(title = "NRC by nominal mapped coverage for indels")

ggsave("../paper/src/figs/nomcov_nrc_indels.pdf")
```

And the same for overall concordance but by *nominal coverage*, for SNPs 

```{r fig.width=7, fig.height=6}
# plot overall concordance for snps by nominal coverage
seq_conc_bysample_plot %>% 
  filter(stat_type == "GCsAF") %>% 
  ggplot(aes(x = `Nominal Mapped Coverage`, y = conc, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  theme_few() + 
  scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) + 
  labs(title = "Overall concordance by nominal mapped coverage for SNPs")

ggsave("../paper/src/figs/nomcov_conc_snps.pdf")
```

and indels

```{r fig.width=7, fig.height=6}
# plot overall concordance for indels by nominal coverage
seq_conc_bysample_plot %>% 
  filter(stat_type == "GCiAF") %>% 
  ggplot(aes(x = `Nominal Mapped Coverage`, y = conc, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  theme_few() + 
  scale_color_tableau() + 
  facet_grid(rows = vars(`Site Type`), cols = vars(super_pop)) + 
  labs(title = "Overall concordance by nominal mapped coverage for indels")

ggsave("../paper/src/figs/nomcov_conc_indels.pdf")
```

### Fitting effective coverage  and NRC
Let's see whether nominal coverage or effective coverage better predict NRC. We use the `rms` package to fit splines with knots at percentile locations on the x axis as by [Harrell](http://biostat.mc.vanderbilt.edu/wiki/pub/Main/BioMod/notes.pdf) and evaluate the spline fits by $r^2$. 

We restrict ourselves to EUR samples sequenced on Illumina machines and consider the NRC for unfiltered SNPs. Modeling the NRC as the response variable for a $k=5$ knot cubic restricted spline with knot locations at the 5%, 27.5%, 50%, 72.5%, and 95% percentiles of the $x$-axis using Harrell's rule of thumb, we can plot the spline fits with 95% confidence intervals along the range of the predictor. 

```{r fig.height=6, fig.width=10}
# filter down to relevant
cov_predict = seq_conc_bysample %>% 
  filter(stat_type == "GCsAF" & site_type == "allsites" & super_pop == "EUR") %>% 
  mutate(
    # replace assay type with Experiment X
    Experiment = case_when(
      assay_type == "ilmn_0.5x" ~ "A",
      assay_type == "ilmn_1x" ~ "B",
      assay_type == "ilmn_1x_repl" ~ "C",
      assay_type == "bgi_1x" ~ "D"
    )
  )
  

# define distribution
ddist = datadist(cov_predict, q.display = c(0, 1))
options(datadist = "ddist")

# fit and plot
model = ols(
  nrc ~ rcs(eff_cov, quantile(eff_cov, c(0, .05, 0.275, 0.5, 0.725, 0.95, 1))), 
  data = cov_predict, 
  x = TRUE, 
  y = TRUE
  )

model_r2 = signif(model$stats[4], 2)

eff_cov_model_plot = ggplot(Predict(model)) + 
  geom_point(data = cov_predict, aes(x = eff_cov, y = nrc, color = Experiment)) + 
  theme_few() + 
  scale_color_tableau() + 
  labs(
    x = "Effective Coverage",
    y = "NRC"
    ) + 
  annotate(
    "text", x = 1, y = .88, label = glue("R ^ 2 == {model_r2}"),
    parse = TRUE
  ) + 
  theme(legend.position = "none")

# model for nominal coverage
model1 =  ols(
  nrc ~ rcs(coverage_deduped_mapped, quantile(coverage_deduped_mapped,  c(0, .05, 0.275, 0.5, 0.725, 0.95, 1))),
  data = cov_predict, 
  x = TRUE, 
  y = TRUE
  )
model1_r2 = signif(model1$stats[4], 2)

nominal_cov_model_plot = ggplot(Predict(model1)) + 
  geom_point(data = cov_predict, aes(x = coverage_deduped_mapped, y = nrc, color = Experiment)) + 
  theme_few() + scale_color_tableau() + 
  labs(
    x = "Nominal Mapped Coverage",
    y = "NRC"
    ) + 
  annotate(
    "text", x = 1.5, y = .88, label = glue("R ^ 2 == {model1_r2}"),
    parse = TRUE
  )

# put these together
patched = (eff_cov_model_plot | nominal_cov_model_plot)
patched + plot_annotation(
  title = "NRC by effective coverage vs. NRC by nominal mapped coverage"
)

ggsave("../paper/src/figs/nrc_coverage_splines.pdf")
```

from which we observe that effective coverage is a better predictor for non-reference concordance. As such, coverage-dependent metrics will be considered relative to effective coverage rather than nominal mapped coverage in the following. 

### Polygenic risk scores and effective coverage

#### Overall pattern 
Polygenic risk scores are calculated for brca and CAD for all individuals; the CAD score comprises ~1.7M variants while the brca score comprises only ~130 variants. If we treat the PRS calculated off the 1KGP3 gold standards as "truth", we can calculate the error of each estimate from imputed data: we can plot this against effective coverage. 

Let's start out (since plots in this section will probably make it into the paper) by recoding some variables like we did for the previous section

```{r}
# replace with human-readable experiment and coverage indications
scores_plot = scores %>% 
  mutate(
    # replace assay type with Experiment X
    Experiment = case_when(
      assay_type == "ilmn_0.5x" ~ "A",
      assay_type == "ilmn_1x" ~ "B",
      assay_type == "ilmn_1x_repl" ~ "C",
      assay_type == "bgi_1x" ~ "D", 
      assay_type == "ilmn_GSA" ~ "E"
    ),
    # coverage columns
    `Nominal Mapped Coverage` = coverage_deduped_mapped,
    `Effective Coverage` = eff_cov,
    # PRS names
    Trait = case_when(
      trait == "brca" ~ "Breast Cancer",
      trait == "cad" ~ "CAD"
    )
  )

scores_plot
```

and we can get started. We can then drill down to sequence scores and calculate absolute and squared error:

```{r}
# filter down to sequence samples only and add absolute error and squared error variables
seq_scores = scores_plot %>% 
  filter(assay_type != "ilmn_GSA") %>% 
  mutate(
    abs_error = abs(error),
    squared_error = error ^ 2,
    
    # the same, but for display
    `Absolute Error` = abs_error,
    `Squared Error` = squared_error
    )
seq_scores
```

#### Effective coverage, absolute error, squared error
We can get an idea of what the "typical" measurement error is in estimating an individual's PRS using a microarray by taking the array scores and calculating the absolute error as well as the squared error in the estimate for each sample, averaging across replicates, and then averaging across cell lines. We can also calculate the standard error of the mean for each of these:

```{r}
# take array scores absolute error and average among replicates, then group by trait and super pop to get population-level error
# we do this for squared error as well
array_score_summ = scores_plot %>% 
  filter(assay_type == "ilmn_GSA") %>% 
  mutate(
    abs_error = abs(error), # compute absolute error
    squared_error = error ^ 2 # compute squared error
    ) %>% 
  group_by(cell_line, Trait, super_pop) %>% 
  summarise(
    mean_abs_error = mean(abs_error), # mean across replicates
    mean_squared_error = mean(squared_error) # mean across replicates
    ) %>% 
  group_by(super_pop, Trait) %>% 
  summarise(
    # absolute error
    population_mean_abs_error = mean(mean_abs_error), # mean across samples
    population_std_abs_error = sqrt(var(mean_abs_error)), # std across samples
    population_mean_abs_error_stderr = population_std_abs_error / sqrt(n()),
    
    # mse
    population_mean_squared_error = mean(mean_squared_error),
    population_std_squared_error = sqrt(var(mean_squared_error)),
    population_mean_squared_error_stderr = population_std_squared_error / sqrt(n())
    )
array_score_summ
```

We can also obtain the average error by sequence experiment group

```{r}
# take MSE on sequence data similarly
seq_score_summ = seq_scores %>% 
  group_by(Experiment, Trait, super_pop, cell_line) %>% 
  summarise(cl_mse = mean(squared_error)) %>% 
  group_by(Experiment, Trait, super_pop) %>% 
  summarise(
    mse = mean(cl_mse),
    mse_se = sqrt(var(cl_mse)) / sqrt(n())
    )
seq_score_summ
```

We can then plot the absolute error of PRS estimated from sequence data and marking the mean absolute error for each super population and trait from array estimates: 

```{r echo = FALSE, fig.width =8 , fig.height = 6}
seq_scores %>% 
  ggplot(aes(x = `Effective Coverage`, y = `Absolute Error`, color = `Experiment`)) + 
  geom_point(alpha = 0.7) + 
  facet_wrap(~Trait + super_pop) + 
  theme_few() + 
  scale_color_tableau() + 
  labs(title = "Absolute error in PRS estimates by effective coverage") + 
  # array means as hlines
  geom_hline(
    data = array_score_summ, 
    aes(yintercept = population_mean_abs_error), 
    linetype = "dashed"
    ) +
  # shaded area about the mean 
  geom_rect(
    data = array_score_summ, 
    aes(xmin = -Inf, xmax = Inf, ymin = population_mean_abs_error - population_mean_abs_error_stderr, ymax = population_mean_abs_error + population_mean_abs_error_stderr), 
    alpha = 0.2, 
    inherit.aes = FALSE
    )

ggsave("../paper/src/figs/abserrs_prs_effcov.pdf")
```

We can also plot the same for squared error instead: with standard error 

```{r echo = FALSE, fig.width = 8, fig.height = 6}
# squared error in PRS estimates w/ standard errors
seq_scores %>% 
  ggplot(aes(x = `Effective Coverage`, y = `Squared Error`, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  facet_wrap(~Trait + super_pop, scales = "free") + 
  theme_few() + 
  scale_color_tableau() + 
  labs(title = "Squared error in PRS estimates by effective coverage") + 
  geom_hline( # horizontal lines for array means 
    data = array_score_summ, 
    aes(yintercept = population_mean_squared_error), 
    linetype = "dashed"
    ) +
  geom_rect( # shade the area within the mean standard errors 
    data = array_score_summ, 
    aes(xmin = -Inf, xmax = Inf, ymin = population_mean_squared_error - population_mean_squared_error_stderr,  ymax = population_mean_squared_error + population_mean_squared_error_stderr), 
    alpha = 0.2, 
    inherit.aes = FALSE) + 
  coord_cartesian(xlim = c(0, 2.5)) + 
  geom_segment(
    data = seq_score_summ, 
    aes(x = 2.5, xend = 2.7, y = mse, yend = mse)
    )

ggsave("../paper/src/figs/sqerrors_prs_effcov.pdf")
```

we can also write out the cohort level means to tables:

```{r}
# array score summary 
array_score_summ_table = array_score_summ %>% 
  select(super_pop, Trait, population_mean_squared_error, population_mean_squared_error_stderr) %>% 
  rename(
    `Super population` = super_pop,
    `MSE` = population_mean_squared_error,
    `SE of MSE` = population_mean_squared_error_stderr
  ) %>% 
  arrange(Trait)

array_score_summ_table %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat( file = "../paper/src/tabs/array-score-mse.tex")

array_score_summ_table
```

and similarly for sequence data

```{r}
# reshape to human readable.
# we also take care to join with the array mse in order to get the ratio 
array_mses_tojoin = array_score_summ %>% 
  select(super_pop, Trait, population_mean_squared_error) %>% 
  rename(array_mse = population_mean_squared_error)

# left join with seq score summs and get the ratio of array / seq 
seq_score_summ_table = seq_score_summ %>% 
  left_join(array_mses_tojoin, by = c("super_pop", "Trait")) %>% 
  mutate(
    `MSE fold decrease` = array_mse / mse
  ) %>% 
  rename(
    `Super population` = super_pop,
    `MSE` = mse,
    `SE of MSE` = mse_se
  ) %>% 
  select(-array_mse) %>% 
  arrange(Trait)

seq_score_summ_table %>% 
  kable("latex", booktabs = TRUE, linesep = "", digits = 4) %>% 
  cat(file = "../paper/src/tabs/seq-score-mse.tex")

seq_score_summ_table
```

### Significance tests

We can perform significance tests for differences in means of the squared error in PRS estimates by trait, population, and comparing the sequence experiment results to the array estimates. 

```{r}
# inelegant variable naming
exp_A_eur_cad = scores_plot %>% filter(Experiment == "A", super_pop == "EUR", trait == "cad") %>% pull(squared_error)
exp_A_eur_brca = scores_plot %>% filter(Experiment == "A", super_pop == "EUR", trait == "brca") %>% pull(squared_error)
exp_B_eur_cad = scores_plot %>% filter(Experiment == "B", super_pop == "EUR", trait =="cad") %>% pull(squared_error)
exp_B_eur_brca = scores_plot %>% filter(Experiment == "B", super_pop == "EUR", trait == "brca") %>% pull(squared_error)
exp_C_eur_cad = scores_plot %>% filter(Experiment == "C", super_pop == "EUR", trait == "cad") %>% pull(squared_error)
exp_C_eur_brca = scores_plot %>% filter(Experiment == "C", super_pop == "EUR", trait == "brca") %>% pull(squared_error)
exp_D_eur_cad = scores_plot %>% filter(Experiment == "D", super_pop == "EUR", trait == "cad") %>% pull(squared_error)
exp_D_eur_brca = scores_plot %>% filter(Experiment == "D", super_pop == "EUR", trait == "brca") %>% pull(squared_error)
exp_E_eur_cad = scores_plot %>% filter(Experiment == "E", super_pop == "EUR", trait == "cad") %>% pull(squared_error)
exp_E_eur_brca = scores_plot %>% filter(Experiment == "E", super_pop == "EUR", trait == "brca") %>% pull(squared_error)

exp_A_afr_cad = scores_plot %>% filter(Experiment == "A", super_pop == "AFR", trait == "cad") %>% pull(squared_error)
exp_A_afr_brca = scores_plot %>% filter(Experiment == "A", super_pop == "AFR", trait == "brca") %>% pull(squared_error)
exp_B_afr_cad = scores_plot %>% filter(Experiment == "B", super_pop == "AFR", trait == "cad") %>% pull(squared_error)
exp_B_afr_brca = scores_plot %>% filter(Experiment == "B", super_pop == "AFR", trait == "brca") %>% pull(squared_error)
# exp_C_afr_cad = scores_plot %>% filter(Experiment == "C", super_pop == "AFR", trait == "cad") %>% pull(squared_error)
# exp_C_afr_brca = scores_plot %>% filter(Experiment == "C", super_pop == "AFR", trait == "brca") %>% pull(squared_error)
exp_D_afr_cad = scores_plot %>% filter(Experiment == "D", super_pop == "AFR", trait == "cad") %>% pull(squared_error)
exp_D_afr_brca = scores_plot %>% filter(Experiment == "D", super_pop == "AFR", trait == "brca") %>% pull(squared_error)
exp_E_afr_cad = scores_plot %>% filter(Experiment == "E", super_pop == "AFR", trait == "cad") %>% pull(squared_error)
exp_E_afr_brca = scores_plot %>% filter(Experiment == "E", super_pop == "AFR", trait == "brca") %>% pull(squared_error)

print("==========EUR=========")
paste("t test p value EUR exp A vs array, cad", t.test(exp_A_eur_cad, exp_E_eur_cad)$p.value)
paste("t test p value EUR exp A vs array, brca", t.test(exp_A_eur_brca, exp_E_eur_brca)$p.value)
paste("t test p value EUR exp B vs array, cad", t.test(exp_B_eur_cad, exp_E_eur_cad)$p.value)
paste("t test p value EUR exp B vs array, brca", t.test(exp_B_eur_brca, exp_E_eur_brca)$p.value)
paste("t test p value EUR exp C vs array, cad", t.test(exp_C_eur_cad, exp_E_eur_cad)$p.value)
paste("t test p value EUR exp C vs array, brca", t.test(exp_C_eur_brca, exp_E_eur_brca)$p.value)
paste("t test p value EUR exp D vs array, cad", t.test(exp_D_eur_cad, exp_E_eur_cad)$p.value)
paste("t test p value EUR exp D vs array, brca", t.test(exp_D_eur_brca, exp_E_eur_brca)$p.value)

print("==========AFR=========")
paste("t test p value AFR exp A vs array, cad", t.test(exp_A_afr_cad, exp_E_afr_cad)$p.value)
paste("t test p value AFR exp A vs array, brca", t.test(exp_A_afr_brca, exp_E_afr_brca)$p.value)
paste("t test p value AFR exp B vs array, cad", t.test(exp_B_afr_cad, exp_E_afr_cad)$p.value)
paste("t test p value AFR exp B vs array, brca", t.test(exp_B_afr_brca, exp_E_afr_brca)$p.value)
paste("t test p value AFR exp D vs array, cad", t.test(exp_D_afr_cad, exp_E_afr_cad)$p.value)
paste("t test p value AFR exp D vs array, brca", t.test(exp_D_afr_brca, exp_E_afr_brca)$p.value)
```


## Head-to-head results

In order to obtain directly comparable results, we chose one sample for each cell line for each assay type for each super population at random out of the replicates. We can filter down to those individuals and then aggregate by assay types.

### Concordance
We filter down the results to the samples in question. Then, as before, since stuff in this section will eventually be put into figures, we must make a version where the column names and values more clear to humans:

```{r}
# just filtered down, no annotations
keep_concordance = concordance %>% 
  filter(sample %in% random_samples$sample)

# for plotting down the line
keep_concordance_plot = keep_concordance %>% 
  mutate(
    # replace assay type with Experiment X
    Experiment = case_when(
      assay_type == "ilmn_0.5x" ~ "A",
      assay_type == "ilmn_1x" ~ "B",
      assay_type == "bgi_1x" ~ "D",
      assay_type == "ilmn_GSA" ~ "E"
    ),
    # unfiltered status 
    `Site Type` = case_when(
      site_type == "allsites" ~ "Unfiltered",
      site_type == "passing" ~ "Filtered",
    ),
    # allele frequency
    `Non-Reference Allele Frequency` = af_midpoint,
    
    `Super population` = super_pop
  )
  
head(keep_concordance)
```

where we pause to note exactly what samples we have here. In order to calculate comparable statistics _across assay types_ for a _given ancestry group_, a cell line must have at least one replicate in each assay type that passed QC upstream. That was not the case for all cell lines: in particular, there were three EUR cell lines which did not satisfy this condition and therefore were excluded.

Specifically, the two TSI cell lines NA20802 and NA20822 did not have a single passing sample, and HG00096 did not have any successful 1x Illumina samples. 

We can tally the number of samples in each of these groups to confirm:

```{r}
keep_concordance %>% 
  select(cell_line, super_pop, assay_type, sample) %>% 
  distinct(sample, .keep_all = TRUE) %>% 
  group_by(super_pop, assay_type) %>% 
  tally()
```

Also, it's important to note that although I'm including the BGI samples here, they can't be as *directly* compared since they aren't the same samples and the $n$ is lower. 

#### Overall results
We can calculate sample-level results and then aggregate by assay type: for each assay type, site type, stat type, and super population, we can take the mean concordance and NRC of the samples in each of those bins. 

##### Coverage and effective coverage for random samples
Let's take a look at what the actual coverages look like for the representative samples:

```{r}
keep_concordance_plot %>% 
  distinct(cell_line, sample, coverage_deduped_mapped, eff_cov, assay_type, super_pop) %>%
  group_by(assay_type, super_pop) %>% 
  summarise(
    mean_eff_cov = mean(eff_cov),
    std_eff_cov = sqrt(var(eff_cov)),
    mean_nominal_cov = mean(coverage_deduped_mapped), 
    std_nominal_cov = sqrt(var(coverage_deduped_mapped))

    )
```

or, if we look at overall experiments (without stratifying into pops),

```{r}
# coverage by experiment 
keep_concordance_plot %>% 
  distinct(cell_line, sample, coverage_deduped_mapped, eff_cov, assay_type, super_pop) %>%
  group_by(assay_type) %>% 
  summarise(
     mean_eff_cov = mean(eff_cov),
    std_eff_cov = sqrt(var(eff_cov)),
    mean_nominal_cov = mean(coverage_deduped_mapped), 
    std_nominal_cov = sqrt(var(coverage_deduped_mapped))

  )
```

and we can also write out all the samples actually kept, along with their metadata
```{r}
keep_concordance_plot %>% 
  distinct(cell_line, sample, pop, super_pop, Experiment) %>% 
  write_csv("../data/out/representative-cohorts.csv")
```


##### Performance
Since we are going to be printing out tables of summary statistics we should also recode the assay type to the Experiment nomenclature.

```{r echo = FALSE}
# aggregate to sample level for kept samples
keep_sample_concordance = keep_concordance %>% 
  group_by(sample, cell_line, stat_type, site_type, super_pop, assay_type) %>% 
  summarise(
    overall_concordance = sum(total_matches) / sum(total_sites),
    nrc = 1 - sum(nrd_numerator) / sum(nrd_denominator),
    total_nrd_numerator = sum(nrd_numerator),
    total_nrd_denominator = sum(nrd_denominator)
    )

# aggregate by assay type, stat type, site type, super pop
# keep_sample_concordance_aggregated = keep_sample_concordance %>% 
#   group_by(super_pop, stat_type, site_type, assay_type) %>% 
#   summarise(
#     mean_concordance = mean(overall_concordance), 
#     mean_nrc = mean(nrc)
#     ) %>% 
#   ungroup()

keep_sample_concordance_aggregated = keep_sample_concordance %>% 
  group_by(super_pop, stat_type, site_type, assay_type) %>% 
  summarise(
    mean_concordance = mean(overall_concordance), 
    mean_nrc = mean(nrc),
    concordance_quantile_25 = quantile(overall_concordance, 1/4),
    concordance_quantile_75 = quantile(overall_concordance, 3/4),
    nrc_quantile_25 = quantile(nrc, 1/4),
    nrc_quantile_75 = quantile(nrc, 3/4),
    
    # display results, which is `mean (25th - 75th)` quantiles
    mean_concordance_display = paste(round(mean_concordance, 4), " (", round(concordance_quantile_25, 4), "-", round(concordance_quantile_75, 4), ")", sep = ""),
    mean_nrc_display = paste(round(mean_nrc, 4), " (", round(nrc_quantile_25, 4), "-", round(nrc_quantile_75, 4), ")", sep = "")
    ) %>% 
  ungroup()

# add experiment variable
keep_sample_concordance_aggregated = keep_sample_concordance_aggregated %>% 
  mutate(
     # replace assay type with Experiment X
    Experiment = case_when(
      assay_type == "ilmn_0.5x" ~ "A",
      assay_type == "ilmn_1x" ~ "B",
      assay_type == "bgi_1x" ~ "D",
      assay_type == "ilmn_GSA" ~ "E"
    ),
    `Super Population` = super_pop
  )

keep_sample_concordance_aggregated
```

which is still grouped by quite a few variables. Comparing assay types directly is a bit clearer if we restrict ourselves to NRC for unfiltered SNPs.

Note that in the course of coming up with these data we also save them in TeX format to disk so that we can directly import them into our manuscript

### Overall results, NRC
```{r}
# show mean NRC across samples for unfiltered SNPs by assay type and super pop
# unfiltered_snps_mean_nrc = keep_sample_concordance_aggregated %>%
#   filter(stat_type == "GCsAF", site_type == "allsites") %>% 
#   select(`Super Population`, Experiment, mean_nrc) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_nrc) %>% 
#   select(`Super Population`, sort(current_vars())) %>% 
#   rename(`E (array)` = E) 

unfiltered_snps_mean_nrc = keep_sample_concordance_aggregated %>%
  filter(stat_type == "GCsAF", site_type == "allsites") %>% 
  select(`Super Population`, Experiment, mean_nrc_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_nrc_display) %>% 
  select(`Super Population`, sort(current_vars())) %>% 
  rename(`E (array)` = E) 

unfiltered_snps_mean_nrc %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  add_header_above(c(" ", "Experiment" = 4)) %>% 
  cat(file = "../paper/src/tabs/unfiltered-snps-mean-nrc.tex")

unfiltered_snps_mean_nrc 
```
from which a sharp contrast can be immediately drawn across assay types and populations --- basically, for this particular set of sites the GSA performs more poorly than all the sequence results. If we consider filtered SNPs instead we find 

```{r}
# do the same but for filtered sites
# filtered_snps_mean_nrc = keep_sample_concordance_aggregated %>% 
#   filter(stat_type == "GCsAF", site_type == "passing") %>% 
#   select(`Super Population`, Experiment, mean_nrc) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_nrc) %>% 
#   select(`Super Population`, sort(current_vars()))

filtered_snps_mean_nrc = keep_sample_concordance_aggregated %>% 
  filter(stat_type == "GCsAF", site_type == "passing") %>% 
  select(`Super Population`, Experiment, mean_nrc_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_nrc_display) %>% 
  select(`Super Population`, sort(current_vars()))


filtered_snps_mean_nrc %>%
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat(file = "../paper/src/tabs/filtered-snps-mean-nrc.tex")

filtered_snps_mean_nrc
```
we find less of a difference; in fact, the GSA in Europeans outperforms the 0.5x Illumina sequence data, while still performing less well on African populations.  

We can do the same for NRC for filtered and unfiltered indels:
```{r}
# show mean NRC across samples for unfiltered indels by assay type and super pop
# unfiltered_indels_mean_nrc = keep_sample_concordance_aggregated %>%
#   filter(stat_type == "GCiAF", site_type == "allsites") %>% 
#   select(`Super Population`, Experiment, mean_nrc) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_nrc) %>% 
#   select(`Super Population`, sort(current_vars()))
unfiltered_indels_mean_nrc = keep_sample_concordance_aggregated %>%
  filter(stat_type == "GCiAF", site_type == "allsites") %>% 
  select(`Super Population`, Experiment, mean_nrc_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_nrc_display) %>% 
  select(`Super Population`, sort(current_vars()))

unfiltered_indels_mean_nrc %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat(file = "../paper/src/tabs/unfiltered-indels-mean-nrc.tex")

unfiltered_indels_mean_nrc
```

for filtered sites:

```{r}
# mean NRC for filtered indels
# filtered_indels_mean_nrc = keep_sample_concordance_aggregated %>%
#   filter(stat_type == "GCiAF", site_type == "passing") %>% 
#   select(`Super Population`, Experiment, mean_nrc) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_nrc) %>% 
#   select(`Super Population`, sort(current_vars()))
filtered_indels_mean_nrc = keep_sample_concordance_aggregated %>%
  filter(stat_type == "GCiAF", site_type == "passing") %>% 
  select(`Super Population`, Experiment, mean_nrc_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_nrc_display) %>% 
  select(`Super Population`, sort(current_vars()))

filtered_indels_mean_nrc %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat(file = "../paper/src/tabs/filtered-indels-mean-nrc.tex")

filtered_indels_mean_nrc
```


### Overall results, overall concordance
Similarly, we can compute the same for overall concordance:

unfiltered:
```{r}
# show mean overall concordance across samples for unfiltered SNPs by assay type and super pop
# unfiltered_snps_mean_oconc = keep_sample_concordance_aggregated %>%
#   filter(stat_type == "GCsAF", site_type == "allsites") %>% 
#   select(`Super Population`, Experiment, mean_concordance) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_concordance)%>% 
#   select(`Super Population`, sort(current_vars()))
unfiltered_snps_mean_oconc = keep_sample_concordance_aggregated %>%
  filter(stat_type == "GCsAF", site_type == "allsites") %>% 
  select(`Super Population`, Experiment, mean_concordance_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_concordance_display)%>% 
  select(`Super Population`, sort(current_vars()))
# print to latex 
unfiltered_snps_mean_oconc %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat(file = "../paper/src/tabs/unfiltered-snps-mean-oconc.tex")

unfiltered_snps_mean_oconc
```

filtered:
```{r}
# show mean overall concordance across samples for filtered SNPs by assay type and super pop
# filtered_snps_mean_oconc = keep_sample_concordance_aggregated %>%
#   filter(stat_type == "GCsAF", site_type == "passing") %>% 
#   select(`Super Population`, Experiment, mean_concordance) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_concordance)%>% 
#   select(`Super Population`, sort(current_vars()))

filtered_snps_mean_oconc = keep_sample_concordance_aggregated %>%
  filter(stat_type == "GCsAF", site_type == "passing") %>% 
  select(`Super Population`, Experiment, mean_concordance_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_concordance_display)%>% 
  select(`Super Population`, sort(current_vars()))

filtered_snps_mean_oconc %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat(file = "../paper/src/tabs/filtered-snps-mean-oconc.tex")

filtered_snps_mean_oconc
```

and the same for indels: unfiltered
```{r}
# show mean overall concordance across samples for unfiltered indels by assay type and super pop
# unfiltered_indels_mean_oconc = keep_sample_concordance_aggregated %>%
#   filter(stat_type == "GCiAF", site_type == "allsites") %>% 
#   select(`Super Population`, Experiment, mean_concordance) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_concordance)%>% 
#   select(`Super Population`, sort(current_vars()))

unfiltered_indels_mean_oconc = keep_sample_concordance_aggregated %>%
  filter(stat_type == "GCiAF", site_type == "allsites") %>% 
  select(`Super Population`, Experiment, mean_concordance_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_concordance_display)%>% 
  select(`Super Population`, sort(current_vars()))

unfiltered_indels_mean_oconc %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat(file = "../paper/src/tabs/unfiltered-indels-mean-oconc.tex")

unfiltered_indels_mean_oconc
```

vs filtered:
```{r}
# show mean overall concordance across samples for filtered indels by assay type and super pop
# filtered_indels_mean_oconc = keep_sample_concordance_aggregated %>%
#   filter(stat_type == "GCiAF", site_type == "passing") %>% 
#   select(`Super Population`, Experiment, mean_concordance) %>% 
#   pivot_wider(names_from = Experiment, values_from = mean_concordance)%>% 
#   select(`Super Population`, sort(current_vars()))

filtered_indels_mean_oconc = keep_sample_concordance_aggregated %>%
  filter(stat_type == "GCiAF", site_type == "passing") %>% 
  select(`Super Population`, Experiment, mean_concordance_display) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_concordance_display)%>% 
  select(`Super Population`, sort(current_vars()))

filtered_indels_mean_oconc %>% 
  kable("latex", booktabs = TRUE, digits = 4) %>% 
  cat(file = "../paper/src/tabs/filtered-indels-mean-oconc.tex")

filtered_indels_mean_oconc
```

We can also compute the mean unfiltered overall concordance by experiment for variants with MAF >1% (for marketing one-pager). 

```{r}
# compute mean concordance metrics at common variants
sample_common_concordance = keep_concordance %>%
  filter(af_midpoint != 0.005) %>% 
  group_by(sample, cell_line, stat_type, site_type, super_pop, assay_type) %>%  
  summarise(
    overall_concordance = sum(total_matches) / sum(total_sites),
    nrc = 1 - sum(nrd_numerator) / sum(nrd_denominator),
    total_nrd_numerator = sum(nrd_numerator),
    total_nrd_denominator = sum(nrd_denominator)
    ) %>% 
  group_by(super_pop, stat_type, site_type, assay_type) %>% 
  summarise(
    mean_concordance = mean(overall_concordance), 
    mean_nrc = mean(nrc)
    ) %>% 
  ungroup() %>%  
  mutate(
        # replace assay type with Experiment X
        Experiment = case_when(
            assay_type == "ilmn_0.5x" ~ "A",
            assay_type == "ilmn_1x" ~ "B",
            assay_type == "bgi_1x" ~ "D",
            assay_type == "ilmn_GSA" ~ "E"
        ),
        `Super Population` = super_pop
    )

# mean concordance at AF >1%
sample_common_concordance %>% 
  filter(stat_type == "GCsAF" & site_type == "allsites") %>% 
  select(`Super Population`, Experiment, mean_concordance) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_concordance) %>% 
  select(`Super Population`, sort(current_vars()))

```

and the same for filtered

```{r}
# mean concordance at AF >1%, filtered
sample_common_concordance %>% 
  filter(stat_type == "GCsAF" & site_type == "passing") %>% 
  select(`Super Population`, Experiment, mean_concordance) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_concordance) %>% 
  select(`Super Population`, sort(current_vars()))

```

and we can also do this for maf >0.05:


```{r}
# concordance at common (af >5%) 
sample_common_concordance_05 = keep_concordance %>%
  filter(af_midpoint > 0.045) %>% 
  group_by(sample, cell_line, stat_type, site_type, super_pop, assay_type) %>%  
  summarise(
    overall_concordance = sum(total_matches) / sum(total_sites),
    nrc = 1 - sum(nrd_numerator) / sum(nrd_denominator),
    total_nrd_numerator = sum(nrd_numerator),
    total_nrd_denominator = sum(nrd_denominator)
    ) %>% 
  group_by(super_pop, stat_type, site_type, assay_type) %>% 
  summarise(
    mean_concordance = mean(overall_concordance), 
    mean_nrc = mean(nrc)
    ) %>% 
  ungroup() %>%  
  mutate(
        # replace assay type with Experiment X
        Experiment = case_when(
            assay_type == "ilmn_0.5x" ~ "A",
            assay_type == "ilmn_1x" ~ "B",
            assay_type == "bgi_1x" ~ "D",
            assay_type == "ilmn_GSA" ~ "E"
        ),
        `Super Population` = super_pop
    )

# mean concordance at AF >1%
sample_common_concordance_05 %>% 
  filter(stat_type == "GCsAF" & site_type == "passing") %>% 
  select(`Super Population`, Experiment, mean_concordance) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_concordance) %>% 
  select(`Super Population`, sort(current_vars()))

```
## Tests for significance

We can perform two-sample paired t-tests for differences in means between the methods, given superpopulation. 

```{r}
# snp data
tdata = keep_concordance %>%
  filter(af_midpoint > 0.045, stat_type == "GCsAF", site_type == "allsites") %>% 
  group_by(sample, cell_line, stat_type, site_type, super_pop, assay_type) %>%  
  summarise(
    overall_concordance = sum(total_matches) / sum(total_sites),
    nrc = 1 - sum(nrd_numerator) / sum(nrd_denominator),
    total_nrd_numerator = sum(nrd_numerator),
    total_nrd_denominator = sum(nrd_denominator)
    )

# compare against array results 
eur_array_nrcs = tdata %>% filter(super_pop == "EUR", assay_type == "ilmn_GSA") %>% ungroup() %>% select(nrc, cell_line)
afr_array_nrcs = tdata %>% filter(super_pop == "AFR", assay_type == "ilmn_GSA") %>% ungroup() %>% select(nrc, cell_line)


print("Two-sample paired t-tests for difference in means between sequence and array NRCs")
for (i in c("ilmn_0.5x", "ilmn_1x", "bgi_1x")) {
    df = tdata %>% filter(super_pop == "EUR", assay_type == i) %>% ungroup() %>% select(nrc, cell_line) %>%
      inner_join(eur_array_nrcs, by = "cell_line")
    print(paste("p value for eur", i, "vs array", t.test(df$nrc.x, df$nrc.y, paired = TRUE)$p.value))
    
    df = tdata %>% filter(super_pop == "AFR", assay_type == i) %>% ungroup() %>% select(nrc, cell_line) %>%
      inner_join(afr_array_nrcs, by = "cell_line")
    print(paste("p value for afr", i, "vs array", t.test(df$nrc.x, df$nrc.y, paired = TRUE)$p.value))
}
```
We can also perform _unpaired_ two-sample t-tests for differences in means between the AFR and EUR results

```{r}
print("two-sample unpaired t-tests for differences in means between populations")
for (i in c("ilmn_0.5x", "ilmn_1x", "bgi_1x")) {
  eur_nrcs = tdata %>% filter(assay_type == i) %>% ungroup() %>% filter(super_pop == "EUR") %>% pull(nrc)
  afr_nrcs = tdata %>% filter(assay_type == i) %>% ungroup() %>% filter(super_pop == "AFR") %>% pull(nrc)
  
  print(paste("p value for", i, "between eur and afr NRCs", t.test(eur_nrcs, afr_nrcs, paired = FALSE)$p.value))
}
```


#### By allele frequency 

##### Aggregating sample metrics for each bin
Here we calculate NRC for each sample for each bin and then average across samples for each bin. 

```{r}
# concordance by nraf 
keep_concordance_plot_bynraf = keep_concordance_plot %>% 
  group_by(Experiment, super_pop, `Site Type`, stat_type, `Non-Reference Allele Frequency`) %>% 
  summarise(
    `Mean NRC` = mean(nr_concordance), 
    n = n(),
    `Mean concordance` = mean(overall_concordance)
    ) 

keep_concordance_plot_bynraf = keep_concordance_plot %>% 
  group_by(Experiment, super_pop, `Site Type`, stat_type, `Non-Reference Allele Frequency`) %>% 
  summarise(
    `Mean NRC` = mean(nr_concordance), 
    n = n(),
    `Mean concordance` = mean(overall_concordance),
    
    # these are replicated to create the display names blow
    mean_nrc = mean(nr_concordance),
    mean_concordance = mean(overall_concordance),
    
    concordance_quantile_25 = quantile(overall_concordance, 1/4),
    concordance_quantile_75 = quantile(overall_concordance, 3/4),
    nrc_quantile_25 = quantile(nr_concordance, 1/4),
    nrc_quantile_75 = quantile(nr_concordance, 3/4),
    
    # display results, which is `mean (25th - 75th)` quantiles
    mean_concordance_display = paste(round(mean_concordance, 3), " (", round(concordance_quantile_25, 3), "-", round(concordance_quantile_75, 3), ")", sep = ""),
    mean_nrc_display = paste(round(mean_nrc, 3), " (", round(nrc_quantile_25, 3), "-", round(nrc_quantile_75, 3), ")", sep = "")
    ) 

keep_concordance_plot_bynraf
```

Since we're plotting here, let's use the new variables we renamed above:

```{r, fig.width=7, fig.height=4}
# for each AF bin, group by assay type and average the NRC for _all_ sites, not just passing
keep_concordance_plot_bynraf %>% 
  ungroup() %>% 
  mutate(
    Experiment = if_else(Experiment == "E", "E (array)", Experiment)
  ) %>% 
  filter(stat_type == "GCsAF" & `Site Type` == "Unfiltered") %>% 
  ggplot(
    aes(x = `Non-Reference Allele Frequency`, y = `Mean NRC`, color = Experiment)
    ) + 
  facet_wrap(~super_pop) + 
  geom_line() + geom_point() + theme_few() + scale_color_tableau() + 
  labs(title = "Average NRC by non-reference allele frequency\nfor unfiltered sites")

ggsave("../paper/src/figs/unfiltered_snp_nrc_by_nraf.pdf")
```
where we observe higher performance for imputed sequence data compared to imputed array data across the allele frequency spectrum for unfiltered variant calls. 

We can do this on a log scale too:

```{r, fig.width=7, fig.height=4}
# for each AF bin, group by assay type and average the NRC for _all_ sites, not just passing
keep_concordance_plot_bynraf %>% 
  ungroup() %>% 
  mutate(
    Experiment = if_else(Experiment == "E", "E (array)", Experiment)
  ) %>% 
  filter(stat_type == "GCsAF" & `Site Type` == "Unfiltered") %>% 
  ggplot(
    aes(x = `Non-Reference Allele Frequency`, y = `Mean NRC`, color = Experiment)
    ) + 
  facet_wrap(~super_pop) + 
  geom_line() + geom_point() + theme_few() + scale_color_tableau() + 
  labs(title = "Average NRC by non-reference allele frequency\nfor unfiltered sites") + 
  scale_x_log10()

ggsave("../paper/src/figs/unfiltered_snp_nrc_by_nraf_logscale.pdf")
```

Filtering down to passing sites, we observe

```{r, fig.width=7, fig.height=4}
# for each AF bin, group by assay type and average the NRC for passing sites
keep_concordance_plot_bynraf %>% 
  ungroup() %>% 
  filter(stat_type == "GCsAF" & `Site Type` == "Filtered") %>% 
  mutate(
    Experiment = if_else(Experiment == "E", "E (array)", Experiment)
  ) %>% 
  ggplot(
    aes(x = `Non-Reference Allele Frequency`, y = `Mean NRC`, color = Experiment)
    ) + 
  facet_wrap(~super_pop) + 
  geom_line() + geom_point() + theme_few() + scale_color_tableau() + 
  labs(title = "Average NRC by non-reference allele frequency for filtered sites")


ggsave("../paper/src/figs/filtered_snp_nrc_by_nraf.pdf")
```

We can also write out the exact numbers here to a latex table for completeness:

```{r}
# values plotted above for unfiltered snps NRC
keep_concordance_plot_bynraf %>%
  filter(stat_type == "GCsAF" & `Site Type` == "Unfiltered") %>%
  ungroup() %>%
  rename(`Super Population` = super_pop) %>%
  select(`Super Population`, `Non-Reference Allele Frequency`, `Mean NRC`, Experiment) %>%
  pivot_wider(names_from = Experiment, values_from = `Mean NRC`) %>%
  kable("latex", booktabs = TRUE, align = c("llrrrr"), digits = 4) %>%
  cat(file = "../paper/src/tabs/unfiltered-snp-nrc-by-nraf.tex")
# 
# keep_concordance_plot_bynraf %>%
#   filter(stat_type == "GCsAF" & `Site Type` == "Unfiltered") %>%
#   ungroup() %>%
#   rename(`Super Population` = super_pop, NRAF =  `Non-Reference Allele Frequency`) %>%
#   select(`Super Population`, NRAF, mean_nrc_display, Experiment) %>%
#   pivot_wider(names_from = Experiment, values_from = mean_nrc_display) %>%
#   kable("latex", booktabs = TRUE, align = c("llrrrr"), digits = 4) %>%
#   cat(file = "../paper/src/tabs/unfiltered-snp-nrc-by-nraf.tex")
```



### Imputation r2s

We have our binned $r^2$s for _biallelic SNPs_ and are interested in looking at the performance across allele frequencies: where the bin breakpoints are 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 and the x-axis is the _minor allele frequency_ in the 1KGP3.

Let's recode the assay type to Experiment notation
```{r}
binned_r2s = binned_r2s %>% 
  mutate(
    # replace assay type with Experiment X
    Experiment = case_when(
      assay_type == "ilmn_0.5x" ~ "A",
      assay_type == "ilmn_1x" ~ "B",
      assay_type == "bgi_1x" ~ "D",
      assay_type == "gsa" ~ "E"
    ),
    `Minor Allele Frequency` = af_midpoint,
  )
```


Looking at the number of variants in each of these bins, we are reminded that the _n_ in each bin differs from cohort to cohort due to the fact that a correlation can only be computed between two vectors if the elements of each are not all the same value. 

```{r}
# for EUR
binned_r2s %>% 
  filter(super_pop == "EUR") %>% 
  ggplot(aes(x = factor(af_midpoint), y = n_snps)) + geom_bar(stat = "identity") + 
  facet_wrap(~assay_type) + 
  theme(axis.text.x = element_text(angle = 90)) + labs(title = "EUR")
```
```{r}
# for AFR
binned_r2s %>% 
  filter(super_pop == "AFR") %>% 
  ggplot(aes(x = factor(af_midpoint), y = n_snps)) + geom_bar(stat = "identity") + 
  facet_wrap(~assay_type) + 
  theme(axis.text.x = element_text(angle = 90)) + labs(title = "AFR")
```

With that in mind, we can proceed to plotting the mean $r^2$ for each cohort for each assay type for each allele frequency bin. On a log scale, we have

```{r echo = FALSE, fig.width = 7, fig.height = 4}
binned_r2s %>% 
  mutate(
    Experiment = if_else(Experiment == "E", "E (array)", Experiment)
  ) %>% 
  ggplot(aes(x = `Minor Allele Frequency`, y = mean_r2, color = Experiment)) + 
  geom_point() + 
  facet_grid(cols = vars(super_pop)) + 
  geom_line() + 
  theme_few() + 
  scale_color_tableau() + 
  scale_x_log10() + 
  labs(
    title = expression(paste("Imputation ", "r"^{2}," vs. minor allele frequency in the 1KGP3")),
    y = expression("Mean "~r^2)
) 

ggsave("../paper/src/figs/r2_logscale.pdf")
```

Interestingly, the GSA performs better than all other assay types for Africans at the lowest AF bin. We can also plot this on a linear scale

```{r echo = FALSE, fig.width = 7, fig.height = 4}
binned_r2s %>% 
  mutate(
    Experiment = if_else(Experiment == "E", "E (array)", Experiment)
  ) %>% 
  ggplot(aes(x = `Minor Allele Frequency`, y = mean_r2, color = Experiment)) + 
  geom_point() + 
  facet_grid(cols = vars(super_pop)) + 
  geom_line() + 
  theme_few() + 
  scale_color_tableau() + 
  labs(
    title = expression(paste("Imputation ", "r"^{2}," vs. minor allele frequency in the 1KGP3")),
    y = expression("Mean "~r^2)
) 

ggsave("../paper/src/figs/r2_linearscale.pdf")
```
 
We can also write out the exact values plotted for each super population to latex tables

```{r}
# AFR
binned_r2s %>% 
  filter(super_pop == "AFR") %>% 
  select(mean_r2, `Minor Allele Frequency`, Experiment) %>% 
  rename(`Mean r^2` = mean_r2) %>% 
  pivot_wider(names_from = Experiment, values_from = `Mean r^2`) %>% 
  kable("latex", booktabs = TRUE, align = c("lrrrr"), digits = 4) %>% 
  cat(file = "../paper/src/tabs/r2-by-maf-afr.tex")

# EUR
binned_r2s %>% 
  filter(super_pop == "EUR") %>% 
  select(mean_r2, `Minor Allele Frequency`, Experiment) %>% 
  rename(`Mean r^2` = mean_r2) %>% 
  pivot_wider(names_from = Experiment, values_from = `Mean r^2`) %>% 
  kable("latex", booktabs = TRUE, align = c("lrrrr"), digits = 4) %>% 
  cat(file = "../paper/src/tabs/r2-by-maf-eur.tex")
```

We can summarize this by taking the average $r^2$s by experiment and population as well: 

```{r}
# we can back-calculate the average for bins 
mean_r2s = binned_r2s %>% 
  select(`Minor Allele Frequency`, mean_r2, super_pop, n_snps, Experiment) %>% 
  group_by(super_pop, Experiment, `Minor Allele Frequency`) %>% 
  mutate(r2xnsnps = mean_r2 * n_snps) %>% 
  group_by(Experiment, super_pop) %>% 
  filter(`Minor Allele Frequency` > 0.045) %>% # anything above a bin midpoint is >5% MAF
  summarise(
    mean_r2 = sum(r2xnsnps) / sum(n_snps)
    ) %>% 
  pivot_wider(names_from = Experiment, values_from = mean_r2) %>% 
  rename(`Super population` = super_pop)
# 
# mean_r2s %>% 
#   kable("latex", booktabs = TRUE, digits = 4) %>% 
#   add_header_above(c(" ", "Experiment" = 4)) %>% 
#   cat(file = "../paper/src/tabs/mean_r2s_common.tex")

mean_r2s
```

 
### Polygenic scores

We can examine the PRS head to head for these cohorts by filtering down to those samples we keep and then getting the MSE within each category:

```{r}
scores_keep = scores %>% filter(sample %in% random_samples$sample)
scores_keep_mse = scores_keep %>% 
  group_by(assay_type, super_pop, trait) %>% 
  summarise(
    mse = mean(squared_error), 
    mean_eff_cov = mean(eff_cov)
    )
head(scores_keep_mse)
```

and see the average squared error 

```{r}
scores_keep_mse %>% ggplot(aes(x = assay_type, fill = super_pop, y = mse)) + geom_bar(stat = "identity", position = "dodge") + facet_wrap(~trait) + theme_few() + scale_fill_tableau() + labs(title = "MSE for PRS by super population and trait")
```

from which we observe that the 1x-target-sequenced Illumina data have consistently lower measurement error compared to imputed GSA data across traits. However, if we pay attention to our nifty new metric $\lambda_{\mathrm{eff}}$, we can see that the target coverages of 0.5x and 1x for Illumina have rather lower mean effective coverages. 

```{r}
scores_keep_mse
```

We can additionally just see how the PRS correlates by experiment

```{r, fig.height = 7, fig.width = 7}
scores_plot %>% 
  ggplot(aes(x = score, y = true_score, color = Experiment)) + 
  facet_grid(cols = vars(super_pop), rows = vars(Trait)) + 
  geom_point(alpha = 0.7) + geom_abline() + theme_few() + scale_color_tableau() + 
  labs(
    title = "PRS estimates vs. true score by population and trait",
    x = "Score",
    y = "True Score"
  )

ggsave("../paper/src/figs/prs-estimated-true-scatterplot.pdf")
```

We can also write out the pearson correlations for each experiment for each super population for each trait

```{r}
# correlation between true and imputed scores
scores_correlation = scores_plot %>% 
  filter(Experiment != "C") %>% 
  group_by(Experiment, super_pop, Trait) %>% 
  summarise(correlation = cor(score, true_score) ^ 2)  %>% 
  pivot_wider(names_from = Experiment, values_from = correlation) %>% 
  arrange(Trait) %>% 
  rename(`Super population` = super_pop)

scores_correlation %>% 
  kable("latex", booktabs = TRUE,  digits = 4) %>% 
  add_header_above(c(" " = 2, "Experiment" = 4)) %>% 
  cat(file = "../paper/src/tabs/prs-corr-with-truth.tex")

scores_correlation
```


<!-- #### PRS precision/accuracy  -->
<!-- We can examine the "inter-rater reliability" of each assay type in estimating a PRS. In other words, since apolygenic risk score is computed as the dot product of a vector of weights and a vector of allelic dosages, the consistency of a PRS estimated on biological replicates is of interest. We can quantify this by calculating the intraclass correlation coefficient (ICC) of PRS estimates derived from multiple samples of the same individual for each assay type.  -->

## Tables as latex

Copy/paste these code snippets and run them in a console to get the latex code:

- Sample breakdown by population as table:

```{r}
concordance %>%
  distinct(cell_line, super_pop, pop) %>% 
  group_by(pop, super_pop) %>% 
  tally() %>% 
  arrange(super_pop, -n) %>% 
  kable("latex", booktabs = TRUE)
```


# Blog post figures

We want to have some figures for our blog post/ one pager. In particular we want the $r^2$ and the PRS figures, the same but with different legends:

```{r echo = FALSE, fig.width = 9, fig.height = 4}
# r^2 by MAF, different legends
binned_r2s %>% 
  mutate(
    # rename labels
    Experiment = case_when(
      Experiment == "A" ~ "Illumina 0.41x effective coverage",
      Experiment == "B" ~ "Illumina 0.72x effective coverage",
      Experiment == "C" ~ "Illumina 0.53x effective coverage",
      Experiment == "D" ~ "BGI 1.2x effective coverage",
      Experiment == "E" ~ "Illumina GSA"
    ),
    # relevel
    Experiment = fct_relevel(
      factor(Experiment), 
      "Illumina 0.41x effective coverage", 
      "Illumina 0.72x effective coverage",
      "Illumina 0.53x effective coverage",
      "BGI 1.2x effective coverage")
    
  ) %>% 
  ggplot(aes(x = `Minor Allele Frequency`, y = mean_r2, color = Experiment)) + 
  geom_point() + 
  facet_grid(cols = vars(super_pop)) + 
  geom_line() + 
  theme_few() + 
  scale_color_tableau() + 
  labs(
    title = expression(paste("Imputation ", "r"^{2}," vs. minor allele frequency in the 1KGP3")),
    y = expression("Mean "~r^2)
  )

ggsave("../paper/src/figs/r2_linearscale_forblog.pdf")
```

and we can also plot this with just the BGI and GSA data. 

```{r echo = FALSE, fig.width = 9, fig.height = 4}
# r^2 by MAF, different legends
binned_r2s %>% 
  mutate(
    # rename labels
    Experiment = case_when(
      Experiment == "A" ~ "Illumina 0.41x effective coverage",
      Experiment == "B" ~ "Illumina 0.72x effective coverage",
      Experiment == "C" ~ "Illumina 0.53x effective coverage",
      Experiment == "D" ~ "BGI 1.2x effective coverage",
      Experiment == "E" ~ "Illumina GSA"
    ),
    # relevel
    Experiment = fct_relevel(
      factor(Experiment), 
      "Illumina 0.41x effective coverage", 
      "Illumina 0.72x effective coverage",
      "Illumina 0.53x effective coverage",
      "BGI 1.2x effective coverage")
    
  ) %>% 
  filter(Experiment %in% c("BGI 1.2x effective coverage", "Illumina GSA")) %>% 
  ggplot(aes(x = `Minor Allele Frequency`, y = mean_r2, color = Experiment)) + 
  geom_point() + 
  facet_grid(cols = vars(super_pop)) + 
  geom_line() + 
  theme_few() + 
  scale_color_tableau() + 
  labs(
    title = expression(paste("Imputation ", "r"^{2}," vs. minor allele frequency in the 1KGP3")),
    y = expression("Mean "~r^2)
  )

ggsave("../paper/src/figs/r2_linearscale_forblog_gsa_bgi.pdf")
```


and the same for PRS

```{r echo = FALSE, fig.width = 10, fig.height = 6}
# squared error in PRS estimates w/ standard errors
seq_scores %>% 
   mutate(
    # rename labels
    Experiment = case_when(
      Experiment == "A" ~ "Illumina 0.41x effective coverage",
      Experiment == "B" ~ "Illumina 0.72x effective coverage",
      Experiment == "C" ~ "Illumina 0.53x effective coverage",
      Experiment == "D" ~ "BGI 1.2x effective coverage",
      Experiment == "E" ~ "Illumina GSA"
    ),
    # relevel
    Experiment = fct_relevel(
      factor(Experiment), 
      "Illumina 0.41x effective coverage", 
      "Illumina 0.72x effective coverage",
      "Illumina 0.53x effective coverage",
      "BGI 1.2x effective coverage")
    
  ) %>% 
  ggplot(aes(x = `Effective Coverage`, y = `Squared Error`, color = Experiment)) + 
  geom_point(alpha = 0.7) + 
  facet_wrap(~Trait + super_pop, scales = "free") + 
  theme_few() + 
  scale_color_tableau() + 
  labs(title = "Squared error in PRS estimates by effective coverage") + 
  geom_hline( # horizontal lines for array means 
    data = array_score_summ, 
    aes(yintercept = population_mean_squared_error), 
    linetype = "dashed",
    show.legend = FALSE
    ) +
  geom_rect( # shade the area within the mean standard errors 
    data = array_score_summ, 
    aes(xmin = -Inf, 
        xmax = Inf, 
        ymin = population_mean_squared_error - population_mean_squared_error_stderr,  
        ymax = population_mean_squared_error + population_mean_squared_error_stderr), 
    alpha = 0.2, 
    inherit.aes = FALSE,
    show.legend = FALSE) + 
  coord_cartesian(xlim = c(0, 2.5)) + 
  geom_segment(
    # reorder the segment labels as well, in order to properly illustrate the mean by experiment 
    data = seq_score_summ %>% ungroup() %>% mutate(
    # rename labels
    Experiment = case_when(
      Experiment == "A" ~ "Illumina 0.41x effective coverage",
      Experiment == "B" ~ "Illumina 0.72x effective coverage",
      Experiment == "C" ~ "Illumina 0.53x effective coverage",
      Experiment == "D" ~ "BGI 1.2x effective coverage",
      Experiment == "E" ~ "Illumina GSA"
    ),
    # relevel
    Experiment = fct_relevel(
      factor(Experiment), 
      "Illumina 0.41x effective coverage", 
      "Illumina 0.72x effective coverage",
      "Illumina 0.53x effective coverage",
      "BGI 1.2x effective coverage")
    
  ), 
    aes(x = 2.5, xend = 2.7, y = mse, yend = mse, color = Experiment),
    show.legend = FALSE,
    inherit.aes = FALSE
    )

ggsave("../paper/src/figs/sqerrors_prs_effcov_forblog.pdf")
```


# Revisions 

Revisions requested included an analysis using different imputation programs for both the array and sequence data. 
We can do this on a subset of the individuals (the randomly selected subset) and using GLIMPSE and IMPUTE5 for the sequence and array data respectively. 

<!-- - Unfiltered SNPs mean NRC for random samples -->

<!-- ```{r} -->
<!-- unfiltered_snps_mean_nrc %>%  -->
<!--   kable("latex", booktabs = TRUE) -->
<!-- ``` -->

<!-- - Filtered SNPs mean NRC for random samples -->

<!-- ```{r} -->
<!-- filtered_snps_mean_nrc %>%  -->
<!--   kable("latex", booktabs = TRUE) -->
<!-- ``` -->

<!-- - Unfiltered indels mean NRC for random samples -->
<!-- ```{r} -->
<!-- unfiltered_indels_mean_nrc %>%  -->
<!--   kabel("latex", booktabs = TRUE) -->
<!-- ``` -->

<!-- - Filtered indels mean NRC for random samples -->
<!-- ```{r} -->
<!-- filtered_indels_mean_nrc %>%  -->
<!--   kabel("latex", booktabs = TRUE) -->
<!-- ``` -->


<!-- It is at this juncture that one stumbles upon the fact that this is a pretty misleading plot and doesn't seem to comport with the previous results. Why is it that if the mean NRC is so much higher for AFRs across the board at the lowest AF bin (and taking also into account the previous fact that the # of sites in that bin is huge)?  -->

<!-- It's because the NRC (for say a given AF bin) is computed at a different set of sites for each sample, since recall that for two genotypes side by side,  -->

<!-- | |0|1|2| -->
<!-- |-|-|-|-| -->
<!-- |0|a|b|c| -->
<!-- |1|d|e|f| -->
<!-- |2|g|h|i| -->

<!-- we define $\mathrm{NRC} = (e + i) / (b + c + d + e + f + g + h + i)$. $a$ is different for each sample since it depends on the calls.  -->

<!-- ##### The right way  -->
<!-- The right way to do this (having actually just now stumbled into this problem) would be to actually recalculate the NRC in each bin across all samples. In other words, since the denominator differs across samples, we can't just take the mean of a fraction; we need to sum the numerators and denominators and do the division at that point.  -->

<!-- Let's backtrack to where we had results for each sample at each allele frequency bin.  -->

<!-- ```{r} -->
<!-- keep_concordance %>% group_by(assay_type, super_pop, site_type, stat_type, af_midpoint) %>%  -->
<!--   summarise( -->
<!--     mean_concordance = sum(total_matches) / sum(total_sites),  -->
<!--     mean_nrc = 1 - sum(nrd_numerator) / sum(nrd_denominator), -->
<!--     nsites = sum(total_sites), # total sites -->
<!--     nsite_nrc_denom = sum(nrd_denominator) # total sites used in each individual sample's NRC calculation -->
<!--     ) -->
<!-- ``` -->


<!-- ##### Aggregating across all samples for each bin -->
<!-- This is an example of one way of aggregating and illustrating the NRC across assays. In this example we compute the NRC for each bin by adding up the numerators for all samples and dividing through by the sum of all denominators.  -->

<!-- We then aggregate across assay type to obtain concordance and NRC by allele frequency bin, variant type, imputation quality, super population, and assay type.  -->

<!-- ```{r} -->
<!-- # aggregate concordance by allele frequency  -->
<!-- keep_concordance_aggregated_af = keep_concordance %>% group_by(assay_type, super_pop, site_type, stat_type, af_midpoint) %>%  -->
<!--   summarise( -->
<!--     mean_concordance = sum(total_matches) / sum(total_sites),  -->
<!--     mean_nrc = 1 - sum(nrd_numerator) / sum(nrd_denominator), -->
<!--     nsites = sum(total_sites), # total sites -->
<!--     total_nrd_numerator = sum(nrd_numerator), -->
<!--     total_nrd_denominator = sum(nrd_denominator) -->
<!--     ) -->
<!-- head(keep_concordance_aggregated_af) -->
<!-- ``` -->

<!-- Let's check the sample size of each bin: -->

<!-- ```{r echo = FALSE} -->
<!-- keep_concordance_aggregated_af %>%  -->
<!--   filter(assay_type == "bgi_1x" & site_type == "allsites") %>%  -->
<!--   ggplot(aes(x = factor(af_midpoint), y = total_nrd_denominator)) + geom_bar(stat = "identity") +  -->
<!--   facet_wrap(~stat_type) +  -->
<!--   theme_few() +  -->
<!--   labs(title = "Number of variants by variant type for NRD denominator") +  -->
<!--   theme(axis.text.x = element_text(angle = 90)) -->
<!-- ``` -->
<!-- Note that this displays a different pattern than if we look at the distribution of _all_ sites by their allele frequency:  -->

<!-- ```{r echo = FALSE} -->
<!-- keep_concordance_aggregated_af %>%  -->
<!--   filter(assay_type == "bgi_1x" & site_type == "allsites") %>%  -->
<!--   ggplot(aes(x = factor(af_midpoint), y = nsites)) + geom_bar(stat = "identity") +  -->
<!--   facet_wrap(~stat_type) +  -->
<!--   theme_few() +  -->
<!--   labs(title = "Number of variants by variant type for all sites") +  -->
<!--   theme(axis.text.x = element_text(angle = 90)) -->
<!-- ``` -->

<!-- With that in mind, we can plot out the mean NRC and overall concordance by non-reference allele frequency on a log scale for SNPs: -->
<!-- ```{r echo = FALSE, fig.height=10, fig.width=7} -->
<!-- keep_concordance_aggregated_af %>%  -->
<!--   filter(stat_type == "GCsAF") %>%  -->
<!--   ggplot(aes(x = af_midpoint, y = mean_nrc, color = assay_type)) +  -->
<!--   facet_wrap(~site_type + super_pop) +  -->
<!--   theme_few() +  -->
<!--   scale_color_tableau() +  -->
<!--   geom_point() +  -->
<!--   geom_line() + -->
<!--   scale_x_log10() +  -->
<!--   labs(title = "NRC by non-reference allele frequency in 1KGP3") -->
<!-- ``` -->

<!-- and if we want to show the sample size with bubble size:  -->
<!-- ```{r echo = FALSE, fig.height=10, fig.width=7} -->
<!-- keep_concordance_aggregated_af %>%  -->
<!--   filter(stat_type == "GCsAF") %>%  -->
<!--   ggplot(aes(x = af_midpoint, y = mean_nrc, color = assay_type, size = total_nrd_denominator)) +  -->
<!--   facet_wrap(~site_type + super_pop) +  -->
<!--   theme_few() +  -->
<!--   scale_color_tableau() +  -->
<!--   geom_point() +  -->
<!-- #  geom_line() +  -->
<!--   scale_x_log10() +  -->
<!--   labs(title = "NRC by non-reference allele frequency in 1KGP3") -->
<!-- ``` -->

<!-- and standard deviations -->

<!-- ```{r echo = FALSE, fig.height = 7, fig.width = 7} -->
<!-- # squared error in PRS estimates w/ standard deviations -->
<!-- seq_scores %>%  -->
<!--   ggplot(aes(x = eff_cov, y = squared_error, color = assay_type)) +  -->
<!--   geom_point() +  -->
<!--   facet_wrap(~trait + super_pop, scales = "free") +  -->
<!--   theme_few() +  -->
<!--   scale_color_tableau() +  -->
<!--   labs(title = "Squared error in PRS estimates by effective coverage\nwith standard deviations") +  -->
<!--   geom_hline(data = array_score_summ, aes(yintercept = population_mean_squared_error), linetype = "dashed") + -->
<!--   geom_rect(data = array_score_summ, aes(xmin = -Inf, xmax = Inf, ymin = population_mean_squared_error - population_std_squared_error,  ymax = population_mean_squared_error + population_std_squared_error), alpha = 0.2, inherit.aes = FALSE) -->
<!-- ``` -->

<!-- Alternatively we can have the error bars be the standard deviation -->
<!-- ```{r echo = FALSE, fig.height = 7, fig.width = 7} -->
<!-- seq_scores %>%  -->
<!--   ggplot(aes(x = eff_cov, y = abs_error, color = assay_type)) +  -->
<!--   geom_point() +  -->
<!--   facet_wrap(~trait + super_pop) +  -->
<!--   theme_few() +  -->
<!--   scale_color_tableau() +  -->
<!--   labs(title = "Absolute error in PRS estimates by effective coverage with\nstandard deviations") +  -->
<!--   geom_hline(data = array_score_summ, aes(yintercept = population_mean_abs_error), linetype = "dashed") + -->
<!--   geom_rect(data = array_score_summ, aes(xmin = -Inf, xmax = Inf, ymin = population_mean_abs_error - population_std_abs_error , ymax = population_mean_abs_error +population_std_abs_error ), alpha = 0.2, inherit.aes = FALSE) -->
<!-- ``` -->








