cancer3x3.Rmd

---
title: "Cancer"
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

Twelve Estrogen Receptor Positive Breast cancer tissues from from patients treated with tamoxifen upon recurrence have been assessed in a proteomics study. Six patients had a good outcome (OR) and the other Six had a poor outcome (PD). The proteomes have been assessed using an LTQ-Orbitrap and the thermo output .RAW files were searched with MaxQuant (version 1.4.1.2) against the human proteome database (FASTA version 2012-09, human canonical proteome).



# Data

We first import the data from peptide.txt file. This is the file containing
your peptide-level intensities. For a MaxQuant search [6],
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.

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


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

peptidesFile <- "https://raw.githubusercontent.com/statOmics/PDA22GTPB/data/quantification/cancer/peptides3vs3.txt"

ecols <- grep(
  "Intensity\\.", 
  names(read.delim(peptidesFile))
  )
```

Next, we read the data and store it in QFeatures object

```{r}
pe <- readQFeatures(
  table = peptidesFile,
  fnames = 1,
  ecol = ecols,
  name = "peptideRaw", sep="\t")
```

The QFeatures object pe currently contains a single assay, named peptideRaw.

We extract the column names from the peptideRaw assay and see that this contains information about the prognosis.

```{r}
colnames(pe[["peptideRaw"]])
```

We rename the colnames by dropping the "Intensity." from the name.

```{r}
(newNames <-  sub(
  pattern = "Intensity\\.", 
  replacement = "",
  colnames(pe[["peptideRaw"]]))
 )
```

```{r}
pe <- renameColname(pe, 
                    i = "peptideRaw",
                    newNames)
pe <- renamePrimary(pe, newNames)
colnames(pe[["peptideRaw"]])
```

In the following code chunk, we add the prognosis of the patients that we can read in the raw file name to the colData.

```{r}
colData(pe)$prognosis <-  
  colnames(pe[["peptideRaw"]]) %>%
  substr(start = 1, stop = 2) %>%
  as.factor
colData(pe)$prognosis
```


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

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


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

Look at the column names of the data to know the variables that you can use for filtering. 

```{r}
pe[["peptideRaw"]] %>% rowData %>% names
```
So we will filter on the "Reverse", "Contaminant" and "nNonZero" column.

## Data exploration

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


# Preprocessing

This section preforms preprocessing for the peptide data. 
This include 

- log transformation, 
- filtering and 
- summarisation of the data.

## Log transform the data

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

## Filtering

1. 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}
pe <- filterFeatures(pe, ~ Proteins %in% smallestUniqueGroups(rowData(pe[["peptideLog"]])$Proteins))
```

2. Remove reverse sequences (decoys) and contaminants

We now remove the contaminants and peptides that map to decoy sequences.

```{r}
pe <- filterFeatures(pe,~Reverse != "+")
pe <- filterFeatures(pe,~ Contaminant != "+")
```

3. Drop peptides that were only identified in one sample

We keep peptides that were observed at last twice.

```{r}
pe <- filterFeatures(pe,~ nNonZero >=2)
nrow(pe[["peptideLog"]])
```

We keep `r nrow(pe[["peptideLog"]])` 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 = "peptideLog", 
                name = "peptideNorm", 
                method = "center.median")
```


## Explore  normalized data

Upon the normalisation the density curves are nicely registered

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

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

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

The first axis in the plot is showing the leading log fold changes
(differences on the log scale) between the samples.


## Summarization to protein level

- By default robust summarization is used:  `fun = MsCoreUtils::robustSummary()`

```{r,warning=FALSE}
pe <- aggregateFeatures(pe,
  i = "peptideNorm",
  fcol = "Proteins",
  na.rm = TRUE,
  name = "protein")
```



```{r}
plotMDS(assay(pe[["protein"]]), col = as.numeric(colData(pe)$prognosis))
```

Note that the samples upon robust summarisation show a separation according to the prognosis.

# 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 `prognosis` of the colData. 
We can specify this model by using a formula with the factor condition as its predictor: 
`formula = ~prognosis`.

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

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

## 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),~prognosis)$plotlist
```

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

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


## Plots

### Volcano-plot


```{r,warning=FALSE}
volcano <- ggplot(rowData(pe[["protein"]])$prognosisPD,
                  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 `r sum(rowData(pe[["protein"]])$prognosisPD$adjPval < 0.05, na.rm = TRUE)` proteins are found to be differentially abundant.

### Heatmap

There were no proteins significantly differentially abundant at the 5% FDR level. 
We will therefore make an heatmap for the top 10 DE proteins. Note, however that they are not statistically significant! 

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

heatmap(assay(pe[["protein"]])[orderProt[1:10],])
```


### Detail plots

We make detail plots for the top 10 proteins. Note, however that their fold changes are not statistically significant at the 5% FDR level. 

```{r, warning=FALSE, message=FALSE}
#for (protName in sigNames)
for (protName in orderProt[1:10])
{
pePlot <- pe[protName, , c("peptideNorm","protein")]
pePlotDf <- data.frame(longFormat(pePlot))
pePlotDf$assay <- factor(pePlotDf$assay,
                        levels = c("peptideNorm", "protein"))
pePlotDf$prognosis <- as.factor(colData(pePlot)[pePlotDf$colname, "prognosis"])

# plotting
p1 <- ggplot(data = pePlotDf,
       aes(x = colname, y = value, group = rowname)) +
    geom_line() + 
    geom_point() +  
    theme(axis.text.x = element_text(angle = 70, hjust = 1, vjust = 0.5)) +
    facet_grid(~assay) + 
    ggtitle(protName)
print(p1)

# plotting 2
p2 <- ggplot(pePlotDf, aes(x = colname, y = value, fill = prognosis)) +
  geom_boxplot(outlier.shape = NA) + 
  geom_point(
    position = position_jitter(width = .1),
    aes(shape = rowname)) +
  scale_shape_manual(values = 1:nrow(pePlotDf)) +
  labs(title = protName, x = "sample", y = "peptide intensity (log2)") + 
  theme(axis.text.x = element_text(angle = 70, hjust = 1, vjust = 0.5)) +
  facet_grid(~assay)
print(p2)
}
```


# 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()
```