---
title: "Introduction to proteomics data analysis - MaxQuant Data Dependent Acquisition spike-in study"
author:
- name: Lieven Clement, Stijn Vandenbulcke, Oliver M. Crook and Christophe Vanderaa
output:
  BiocStyle::html_document:
    self_contained: yes
    toc: true
    toc_float: true
    toc_depth: 3
    code_folding: show
bibliography: msqrob2.bib
vignette: >
  %\VignetteIndexEntry{A. label-free workflow with two group design}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Introduction

Label-Free Quantitative mass spectrometry based workflows for differential 
expression (DE) analysis of proteins is often challenging due to 
peptide-specific effects and context-sensitive missingness of peptide 
intensities.

- `msqrob2` provides  peptide-based workflows that can assess for DE directly 
from peptide intensities and outperform summarisation methods which first 
aggregate MS1 peptide intensities to protein intensities before DE analysis.
However, they are computationally expensive, often hard to understand for the 
non-specialised end-user, and they do not provide protein summaries, which are 
important for visualisation or downstream processing.

- `msqrob2` therefore also proposes a novel summarisation strategy, which 
estimates MSqRob's model parameters in a two-stage procedure circumventing the 
drawbacks of peptide-based workflows.

- the summarisation based workflow in `msqrob2` maintains MSqRob's superior 
performance, while providing useful protein expression summaries for plotting 
and downstream analysis.
Summarising peptide to protein intensities considerably reduces the 
computational complexity, the memory footprint and the model complexity.
Moreover, it renders the analysis framework to become modular, providing users 
the flexibility to develop workflows tailored towards specific applications.

In this vignette we will demonstrate how to perform  `msqrob`'s summarisation 
based workflow  starting from a Maxquant search on a subset of the cptac 
spike-in study.

