---
title: "Platelet proteomics analysis"
date: "`r format(Sys.time(), '%d %B, %Y')`"
author: "<br><br>Rodolphe Poupardin"
output:
  html_document:
    toc: true
    code_folding: hide
    toc_float:
      smooth_scroll: false
      
---


This is the pipeline used for the label free platelets proteomics data.

# Libraries and functions used

```{r message=FALSE, warning=FALSE, paged.print=FALSE}
library("DEP")
library("dplyr")
library("SummarizedExperiment")
library("openxlsx2")
library("ggplot2")
library("ggrepel")
library("ggiraph")
library("reactable")
library("crosstalk")
library("GOplot")
library("pathview")
library("clusterProfiler")
library("topGO")
library("GSEABase")
library("ReactomePA")
library("RColorBrewer")
library("org.Hs.eg.db")
library("STRINGdb")



```


# Data preparation 

Processed Maxquant files *"proteinGroups_Drerup.txt"* were used for the downstream analysis. The dataset has unique Uniprot identifiers, however those are not immediately informative. The associated gene names are informative, however these are not always unique. For further analysis these proteins must get unique names. Additionally, some proteins do not have an annotated gene name and for those we will use the Uniprot ID.

```{r, message= FALSE, warning= FALSE}


data <- read.delim("proteinGroups_AD_patients.txt", h = TRUE, stringsAsFactors = FALSE)


## We filter for contaminant proteins and decoy database hits, which are indicated by "+" in the columns "Potential.contaminants" and "Reverse", respectively. 


data <- data[data$Reverse != "+",]



# Make unique names using the annotation in the "Gene.names" column as primary names and the annotation in "Protein.IDs" as name for those that do not have an gene name.
data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")


### And the first batch of data
data_OLD <- read.delim("proteinGroups_Drerup.txt", h = TRUE, stringsAsFactors = FALSE)




## We filter for contaminant proteins and decoy database hits, which are indicated by "+" in the columns "Potential.contaminants" and "Reverse", respectively. 
data_OLD <- data_OLD[data_OLD$Reverse != "+" & data_OLD$Potential.contaminant != "+",]

# Make unique names using the annotation in the "Gene.names" column as primary names and the annotation in "Protein.IDs" as name for those that do not have an gene name.
data_OLD_unique <- make_unique(data_OLD, "Gene.names", "Protein.IDs", delim = ";")

# Remove bad samples/outliers (based on previous analysis)
data_OLD_unique$LFQ.intensity.CO6 <- NULL
data_OLD_unique$LFQ.intensity.CO7  <- NULL
data_OLD_unique$LFQ.intensity.CO8 <- NULL
data_OLD_unique$LFQ.intensity.CY4 <- NULL

```


# Merge both table

```{r}
#data_unique <- merge(data_unique, data_OLD_unique, by = "Protein.IDs", all.x = TRUE, all.y = TRUE)
data_unique <- merge(data_unique, data_OLD_unique, by = "Protein.IDs")

data_unique$Gene.names.x <- ifelse(is.na(data_unique$Gene.names.x) == TRUE, data_unique$Gene.names.y, data_unique$Gene.names.x)

data_unique <- make_unique(data_unique, "Gene.names.x", "Protein.IDs", delim = ";")
```

# Generate a SummarizedExperiment object

Many Bioconductor packages use SummarizedExperiment objects as input and/or output. This class of objects contains and coordinates the actual (assay) data, information on the samples as well as feature annotation. 

```{r}
LFQ_columns <- grep("LFQ.", colnames(data_unique)) 

experimental_design <- data.frame(label = colnames(data_unique[,LFQ_columns]))


experimental_design$condition <- experimental_design$label
experimental_design$condition <- gsub("LFQ.intensity.CO.*", "Old", experimental_design$condition)
experimental_design$condition <- gsub("LFQ.intensity.CY.*", "Young", experimental_design$condition)

experimental_design$condition <- gsub("LFQ.intensity.AD.*", "AD", experimental_design$condition)



experimental_design$replicate <- gsub("LFQ.intensity.|LFQ.intensity.", "", experimental_design$label)

experimental_design$replicate <- gsub("CO_|AD_|AD|CO", "", experimental_design$replicate)
experimental_design$replicate <- gsub(".2020|_2020", "", experimental_design$replicate)


#convert to character
experimental_design$label <- as.character(experimental_design$label)

experimental_design$label <- gsub("\\.[xy]|LFQ.intensity.", "", experimental_design$label)

# get extra inbformation

extra_info_participants <- read.delim("Feb_2023/APPA-participants_02.2023.csv", h = T, sep = ",")

experimental_design <- merge(experimental_design, extra_info_participants, by.x = 1, by.y = 1, all.x = T)

# Make SummarizedExperiment
data_se <- make_se(data_unique, LFQ_columns, experimental_design)

reactable(experimental_design, minRows = 19, pagination = FALSE, style = list(fontFamily = "Arial", fontSize = "14px"))
  

```

