## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
library(prolfqua)
library(dplyr)
library(ggplot2)

## ----simulate_two_factor_subgroups--------------------------------------------
dd <- sim_lfq_data_2factor_config(
  Nprot = 80,
  with_missing = TRUE,
  weight_missing = 1,
  PEPTIDE = TRUE,
  seed = 7
)

cfg_subgroup <- dd$config$clone(deep = TRUE)
data_subgroup <- dd$data |>
  dplyr::mutate(
    Condition = dplyr::recode(Treatment, A = "C", B = "T"),
    Genotype = dplyr::recode(Background, Z = "WT", X = "KO"),
    Subgroup = paste(Condition, Genotype, sep = "_")
  )

cfg_subgroup$factors <- c(Subgroup = "Subgroup")
cfg_subgroup$factor_depth <- 1

lfq_peptide_2f <- LFQData$new(data_subgroup, cfg_subgroup)
lfq_peptide_2f <- lfq_peptide_2f$get_Transformer()$log2()$lfq

lfq_protein_2f <- lfq_peptide_2f$get_Aggregator()$aggregate()

lfq_peptide_2f$hierarchy_keys()
lfq_protein_2f$hierarchy_keys()
lfq_protein_2f$nr_children_col()

## ----contrasts_two_factor-----------------------------------------------------
contrasts_2f <- c(
  T_C_gv_WT = "SubgroupT_WT - SubgroupC_WT",
  T_C_gv_KO = "SubgroupT_KO - SubgroupC_KO",
  WT_KO_comp = "T_C_gv_WT - T_C_gv_KO"
)

## ----facades_protein_two_factor-----------------------------------------------
fa_lm_2f <- build_contrast_analysis(
  lfq_protein_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "lm"
)

fa_limma_2f <- build_contrast_analysis(
  lfq_protein_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "limma"
)

fa_limma_impute_2f <- build_contrast_analysis(
  lfq_protein_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "limma_impute"
)

fa_deqms_2f <- build_contrast_analysis(
  lfq_protein_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "deqms"
)

fa_lm_impute_2f <- build_contrast_analysis(
  lfq_protein_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "lm_impute"
)

fa_firth_protein_2f <- build_contrast_analysis(
  lfq_protein_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "firth"
)

## ----collect_protein_two_factor-----------------------------------------------
# Proteins missing in the baseline lm facade (used to flag rescued proteins)
lm_missing_ids_2f <- fa_lm_2f$get_missing() |>
  dplyr::select(protein_Id, contrast) |>
  dplyr::mutate(rescued = TRUE)

results_protein_2f <- bind_rows(
  fa_lm_2f$get_contrasts(),
  fa_limma_2f$get_contrasts(),
  fa_limma_impute_2f$get_contrasts(),
  fa_lm_impute_2f$get_contrasts(),
  fa_deqms_2f$get_contrasts(),
  fa_firth_protein_2f$get_contrasts()
) |>
  dplyr::select(dplyr::any_of(c(
    "modelName", "estimate_type", "protein_Id", "contrast", "avgAbd", "diff", "FDR",
    "statistic", "std.error", "df", "p.value", "conf.low", "conf.high",
    "sigma"
  ))) |>
  dplyr::left_join(lm_missing_ids_2f, by = c("protein_Id", "contrast")) |>
  dplyr::mutate(
    rescued = dplyr::coalesce(rescued, FALSE),
    significant = FDR < 0.1 & abs(diff) > 0.5
  )

results_protein_2f |>
  dplyr::count(modelName, name = "n_results")

## ----collect_protein_model_names_2f-------------------------------------------
results_protein_2f |>
  dplyr::count(modelName, contrast, estimate_type, name = "n_results")

## ----volcano_standard_2f, fig.height=10, fig.width=12, fig.cap="Volcano plots for the standard protein-level facades (lm, limma, deqms)."----
standard_facades_2f <- c("lm", "limma", "deqms")
results_standard_2f <- results_protein_2f |>
  dplyr::filter(modelName %in% standard_facades_2f)

ggplot(results_standard_2f, aes(x = diff, y = -log10(p.value), color = significant)) +
  geom_point(alpha = 0.6, size = 1.5) +
  facet_grid(contrast ~ modelName, scales = "free_y") +
  geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "grey60") +
  geom_hline(yintercept = -log10(0.1), linetype = "dashed", color = "grey60") +
  scale_color_manual(values = c(`TRUE` = "firebrick", `FALSE` = "grey70")) +
  labs(x = "log2 fold change", y = "-log10(p.value)", color = "FDR < 0.1\nand |diff| > 0.5") +
  theme_minimal(base_size = 12)

## ----volcano_impute_2f, fig.height=10, fig.width=14, fig.cap="Volcano plots for the imputation and missingness facades. Large red diamonds mark proteins rescued (missing in plain lm)."----
impute_facades_2f <- c("lm_impute", "limma_impute", "firth")
results_impute_2f <- results_protein_2f |>
  dplyr::filter(modelName %in% impute_facades_2f) |>
  dplyr::mutate(modelName = factor(modelName, levels = impute_facades_2f))

