## ----my_setup, echo=FALSE, warning=FALSE, include=FALSE-----------------------
knitr::opts_chunk$set(echo = TRUE, fig.align="center")
devtools::load_all()

## ----logo, echo=FALSE, warning=FALSE------------------------------------------
htmltools::img(
  src = knitr::image_uri("MiDAS_logo.png"),
  alt = 'logo',
  style = 'position:absolute; top:0; right:0; padding:10px; width:250px'
  )

## ----get_data, echo=TRUE, warning=FALSE---------------------------------------
dat_pheno <-
  read.table(
  file = system.file("extdata", "MiDAS_tut_pheno.txt", package = "midasHLA"),
  header = TRUE
  )
dat_HLA <-
  readHlaCalls(
  file = system.file("extdata", "MiDAS_tut_HLA.txt", package = "midasHLA"),
  resolution = 4
  )

## ----show_imported_data, echo=FALSE, warning=FALSE----------------------------
dat_HLA %>%
  head(20) %>%
  kable(caption = "HLA data as imported by MiDAS") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
  scroll_box(width = "100%", height = "300px")

## ----get_frequencies, echo=TRUE, warning=FALSE--------------------------------
freq_HLA <- getHlaFrequencies(hla_calls = dat_HLA, compare = TRUE) %>%
  filter(Freq > 0.01)

## ----show_frequencies, echo=FALSE, warning=FALSE------------------------------
freq_HLA %>%
  head() %>%
  kable(caption = "HLA frequencies, compared to published references") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

## ----plot_frequencies, echo=TRUE, warning=FALSE, fig.height=7, fig.width=7----
freq_HLA_long <- tidyr::gather(
  data = freq_HLA,
  key = "population",
  value = "freq",
  "Freq",
  "USA NMDP European Caucasian",
  "USA NMDP Chinese",
  "USA NMDP African American pop 2",
  factor_key = TRUE
) %>% 
  filter(grepl("^A", allele))

plot_HLAfreq <-
  ggplot2::ggplot(data = freq_HLA_long, ggplot2::aes(x = allele, y = freq, fill = population)) +
  ggplot2::geom_bar(
    stat = "identity",
    position = ggplot2::position_dodge(0.7),
    width = 0.7,
    colour = "black"
  ) +
  ggplot2::coord_flip() +
  ggplot2::scale_y_continuous(labels = formattable::percent)

plot_HLAfreq

## ----prep_HLA, echo=TRUE, warning=FALSE---------------------------------------
HLA <- prepareMiDAS(
  hla_calls = dat_HLA,
  colData = dat_pheno,
  experiment = "hla_alleles"
)

## ----HWE_filtering, echo=TRUE, warning=FALSE----------------------------------
HLA <- HWETest(
  object = HLA,
  experiment = "hla_alleles",
  HWE_cutoff = 0.05 / 447,
  as.MiDAS = TRUE
)

## ----run_HLA_allelic, echo=TRUE, warning=FALSE, message=FALSE-----------------
HLA_model <- glm(disease ~ term, data = HLA, family = binomial())
HLA_results <- runMiDAS(
  object = HLA_model, 
  experiment = "hla_alleles", 
  inheritance_model = "dominant",
  lower_frequency_cutoff = 0.02, 
  upper_frequency_cutoff = 0.98, 
  correction = "bonferroni", 
  exponentiate = TRUE
)

kableResults(HLA_results)

## ----run_HLA_allelic_cond, echo=TRUE, warning=FALSE, message=FALSE------------
HLA_results_cond <- runMiDAS(
  object = HLA_model, 
  experiment = "hla_alleles", 
  inheritance_model = "dominant", 
  conditional = TRUE,
  lower_frequency_cutoff = 0.02, 
  upper_frequency_cutoff = 0.98, 
  correction = "bonferroni", 
  exponentiate = TRUE
)

kableResults(HLA_results_cond, scroll_box_height = "200px")

## ----prepare_HLA_AA, echo=TRUE, warning=FALSE---------------------------------
HLA_AA <- prepareMiDAS(
  hla_calls = dat_HLA,
  colData = dat_pheno,
  experiment = "hla_aa"
)

## ----HLA_AA_to_df, echo=TRUE, warning=FALSE-----------------------------------
dat_HLA_AA <- HLA_AA[["hla_aa"]] %>% 
  assay() %>% 
  t() %>% 
  as.data.frame() %>% 
  select(starts_with("B_97_")) %>% 
  head()

## ----show_HLA_AA, echo=FALSE, warning=FALSE-----------------------------------
kable(dat_HLA_AA, caption = "HLA amino acid data as inferred by MiDAS") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

## ----run_HLA_AA_omnibus, echo=TRUE, warning=FALSE, message=FALSE--------------
HLA_AA_model <- glm(disease ~ term, data = HLA_AA, family = binomial())
HLA_AA_omnibus_results <- runMiDAS(
  HLA_AA_model,
  experiment = "hla_aa",
  inheritance_model = "dominant",
  conditional = FALSE,
  omnibus = TRUE,
  lower_frequency_cutoff = 0.02,
  upper_frequency_cutoff = 0.98,
  correction = "bonferroni"
)