# Filter on missing values

The dataset contains proteins which are not quantified in all replicates. Some proteins are even only quantified in a single replicate.

```{r}
plot_frequency(data_se)
```

This leaves our dataset with missing values, which need to be imputed. However, this should not be done for proteins that contain too many missing values. Therefore, we first filter out proteins that contain too many missing values. This is done by setting the threshold for the allowed number of missing values per condition in the filter_missval function.

```{r}
# Filter for proteins that have a maximum of four missing values in one condition
data_filt2 <- filter_missval(data_se, thr = 4)
```

After filtering, the number of identified proteins per sample can be plotted as well as the overlap in identifications between samples.

```{r}
# Plot a barplot of the number of identified proteins per samples
plot_numbers(data_filt2)
```

Plot a barplot of the protein identification overlap between samples
```{r}
# Plot a barplot of the protein identification overlap between samples
plot_coverage(data_filt2)
```

# Normalization

mass spectrometry (MS) data remain inherently biased as a result of reasons ranging from sample handling to differences caused by the instrumentation. Normalization is the process that aims to account for the bias and make samples more comparable. The selection of a proper normalization method is a pivotal task for the reliability of the downstream analysis and results. The data is background corrected and normalized by variance stabilizing transformation (vsn).

```{r, message=FALSE}
# Normalize the data
data_norm <- normalize_vsn(data_filt2)
```

The normalization can be inspected by checking the distributions of the samples before and after normalization.

```{r, fig.height= 8, fig.width = 7, fig.align = "center"}
# Visualize normalization by boxplots for all samples before and after normalization
plot_normalization(data_filt2, data_norm)
```

Standard deviation and mean are calculated row-wise from the expression matrix (in) x. The scatterplot of these versus each other allows you to visually verify whether there is a dependence of the standard deviation (or variance) on the mean. The red line depicts the running median estimator (window-width 10%). If there is no variance-mean dependence, then the line should be approximately horizontal which seems to be our case here.

```{r}
meanSdPlot(data_norm)
```


# Impute data for missing values

The remaining missing values in the dataset need to be imputed. The data can be missing at random (MAR), for example if proteins are quantified in some replicates but not in others. Data can also be missing not at random (MNAR), for example if proteins are not quantified in specific conditions (e.g. in the control samples). MNAR can indicate that proteins are below the detection limit in specific samples, which could be very well the case in proteomics experiments. 
To explore the pattern of missing values in the data, a heatmap is plotted indicating whether values are missing (0) or not (1). Only proteins with at least one missing value are visualized.

```{r, fig.height= 8, fig.width = 7, fig.align = "center"}
# Plot a heatmap of proteins with missing values
plot_missval(data_filt2)
```

This heatmap indicates that missing values are highly biased to specific samples. To check whether missing values are biased to lower intense proteins, the densities and cumulative fractions are plotted for proteins with and without missing values.

```{r}
# Plot intensity distributions and cumulative fraction of proteins with and without missing values
plot_detect(data_filt2)
```

Indeed the proteins with missing values have on average low intensities. This data (MNAR and close to the detection limit) should be imputed by a left-censored imputation method, such as the quantile regression-based left-censored function (“QRILC”) or random draws from a left-shifted distribution (“MinProb” and “man”). Here, the "QRILC" function was used.


```{r}
data_imp <- impute(data_norm, fun = "QRILC") 
```

The effect of the imputation on the distributions can be visualized.
```{r}
# Plot intensity distributions before and after imputation
plot_imputation(data_norm, data_imp)
```

# Differential enrichment analysis

Protein-wise linear models combined with empirical Bayes statistics are used for the differential enrichment analysis (or differential expression analysis). The test_diff function introduced here uses limma and automatically generates the contrasts to be tested. For the contrasts generation, the control sample has to be specified.

```{r}
# Differential enrichment analysis  based on linear models and empherical Bayes statistics

# Test every sample versus control 
## Chose 
data_diff <- test_diff(data_imp, type = "control", control = "Old")
# OR
data_diff <- test_diff(data_imp, type = "control", control = "Young")


data_diff <- test_diff(data_imp, type = "all", design_formula = formula(~ 0 + condition + stress))
```

Finally, significant proteins are defined by user-defined cutoffs using add_rejections.