Examples on our peptide-based workflows and on the analysis of more complex 
designs can be found on our companion website 
[msqrob2 book](https://statomics.github.io/msqrob2book).

Technical details on our methods can be found in [@goeminne2016], 
[@goeminne2020], [@sticker2020] and [@vandenbulcke2025].

# Background

This case-study is a subset of the data of the 6th study of the Clinical
Proteomic Technology Assessment for Cancer (CPTAC).
In this experiment, the authors spiked the Sigma Universal Protein Standard
mixture 1 (UPS1) containing 48 different human proteins in a protein background
of 60 ng/$\mu$L Saccharomyces cerevisiae strain BY4741.
Two different spike-in concentrations were used:
6A (0.25 fmol UPS1 proteins/$\mu$L) and 6B (0.74 fmol UPS1 proteins/$\mu$L) [5].
We limited ourselves to the data of LTQ-Orbitrap W at site 56.
The data were searched with MaxQuant version 1.5.2.8, and
detailed search settings were described in Goeminne et al. (2016) [1].
Three replicates are available for each concentration.

# Load packages

We load the `msqrob2` package, along with additional packages for
data manipulation and visualisation.


```{r load_libraries}
suppressPackageStartupMessages({
library("QFeatures")  
library("dplyr") 
library("tidyr")
library("ggplot2")
library("msqrob2")    
library("stringr")
library("ExploreModelMatrix")
library("kableExtra")
library("ComplexHeatmap")
library("scater")
})
```

# Data

## Peptide table

We first import the data from the peptide.txt file. This is the file containing
your peptide-level intensities searched by MaxQuant search,
this peptide.txt file can be found by default in the
"path_to_raw_files/combined/txt/" folder from the MaxQuant output,
with "path_to_raw_files" the folder where the raw files were saved.
In this vignette, we use a MaxQuant peptide file which is a subset
of the cptac study. 
This data is available in 
[msqrob2data](https://github.com/statOmics/msqrob2data) GitHub repository.
To import the data we use the `QFeatures` package.

We generate the object peptideFile with the path to the peptide.txt file.
This can be a path to the data on the web or a path to data on your local 
computer. 

```{r file_location}
peptideFile <- "https://github.com/statOmics/msqrob2data/raw/refs/heads/main/dda/cptacAvsB_lab3/peptides.txt"
```


We will import the data using the `fread` function of the data.table package. We set the `check.names` argument equal to  `TRUE` so that all column names are retrieved using the `$` operator.  

```{r import_data}
peptides <- data.table::fread(peptideFile, check.names = TRUE)
head(peptides)
```

Each row in the peptide data table contains information about one peptide 
(the table above shows the first 6 rows). The columns contains various 
descriptors about the peptide, such as its sequence, its charge, the amino acid 
composition, etc. Some of these columns (those starting with `Intensity `) 
contain the quantification values for each sample. The table format where the 
quantitative values for each sample are contained in a separate column is 
depicted as the "wide format", as opposed to the "long format" 
(e.g., the [PSM table]).

```{r extract_quantCols}
quantCols <- grep("Intensity[.]", names(peptides), value = TRUE)
```

Unlike for real experiments, we known the ground truth for this dataset: UPS 
proteins were spiked in at different concentrations and are thus differentially 
abundant (DA, spiked in). 
Yeast proteins consist of the background proteome that is the same in all 
samples and are thus non-DA.  

```{r add_species_info}
peptides <- peptides |> 
  mutate(species = grepl(pattern = "UPS",Proteins) |> 
           as.factor() |>
           recode("TRUE"="ups","FALSE" = "yeast"))
```

## Sample annotation 

Each row in the annotation table contains information about one sample. The 
columns contain various descriptors about the sample, such as the name of the 
sample or the MS run, the treatment (here the spike-in condition), or any other 
biological or technical information that may impact the data quality or the 
quantification. Without an annotation table, no analysis can be performed. The 
sample annotations are generated by the researcher. In this example, the 
annotations are extracted from the sample names, although reporting a detailed 
design of experiments in a table is better practice.

1. We first generate a variable quantCols, which contain the names of the 
quantification columns. 
2. Next we generate variable `sampleId`, to pinpoint the different samples in 
the dataset. Since all quantification columns start with the generic string 
`Intensity 6` we replace the string with an empty string `""` so we keep the 
last 3 characters to refer to the sample. 
3. We construct a variable condition, by using the first letter of the 
sampleId, which refers to the spikein condition. 
4. We construct a variable rep, by using the third letter of the 
sampleId.


```{r make_annotation}
(annot <- data.frame(quantCols = quantCols) |> #1.
  mutate(sampleId = gsub("Intensity[.]6","", quantCols), #2. 
         condition = substr(sampleId,1,1), #3.
         rep = substr(sampleId, 3, 3)) #4.
)
```

## Convert to QFeatures


The `readQFeatures()` enables a seamless conversion of tabular data into a 
`QFeatures` object. We provide the peptide table and the sample annotation table. 
The function will use the `quantCols` column in the sample annotation table to 
understand which columns in `peptides` contain the quantitative values, and 
automatically link the corresponding sample annotation with the quantitative 
values. We also tell the function to use the `Sequence` column as peptide 
identifier, which will be used as rownames. See `?readQFeatures()` for more 
details.

```{r convert_to_QFeatures}
(qf <- readQFeatures(assayData = peptides, 
                     colData = annot,
                     name = "peptides",
                     fnames = "Sequence")
)
```

# Data preprocessing 

Since we have a `QFeatures` object, we can directly make use of `QFeatures’` 
data preprocessing functionality. A major advantage of QFeatures is it provides 
the power to build highly modular workflows, where each step is carried out by a 
dedicated function with a large choice of available methods and parameters. This 
means that users can adapt the workflow to their specific use case and their 
specific needs.

## Encoding missing values

The first preprocessing step is to correctly encode the missing values. 
It is important that missing values are encoded using `NA`. For instance, 
non-observed values should not be encoded with a zero because true zeros 
(the proteomic feature is absent from the sample) cannot be distinguished from 
technical zeros (the feature was missed by the instrument or could not be 
identified). We therefore replace any zero in the quantitative data with an `NA`.

```{r convert0_to_NA}
qf <- zeroIsNA(qf, "peptides") # convert 0 to NA
```

Note that `msqrob2` can handle missing data without having to rely on hard-to-verify imputation assumptions, which is our general recommendation. However, `msqrob2` does not prevent users from using imputation, which can be performed with `impute()` from the QFeatures package. 

## $\log_2$ transformation

We perform log2-transformation with `logTransform()` from the `QFeatures` 
package.

```{r log2_transformation}
qf <- logTransform(qf, base = 2, i = "peptides", name = "peptides_log")
```


## Peptide-Filtering

Filtering removes low-quality and unreliable peptides that would otherwise introduce noise and artefacts in the data. 

### Remove failed protein inference
We remove peptides that could not be mapped to a protein. We also remove peptides that cannot be uniquely mapped to a single proteins because its sequence is shared across multiple proteins. This often occurs for homologs or for proteins that share similar functional domains. Shared peptides are often mapped to protein groups, which are the set of proteins that share the given peptides. Protein groups are encoded by appending the constituting protein identifiers, separated by a `";"`.

```{r failed_protein_inference_filtering}
qf <- filterFeatures(
  qf, ~ Proteins != "" & ## Remove failed protein inference
    !grepl(";", Proteins)) ## Remove protein groups
```

### Remove reverse sequences and contaminants

We now remove the peptides that map to decoy sequences or to
contaminants proteins. In our `peptide.txt` 
file decoys and contaminants are indicated with a "+" character in the 
column `Reverse` and `Potential.contaminant`, respectively. 
Note, that depending on the version of MaxQuant the name of the contaminant variable differs (e.g. `Contaminant` or `Potential contaminant`). 
Note, that while reading in the `peptides` table spaces in the column names are replaced by a dot. 


```{r filtering_decoys_contaminants}
qf <- filterFeatures(qf, ~ Reverse != "+" & ## Remove decoys
                          Potential.contaminant != "+") ## Remove contaminants
```

### Remove highly missing peptides

The data are characterised by a high proportion of missing values as
reported by `nNA()` from the `QFeatures` package. The function
computes the number and the proportion of missing values across
peptides, samples or for the whole data set and return a list of three
tables called `nNArows`, `nNAcols`, and `nNA`, respectively.

```{r assess_NA_features}
nNaRes <- nNA(qf, "peptides")
range(nNaRes$nNAcols$pNA)
```

The samples contain between `r round(min(nNaRes$nNAcols$pNA)*100, 1)`%
and `r round(max(nNaRes$nNAcols$pNA)*100, 1)`% missing values. The
missingness within each peptide is more variable, with most peptides
found accross all samples (`nNA` is 0) and other that could not be
quantified in any sample (`nNA` is `r ncol(qf[[1]])`), as depicted
by the histogram below.

```{r plot_NA_features}
data.frame(nNaRes$nNArows) |> 
  ggplot() +
  aes(x = nNA) +
  geom_histogram() +
  ggtitle("Number of missing values for each peptide") +
  theme_bw()
```


We keep peptides that were observed at last 3 times out of the $n
= 6$ samples, so that we can estimate the peptide characteristics. We tolerate 
the following proportion of NAs:
$\text{pNA} = \frac{(n - 3)}{n} = 0.5$, so we keep peptides that
are observed in at least 50% of the samples, which corresponds to one treatment 
condition. 
This is an arbitrary value that may need to be adjusted depending on 
the experiment and the data set.

```{r filter_features_with_many_NA}
nObs <- 3
n <- ncol(qf[["peptides_log"]])
(qf <- filterNA(qf, i = "peptides_log", pNA = (n - nObs) / n))
```

We keep `r nrow(qf[["peptides_log"]])` peptides upon filtering.

## Normalisation{#sec-norm}

The most common objective of MS-based proteomics experiments is to
understand the biological changes in protein abundance between
experimental conditions. However, changes in measurements between
groups can be caused due to technical factors. For instance, there are
systematic fluctuations from run-to-run that shift the measured
intensity distribution. We can this explore as follows:

1. We extract the sets containing the log transformed data. This is 
performed using `QFeatures`' 3-way subsetting.
2. We use `longForm()` to convert the `QFeatures` object into a long
table, including all `colData` for filtering and
colouring.
3. We visualise the density of the quantitative values within each
sample. We colour each sample based on its spike-in condition.

```{r normalisation_needed}
qf[, , "peptides_log"] |> #1.
  longForm(colvars = colnames(colData(qf))) |>  #2. 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + #3.
  aes(x = value,
      colour = condition,
      group = colname) +
  geom_density() +
  theme_minimal()
```

Even in this clean synthetic data set (same yeast background with 
some spiked UPS proteins), the marginal precursor intensity
distributions across samples are not well aligned. Ignoring this
effect will increase the noise and reduce the statistical power of the
experiment, and may also, in case of unbalanced designs, introduce
confounding effects that will bias the results.

There are many ways to perform the normalisation, e.g.  
median centering is a popular choice. 
If we subtract the sample median at the log scale from each precursor 
log2 intensity, this basically boils down to calculating log-ratio's between 
the precursor intensity and its sample median. 
Here, we choose to work with offsets, that are chosen in such a way so that 
the distributions are centered around the location of a reference sample. 
E.g. the median of the sample medians. 

Note that users can also use all normalisation functions that are provided in 
QFeatures, i.e. by replacing the sweep function by `QFeatures::normalize`.


```{r normalize_peptides}
qf <- sweep( #Subtract log2 norm factor column-by-column (MARGIN = 2)
  qf,
  MARGIN = 2,
  STATS = nfLogMedian(qf,"peptides_log"),
  i = "peptides_log",
  name = "peptides_norm"
)
```


We explore the effect of the global normalisation in the subsequent plot. 

Formally, the function applies the following operation on each sample
$i$ across all precursors $p$:

$$
y_{ip}^{\text{norm}} = y_{ip} - \log_2(nf_i)
$$

with $y_{ip}$ the log2-transformed intensities and $\log_2(nf_i)$ the log2-transformed 
norm factor. Upon normalisation, we can see that the distribution of the 
$y_{ip}^{\text{norm}}$ nicely overlap (using the same code as above)

```{r assess_normalisation}
qf[, , "peptides_norm"] |>
  longForm(colvars = colnames(colData(qf))) |> 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + 
  aes(x = value,
      colour = condition,
      group = colname) +
  geom_density() +
  labs(subtitle = "Normalised log2 precursor intensities") +
  theme_minimal()
```

## Summarization to protein level

Here, we summarise the peptide-level data into protein intensities using
`aggregateFeatures()`, which streamlines summarisation. It requires the name
of a `rowData` column to group the peptides into proteins (or
protein groups), here `Proteins`. We can provide the summarisation
approach through the `fun` argument. 
We use the standard summarisation in aggregateFeatures, which is a
robust summarisation method.

```{r summarize_peptides_to_proteins, warning=FALSE}
qf <- aggregateFeatures(qf,
    i = "peptides_norm", 
    fcol = "Proteins", 
    na.rm = TRUE,
    name = "proteins"
)
```

Note, that robust summarisation can take a long time for large datasets. 
In that case we suggest to use `fun = MsCoreUtils::medianPolish, na.rm = TRUE`. 
Note that the optional argument na.rm = TRUE is needed when using `medianPolish` 
because the data contains missing values as we did not adopt imputation and 
`medianPolish` by default cannot handle missing values. 

# Data exploration and QC

Data exploration aims to highlight the main sources of variation in
the data prior to data modelling and can pinpoint to outlying or
off-behaving samples. 

## Marginal distribution at precursor and protein level 

```{r qc_normalisation_peptides}
qf[, , "peptides_norm"] |>
  longForm(colvars = colnames(colData(qf))) |> 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + 
  aes(x = value,
      colour = condition,
      group = colname) +
  geom_density() +
  theme_minimal() +
   labs(subtitle = "Normalised log2 precursor intensities")
```

```{r qc_normalisation_peptides_boxplot}
qf[, , "peptides_norm"] |>
  longForm(colvars = colnames(colData(qf))) |> 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + 
  aes(x = sampleId,
      y = value,
      colour = condition,
      group = colname) +
  xlab("sample") +
  geom_boxplot() +
  theme_minimal() +
   labs(subtitle = "Normalised log2 precursor intensities")
```

We do the same for the protein-level values. 

```{r qc_normalisation_proteins}
qf[, , "proteins"] |>
  longForm(colvars = colnames(colData(qf))) |> 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + 
  aes(x = value,
      colour = condition,
      group = colname) +
  geom_density() +
  theme_minimal() +
   labs(subtitle = "Normalised log2 protein intensities")
```

```{r qc_normalisation_proteins_boxplot}
qf[, , "proteins"] |>
  longForm(colvars = colnames(colData(qf))) |> 
  data.frame() |> 
  filter(!is.na(value)) |>
  ggplot() + 
  aes(x = sampleId,
      y = value,
      colour = condition,
      group = colname) +
  xlab("sample") +
  geom_boxplot() +
  theme_minimal() +
   labs(subtitle = "Normalised log2 protein intensities")
```

## Identifications per sample 

```{r qc_identifications}
qf[,,"peptides_norm"] |> 
  longForm(colvars = colnames(colData(qf)), 
           rowvars= c("Sequence", 
                      "Proteins")) |>
  data.frame() |>
  filter(!is.na(value)) |>
  group_by(condition, sampleId) |>
  summarise(Precursors = length(unique(Sequence)),
            `Protein Groups` = length(unique(Proteins))) |> 
  pivot_longer(-(1:2),
               names_to = "Feature",
               values_to = "IDs") |> 
  ggplot(aes(x = sampleId, y = IDs, fill = condition)) +
  geom_col() +
  #scale_fill_observable() +
  facet_wrap(~Feature,
             scales = "free_y") +
  labs(title = "Peptide and protein group identificiations per sample",
       x = "Sample",
       y = "Identifications") +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 90))
```

## Dimensionality reduction plot 

A common approach for data exploration is to
perform dimension reduction, such as Multi Dimensional Scaling (MDS).
We will first extract the set to explore along the sample annotations
(used for plot colouring).

```{r qc_mds_proteins}
getWithColData(qf, "proteins") |> 
  as("SingleCellExperiment") |> 
  runMDS(exprs_values = 1) |> 
  plotMDS(colour_by = "condition", shape_by = "rep", point_size = 3)
```


This plot reveals interesting information. First, we see that the
samples are nicely separated according to their spike-in condition.
The also seem to line up according to `rep`. 

## Correlation matrix

```{r qc_correlation_proteins}
corMat <- qf[["proteins"]] |> 
  assay() |>
  cor(method = "spearman", use = "pairwise.complete.obs") 
 
colnames(corMat) <- qf$sampleId
rownames(corMat) <- qf$sampleId
corMat |> 
  ggcorrplot::ggcorrplot() +
  scale_fill_viridis_c() +
  labs(title = "Correlation matrix",
       fill = "Correlation") +
  theme(axis.text.x = element_text(angle = 90))
```


# Data Modeling (Robust Regression){#sec-modelling}

## Sources of variation 

For this experiment, we can only model a single source of variation: spike-in 
concentration. All remaining sources of variability, e.g. pipetting errors, 
MS run to run variability, etc will be lumped in the residual variability 
(error term).

Spike-in condition effects: we model the source of variation induced by the 
experimental treatment of interest as a fixed effect, which we consider 
non-random, i.e. the treatment effect is assumed to be the same in repeated 
experiments, but it is unknown and has to be estimated.

When modelling a typical label-free experiment at the protein level, the model 
boils down to the following linear model:

$$ y_i = \mathbf{x}_i^T\boldsymbol{\beta} +\epsilon_i$$
with the $\log_2$-normalised and summarised protein intensities in sample; 
$\mathbf{x}_i$ a vector with the covariate pattern for the sample encoding the 
intercept, spike-in condition, or other experimental factors;  
$\boldsymbol{\beta}$ the vector of parameters that model the association between 
the covariates and the outcome; and $\epsilon_i$ the residuals reflecting 
variation that is not captured by the fixed effects. 
Note that allows for a flexible parameterisation of the treatment 
beyond a single covariate, i.e. including a 1 for the intercept, continuous and 
categorical variables as well as their interactions. We assume the residuals to 
be independent and identically distributed (i.i.d) according to a normal 
distribution with zero mean and constant variance, 
i.e. $\epsilon_i \sim N(0, \sigma^2_\epsilon)$, that can differ from protein to 
protein.

In R, this model is encoded using the following simple formula:
```{r define_model}
model <- ~ condition
```

## Model Estimation 

We estimate the model with `msqrob()`. The function takes the QFeatures object, 
extracts the quantitative values from the `"proteins"` set generated during 
summarisation, and fits a simple linear model with `condition` as covariate, 
which are automatically retrieved from `colData(qf)`.
Note, that `msqrob()` can also take a `SummarizedExperiment` of 
precursor, peptide or protein intensities as its input. 

```{r fit_models, warning=FALSE}
  qf <- msqrob(
    qf,  
    i = "proteins",
    formula = model,
    robust = TRUE)
```
  

We enable M-estimation (`robust = TRUE`) for improved robustness against outliers.

The fitting results are available in the `msqrobModels` column of the rowData. More specifically, the modelling output is stored in the `rowData` as a `statModel` object, one model per row (protein). We will see in a later section how to perform statistical inference on the estimated parameters.

```{r print_models}
models <- rowData(qf[["proteins"]])[["msqrobModels"]]
models[1:3]
```

## Statistical inference 

We can now convert the research question “do the spike-in conditions affect the protein intensities?” into a statistical hypothesis.

In other words, we need to translate this question in a null and alternative hypothesis on a single model parameter or on a linear combination of model parameters, which is also referred to with a contrast.

To aid defining contrasts, we will visualise the experimental design using the `ExploreModelMatrix` package.


```{r explore_model_matrix}  
vd <- ExploreModelMatrix::VisualizeDesign(
    sampleData =  colData(qf),
    designFormula = model,
    textSizeFitted = 4
)
vd$plotlist
```


This plot shows that the average protein $\log_2$ intensity for conditionA is modelled by `(Intercept)` and that for the `conditionB` group by `(Intercept) + conditionB` .
 
Hence, to assess differential abundance of proteins between the two spike-in conditions, we will have to assess the following null hypothesis:  $\log_2 FC_{B-A} = (Intercept + conditionB) - Intercept = conditionB = 0$. 

Based on this null hypothesis we now generate a contrast matrix. 

```{r make_contrasts}
(L <- makeContrast(
  "conditionB = 0",
  parameterNames = colnames(vd$designmatrix)
))
```


We can now test our null hypothesis using `hypothesisTest()` which
takes the `QFeatures` object with the fitted model and the contrast we
just built. `msqrob2` automatically applies the hypothesis
testing to all features in the assay.

```{r hypthesis_tests}
qf <- hypothesisTest(qf, i = "proteins", contrast = L)
```

The results tables for the contrast are stored in the row data of the 
`proteins` assay in `qf` under the colomnames of the contrast matrix `L`.

Users can directly access the results tables via de colData. However, `msqrob2` 
also provides the `msqrobCollect` function. This function will automatically 
retrieve the results tables for all contrasts in contrast matrix `L` and combine 
them by default in one  table. This table contains two additional columns: 
`contrast` with the name of the contrast and `feature` with the name of the 
modelled `feature`, here protein. 

```{r collect_results}
inferences <- 
  msqrobCollect(qf[["proteins"]], L) 
head(inferences)
```


### Results table for significant proteins

We will return a table of proteins for which the contrast is significant at the 5% FDR level. 

1. We set the significance level 
3. We filter the significant results for the contrast from the table
4. We print the table 

```{r significant_table, results='asis'}
alpha <- 0.05 #1.
sigList <- inferences |> 
    filter(adjPval < alpha) #3.
kable(sigList) |>
      kable_styling(full_width = FALSE) |>
      scroll_box(height = "250px") #4.
```


### Volcanoplots 

A volcano plot is a common visualisation that provides an overview of
the hypothesis testing results, plotting the $-\log_{10}$ p-value^[Note, 
that upon this transformation a value of 1 represents a p-value of 0.1, 2 
a p-value of 0.01, 3 a p-value of 0.001, etc.] as a function of
the estimated log fold change. Volcano plots can be used to highlight
the most interesting proteins that have large fold changes and/or are
highly significant. We can use the table above directly to build a
volcano plot using `ggplot2` functionality. We also highlight which
proteins are UPS standards, known to be differentially abundant by
experimental design.

Experienced users can make the plot themselves. However, in msqrob2 we also 
provide the plotVolcano function that generates volcanoplots based on `msqrob2`
inference tables generated with the `hypothesisTest` function. 


```{r volcanoplot}
inferences |> 
  plotVolcano()
```

### Heatmaps

We construct heatmaps for the significant proteins for each contrast. 

1. We get the names of the significant features.
2. We extract the quant data for the significant features and scale it. 
3. We extract annotation
4. We make the heatmap using the quants data. We do not cluster columns (samples) to keep them together according to the design.


```{r heatmap}
sig <- sigList |> pull(feature) #1.


quants <- assay(qf,"proteins")[sig,] |> 
                     t() |> 
                     scale() |>
                     t() #2.
 
colnames(quants) <- qf$sampleId #3.
annotations <- columnAnnotation(condition = qf$condition) #3.

set.seed(1234) ## annotation colours are randomly generated by default
Heatmap(
  quants,
  name = "log2 intensity",
  top_annotation = annotations
  )  #4.
```

Note, that `r nrow(sigList)` DA proteins are returned and that they are all UPS proteins! 

A typical workflow for the analysis of a proteomics experiment will end here. 


# Assess performance

Note, that the evaluations in this section are typically not possible on real experimental data. 
Indeed, we are using a spike-in dataset so we know the ground truth: all human UPS proteins are differentially abundant (DA) between the conditions and the yeast proteins are non-DA.

We use a nominal FDR level of 0.05

```{r define_alpha}
alpha <- 0.05
```

## Real Fold changes 

As this is a spike-in study with known ground truth, we can also plot the log2 fold change distributions against the expected values, in this case 0 for the yeast proteins and the difference of the log concentration for the spiked-in UPS standards, i.e. $logFC = \log_2 (0.74 fmol/0.25 fmol) = \log_2 2.96 = 1.57$. 

```{r real_logFC}
realLogFC <- log2(0.74/0.25)
```

We collect the inference table again and now add the species information.  

```{r collect_results_species}
inferences <- 
  msqrobCollect(qf[["proteins"]], L)  |> 
  mutate(species =  rowData(qf[["proteins"]])$species)
```

We now evaluate if the estimated fold changes correspond to the real fold changes.

1. We plot the log fold changes in function of the spikein condition and color according to spikein condition.
2. We add a boxplot layer.
3. We use custom colors. 
4. We add a horizontal line at 0, which corresponds to the known log2 fold change for yeast proteins
5. We add a horizontal lines for known log2 fold changes of the spiked UPS proteins. 

```{r assess_logFC}
logFC <- inferences |> 
  filter(!is.na(logFC)) |>
  ggplot(aes(x = species, y = logFC, color = species)) + #1.
  geom_boxplot() + #2.
  theme_bw() + 
  scale_color_manual( 
    values = c("grey20", "firebrick"), #3. 
    name = "",
    labels = c("yeast", "ups")
    ) + 
  geom_hline(yintercept = 0, color="grey20") + # 4. 
  geom_hline(yintercept = realLogFC, color="firebrick") #5. 
logFC
```

Note, that the msqrob2 estimates are nicely distributed around the true fold changes. 


## True and false positives

All human UPS proteins are differentially abundant (DA) between the conditions and the yeast proteins are non-DA.

```{r assess_TP_FP_FDP}
inferences |>
    filter(adjPval < alpha) |>
    summarise("TP" = sum(species == "ups"),
              "FP" = sum(species != "ups"),
              "FDP" = mean(species != "ups"))
```

Note, that at the 5% FDR level 20 true positives (TP, UPS proteins) are returned as DA, 0 false positives (FP, yeast proteins) are returned DA, resulting in a false discovery proportion (FDP) of 0, confirming that msqrob2 is providing strong FDR control for this dataset. 

## TPR - FDP curves

We generate the TPR-FDP curves to assess the performance to prioritise differentially abundant proteins. Again, these curves are built using the ground truth information about which proteins are differentially abundant (spiked in) and which proteins are constant across samples. We create two functions to compute the TPR and the FDP.

```{r FDR_TPR_functions}
computeFDP <- function(pval, tp) {
    ord <- order(pval)
    fdp <- cumsum(!tp[ord]) / 1:length(tp)
    fdp[order(ord)]
}
computeTPR <- function(pval, tp, nTP = NULL) {
    if (is.null(nTP)) nTP <- sum(tp)
    ord <- order(pval)
    tpr <- cumsum(tp[ord]) / nTP
    tpr[order(ord)]
}
```

We apply these functions and compute the corresponding metric using the statistical inference results and the ground truth information.

```{r calculate_performance}
performance <- inferences |> group_by(contrast) |> 
    na.exclude() |>
    mutate(tpr = computeTPR(pval, species == "ups"),
           fdp = computeFDP(pval, species == "ups")) |>
    arrange(fdp)
```

We also highlight the observed FDP at a $alpha = 5\%$ FDR threshold.

```{r calculate_workpoints}
workPoints <- performance |> 
  group_by(contrast) |>
    filter(adjPval < 0.05) |>
    slice_max(pval)
```

```{r plot_tpr_fdp}
ggplot(performance) +
    aes(
        y = fdp,
        x = tpr,
    ) +
    geom_line() +
    geom_point(data = workPoints, size = 3) +
    geom_hline(yintercept = 0.05, linetype = 2) +
    coord_flip(ylim = c(0, 0.2)) +
    theme_minimal()
```

# References

Note, that more examples can be found in our [msqrob2 book](https://statomics.github.io/msqrob2book).

```{r sessionInfo}
sessionInfo()
```

