## ----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, weight_missing = 0.5)
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

## ----computeContrasts---------------------------------------------------------
contr_spec <- c(
  "AvsC" = "group_A - group_Ctrl",
  "BvsC" = "group_B - group_Ctrl"
)
fa <- build_contrast_analysis(transformed, "~ group_", contr_spec, method = "lm_impute")
contrdf <- fa$get_contrasts()

## ----estimateType-------------------------------------------------------------
dplyr::count(contrdf, contrast, estimate_type)

## ----plotVolcanos, fig.cap="Volcano plots of the two contrasts."--------------
plotter <- fa$get_Plotter()
plotter$volcano()$FDR

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

## ----createSigHeatmap---------------------------------------------------------
sig <- dplyr::filter(contrdf, FDR < 0.1)
signif <- transformed$get_copy()
signif$set_data(signif$data_long() |> dplyr::filter(protein_Id %in% sig$protein_Id))
hmSig <- signif$get_Plotter()$heatmap()

## ----showSigHeatmap, fig.cap="Heatmap of the proteins with FDR < 0.1 in at least one contrast"----
hmSig

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

## ----missingContrasts---------------------------------------------------------
fa$get_missing()

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