## ----load data, include = FALSE-----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  comment = "#>",
  fig.width = 10
)

## ----install bioconductor, eval = FALSE---------------------------------------
# if (!require("BiocManager", quietly = TRUE)){
#   install.packages("BiocManager")
# }
# 
# BiocManager::install("sigvar")

## ----install github, eval=FALSE-----------------------------------------------
# if (!require("devtools", quietly = TRUE)){
#   install.packages("devtools")
# }
# 
# devtools::install_github("MaikeMorrison/sigvar",
#                          dependencies = TRUE,
#                          build_vignettes = TRUE)

## ----load vignette, eval = FALSE----------------------------------------------
# vignette("sigvar_tutorial", package = "sigvar")
# 
# # or
# 
# browseVignettes("sigvar")

## ----data schematic, out.width="500%", echo = FALSE---------------------------
knitr::include_graphics("../man/figures/schematic_data_structure.png")

## ----view cosmic SBS, echo = FALSE--------------------------------------------
library(dplyr)
data(COSMIC3.3.1_SBS, package = "sigvar")

knitr::kable(COSMIC3.3.1_SBS %>% mutate(across(SBS1:SBS95, function(col) round(col, 3)))) %>%
  kableExtra::scroll_box(width = "600px", height = "300px")

## ----check cossim matrix------------------------------------------------------
# Exclude the first column, "Type", which defines each  mutational type
cosmic_sbs_sim <- sigvar::cossim(as.matrix(COSMIC3.3.1_SBS[, -1]))

# View the first 10 rows and 5 columns of the similarity matrix
cosmic_sbs_sim[1:10, 1:5]

# Confirm that the diagonal elements of the matrix are all equal to 1
all(diag(cosmic_sbs_sim) == 1)

# Confirm that all elements are non-negative
all(cosmic_sbs_sim >= 0)

# Confirm that no elements exceed 1
all(cosmic_sbs_sim <= 1)

## ----demonstrate reordering cossim matrix-------------------------------------
# Suppose this is the list of signatures included
# in your relative activity matrix:
sigs <- c("SBS95", "SBS10a", "SBS5", "SBS1", "SBS17b")

# Re-order the rows and columns of the similarity matrix
sim_reordered <- cosmic_sbs_sim[sigs, sigs]

sim_reordered

## ----setup--------------------------------------------------------------------
library(sigvar)
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)

# load data from the sigvar R package
data(ESCC_sig_activity, package = "sigvar")
data(ESCC_sig_similarity, package = "sigvar")

## ----explore data, eval=FALSE-------------------------------------------------
# # open the data set in a new window
# View(ESCC_sig_activity)
# 
# # view the structure of the data set
# str(ESCC_sig_activity)

## ----show ESCC_sig_activity, echo=FALSE---------------------------------------
knitr::kable(ESCC_sig_activity[1:40, 1:20]) %>%
  kableExtra::scroll_box(width = "600px", height = "300px")

## ----show sim matrix for ESCC, echo = FALSE-----------------------------------
knitr::kable(ESCC_sig_similarity[1:20, 1:20]) %>%
  kableExtra::scroll_box(width = "600px", height = "300px")

## ----cossim heat map, fig.height = 7, echo = FALSE----------------------------
ggplot(ESCC_sig_similarity %>%
  data.frame() %>%
  mutate(name2 = rownames(ESCC_sig_similarity)) %>%
  pivot_longer(
    cols = 1:43,
    values_to = "Similarity"
  ) %>%
  mutate(across(
    c(name, name2),
    function(col) {
      factor(col,
        ordered = TRUE,
        levels = colnames(ESCC_sig_similarity)
      )
    }
  ))) +
  geom_raster(aes(
    x = name,
    y = name2,
    fill = Similarity
  )) +
  theme_minimal() +
  scale_fill_viridis_c() +
  theme(
    axis.text.y = element_text(size = 6),
    axis.text.x = element_text(
      size = 6,
      angle = -90,
      hjust = 0
    ),
    axis.title = element_blank()
  )

