---
title: "Introduction to proteomics data analysis"
author:
- name: Lieven Clement
output:
  BiocStyle::html_document
---


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

# Data

We first import the peptides.txt file. This is the file that contains your peptide-level intensities. For a MaxQuant search [6], this peptides.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 raw files were saved. In this tutorial, we will use a MaxQuant peptides file from MaxQuant that can be found on the pdaData repository.
We will use the MSnbase package to import the data.

We generate the object peptideFile with the path to the peptides.txt file.
In this file we will fetch the data from the github repository linked to the pda course: https://statomics.github.io/SGA2019.
You can also replace the peptideFile with a string that points to the path of a file on your local hard drive.
With the grepEcols function we find the columns that are containing the expression data of the peptides in the peptides.txt file.

```{r, warning=FALSE, message=FALSE}
library(MSqRob)
library(MSnbase)
library(tidyverse)
library(limma)
peptidesFile <- "https://raw.githubusercontent.com/statOmics/SGA2019/data/quantification/cptacAvsB_lab3/peptides.txt"
proteinGroupsFile <-"https://raw.githubusercontent.com/statOmics/SGA2019/data/quantification/cptacAvsB_lab3/proteinGroups.txt"
ecols<-grepEcols(peptidesFile, "Intensity ", split = "\t")
```

Next we import the peptide intensities

```{r}
pepData<-readMSnSet2(peptidesFile,ecol=ecols,fnames="Sequence",sep="\t")
pepData
```

The pepData object is an MSnSet, a container for the data, features information and experimental annotation.
They can be accessed using the accessor functions 'exprs' (matrix of intensities, features in rows, samples in columns), 'fData' (properties for each feature, peptide or protein, on the rows) and 'pData' (properties for the samples on the columns).

We will make use from data wrangling functionalities from the tidyverse package.
The %>% operator allows us to pipe the output of one function to the next function.

```{r}
head(exprs(pepData))
pepData %>% exprs %>% head
pepData %>% sampleNames
```

The sample names are rather long and contain information on the spike-in concentration and the repeat.
We will remove "Intensity.6" from the string

```{r}
sampleNames(pepData) <- pepData %>% sampleNames %>% str_replace(., pattern="Intensity.6", replacement="")
```

Next we will add information on the proteins to the feature data.
The feature data of the object pepData contains information for each peptide in the experiment.
This info on the proteins can be found in the proteingroups.txt file of maxquant.
The peptides.txt file contains many data on each feature (peptide).
We will only retain specified columns Protein, Sequence, Reverse and Potential.Contaminant, indicating the protein, peptide sequence, if the peptide is a decoy and if it is a contaminant, respectively.

In older and more recent versions of maxquant the names of the fields might vary. We can therefore explore the names we want to select with the shiny app that can be launched from the selectFeatureData from the MSnbase package, i.e. by running the function without an fcol argument that specifies the columns you would like to select.

The code can be found below, but is not executed by this Rmd script (note, that the flag eval=FALSE has been used as option for the R chunk).

```{r,eval=FALSE}
pepData<-selectFeatureData(pepData)
```

If we know the columns, we can specify them as a vector in the fcol slot.

```{r}
pepData<-selectFeatureData(pepData,fcol=c("Proteins","Potential.contaminant","Reverse","Sequence"))
```

Calculate how many non zero intensities we have per peptide. This will be useful for filtering.
```{r}
fData(pepData)$nNonZero<- rowSums(exprs(pepData) > 0)
```

Peptides with zero intensities are missing peptides and should be represent with a `NA` value instead of `0`.
```{r}
exprs(pepData)[exprs(pepData)==0]<-NA
```


Next we create the data on the experimental layout.
We can do this based on the samplenames.
The information on the spike in condition is available as the first letter of the sample name.
We extract this information and recast the character vector in a factor.

```{r}
pData(pepData)$condition <- pepData %>% sampleNames %>% substr(1,1) %>% as.factor
```

## Data exploration

We can inspect the missingness in our data with the `plotNA()` function provided with `MSnbase`.
`r format(mean(is.na(exprs(pepData)))*100,digits=2)`% of all peptide intensities are missing and for some peptides we don't even measure a signal in any sample. The missingness is similar across samples.
Note, that we plot the peptide data, so the label protein in the plot refers to peptides.
```{r, cache= TRUE}
plotNA(pepData)
```

When we look at density plots of the intensities, we see that most intensities are low and that there is a long tail to the right.
It is more useful to explore the data upon log transformation.

```{r}
plotDensities(exprs(pepData))
```

```{r}
nrow(pepData)
```

There are `r nrow(pepData)` peptides before preprocessing.

# Preprocessing

We will log transform, normalize, filter and summarize the data.


## Log transform the data

