## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()

## ----env, include=FALSE, echo=FALSE, cache=FALSE------------------------------
library("knitr")
opts_chunk$set(stop_on_error = 1L)
suppressPackageStartupMessages(library("MatrixQCvis"))

## ----install, eval=FALSE------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("MatrixQCvis")

## ----load_MatrixQCvis, eval=TRUE----------------------------------------------
library(MatrixQCvis)

## ---- prepare_se, eval=TRUE, echo=TRUE----------------------------------------
## create synthetic assay using the generate_synthetic_data function from the
## proDA package
library(proDA)
n_samples <- 40
n_feat <- 1000
data <- generate_synthetic_data(n_feat, n_conditions = n_samples / 10, 
    n_replicates = n_samples / 4, frac_changed = 0.1)
a <- data$Y
colnames(a) <- gsub(colnames(a), pattern = "Condition", replacement = "Sample")

## add some treatment-specific effects
set.seed(1)
a[, 1:5] <- a[, 1:5] + rnorm(5000, mean = 1.0, sd = 0.5)
a[, 11:15] <- a[, 11:15] + rnorm(5000, mean = 0.8, sd = 0.5)
a[, 21:25] <- a[, 21:25] + rnorm(5000, mean = 1.2, sd = 0.5)
a[, 31:35] <- a[, 31:35] + rnorm(5000, mean = 0.7, sd = 0.5)

## create information on the samples
type_sample <- gsub(data$groups, pattern = "Condition", replacement = "Type")
trmt_sample <- paste(
   c(rep("1", 10), rep("2", 10), rep("3", 10), rep("4", 10)),
   c(rep("A", 5), rep("B", 5)), sep = "_")
cD <- data.frame(name = colnames(a), type = type_sample, 
                     treatment = trmt_sample)

## create information on the proteins
rD <- data.frame(spectra = rownames(a))

## create se
se <- SummarizedExperiment(assay = a, rowData = rD, colData = cD)

## ---- eval=FALSE--------------------------------------------------------------
#  qc <- shinyQC(se)

## ----hist_sample, eval=TRUE, echo=FALSE, crop=NULL----------------------------
tbl <- hist_sample_num(se = se, category = "type")
partial_bundle(hist_sample(tbl, category = "type"))

## ----mosaic, eval=TRUE, echo=FALSE, crop=NULL---------------------------------
mosaic(se = se, f1 = "type", f2 = "treatment")

## ----samplesMeasuredMissing, eval=TRUE, echo=FALSE, crop=NULL-----------------
tbl <- samplesMeasuredMissing(se) 
partial_bundle(barplotSamplesMeasuredMissing(tbl, measured = TRUE))

## ----histFeature_measured, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, crop=NULL----
partial_bundle(histFeature(assay(se), measured = TRUE))

## ----histFeature_missing, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, crop=NULL----
partial_bundle(histFeature(assay(se), measured = FALSE))

## ----histFeatureCategory, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, crop=NULL----
partial_bundle(histFeatureCategory(se, measured = TRUE, category = "type"))

## ----upsetmeasured, eval=TRUE, echo=FALSE, crop=NULL--------------------------
upsetCategory(se, category = "type", measured = TRUE) 

## ----upsetmissing, eval=TRUE, echo=FALSE, crop=NULL---------------------------
upsetCategory(se, category = "type", measured = FALSE)

## ----boxplot, eval=TRUE, echo=FALSE, warning=FALSE, crop=NULL-----------------
createBoxplot(se, orderCategory = "name", title = "raw", log2 = FALSE, 
               violin = TRUE)

## ----drift, eval=TRUE, echo=FALSE, message=FALSE, warning=FALSE, crop=NULL----
partial_bundle(driftPlot(se, aggregation = "median", category = "type",
          orderCategory = "type", level = "Type_1", method = "loess"))

## ----cv, eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE, crop=NULL-------
a <- assay(se)
a_n <- normalizeAssay(a, "sum")
se_n <- MatrixQCvis:::updateSE(se = se, assay = a_n)
a_b <- batchCorrectionAssay(se_n, method = "none")
a_t <- transformAssay(a_b, "vsn")
se_t <- MatrixQCvis:::updateSE(se = se, assay = a_t)

a_i <- imputeAssay(a_t, "MinDet")

## calculate cv values
cv_r <- cv(a, "raw")
cv_n <- cv(a_n, "normalized")
cv_b <- cv(a_b, "batch corrected")
cv_t <- cv(a_t, "transformed")
cv_i <- cv(a_i, "imputed")
 
df <- data.frame(cv_r, cv_n, cv_b, cv_t, cv_i)
plotCV(df)

## ----meansdplot, eval=TRUE, echo=FALSE, warning=FALSE, crop=NULL--------------
vsn::meanSdPlot(assay(se_t))

## ----maplot, eval=TRUE, echo = FALSE, show=FALSE, crop=NULL-------------------
tbl <- MAvalues(se, group = "all")
MAplot(tbl, group = "all", plot = "Sample_1-1")

## ----ecdf, eval=TRUE, echo=FALSE, crop=NULL-----------------------------------
ECDF(se, "Sample_1-1", group = "all")

## ----distSampe, eval=TRUE, echo=FALSE, crop=NULL------------------------------
a <- assay(se)
dist_mat <- distShiny(a)
distSample(dist_mat, se, "type")
partial_bundle(sumDistSample(dist_mat, title = "raw"))

## ----featurePlot, eval=TRUE, echo=FALSE, crop=NULL----------------------------
l_assays <- list(raw = a, normalized = a_n, `batch.corrected` = a_b,
   transformed = a_t, imputed = a_i)
df_feature <- createDfFeature(l_assays, feature = rownames(se)[1])
featurePlot(df_feature)

## ----cvFeaturePlot, eval=TRUE, echo=FALSE, crop=NULL--------------------------
partial_bundle(cvFeaturePlot(l_assays, lines = FALSE))

## ----colData, echo=FALSE, eval=TRUE-------------------------------------------
colData(se)

## ----modelMatrix, echo=FALSE, eval=TRUE---------------------------------------
mM <- model.matrix(formula("~ treatment + 0"), data = colData(se))
head(mM)

## ----contrasts, echo=FALSE, eval=TRUE, warning=FALSE, message=FALSE-----------
library(limma)
cM <- makeContrasts(contrasts = "treatment1_B-treatment1_A", levels = mM)
head(cM)

## ----topDE, echo=FALSE, eval=TRUE, warning=FALSE, message=FALSE---------------
a <- assay(se)
a_n <- a %>% normalizeAssay(., "none")
a_t <- a_n %>% transformAssay(., "none")
a_i <- a_t %>% imputeAssay(., "MinDet")

fit <- lmFit(a_i, design = mM)
fit <- contrasts.fit(fit, contrasts = cM)
fit <- eBayes(fit, trend = TRUE, robust = TRUE)
tT <- topTable(fit, number = Inf, adjust.method = "fdr", p.value = 0.05)
rmarkdown::paged_table(tT)

## ----volcano, eval=TRUE, echo=FALSE, crop=NULL--------------------------------
partial_bundle(volcanoPlot(tT, type = "ttest"))

## ----session,eval=TRUE, echo=FALSE--------------------------------------------
sessionInfo()