kableResults(HLA_AA_omnibus_results)

## ----run_HLA_AA_DQB1_9, echo=TRUE, warning=FALSE------------------------------
HLA_AA_DQB1_9_results <- runMiDAS(
  HLA_AA_model,
  experiment = "hla_aa",
  inheritance_model = "dominant",
  omnibus_groups_filter = "DQB1_9",
  lower_frequency_cutoff = 0.02,
  upper_frequency_cutoff = 0.98,
  correction = "bonferroni",
  exponentiate = TRUE
)

kableResults(HLA_AA_DQB1_9_results, scroll_box_height = "250px")

## ----alleles_for_aa, echo=TRUE, warning=FALSE---------------------------------
HLA_AA_DQB1_9_alleles <- getAllelesForAA(HLA_AA,"DQB1_9")

## ----get_HLA_AA_DQB1_9_alleles, echo=FALSE, warning=FALSE---------------------
kable(HLA_AA_DQB1_9_alleles, escape = FALSE) %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

## ----prep_KIR_ligand, echo=TRUE, warning=FALSE--------------------------------
NKlig <- prepareMiDAS(
  hla_calls = dat_HLA,
  colData = dat_pheno,
  experiment = c("hla_alleles", "hla_NK_ligands")
)

## ----run_KIR_ligand, echo=TRUE, warning=FALSE, message=FALSE------------------
NKlig_model <- glm(disease ~ term, data = NKlig, family = binomial())
NKlig_results <- runMiDAS(
  NKlig_model,
  experiment = "hla_NK_ligands",
  inheritance_model = "dominant",
  correction = "bonferroni",
  exponentiate = TRUE
)

kableResults(NKlig_results)

## ----get_KIR, echo=TRUE, warning=FALSE----------------------------------------
dat_KIR <- readKirCalls(
  file = system.file("extdata", "MiDAS_tut_KIR.txt", package = "midasHLA")
)

## ----show_KIR, echo=FALSE, warning=FALSE--------------------------------------
kable(dat_KIR, caption = "KIR data as imported by MiDAS") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  scroll_box(height = "400px")

## ----KIR_frequency, echo=TRUE, warning=FALSE----------------------------------
kir_freq <- getKIRFrequencies(dat_KIR)

## ----show_KIR_frequency, echo=FALSE, warning=FALSE----------------------------
kable(kir_freq , caption = "KIR data as imported by MiDAS") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
  scroll_box(height = "400px")

## ----run_KIR_gene, echo=TRUE, warning=FALSE, message=FALSE--------------------
HLAKIR <- prepareMiDAS(
  hla_calls = dat_HLA,
  kir_calls = dat_KIR,
  colData = dat_pheno,
  experiment = c("hla_NK_ligands","kir_genes", "hla_kir_interactions")
)
KIR_model <- glm(disease ~ term, data = HLAKIR, family = binomial())
KIR_results <- runMiDAS(
  KIR_model,
  experiment = "kir_genes",
  lower_frequency_cutoff = 0.02,
  upper_frequency_cutoff = 0.98,
  exponentiate = TRUE
)

kableResults(KIR_results)

## ----HLA-KIR_interactions_frequency, echo=TRUE, warning=FALSE-----------------
hlakir_freq <- getFrequencies(HLAKIR, experiment = "hla_kir_interactions")

## ----show_HLA-KIR_interactions_frequency, echo=FALSE, warning=FALSE-----------
kable(hlakir_freq, caption = "HLA-KIR interaction frequencies") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
  scroll_box(height = "400px")

## ----run_HLA-KIR_interactions, echo=TRUE, warning=FALSE, message=FALSE--------
HLAKIR_model <- glm(disease ~ term, data = HLAKIR, family = binomial())
HLA_KIR_results <- runMiDAS(
  HLAKIR_model,
  experiment = "hla_kir_interactions",
  lower_frequency_cutoff = 0.02,
  upper_frequency_cutoff = 0.98,
  exponentiate = TRUE
)

kableResults(HLA_KIR_results)

## ----run_HLAhet, echo=TRUE, warning=FALSE, message=FALSE----------------------

HLA_het <- prepareMiDAS(
  hla_calls = dat_HLA, 
  colData = dat_pheno, 
  experiment = c("hla_het","hla_divergence")
)

HLA_het_model <- glm(outcome ~ term, data=HLA_het, family=binomial())

HLA_het_results <- runMiDAS(HLA_het_model, 
  experiment = "hla_het", 
  exponentiate = TRUE
)

kableResults(HLA_het_results)

HLA_div_results <- runMiDAS(HLA_het_model, 
  experiment = "hla_divergence", 
  exponentiate = TRUE
)

kableResults(HLA_div_results, scroll_box_height = "250px")