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

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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

## ----warning=F,message=F------------------------------------------------------
library(PolySTest)
library(SummarizedExperiment)
# Load the data file "LiverAllProteins.csv" from the PolySTest package
filename <- system.file("extdata", "LiverAllProteins.csv", 
                        package = "PolySTest")

dat <- read.csv(filename, row.names=1, stringsAsFactors = FALSE)
dat <- dat[1:200,] # Use a subset of the data for this example

## -----------------------------------------------------------------------------
# Normalization (example)
dat <- t(t(dat) - colMedians(as.matrix(dat), na.rm = TRUE))

## -----------------------------------------------------------------------------
sampleMetadata <- data.frame(Condition = rep(paste("Condition", 1:4), each=3),
                             Replicate = rep(1:3, each=4))

fulldata <- SummarizedExperiment(assays = list(quant = dat), 
                                 colData = sampleMetadata)
rowData(fulldata) <- rownames(dat)
metadata(fulldata) <- list(NumReps = 3, NumCond = 4)

assay(fulldata, "quant") <- dat

fulldata <- update_conditions_with_lcs(fulldata)

## -----------------------------------------------------------------------------
conditions <- unique(colData(fulldata)$Condition)

allComps <- create_pairwise_comparisons(conditions, 1)

## -----------------------------------------------------------------------------
fulldata_unpaired <- PolySTest_unpaired(fulldata, allComps)


## ----fig.width=7, fig.height=4------------------------------------------------
# Define comparisons to visualize from available ones
compNames <- metadata(fulldata_unpaired)$compNames
print(compNames)
comps <- compNames[1]

# Plotting the results
plotPvalueDistr(fulldata_unpaired,  comps, c("limma", "Miss Test"))

## ----fig.width=7, fig.height=4------------------------------------------------
# Select proteins with FDR < 0.01

sigProts <- which(rowData(fulldata_unpaired)[, paste0("FDR_PolySTest_", 
                                                      comps[1])] < 0.01)

# Volcano plot
plotVolcano(fulldata_unpaired, compNames = comps, sel_prots = sigProts, 
            testNames = c("PolySTest", "limma", "Miss Test"))

## ----fig.width=7, fig.height=4------------------------------------------------
plotUpset(fulldata_unpaired, qlim=0.01)

## ----fig.width=8, fig.height=6------------------------------------------------
plotExpression(fulldata_unpaired, sel_prots = sigProts)

## ----fig.width=7, fig.height=4------------------------------------------------
plotRegNumber(fulldata_unpaired)

## ----fig.width=6, fig.height=6------------------------------------------------
plotHeatmaply(fulldata_unpaired, sel_prots = sigProts, heatmap_scale = "row")

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