cptac_maxLFQ.Rmd

---
title: "Introduction to proteomics data analysis: maxLFQ summarization"
author: "Lieven Clement"
date: "statOmics, Ghent University (https://statomics.github.io)"
output:
    html_document:
      code_download: true
      theme: flatly
      toc: true
      toc_float: true
      highlight: tango
      number_sections: true
    pdf_document:
      toc: true
      number_sections: true
linkcolor: blue
urlcolor: blue
citecolor: blue

bibliography: msqrob2.bib

---

<a rel="license" href="https://creativecommons.org/licenses/by-nc-sa/4.0"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png" /></a>

This is part of the online course [Proteomics Data Analysis (PDA)](https://statomics.github.io/PDA/)

# 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.


- NOTE THAT maxLFQ SUMMARISATION IS SUBOPTIMAL! 
- THIS IS FOR DIDACTICAL PURPOSES ONLY.

# Data

We first import the data from proteinGroups.txt file. This is the file containing
maxLFQ summarized protein-level intensities. For a MaxQuant search [6],
this proteinGroups.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 proteinRaws file which is a subset
of the cptac study. 
To import the data we use the `QFeatures` package.

We generate the object proteinRawFile with the path to the proteinGroups.txt file.
Using the `grepEcols` function, we find the columns that contain the LFQ expression
data of the proteinRaws in the proteinGroups.txt file.



```{r, warning=FALSE, message=FALSE}
library(tidyverse)
library(limma)
library(QFeatures)
library(msqrob2)
library(plotly)
library(gridExtra)

proteinsFile <- "https://raw.githubusercontent.com/statOmics/PDA/data/quantification/cptacAvsB_lab3/proteinGroups.txt"

ecols <- grep("LFQ\\.intensity\\.", names(read.delim(proteinsFile)))

pe <- readQFeatures(
    table = proteinsFile, fnames = 1, ecol = ecols,
    name = "proteinRaw", sep = "\t"
)
```

In the following code chunk, we can extract the spikein condition from the raw file name.

```{r}
cond <- which(
  strsplit(colnames(pe)[[1]][1], split = "")[[1]] == "A") # find where condition is stored

colData(pe)$condition <- substr(colnames(pe), cond, cond) %>%
  unlist %>%  
  as.factor
```


We calculate how many non zero intensities we have per protein and this
will be useful for filtering.

```{r}
rowData(pe[["proteinRaw"]])$nNonZero <- rowSums(assay(pe[["proteinRaw"]]) > 0)
```


Proteins with zero intensities are missing and should be represent
with a `NA` value rather than `0`.
```{r}
pe <- zeroIsNA(pe, "proteinRaw") # convert 0 to NA
```


## Data exploration

`r format(mean(is.na(assay(pe[["proteinRaw"]])))*100,digits=2)`% of all peptide
intensities are missing and for some proteins we do not even measure a signal
in any sample.


# Preprocessing

This section preforms preprocessing for the peptide data. 
This include 

- log transformation, 
- filtering 

## Log transform the data

```{r}
pe <- logTransform(pe, base = 2, i = "proteinRaw", name = "proteinLog")
```

## Filtering

1. Remove reverse sequences (decoys) and contaminants

We now remove the contaminants and proteins that map to decoys.

```{r}
pe[["proteinLog"]] <- pe[["proteinLog"]][rowData(pe[["proteinLog"]])$Reverse != "+", ]
pe[["proteinLog"]] <- pe[["proteinLog"]][rowData(pe[["proteinLog"]])$
    Potential.contaminant != "+", ]
```
We keep `r nrow(pe[["proteinLog"]])` peptides upon filtering.


## Normalize the data using median centering 

We normalize the data by substracting the sample median from every intensity for peptide $p$  in a sample $i$: 

$$y_{ip}^\text{norm} = y_{ip} - \hat\mu_i$$ 

with $\hat\mu_i$ the median intensity over all observed peptides in sample $i$.

```{r}
pe <- normalize(pe, 
                i = "proteinLog", 
                name = "protein", 
                method = "center.median")
```


## Explore  normalized data

Upon the normalisation the density curves are nicely registered

```{r}
pe[["protein"]] %>% 
  assay %>%
  as.data.frame() %>%
  gather(sample, intensity) %>% 
  mutate(condition = colData(pe)[sample,"condition"]) %>%
  ggplot(aes(x = intensity,group = sample,color = condition)) + 
    geom_density()
```

We can visualize our data using a Multi Dimensional Scaling plot,
eg. as provided by the `limma` package.

```{r}
pe[["protein"]] %>% 
  assay %>%
  limma::plotMDS(col = as.numeric(colData(pe)$condition))
```

Note that the samples show a clear separation according to the spike-in condition in the second dimension of the MDS plot.

# Data Analysis

## Estimation

We model the protein level expression values using `msqrob`.
By default `msqrob2` estimates the model parameters using robust regression.

We will model the data with a different group mean. 
The group is incoded in the variable `condition` of the colData. 
We can specify this model by using a formula with the factor condition as its predictor: 
`formula = ~condition`.

Note, that a formula always starts with a symbol '~'.

```{r, warning=FALSE}
pe <- msqrob(object = pe, i = "protein", formula = ~condition)
```

## Inference

First, we extract the parameter names of the model by looking at the first model. 
The models are stored in the row data of the assay under the default name msqrobModels. 

```{r}
getCoef(rowData(pe[["protein"]])$msqrobModels[[1]])
```

We can also explore the design of the model that we specified using the the package `ExploreModelMatrix` 

```{r}
library(ExploreModelMatrix)
VisualizeDesign(colData(pe),~condition)$plotlist[[1]]
```

Spike-in condition `A` is the reference class. So the mean log2 expression
for samples from condition A is '(Intercept).
The mean log2 expression for samples from condition B is '(Intercept)+conditionB'.
Hence, the average log2 fold change between condition b and
condition a is modelled using the parameter 'conditionB'.
Thus, we assess the contrast 'conditionB = 0' with our statistical test.

```{r}
L <- makeContrast("conditionB=0", parameterNames = c("conditionB"))
pe <- hypothesisTest(object = pe, i = "protein", contrast = L)
```


## Plots

### Volcano-plot


```{r,warning=FALSE}
volcano <- ggplot(rowData(pe[["protein"]])$conditionB,
                  aes(x = logFC, y = -log10(pval), color = adjPval < 0.05)) +
  geom_point(cex = 2.5) +
  scale_color_manual(values = alpha(c("black", "red"), 0.5)) + theme_minimal()
volcano
```

Note, that only `r sum(rowData(pe[["protein"]])$conditionB$adjPval < 0.05, na.rm = TRUE)` proteins are found to be differentially abundant.

### Heatmap

We first select the names of the proteins that were declared signficant.

```{r}
sigNames <- rowData(pe[["protein"]])$conditionB %>%
  rownames_to_column("protein") %>%
  filter(adjPval<0.05) %>%
  pull(protein)
heatmap(assay(pe[["protein"]])[sigNames, ])
```

The majority of the proteins are indeed UPS proteins. 
1 yeast protein is returned. 
Note, that the yeast protein indeed shows evidence for differential abundance. 

### Boxplots

We make boxplot of the log2 FC and stratify according to the whether a protein is spiked or not.

```{r}
rowData(pe[["protein"]])$conditionB %>%
  rownames_to_column(var = "protein") %>%
  ggplot(aes(x=grepl("UPS",protein),y=logFC)) +
  geom_boxplot() +
  xlab("UPS") +
  geom_segment(
    x = 1.5,
    xend = 2.5,
    y = log2(0.74/0.25),
    yend = log2(0.74/0.25),
    colour="red") +
  geom_segment(
    x = 0.5,
    xend = 1.5,
    y = 0,
    yend = 0,
    colour="red") +
  annotate(
    "text",
    x = c(1,2),
    y = c(0,log2(0.74/0.25))+.1,
    label = c(
      "log2 FC Ecoli = 0",
      paste0("log2 FC UPS = ",round(log2(0.74/0.25),2))
      ),
    colour = "red")
```

What do you observe?

# Session Info

With respect to reproducibility, it is highly recommended to include a session info in your script so that readers of your output can see your particular setup of R. 

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