## ----plot sig activity--------------------------------------------------------
data(all_sig_pal, package = "sigvar")

plot_signature_prop(
  relab_matrix = ESCC_sig_activity,
  K = 43,
  group = "Country",
  arrange = TRUE
) +
  scale_color_manual(values = all_sig_pal) +
  scale_fill_manual(values = all_sig_pal)

# all_sig_pal is a vector of hex color codes, each corresponding to
# one signature used in our paper.
# Below we show the first 6 entries as an example.
all_sig_pal[1:6]

## ----dot plots, fig.width=8---------------------------------------------------
plot_dots(
  sig_activity = ESCC_sig_activity,
  K = 43, group = "Country",
  facet = "Incidence_Level",
  pivot = TRUE
)

## ----compute sigvar-----------------------------------------------------------
sigvar <- sigvar(
  sig_activity = ESCC_sig_activity,
  K = 43,
  S = ESCC_sig_similarity,
  group = "Country"
)

knitr::kable(sigvar)

## ----sigvar by country--------------------------------------------------------
sigvar_incidence <- ESCC_sig_activity %>%
  transmute(
    Country = as.character(Country),
    Incidence_Level
  ) %>%
  distinct() %>%
  right_join(sigvar)

knitr::kable(sigvar_incidence)

ggplot(
  sigvar_incidence,
  aes(
    x = mean_within_sample_diversity,
    y = across_sample_heterogeneity,
    color = Incidence_Level
  )
) +
  geom_point(size = 4) +
  ggrepel::geom_text_repel(aes(label = Country)) +
  theme_bw()

ggplot(
  sigvar_incidence %>%
    pivot_longer(
      cols = c(
        "across_sample_heterogeneity",
        "mean_within_sample_diversity"
      ),
      names_to = "Statistic",
      values_to = "sigvar"
    ),
  aes(x = Incidence_Level, y = sigvar)
) +
  geom_boxplot(aes(color = Incidence_Level)) +
  geom_point(aes(color = Incidence_Level), size = 4, alpha = 0.6) +
  facet_wrap(~Statistic, scales = "free_y") +
  ggpubr::stat_compare_means(method = "t.test") +
  theme_bw()

## ----bootstrapping------------------------------------------------------------
comparison_boot <- sigboot(
  sig_activity = ESCC_sig_activity %>%
    filter(Country %in% c("Japan", "Kenya")),
  K = 43,
  group = "Country",
  S = ESCC_sig_similarity,
  n_replicates = 500,
  seed = 1
)

comparison_boot$bootstrap_distribution_plot

## ----view p values------------------------------------------------------------
comparison_boot$P_values

## ----load TP53 spectrum-------------------------------------------------------
ref_genome <- "BSgenome.Hsapiens.UCSC.hg38"
library(ref_genome, character.only = TRUE)
TP53_spectrum <- get_SBS96_spectrum(transcript = "ENST00000269305.9",ref_genome=ref_genome)

# look at first 10 entries
TP53_spectrum[1:10]

## ----TP53 drivers-------------------------------------------------------------
data(TP53_drivers_intogen_LUAD, package = "sigvar")
ref_genome <- "BSgenome.Hsapiens.UCSC.hg38"
library(ref_genome, character.only = TRUE)
TP53_driver_spectrum <- get_SBS96_driver_spectrum(TP53_drivers_intogen_LUAD,ref_genome = ref_genome)

# look at first 10 entries
TP53_driver_spectrum[1:10]

## ----plot TP53 SBS spectrum---------------------------------------------------
plot_SBS_spectrum(data.frame(TP53_driver_spectrum) %>% 
                    `colnames<-`(c("ref", "TP53")) %>% 
                    select(TP53))

## ----ratio of spectra---------------------------------------------------------
TP53_driver_prop <- TP53_driver_spectrum / TP53_spectrum

TP53_driver_prop[1:10]

## ----session info-------------------------------------------------------------
sessionInfo()