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

## ----createConfig-------------------------------------------------------------
library(prolfqua)

## ----LoadDataAndConfigure-----------------------------------------------------
xx <- prolfqua::sim_lfq_data_protein_config(Nprot = 100)
xx

## ----removeSmallIntensities---------------------------------------------------
lfqdata <- prolfqua::LFQData$new(xx$data, xx$config)
lfqdata$remove_small_intensities()
lfqdata$factors()

## ----showWide, eval = TRUE----------------------------------------------------
lfqdata$data_wide()$data[1:3,1:7]

## ----getPlotter---------------------------------------------------------------
lfqplotter <- lfqdata$get_Plotter()
density_nn <- lfqplotter$intensity_distribution_density()

## ----makeMissingHeatmap, fig.cap="Heatmap of missing protein intensities, black - missing, white - present"----
lfqplotter$na_heatmap()

## ----missignessPerGroup, fig.cap="# of proteins with 0,1,...N missing values"----
lfqdata$get_Summariser()$plot_missingness_per_group()

## ----missignessHistogram, fig.cap="Intensity distribution of proteins depending on # of missing values"----
lfqplotter$missigness_histogram()

## ----PlotCVDistributions, fig.cap="Violin plots of CVs in the different groups and among all groups"----
stats <- lfqdata$get_Stats()
stats$violin()
prolfqua::table_facade( stats$stats_quantiles()$wide, paste0("quantile of ",stats$stat ))

## ----plotCVsplitbyMedianIntensity, fig.cap="Distribution of CV's for top 50% and bottom 50% proteins by intensity."----
stats$density_median()


## ----normalizedata------------------------------------------------------------
lt <- lfqdata$get_Transformer()
transformed <- lt$log2()$robscale()$lfq
transformed$is_transformed()


## ----genplotNorm, fig.cap="Normalized intensities."---------------------------
pl <- transformed$get_Plotter()
density_norm <- pl$intensity_distribution_density()

## ----showIntensityDistributions, fig.cap="Distribution of intensities before and after normalization."----
gridExtra::grid.arrange(density_nn, density_norm)

## ----plotScatterMatrix, fig.cap = "Scatterplot matrix"------------------------
pl$pairs_smooth()

## ----createHeatmap------------------------------------------------------------
p <- pl$heatmap_cor()

## ----plotHeatmap, fig.cap="Heatmap, Rows - proteins, Columns - samples", fig.align=5, fig.height=5----
p

## ----lookatfactors------------------------------------------------------------
transformed$factors()


## ----defineModelAndContrasts--------------------------------------------------
formula_Condition <-  strategy_lm("transformedIntensity ~ group_")

# specify model definition
modelName  <- "Model"
contr_spec <- c("AvsC" = "group_A - group_Ctrl",
"BvsC" = "group_B - group_Ctrl")


## ----buildModel---------------------------------------------------------------
mod <- prolfqua::build_model(
  transformed$data_long(),
  formula_Condition,
  subject_id = transformed$hierarchy_keys() )


## ----showANOVA, fig.cap="Distribution of adjusted p-values (FDR)"-------------
mod$anova_histogram("FDR")

## ----filterDataForFDR---------------------------------------------------------
aovtable <- mod$get_anova()
head(aovtable)
dim(aovtable)
xx <- aovtable |> dplyr::filter(FDR < 0.2)
signif <- transformed$get_copy()
signif$set_data(signif$data_long() |> dplyr::filter(protein_Id %in% xx$protein_Id))
hmSig <- signif$get_Plotter()$heatmap()


## ----showSigHeatmap, fig.cap="Heatmap for proteins with FDR < 0.2 in the analysis of variance"----
hmSig


## ----computeContrasts---------------------------------------------------------
contr <- prolfqua::Contrasts$new(mod, contr_spec)
v1 <- contr$get_Plotter()$volcano()

## ----computeContrastsModerated------------------------------------------------
contr <- prolfqua::ContrastsModerated$new(contr)
contrdf <- contr$get_contrasts()

## ----plotVolcanos, fig.cap="Volcano plot, Left panel - no moderation, Right panel - with moderation."----
plotter <- contr$get_Plotter()
v2 <- plotter$volcano()
gridExtra::grid.arrange(v1$FDR,v2$FDR, ncol = 1)

## ----showMAplot, fig.cap="MA plot showing the dependency of mean abundance with respect to the difference"----
plotter$ma_plotly()

## ----checkProteinsInOutput----------------------------------------------------
myProteinIDS <- head(unique(contrdf$protein_Id), 3)
dplyr::filter(contrdf, protein_Id %in% myProteinIDS)

## ----computeMissing, eval=TRUE------------------------------------------------
mC <- prolfqua::ContrastsLMImputeFacade$new(transformed, "~ group_", contr_spec)$contrast
colnames(mC$get_contrasts())


## ----mergeResults, fig.cap="Volcano plots for the two contrasts with missing value imputation from the lm_impute model."----

merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC)$merged
plotter <- merged$get_Plotter()
tmp <- plotter$volcano()
tmp$FDR


## ----mergedMore---------------------------------------------------------------
merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC)

moreProt <- transformed$get_copy()
moreProt$set_data(moreProt$data_long() |> dplyr::filter(protein_Id %in% merged$more$contrast_result$protein_Id))
moreProt$get_Plotter()$raster()

# here we do not get anything because there is nothing imputed!

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