```{r}
pepData <- log(pepData, base = 2)
plotDensities(exprs(pepData))
```

When we look at the densities we see small shifts in location and shape of the density distributions.
So normalisation will be required.
But, we will first perform filtering.

## Filtering

### Handling overlapping protein groups
In our approach a peptide can map to multiple proteins, as long as there is none of these proteins present in a smaller subgroup.

```{r}
pepData<-pepData[fData(pepData)$Proteins %in% smallestUniqueGroups(fData(pepData)$Proteins),]
```

### Remove reverse sequences (decoys) and contaminants

We now remove the contaminants, peptides that map to decoy sequences and proteins, which were only identified by peptides with modifications.

```{r}
pepData <- pepData[fData(pepData)$Reverse!="+",]
pepData <- pepData[fData(pepData)$Potential.contaminant!="+",]
```

### Remove peptides of proteins that were only identified with modified peptides

Proteins for which all peptides are carrying modifications (PTMs) can be considered as unreliable.
We will filter out these proteins.
This information is included in the Only.identified.by.site column of proteinGroups.txt file of maxquant. The variable is containing a "+" character if the protein is only identified by modified peptides and an empty string "" if this is not the case.
Sometimes an NA value is also present. We will replace these by the empty character string.
The information of the peptides.txt file can be linked to the proteinGroups.txt file by using the Proteins column from the peptides.txt file and the Protein.IDs column in the proteinGroups.txt file.
If NAs are included in the

```{r}
proteinGroups<-read.table(proteinGroupsFile, sep = "\t", header = TRUE, quote = "", comment.char = "")
only_site <- proteinGroups$Only.identified.by.site
only_site[is.na(only_site)] <- ""
proteinsOnlyIdentifiedWithPTMs <- proteinGroups$Protein.IDs[only_site=="+"]
pepData<-pepData[!(fData(pepData)$Proteins %in% proteinsOnlyIdentifiedWithPTMs),]
```

We will now remove the proteinGroups object from the R session to free memory.

```{r}
rm(proteinGroups,only_site)
gc()
```

### Drop peptides that were only identified in one sample

We want to keep peptide that were at least observed twice.

```{r}
pepData<-pepData[fData(pepData)$nNonZero>=2,]
nrow(pepData)
```

We keep `r nrow(pepData)` peptides upon filtering.  


## Quantile normalize the data
```{r}
pepData <- normalise(pepData, "quantiles")
```


## Explore quantile normalized data

Upon normalisation the density curves for all samples coincide.

```{r}
plotDensities(exprs(pepData))
```

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

```{r}
plotMDS(exprs(pepData),col=as.double(pData(pepData)$condition))
```

The first axis in the plot is showing the leading log fold changes (differences on the log scale) between the samples.
We notice that the leading differences (log FC) in the peptide data seems to be driven by technical variability.
Indeed the samples do not seem to be clearly separated according to the spike in condition.
Because there are missing values for the peptide level data we use the `na.rm=TRUE` argument to summarize the data based on the observed peptides for every protein.

## Summarization to protein level
```{r,warning=FALSE}
protData<-combineFeatures(pepData,fcol="Proteins",method="robust",na.rm=TRUE)
```

If you summarize using robust summarization please refer to Sticker et al. (2019) https://www.biorxiv.org/content/10.1101/668863v1

We notice that the leading differences (log FC) in the protein data is still according to technical variation. On the second dimension, however, we also observe a clear separation according to the spike-in condition. Hence, the summarization that accounts for peptide specific effects makes the effects due to the spike-in condition more prominent!

```{r}
plotMDS(exprs(protData),col=as.double(pData(pepData)$condition))
```


# Data Analysis

## Estimation
MSqRob is currently working with a format where we have one dataframe for each protein.
This will be changed in the next release.
Therefore we first have to reorganise the data.

Next the models are fitted.
This is done using the fit.model function.
We only have to model the data using the factor condition from the pheno data of the protein level MSnSet.
The name of the factor variable is specified in the fixed argument (if multiple predictors have to be incorporated in the model, a vector of variable names has to be provided in this argument.).
The argument shrinkage is used to specify if ridge regression has to be adopted.
For the sake of speed we do not do this in the tutorial.
The shrinkage has to be specified for each variable in the fixed effects.
We also have to indicate this for the intercept (which we never shrink).
So we specify it at c(0,0) to indicate that the intercept (first 0) and the parameters for the factor condition (second 0) are not penalized.
We set the robust_var function equal to FALSE, this functionality will be removed from the package in the next release.


```{r}
protMSqRob <- MSnSet2protdata(protData, "Proteins")
models <- fit.model(protdata=protMSqRob, response="quant_value", fixed="condition",shrinkage.fixed=c(0,0),robust_var=FALSE)
```

