library(data.table)
library(dplyr)
library(tidyr)
library(ggplot2)

# Question 2.4
# Can you create a different QQ plot for each of the tests that have been run?
dt <- fread("ukbiobank/pilot-traits_uk-biobank_gene_uk-biobank.palmer.pilot.AFib.JULY23Freeze.ALL.EUR.28671.373704.SAIGE.gene.20240110.txt.gz")
dt <- dt %>% filter(Group != "Cauchy")
dt <- dt %>% 
  pivot_longer(cols = c(Pvalue, Pvalue_Burden, Pvalue_SKAT),
               names_to = "class", values_to = "Pvalue") %>%
  mutate(
    class = recode(class, 
                   "Pvalue" = "SKAT-O",
                   "Pvalue_Burden" = "Burden",
                   "Pvalue_SKAT" = "SKAT"),
    Group = recode(Group,
                   "damaging_missense_or_protein_altering" = "Damaging missense",
                   "other_missense_or_protein_altering" = "Other missense",
                   "pLoF" = "pLoF",
                   "pLoF;damaging_missense_or_protein_altering" = "pLoF, damaging missense",
                   "pLoF;damaging_missense_or_protein_altering;other_missense_or_protein_altering;synonymous" = "pLoF, damaging missense,\nother missense, synonymous",
                   "synonymous" = "Synonymous",
    )
  ) %>%
  group_by(class, Group, max_MAF) %>% 
  setorder("Pvalue") %>%
  mutate(Pvalue_exp = seq(1,n())/(n()+1))

ggplot(dt %>% setorder("class"), aes(x=-log10(Pvalue_exp), y=-log10(Pvalue), col=class)) + 
  geom_point() + facet_wrap(~max_MAF+Group, scales = 'free') + 
  theme_bw() + geom_abline(intercept=0,slope=1, col='grey')
# Note: you'll need to wait a little while for the plots to arrive!