```{r}
# Denote significant proteins based on user defined cutoffs
dep <- add_rejections(data_diff, alpha = 1, lfc = log2(1))
```

# PCA plot

The PCA plot can be used to get a high-level overview of the data. This can be very useful to observe batch effects, such as clear differences between replicates.

```{r, fig.height= 6, fig.width = 10, fig.align = "center", message=FALSE, warning=FALSE}
# Plot the first and second principal components
plot_pca(dep, x = 1, y = 2, n = 100, point_size = 4, indicate = c("condition"), label = TRUE)

```

Most of the variance (PC1) is caused by some outliers. These outliers are the samples with more missing values. 
> Recommendation : remove the 5 samples out of the analysis

# Remove outliers and reconduct previous steps

Five outliers removed from analysis:<br>
- AD10<br>
- AD3<br>
- AD7<br>
- CO7<br>
- CO8<br>

```{r, fig.height= 6, fig.width = 10, fig.align = "center", message=FALSE, warning=FALSE}
# Remove bad samples/outliers (based on previous analysis)
data_unique$LFQ.intensity.AD10 <- NULL
data_unique$LFQ.intensity.AD3 <- NULL
data_unique$LFQ.intensity.AD7 <- NULL
data_unique$LFQ.intensity.CO7 <- NULL
data_unique$LFQ.intensity.CO8 <- NULL


LFQ_columns <- grep("LFQ.", colnames(data_unique))


LFQ_columns <- grep("LFQ.", colnames(data_OLD_unique))



experimental_design <-
  data.frame(label = colnames(data_unique[, LFQ_columns]))

experimental_design <-
  data.frame(label = colnames(data_OLD_unique[, LFQ_columns]))


experimental_design$condition <- experimental_design$label
experimental_design$condition <-
  gsub("LFQ.intensity.CO.*", "Old", experimental_design$condition)
experimental_design$condition <-
  gsub("LFQ.intensity.CY.*", "Young", experimental_design$condition)

experimental_design$condition <-
  gsub("LFQ.intensity.AD.*", "AD", experimental_design$condition)


experimental_design$replicate <-
  gsub("LFQ.intensity.|LFQ.intensity.",
       "",
       experimental_design$label)

experimental_design$replicate <-
  gsub("CO_|AD_|AD|CO", "", experimental_design$replicate)
experimental_design$replicate <-
  gsub(".2020|_2020", "", experimental_design$replicate)

#convert to character
experimental_design$label <- as.character(experimental_design$label)
experimental_design$label <- gsub("\\.[xy]|LFQ.intensity.", "", experimental_design$label)

# get extra information

extra_info_participants <- read.delim("APPA-participants.csv", h = T, sep = ",")

experimental_design <- merge(experimental_design, extra_info_participants, by.x = 1, by.y = 1, all.x = T)

experimental_design$label <- colnames(data_unique[, LFQ_columns])
experimental_design$label <- colnames(data_OLD_unique[, LFQ_columns])

# Make SummarizedExperiment
data_se <- make_se(data_unique, LFQ_columns, experimental_design)
data_se <- make_se(data_OLD_unique, LFQ_columns, experimental_design)


# Filter for proteins that have a maximum of three missing values in one condition
data_filt2 <- filter_missval(data_se, thr = 4)

# Normalize the data
data_norm <- normalize_vsn(data_filt2)

data_imp <- impute(data_norm, fun = "QRILC")

data_diff <- test_diff(data_imp, type = "control", control = "Old")
data_diff <- test_diff(data_imp, type = "control", control = "Young")

dep <- add_rejections(data_diff, alpha = 1, lfc = log2(1))

# Plot the first and second principal components
PCA_PLOT <- plot_pca(
  dep,
  x = 1,
  y = 2,
  n = 197,
  point_size = 4,
  indicate = c("condition", "Gender"),
  label = F,
)


# Other PCA plot

var <- apply(assay(dep), 1, sd)
df <- assay(dep)[order(var, decreasing = TRUE)[seq_len(197)],]

pca <- prcomp(t(df), scale = FALSE)
pca_df <- pca$x %>% data.frame() %>% tibble::rownames_to_column() %>% left_join(., data.frame(colData(dep)), by = c(rowname = "ID"))
percent <- round(100 * pca$sdev^2/sum(pca$sdev^2), 1)


ggplot(pca_df, aes(x = PC1, y = PC2, colour = condition, shape = Gender, size  = (Age), label = condition)) +
  geom_point() +
  xlab(paste0("PC1: ", percent[1], "% variance")) +
  scale_colour_manual(values = c("#00BFC4", "#7CAE00")) +
  ylab(paste0("PC2: ", percent[2], "% variance")) +
  coord_fixed() +
  theme_minimal() +
  labs(condition = "Condition", Gender = "Sex", size="Age")


```

