## ----setup-knitr, include = FALSE---------------------------------------------
knitr::opts_chunk$set(comment = "#>", fig.width = 6.5, fig.height = 4,
                      dpi = 100, message = FALSE, warning = FALSE)

library(dplyr)
library(tidyr)
library(ggplot2)
library(MRManalyzeR)

## ----dataset------------------------------------------------------------------
xlsx  <- system.file("extdata", "example_data.xlsx",
                     package = "MRManalyzeR")
fdata <- openxlsx::read.xlsx(xlsx, sheet = "feature_metadata")
meta  <- openxlsx::read.xlsx(xlsx, sheet = "sample_metadata")

c(injections = nrow(meta), compounds = nrow(fdata))

## ----dataset-samples----------------------------------------------------------
DT::datatable(as.data.frame(table(Sample_type = meta$Sample_type)),
   rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

## ----dataset-features---------------------------------------------------------
DT::datatable(fdata[fdata$Report == "YES",
         c("Compound", "Enzymatic_pathway", "PUFA", "RT", "Report")],
   rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

## ----read---------------------------------------------------------------------
m <- readTargetLynx(xlsx, datatype = "Area", snr = 3)
dim(m)

## ----read-table---------------------------------------------------------------
DT::datatable(round(m[, 1:6], 1), filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

## ----dims-assembled-----------------------------------------------------------
de_real <- assembleDataset(m, fdata, meta)
dim(as.data.frame(de_real$data))

## ----dims-names---------------------------------------------------------------
rep18 <- fdata[fdata$Report == "YES", c("Processing_name", "Compound")]
rep18[rep18$Processing_name != rep18$Compound, ]

## ----syn-data-----------------------------------------------------------------
sxlsx <- system.file("extdata", "example_synthetic.xlsx",
                     package = "MRManalyzeR")
sfeat <- openxlsx::read.xlsx(sxlsx, sheet = "feature_metadata")
smeta <- openxlsx::read.xlsx(sxlsx, sheet = "sample_metadata")

table(smeta$Chrom_Batch, smeta$Sample_type)

## ----step-assemble------------------------------------------------------------
sm <- readSkyline(sxlsx, "skyline_data", sfeat, signal_filter = "LOD")
de <- assembleDataset(sm, sfeat, smeta)
de

## ----step-blanks-demo---------------------------------------------------------
bl <- assembleDataset(
  data.frame(well_above = c(100, 120, 110, 10, 12),
             near_blank = c(30, 28, 26, 25, 24),
             row.names  = c("S1", "S2", "S3", "B1", "B2")),
  data.frame(Compound = c("well_above", "near_blank"), Report = "YES"),
  data.frame(Name        = c("S1", "S2", "S3", "B1", "B2"),
             Sample_type = c("Sample", "Sample", "Sample", "Blank", "Blank")))

before <- as.data.frame(bl$data)
bl     <- filterBlanks(bl, blank_filter = 3, blank_head = "Sample_type",
                       blank_name = "Blank")
cbind(before, after = as.data.frame(bl$data))

## ----step-blanks-record-------------------------------------------------------
as.data.frame(bl$variable_meta)[, c("blank_threshold", "n_masked_blank")]

## ----step-blanks--------------------------------------------------------------
blank_rows <- rownames(de$data)[
  as.data.frame(de$sample_meta)$Sample_type == "Blank"]
c(missing_total = sum(is.na(de$data)),
  missing_in_blanks = sum(is.na(as.data.frame(de$data)[blank_rows, ])),
  blank_cells = length(blank_rows) * ncol(de$data))

## ----step-filter-demo---------------------------------------------------------
demo <- assembleDataset(
  data.frame(everywhere    = c(10, 12, 11, 9, 10, 11),
             group_specific = c(NA, NA, NA, 20, 22, 21),
             sparse         = c(5, NA, NA, NA, NA, 6),
             row.names = paste0("S", 1:6)),
  # assembleDataset() keeps only rows flagged Report == "YES"
  data.frame(Compound = c("everywhere", "group_specific", "sparse"),
             Report   = "YES"),
  data.frame(Name        = paste0("S", 1:6),
             Sample_type = rep("Sample", 6),
             Treatment   = rep(c("PBS", "HDM"), each = 3)))

across <- filterFeatures(demo, min_frac = 0.6, method = "across")
within <- filterFeatures(demo, min_frac = 0.6, method = "within",
                         group_col = "Treatment")

list(across = colnames(across$data), within = colnames(within$data))

## ----step-filter-features-----------------------------------------------------
de <- filterFeatures(de, min_frac = 0.5, method = "within",
                     group_col = "Treatment")

## ----step-filter-samples------------------------------------------------------
de <- filterSamples(de, max_na = 0.8)
summary(as.data.frame(de$sample_meta)$na_frac)

## ----step-batch---------------------------------------------------------------
batch     <- stats::setNames(de$sample_meta$Chrom_Batch,
                             rownames(as.data.frame(de$data)))
de_uncorr <- de

de <- correctBatch(de, qc_label = "QC",
                    factor_name = "Sample_type",
                    batch_head  = "Chrom_Batch")
de_bc     <- de

## ----step-batch-medians-------------------------------------------------------
qc     <- as.data.frame(de_uncorr$sample_meta)$Sample_type == "QC"
j      <- which.min(colSums(is.na(as.data.frame(de_uncorr$data))))
zeroed <- function(v) median(replace(v, is.na(v), 0))

rbind(
  before = tapply(as.data.frame(de_uncorr$data)[qc, j], batch[qc], zeroed),
  after  = tapply(as.data.frame(de_bc$data)[qc, j],     batch[qc], zeroed))

## ----step-batch-medians-all---------------------------------------------------
rbind(
  before = tapply(as.data.frame(de_uncorr$data)[[j]], batch, zeroed),
  after  = tapply(as.data.frame(de_bc$data)[[j]],     batch, zeroed))

## ----step-batch-failed--------------------------------------------------------
vm <- as.data.frame(de_bc$variable_meta)
vm[nzchar(vm$bc_failed_batches), c("Compound", "bc_failed_batches")]

## ----step-normalise-----------------------------------------------------------
de <- normaliseMatrix(de, column = "protein_ug")
round(as.data.frame(de$data)[1:3, 1:3], 3)

## ----step-adjust--------------------------------------------------------------
de <- adjustConcentration(de,
                           sample_vol_col   = "sample_volume_uL",
                           cal_vol_col      = "cal_vol_uL",
                           sample_IS_col    = "sample_IS_vol_uL",
                           cal_IS_col       = "cal_IS_vol_uL",
                           starting_vol_col = "starting_vol_uL")
round(as.data.frame(de$data)[1:3, 1:3], 5)

## ----step-impute--------------------------------------------------------------
de <- imputeMissing(de, scalar = 0.2, blank_head = "Sample_type",
                     blank_name = "Blank")
sum(is.na(de$data))

## ----step-impute-record-------------------------------------------------------
vm <- as.data.frame(de$variable_meta)
head(vm[order(-vm$frac_imputed),
        c("impute_max", "n_imputed", "frac_imputed")])

## ----step-process-------------------------------------------------------------
out <- processDataset(
  sfeat, smeta,
  xlsx_path        = sxlsx,
  data_source      = "skyline",
  data_tab_names   = "skyline_data",
  signal_filter    = "LOD",
  blank_filter     = 3,
  normalize        = "protein_ug",
  adjust_conc      = TRUE,
  starting_vol_col = "starting_vol_uL",
  replace_MVs      = 0.2,
  batch_correction = TRUE,
  bc_qc_label      = "QC",
  bc_header        = "Chrom_Batch",
  blank_head       = "Sample_type")

dim(out[[1]]$data)

## ----step-process-measured----------------------------------------------------
sum(is.na(out$measured$data))                         # gaps as measured
sum(as.data.frame(out[[1]]$variable_meta)$n_imputed)  # filled in out[[1]]

## ----syn-load-----------------------------------------------------------------
syn <- out$measured
dim(as.data.frame(syn$data))
sum(is.na(as.data.frame(syn$data)))

## ----qc-cv--------------------------------------------------------------------
syn <- addCVMetrics(syn, qc_label = "QC")

cv <- as.data.frame(syn$variable_meta)[
  , c("Compound", "Enzymatic_pathway", "CV_QC", "CV_sample",
      "CV_sample_vs_QC")]

DT::datatable(cv[order(cv$CV_QC), ], rownames = FALSE, filter = "top",
              options = list(pageLength = 10, scrollX = TRUE))

## ----qc-drift-----------------------------------------------------------------
plotDrift(syn, "Analyte_02", value_label = "response / ug protein")

## ----qc-qq, fig.width = 7, fig.height = 3.6-----------------------------------
qq_grp <- function(level)
  plotQQ(subsetDataset(syn, conditions = list(Treatment = level,
                                              Sample_type = "Sample")),
         "Analyte_02", transform = "log2") +
    ggplot2::labs(title = level)

cowplot::plot_grid(qq_grp("PBS"), qq_grp("HDM"), nrow = 1)

## ----qc-pca-------------------------------------------------------------------
syn_bio <- subsetDataset(syn, conditions = list(Sample_type = c("Sample", "QC")))

pca <- mrmRunPCA(syn_bio, impute = "gaussian", impute_frac = 0.2,
               transform = "log2", feat_na_max = 0.6,
               sample_na_max = 0.6)
pca$steps

## ----qc-pca-type, fig.height = 4.2--------------------------------------------
mrmPlotPCA(pca, syn_bio, colour_by = "Sample_type")

## ----qc-pca-batch, fig.height = 3.4-------------------------------------------
keep <- list(Sample_type = c("Sample", "QC"))
before <- subsetDataset(de_uncorr, conditions = keep)
after  <- subsetDataset(de_bc,     conditions = keep)

pca_before <- mrmRunPCA(before, impute = "gaussian", impute_frac = 0.2,
                        transform = "log2")
pca_after  <- mrmRunPCA(after,  impute = "gaussian", impute_frac = 0.2,
                        transform = "log2")

cowplot::plot_grid(
  mrmPlotPCA(pca_before, before, colour_by = "Chrom_Batch") +
    ggplot2::labs(title = "Before batch correction"),
  mrmPlotPCA(pca_after, after, colour_by = "Chrom_Batch") +
    ggplot2::labs(title = "After batch correction"),
  nrow = 1)

## ----qc-loadings, fig.width = 7, fig.height = 5-------------------------------
plotLoadings(pca, syn_bio, colour_by = "Enzymatic_pathway")

## ----qc-features--------------------------------------------------------------
DT::datatable(summariseFeatures(syn, out[[2]]), rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

## ----stats-subset-------------------------------------------------------------
samples <- subsetDataset(syn,
                          conditions = list(Sample_type = "Sample"))

table(samples$sample_meta$Treatment, samples$sample_meta$Type)

## ----stats-summary------------------------------------------------------------
gs <- summariseGroups(samples, "Treatment")
DT::datatable(cbind(gs[, 1:3], round(gs[, 4:6], 3)), rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

## ----stats-run----------------------------------------------------------------
params <- list(
  sig_threshold = 0.05,
  comparisons = list(enabled = TRUE, entries = list(
    list(name    = "PBS_vs_HDM",
         method  = "welch",
         compare = list(factor = "Treatment",
                        levels = c("PBS", "HDM"))))),
  ion_ratios = list(enabled = TRUE, ratios = list(
    list(num = "Analyte_02", den = "Analyte_05",
         label = "A02/A05"))))

st <- runStats(samples, params)

## ----stats-truth--------------------------------------------------------------
tt <- st$stats[st$stats$method == "welch",
               c("feature", "n_g1", "n_g2", "log2FC", "p_value", "p_adj")]
truth <- as.data.frame(syn$variable_meta)[, c("Compound",
                                              "Simulated_log2FC")]
chk <- merge(tt, truth, by.x = "feature", by.y = "Compound")
chk <- chk[order(-abs(chk$Simulated_log2FC), chk$p_value), ]

chk$log2FC           <- round(chk$log2FC, 3)
chk$Simulated_log2FC <- round(chk$Simulated_log2FC, 3)
chk$p_value          <- signif(chk$p_value, 3)
chk$p_adj            <- signif(chk$p_adj, 3)

DT::datatable(chk, rownames = FALSE, filter = "top",
              options = list(pageLength = 8, scrollX = TRUE))

## ----stats-box, fig.width = 5, fig.height = 4---------------------------------
plotBoxplot(samples, "Analyte_02", "Treatment",
             value_label = "response / ug protein")

## ----stats-bar, fig.width = 5, fig.height = 4---------------------------------
plotBoxplot(samples, "Analyte_02", "Treatment", type = "bar",
             error_bar = "SE", value_label = "response / ug protein")

## ----stats-volcano, fig.width = 5.5, fig.height = 5---------------------------
plotVolcano(st, "PBS_vs_HDM", use_adjusted = TRUE, sig_threshold = 0.05)

## ----stats-heatmap, fig.width = 7.5, fig.height = 6---------------------------
mrmPlotHeatmap(samples,
             color_samples_by  = "Treatment",
             group_features_by = "Enzymatic_pathway",
             transform         = "log2",
             title             = "All analytes")

## ----stats-group-heatmap, fig.width = 7.5, fig.height = 6---------------------
plotGroupHeatmap(samples,
                   color_samples_by  = "Treatment",
                   group_features_by = "Enzymatic_pathway",
                   transform         = "log2",
                   title             = "All analytes, group means")

## ----stats-corr, fig.width = 6.5, fig.height = 6------------------------------
cp <- list(correlations = list(enabled = TRUE, entries = list(
  list(name = "pairs", methods = "spearman",
       subsets = list(list())))))

plotCorrelation(runStats(samples, cp), samples,
                 group_by = "Enzymatic_pathway")

## ----stats-ratios-------------------------------------------------------------
ir <- st$ion_ratios[!is.na(st$ion_ratios$comparison), ]
DT::datatable(cbind(ir[, c("ratio", "comparison", "sample", "group")],
         round(ir[, c("value", "p_value")], 4)),
   rownames = FALSE, filter = "top",
              options = list(pageLength = 5, scrollX = TRUE))

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