## Inference
Many biologists have problems with the reference coding.
In MSqRob we have opted to formulate contrasts using all levels of a factor.
Internally, the contrasts are than recasted according to the factor level that is the reference class.

```{r}
L <- makeContrast("conditionB - conditionA",  levels=c("conditionA","conditionB"))
tests <- test.contrast_adjust(models, L)
nSig <- sum(tests$signif,na.rm=TRUE)
head(tests,nSig)
```

There are `r nSig` proteins with a significant effect at the 5% FDR level.

## Plots

### Volcano-plot


```{r,warning=FALSE}
volcano<- ggplot(tests,aes(x=estimate,y=-log10(pval),color=signif)) + geom_point() + scale_color_manual(values=c("black","red"))
volcano
```

An interactive volcano plot can be obtained by using the plotly library.
We alter the object tests to add the protein name to the data frame in ggplot, this will allow us to see the protein name in the interactive plot.

```{r,eval=FALSE}
library(plotly)
volcano<- ggplot(tests %>% rownames_to_column("protein"),aes(x=estimate,y=-log10(pval),color=signif,protein=protein)) + geom_point() + scale_color_manual(values=c("black","red"))
ggplotly(volcano)
```

### Heatmap

We first select the names of the significant proteins.

```{r}
sigNames<- tests %>% rownames_to_column("protein") %>% filter(signif) %>% pull(protein)
heatmap(exprs(protData)[sigNames,])
```

### Detail plot for the top 5 proteins.

We first extract the normalized protein expression values for a particular protein and we create a dataframe that contains the protein expression values and the relevant variables of the experimental design.  

```{r, warning=FALSE, message=FALSE}
for (protName in rownames(tests)[1:5])
{
plotDat<-pepData[fData(pepData)$Proteins==protName] %>% exprs
plotDataStack<-data.frame(quant_value=c(plotDat),sequence=rep(rownames(plotDat),rep=ncol(plotDat)),sample=rep(colnames(plotDat),each=nrow(plotDat)),condition=rep(pData(pepData)$condition,each=nrow(plotDat)))
plot1 <- ggplot(plotDataStack, aes(x=sample, y=quant_value,fill=condition))
print(plot1 +geom_boxplot(outlier.shape=NA) + geom_point(position=position_jitter(width=.1),aes(shape=sequence)) + scale_shape_manual(values=1:nrow(plotDat)) +labs(title = protName, x="sample", y="Peptide intensity (log2)"))
}
```

### Sensitivity FDP plot

Because we know the ground truth for the cptac study, i.e. we know that only the spike-in proteins (UPS) are differentially expressed, we can evalute the fold changes.
Yeast proteins should be not differentially expressed and their log fold changes should be centered around 0. These of UPS proteins are spiked at differt concentrations and their log2 fold changes should be centered around log2(concB/concA)=log2(0.74/0.25)=`r log2(0.74/0.25)`.

We will first create a novel factor in the results table that is called `spike` to indicate which proteins are spiked.

```{r,warning=FALSE}
tests$spike <- tests %>% rownames %>% grepl(.,pattern="UPS")
ggplot(tests, aes(x=spike, y=estimate)) +
  geom_boxplot() +
  geom_hline(yintercept=c(0,log(0.74/0.25,base=2)),color="red")+
  ylab("log2 FC")
```

### Sensitivity FDP plot

Because we know the ground truth for the cptac study, i.e. we know that only the spike-in proteins (UPS) are differentially expressed, we can calculate

- the sensitivity or true positive rate (TPR), the proportion of actual positives that are correctly identified, in the protein list that we return
$$TPR=\frac{TP}{\text{#actual positives}},$$
here TP are the true positives in the list. The TPR is thus the fraction of ups proteins that we can recall.

- false discovery proportion (FPD): fraction of false positives in the protein list that we return:
$$FPD=\frac{FP}{FP+TP},$$
with FP the false positives. In our case the yeast proteins that are in our list.

Instead of only calculating that for the protein list that is returned for the chosen FDR level, we can do this for all possible FDR cutoffs so that we get an overview of the quality of the ranking of the proteins in the protein list.

```{r}
tests<-tests %>% rownames_to_column("protein") %>%
mutate(
FDP=cumsum(!spike)/(1:length(spike)),
TPR=cumsum(spike)/sum(spike)) %>%
column_to_rownames("protein")

ggplot(tests,aes(x=FDP,y=TPR)) + geom_path() +geom_vline(xintercept=0.05,lty=2) + geom_point(data=tests[sum(tests$signif,na.rm=TRUE),],aes(x=FDP,y=TPR),cex=2)
```

We observe that there is a good FDR control: the FDP at the 5% FDR level is indeed close to 0.05.
