Tables summarising the GWAS results

library(tidyverse)
library(kableExtra)
library(DT)

kable_table <- function(df) {
  kable(df, "html") %>%
  kable_styling() %>%
  scroll_box(height = "300px")
}

my_data_table <- function(df){ # Make html tables:
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons = 
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'GWAS_sig_loci')),
      pageLength = 50
    )
  )
}

db <- DBI::dbConnect(RSQLite::SQLite(), "data/derived/annotations.sqlite3")


# Load the SNP-clumped GWAS results, created with GEMMA and clumped with PLINK
xx <- bind_rows(
  read_tsv("data/derived/female_early_SNP_clumps.txt"),
  read_tsv("data/derived/female_late_SNP_clumps.txt"),
  read_tsv("data/derived/male_early_SNP_clumps.txt"),
  read_tsv("data/derived/male_late_SNP_clumps.txt")) %>%
  distinct(index_SNP, .keep_all = TRUE) %>%
  mutate(FBIDs = str_remove_all(FBIDs, ", NA"),
         FBIDs = replace(FBIDs, is.na(FBIDs), "")) %>%
  mutate(genes = str_remove_all(genes, ", NA"),
         genes = replace(genes, is.na(genes), "")) %>%
  mutate(`Index variant log10 p` = round(-1 * log10(index_SNP_p_value), 2)) %>%
  select(index_SNP, other_snps, FBIDs, genes, `Index variant log10 p`) 

# Because we clumped separately for the 4 traits, some of the "index SNPs" appear inside the ranges
# tagged by other index SNPs, so there is some redundancy. The following ugly code sorts this out 
# by putting index SNPs that appear inside the range of another index SNP into the 'other SNPs' 
# column for that index SNP
ranges <- xx$other_snps %>% 
  str_extract_all("[_][:digit:]+") %>% 
  map(~ as.numeric(str_remove_all(.x, "_"))) %>% 
    map(~ c(min(.x), max(.x)))

ranges <- map_df(1:nrow(xx), 
                 ~ data.frame(index = xx$index_SNP[.x], 
                              min =ranges[[.x]][1], 
                              max = ranges[[.x]][2]))%>%
  filter(!is.na(min)) %>% 
  mutate(chr = substr(index, 1,2))

to_modify <- lapply(1:nrow(xx), function(i){
  chr <- substr(xx$index_SNP[i],1,2) # get chr and pos
  pos <- as.numeric(str_remove_all(str_extract_all(xx$index_SNP[i], "_[:digit:]+"), "_"))
  foc_ranges <- ranges[ranges$chr == chr & 
                         ranges$min <= pos & 
                         ranges$max >= pos &
                         ranges$index != xx$index_SNP[i], ]
  if(nrow(foc_ranges) > 0) {
    return(data.frame(index_SNP = foc_ranges$index, 
                      to_add = xx$index_SNP[i]))
  }
  return(NULL)
}) %>% bind_rows()

xx <- xx %>% 
  filter(!(index_SNP %in% to_modify$to_add)) %>% 
  left_join(to_modify) %>% 
  mutate(to_add = replace(to_add, is.na(to_add), "NA"))

xx <- xx %>% 
  mutate(other_snps = str_split(str_c(xx$other_snps, to_add, sep = ", "), ", "),
              other_snps = map_chr(other_snps, 
                                   ~ paste0(.x[.x != "NA"], collapse = ", "))) %>%
  select(-to_add) %>%
  mutate(other_snps = str_replace_all(other_snps, "NA", ""))

# Now get the mashr-ED adjusted effect sizes for the index SNPs, 
# and their mixture assignment probs.
univariate_lmm_results <- 
  tbl(db, "univariate_lmm_results") %>%
  filter(SNP %in% !! xx$index_SNP) %>%
  select(-contains("canonical"), 
         -contains("raw")) %>% 
  inner_join(tbl(db, "variants") %>% 
               select(-chr, -position), 
             by = "SNP") %>%
  left_join(
    tbl(db, "genes") %>% 
      select(FBID, gene_name), by = "FBID") %>%
  collect(n=Inf) %>%
  arrange(LFSR_female_early_mashr_ED) %>%
  mutate_at(vars(contains("LFSR")), ~ -log10(.x)) %>%
  mutate_if(is.numeric, ~ round(.x, 2)) %>%
  rename_all(~ gsub("beta_", "", .x)) %>%
  rename_all(~ gsub("_mashr_ED", "", .x)) %>%
  arrange(-P_equal_effects) %>%
  select(SNP, starts_with("female"), 
         starts_with("male"), starts_with("P_")) %>%
  distinct()

GWAS_table <- xx %>%
  left_join(univariate_lmm_results, by = c("index_SNP" = "SNP")) %>%
  arrange(index_SNP) %>%
  rename(`Index variant` = index_SNP,
         `Other variants` = other_snps,
         `Genes` = genes, 
         `Female early` = female_early,
         `Female late` = female_late,
         `Male early` = male_early,
         `Male late` = male_late) %>%
  rename_all(~ str_replace_all(.x, "_", " "))

GWAS_table %>%
  my_data_table()