ggplot(results_impute_2f, aes(x = diff, y = -log10(p.value))) +
  geom_point(data = dplyr::filter(results_impute_2f, !rescued),
             aes(color = significant), alpha = 0.5, size = 1.2) +
  geom_point(data = dplyr::filter(results_impute_2f, rescued),
             color = "red", shape = 18, size = 5, alpha = 1) +
  facet_grid(contrast ~ modelName, scales = "free_y") +
  geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "grey60") +
  geom_hline(yintercept = -log10(0.1), linetype = "dashed", color = "grey60") +
  scale_color_manual(values = c(`TRUE` = "firebrick", `FALSE` = "grey70")) +
  labs(x = "log2 fold change", y = "-log10(p.value)", color = "FDR < 0.1\nand |diff| > 0.5") +
  theme_minimal(base_size = 12)

## ----top_hits_protein_2f------------------------------------------------------
results_protein_2f |>
  dplyr::group_by(modelName, contrast) |>
  dplyr::slice_min(order_by = p.value, n = 5, with_ties = FALSE) |>
  dplyr::ungroup() |>
  dplyr::select(modelName, estimate_type, contrast, protein_Id, diff, p.value, FDR)

## ----missing_protein_facades_2f-----------------------------------------------
missing_all_2f <- dplyr::bind_rows(
  fa_lm_2f$get_missing() |> dplyr::mutate(facade = "lm"),
  fa_limma_2f$get_missing() |> dplyr::mutate(facade = "limma"),
  fa_limma_impute_2f$get_missing() |> dplyr::mutate(facade = "limma_impute"),
  fa_lm_impute_2f$get_missing() |> dplyr::mutate(facade = "lm_impute"),
  fa_deqms_2f$get_missing() |> dplyr::mutate(facade = "deqms"),
  fa_firth_protein_2f$get_missing() |> dplyr::mutate(facade = "firth")
)

missing_all_2f |>
  dplyr::count(facade, contrast, name = "n_missing") |>
  knitr::kable(caption = "Number of missing protein x contrast pairs per facade")

## ----missing_intensities_2f---------------------------------------------------
missing_proteins_2f <- unique(missing_all_2f$protein_Id)

if (length(missing_proteins_2f) > 0) {
  lfq_protein_2f$data_long() |>
    dplyr::filter(protein_Id %in% missing_proteins_2f) |>
    dplyr::select(protein_Id, sampleName,
                  !!rlang::sym(lfq_protein_2f$response())) |>
    tidyr::pivot_wider(names_from = sampleName,
                       values_from = !!rlang::sym(lfq_protein_2f$response())) |>
    knitr::kable(digits = 2, caption = "Per-sample intensities of proteins that could not be estimated")
}

## ----rescued_estimates_2f-----------------------------------------------------
lm_missing_proteins_2f <- fa_lm_2f$get_missing()$protein_Id |> unique()

if (length(lm_missing_proteins_2f) > 0) {
  rescued_2f <- results_protein_2f |>
    dplyr::filter(
      protein_Id %in% lm_missing_proteins_2f,
      modelName %in% c("lm_impute", "limma_impute")
    ) |>
    dplyr::arrange(protein_Id, contrast, modelName)

  rescued_2f |>
    knitr::kable(
      digits = 3,
      caption = "Contrast estimates from lm_impute and limma_impute for proteins that plain lm could not estimate"
    )
} else {
  cat("All proteins were estimable by plain lm — no rescued estimates to show.")
}

## ----facades_peptide_two_factor-----------------------------------------------
fa_lmer_2f <- build_contrast_analysis(
  lfq_peptide_2f,
  "~ Subgroup + (1 | peptide_Id) + (1 | sampleName)",
  contrasts_2f,
  method = "lmer_nested"
)

fa_ropeca_2f <- build_contrast_analysis(
  lfq_peptide_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "ropeca_nested"
)

fa_firth_peptide_2f <- build_contrast_analysis(
  lfq_peptide_2f,
  "~ Subgroup",
  contrasts_2f,
  method = "firth_nested"
)

## ----collect_peptide_two_factor-----------------------------------------------
results_peptide_2f <- bind_rows(
  fa_lmer_2f$get_contrasts(),
  fa_ropeca_2f$get_contrasts(),
  fa_firth_peptide_2f$get_contrasts()
) |>
  dplyr::select(dplyr::any_of(c(
    "modelName", "estimate_type", "protein_Id", "contrast", "avgAbd", "diff", "FDR",
    "statistic", "std.error", "df", "p.value", "conf.low", "conf.high",
    "sigma"
  ))) |>
  dplyr::mutate(
    significant = FDR < 0.1 & abs(diff) > 0.5
  )

results_peptide_2f |>
  dplyr::count(modelName, name = "n_results")

## ----volcano_peptide_2f, fig.height=10, fig.width=10, fig.cap="Volcano plots for the peptide-level facades. Rows are contrasts, columns are backends."----
ggplot(results_peptide_2f, aes(x = diff, y = -log10(p.value), color = significant)) +
  geom_point(alpha = 0.6, size = 1.2) +
  facet_grid(contrast ~ modelName, scales = "free_y") +
  geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "grey60") +
  geom_hline(yintercept = -log10(0.1), linetype = "dashed", color = "grey60") +
  scale_color_manual(values = c(`TRUE` = "firebrick", `FALSE` = "grey70")) +
  labs(
    x = "log2 fold change",
    y = "-log10(p.value)",
    color = "FDR < 0.1\nand |diff| > 0.5"
  ) +
  theme_minimal(base_size = 12)

## ----session_info-------------------------------------------------------------
sessionInfo()