## ----setup, include=FALSE-----------------------------------------------------
library(quantbayes)
library(ggplot2)
library(dplyr)
library(tidyr)
library(purrr)

# widen default chunk output
knitr::opts_chunk$set(
fig.width = 8,
fig.height = 5,
dpi = 120,
out.width = "100%"
)

#CSS to prevent html_vignette from squashing figures
knitr::asis_output("
")

## -----------------------------------------------------------------------------
data(core_test_data)
head(core_test_data)

## -----------------------------------------------------------------------------
x <- as.matrix(core_test_data[, -1])
rownames(x) <- core_test_data[[1]]
dim(x)

## -----------------------------------------------------------------------------
res <- quant_es_core(x)

res$global
head(res$variants)

## -----------------------------------------------------------------------------
plots <- quant_es_plots(res, x)

## -----------------------------------------------------------------------------
plots$p_global

## -----------------------------------------------------------------------------
plots$p_overlay

## -----------------------------------------------------------------------------
plots$p_matrix

## -----------------------------------------------------------------------------
plots$p_p_hat

## -----------------------------------------------------------------------------
plots$p_theta_ci

## -----------------------------------------------------------------------------
plots$p_combined

## -----------------------------------------------------------------------------
highlight_demo <- list(
  list(id = "X-119469998-CT-C_XR", colour = "#ee4035", size = 4),
  list(id = res$variants$variant_id[1], colour = "#2f4356", size = 4)
)

plots2 <- quant_es_plots(res, x, highlight_points = highlight_demo)
plots2$p_overlay

## -----------------------------------------------------------------------------
# Extract evidence columns
ev <- core_test_data[, -1]

# Keep only columns containing valid 0, 1, NA entries
is_binary_col <- function(x) all(x %in% c(0, 1, NA))
ev <- ev[, vapply(ev, is_binary_col, logical(1))]

tmpfile <- tempfile(fileext = ".tsv")

write.table(
  ev,
  tmpfile,
  sep = "\t",
  quote = FALSE,
  row.names = FALSE
)

res_file <- quant_es_from_binary_table(tmpfile, variant_col = NULL)
res_file$global

## -----------------------------------------------------------------------------
outdir <- tempdir()
if (!dir.exists(outdir)) dir.create(outdir)

write.csv(res$variants,
          file.path(outdir, "variants_results.csv"),
          row.names = FALSE)

write.csv(as.data.frame(res$global),
          file.path(outdir, "global_summary.csv"),
          row.names = FALSE)

## -----------------------------------------------------------------------------
ggsave(file.path(outdir, "overlay.png"),
       plots$p_overlay,
       width = 7,
       height = 4,
       dpi = 120)

ggsave(file.path(outdir, "matrix.png"),
       plots$p_matrix,
       width = 7,
       height = 6,
       dpi = 120)

## -----------------------------------------------------------------------------
plots$p_overlay + theme(legend.position = "none")

## -----------------------------------------------------------------------------
plots$p_overlay + theme_minimal()

## -----------------------------------------------------------------------------
plots$p_overlay + theme_void()

## -----------------------------------------------------------------------------
pal10 <- colorRampPalette(c("black", "grey"))(10)
pal20 <- colorRampPalette(c("skyblue", "navy"))(20)

plots_custom <- quant_es_plots(
  res,
  x,
  palette10 = pal10,
  palette20 = pal20
)

plots_custom$p_overlay

## -----------------------------------------------------------------------------
# Two highlighted candidates
swiss_red <- "#ee4035"
federal_blue <- "#2f4356"

# provide the points to highlight
highlight_flagship <- list(
  list(id = core_test_data[[1]][1], colour = swiss_red, size = 4),
  list(id = "6-32664815-C-G_AR", colour = federal_blue, size = 4)
)

plots_flagship <- quant_es_plots(
  res,
  x,
  highlight_points = highlight_flagship
)

# add custom ggplot settings
flagship_plot <- plots_flagship$p_overlay +
  ggplot2::guides(
    fill = ggplot2::guide_legend(title = "highlighted variants"),
    size = "none"
  ) +
  ggplot2::labs(
    title = "Posterior theta distribution",
    subtitle = "Top 10 CrI estimates with evidence available\nand two highlighted variants"
  ) +
  ggplot2::theme(
    legend.position = "right",
    legend.title = ggplot2::element_text(size = 10),
    legend.text = ggplot2::element_text(size = 9)
  )

flagship_plot

# Save flagship plot
outdir <- tempdir()
if (!dir.exists(outdir)) dir.create(outdir)

ggplot2::ggsave(
  filename = file.path(outdir, "flagship_plot.pdf"),
  plot = flagship_plot,
  width = 8,
  height = 4
)

## -----------------------------------------------------------------------------
res_df <- as.data.frame(res$variants)
global_df <- as.data.frame(res$global)

res_df <- res_df[order(res_df$theta_mean, decreasing = TRUE), ]

head(res_df)

## -----------------------------------------------------------------------------
# choose the top variant
v <- res_df[1, ]
g <- global_df

variant_line <- sprintf(
  "Posterior evidence sufficiency: %.3f (credible interval %.3f to %.3f, percentile %.2f)",
  v$theta_mean,
  v$theta_lower,
  v$theta_upper,
  v$percentile
)

global_line <- sprintf(
  "Global evidence sufficiency: %.2f (credible interval %.2f to %.2f)",
  g$mean_theta,
  g$lower_theta,
  g$upper_theta
)

variant_line
global_line