After removing the outliers that had a lower number of proteins detected, most of the variance is now explained by the comparison AD / Old <br>


# 3D PCA

```{r setup}
library(knitr)
library(rgl)
knit_hooks$set(webgl = hook_webgl)
```

```{r testgl, webgl=TRUE}
x <- sort(rnorm(1000))
y <- rnorm(1000)
z <- rnorm(1000) + atan2(x,y)
plot3d(x, y, z, col=rainbow(1000))
```


# Correlation matrix

A correlation matrix can be plotted as a heatmap, to visualize the Pearson correlations between the different samples.

```{r}
# Plot the Pearson correlation matrix
plot_cor(dep, significant = FALSE, lower = 0, pal = "Greens")
```



Samples cluster nicely together according to their age category


# Heatmap of all significant proteins


```{r, fig.height= 10, fig.width= 7}
# Plot a heatmap of all significant proteins with the data centered per protein
plot_heatmap(dep, type = "centered", kmeans = TRUE, 
             k = 2, col_limit = 4, show_row_names = FALSE,
             indicate = c("condition"))

```

# Volcano Plots

```{r, fig.height= 8, fig.width= 8, warning = FALSE, message= FALSE}
data_results <- get_results(dep)
results = mutate(data_results,
                 sig = ifelse(data_results$AD_vs_Old_p.val < 0.05, "Pval < 0.05", "Not Sig"))
results$size = (-log10(results$AD_vs_Old_p.val)) * (abs(results$AD_vs_Old_ratio))

p1 = ggplot(results, aes(AD_vs_Old_ratio,-log10(AD_vs_Old_p.val))) +
  geom_point_interactive(aes(
    tooltip = name,
    data_id = name,
    col = sig,
    size = size
  ), alpha = 0.4) +
  scale_color_manual(values = c("grey60", "red")) +
  #ylim(0, 20) +
  #xlim(-7, 7) +
  ggtitle("Old vs young") +
  theme(legend.position = "bottom") +
  theme_bw() +
  guides(alpha = "none", size = "none") +
  labs(col = "",
       size = "",
       x = "Log2 fold change AD vs Old",
       y = "-Log10(P-value)") +
  theme(legend.position = "bottom")


p2 <- p1 +
  geom_text_repel(
    data = dplyr::filter(results, sig == "Pval < 0.05" &
                           abs(AD_vs_Old_ratio) > 2),
    aes(label = name),
    col = "grey10",
    alpha = 0.9,
    max.overlaps = 30
  )

girafe(
  ggobj = p2,
  width_svg = 8,
  height_svg = 8,
  options = list(opts_hover(css = "stroke-width:10;"))
)


```





#### Figures for Manuscript 2023

# Volcano Plots AD vs OLD

```{r, fig.height= 8, fig.width= 8, warning = FALSE, message= FALSE}

data_diff <- test_diff(data_imp, type = "control", control = "Old")
#data_diff <- test_diff(data_imp, type = "control", control = "Young")

dep <- add_rejections(data_diff, alpha = 1, lfc = log2(1))

data_results <- get_results(dep)
results = mutate(data_results,
                 sig = ifelse(data_results$AD_vs_Old_p.val < 0.05, "Pval < 0.05", "Not Sig"))
results$size = (-log10(results$AD_vs_Old_p.val)) * (abs(results$AD_vs_Old_ratio))

p1 = ggplot(results, aes(AD_vs_Old_ratio,-log10(AD_vs_Old_p.val))) +
  geom_point_interactive(aes(
    tooltip = name,
    data_id = name,
    col = sig,
    size = size
  ), alpha = 0.4) +
  scale_color_manual(values = c("grey60", "red")) +
  #ylim(0, 20) +
  #xlim(-7, 7) +
  ggtitle("AD vs Old") +
  theme(legend.position = "bottom") +
  theme_bw() +
  guides(alpha = "none", size = "none") +
  labs(col = "",
       size = "",
       x = "Log2 fold change AD vs Old",
       y = "-Log10(P-value)") +
  theme(legend.position = "bottom")


p2 <- p1 +
  geom_text_repel(
    data = dplyr::filter(results, sig == "Pval < 0.05" &
                           abs(AD_vs_Old_ratio) > 2),
    aes(label = name),
    col = "grey10",
    alpha = 0.9,
    max.overlaps = 30
  )

girafe(
  ggobj = p2,
  width_svg = 8,
  height_svg = 8,
  options = list(opts_hover(css = "stroke-width:10;"))
)


```





