---
title: "Proteomics Analysis"
output:
  html_document:
    df_print: paged
editor_options:
  chunk_output_type: inline
---



```{r}
# Attach the different packages
suppressPackageStartupMessages(library(fontawesome)) # package to add icons in the html report
suppressPackageStartupMessages(library(tidyverse)) # collection of R packages for data science.
suppressPackageStartupMessages(library(DT)) # datatable function
suppressPackageStartupMessages(library(seqinr)) # trimSpace function
suppressPackageStartupMessages(library(edgeR)) # TMM normalisation
suppressPackageStartupMessages(library(timecourse))  # mb.long function

```


```{r config, echo=FALSE, warning=FALSE, message=FALSE}
# Define the path of the project
project.path <- "PROJECT PATH"

# Path to the Proteome Discoverer results
txt.path <- file.path(project.path, "PD") 

```

```{r Functions}
# Function for the Sample Loading Normalization
SL <- function(intensities){
    # Normalize for each TMT run
    target <- mean(colSums(intensities, na.rm = TRUE), na.rm = TRUE)
    
    # do the sample loading normalization
    # there is a different correction factor for each column
    norm_factors<- target / colSums(intensities, na.rm = TRUE)
    
    # SL Normalization using the correction factor
    SL <- sweep(intensities, 2, norm_factors, FUN = "*")
    
    print(data.frame(Intensity.Sum.perChannel = trimSpace(format(round(colSums(SL), digits = 0), big.mark = ","))))
    
    return(SL)
}

# Function to extract the Gene name from the description column
get_Gene_name_from_Description <- function(x){
  # Look for the "GN=" string and get the following string
   if(grepl("GN", x)){
        desc <- strsplit(x, split=" ")[[1]]
          n <- desc[grep("GN=",desc)[1]]
          n <- sub("GN=", "", n)
   } else { 
       n <- "" # If no "GN=" found
   }
      return(n)

}

# Function to extract the Protein description from the description column
get_ProteinDesc_from_Description <- function(x){
  # Look for the "OS=" string and get the following string
   if(grepl("OS", x)){
        n <- sub("(.*) OS=.*", "\\1", x)
   } else {
      n <- "" # If no "OS=" found
   }
   return(n)
}


# Function to get the Limma results and add a status column (used in the volcano plots)
get_results <- function(limma_Results, d, Protein.IDs, Gene.names) {

        results <- d
  
        # Add the column of Accessions
        results$Accession <- Protein.IDs

        # get the log Fold-Change, the P values and FDR from the limma test
        results$logFC <- limma_Results$logFC
        results$p_value <- limma_Results$P.Value
        results$fdr <- limma_Results$adj.P.Val

        # add a DE candidate status column
        results$candidate <- "FDR > 0.05"
        results$candidate[which(results$fdr <= 0.05 & results$fdr > 0.01)] <- "FDR <= 0.05"
        results$candidate[which(results$fdr <= 0.01)] <- "FDR <= 0.01"
        results$candidate <- factor(results$candidate, levels = c( "FDR <= 0.01", "FDR <= 0.05", "FDR > 0.05"))

        #Add gene names and protein description
        results <- results %>% dplyr::left_join(Gene.names, by = "Accession")
        
        return(results)
}


```



# Labels

```{r, fig.width=15, fig.height=10}

# Get the TMT Labels information
Labels <- read.table("../Labels.txt", header=TRUE, sep="\t", stringsAsFactors = FALSE, comment="", quote="\"")
Labels <- Labels %>% unite(TMT_Label, c(TMT, Label))
Labels <- Labels %>% separate(Biological.name, c("Day", "Replicate"), sep="_", remove=FALSE) 
datatable(Labels)

```


```{r}
# Define the extension of the file names from Proteome Discoverer
Extension <- "VT_210305_Lumos_TMT_B"

# Read the Protein text file from Proteome Discoverer and create a dataframe
Proteins <- read.table(file.path(txt.path, paste(Extension, "Proteins.txt", sep="_")), header=TRUE, sep="\t", stringsAsFactors = FALSE, comment="", quote="\"")

 # Remove contaminants and keep proteins with a high confidence
Proteins <- Proteins %>% filter(Contaminant == "False" &  (Protein.FDR.Confidence.Combined == "High"| Protein.FDR.Confidence.Combined == "Medium") )

# Group the proteins into protein groups
Proteins <- Proteins %>% filter(Protein.Group.IDs != "") %>%
                        group_by(Protein.Group.IDs) %>% 
                        mutate(Protein.Group.Accessions = paste0(Accession, collapse = ";")) %>% 
                        ungroup() %>% filter(Master == "IsMasterProtein") %>% 
                          mutate(Gene.names = pmap_chr(list(x = Description), get_Gene_name_from_Description),
                                 Protein.Desc = pmap_chr(list(x = Description), get_ProteinDesc_from_Description))




```


```{r}
# Keep only the Abundance columns
TMT <- Proteins %>% dplyr::select( matches("Abundance")) 

# Add the accession information
TMT$Accession <- Proteins$Accession 

# Add the Master information
TMT$Master <- Proteins$Master 

# Keep only the proteins that are considered as Master proteins by Proteome Discoverer
TMT <- TMT %>% filter(Master == "IsMasterProtein")

```



# Normalisation 

## Sample Loading Normalisation 

<div id="info">
<div id= "info-icon">`r fa("info-circle", fill = "steelblue", height="18px")`</div>
<div id= "info-text">

The Sample Loading normalisation step scales each TMT channel so that its sum of reporter ion signals equals the average grand total across samples.

</div>
</div>


```{r}
# Get the abundances in a dataframe
corrected.int.only <- TMT %>% dplyr::select(-Master)

# Add the accession names as row names
row.names(corrected.int.only) <- NULL
corrected.int.only <- corrected.int.only %>% column_to_rownames("Accession")

# Remove the missing values
corrected.int.only[ corrected.int.only == 0 ] <- NA    
corrected.int.only <- corrected.int.only %>% na.omit()   

# Sample Loading Normalisation
SL.corrected.int <- SL(corrected.int.only)

# Add the Accessions
row.names(SL.corrected.int) <- row.names(corrected.int.only)

```


## TMM normalisation on SL normalized data

<div id="info">
<div id= "info-icon">`r fa("info-circle", fill = "steelblue", height="18px")`</div>
<div id= "info-text">

Trimmed M-Mean normalisation assumes that the samples mostly have a lot of common proteins that are not differentially expressed. High and low abundance proteins are removed, and larger and smaller fold-changes are also trimmed. Then, the remaining proteins should be similar in expression. TMM computes single global scaling factors that move those presumed unchanged proteins to more similar values. The TMM routine does include the equivalent of an SL norm step. (Phillip Wilmarth)

</div>
</div>

```{r}

# TMM Normalisation from the edge package
tmm_sl_normFactors <- calcNormFactors(SL.corrected.int, method="TMM")
TMM.SL.corrected.int <- sweep(SL.corrected.int, 2, tmm_sl_normFactors, FUN = "/") # this is data after SL on original scale

# Add the row names as a dataframe column
TMM.SL.corrected.int.2Write <- TMM.SL.corrected.int %>% rownames_to_column("Accession")
write.table(TMM.SL.corrected.int.2Write, file.path("TMM_SL_normalised_Intensities.txt"), sep="\t", row.names = FALSE)

```


```{r}
# EITHER
# Get the normalized data from the text file previously written to keep track (if new .rmd file for example)
# Norm_TMTs <- read.table(file.path(txt.path, "TMM_SL_normalised_Intensities.txt"), header=TRUE, sep="\t", stringsAsFactors = FALSE, comment="", quote="\"")

# OR keep going with the previous dataframe
Norm_TMTs <- TMM.SL.corrected.int.2Write

```

```{r}
# # 
# x = 2
# y = 6.5
# 
# xlim = c(-x, x)
# ylim= c(0, y)
```


```{r}
# Add human readable condition names and transform in log2
# --------------------------------------------------------

# Pivot the dataframe in a longer format
TMT.pivot <- Norm_TMTs %>% pivot_longer(cols = c(2:10), names_to = "Sample", values_to = "Abundances" ) %>% separate(Sample, c("A", "F", "Label", "S", "TMT"), sep="[.]", remove=FALSE) %>% dplyr::select(-A,  -S, -F)  

# Transform the raw abundances in log2
TMT.pivot$Abundances <- log2(TMT.pivot$Abundances)

# Remove the string "TMT" from the TMT column
TMT.pivot$TMT = sub("TMT", "", TMT.pivot$TMT)

#Add a variable to do the bridge with the normalized reporter abundances of the dataframe "TMTs" and the "Labels" dataframe
TMT.pivot <- TMT.pivot %>% unite(TMT_Label, c(TMT, Label)) %>% dplyr::left_join(Labels, by = "TMT_Label")

# Pivot again the dataframe in a wider format
TMTs <- TMT.pivot %>% pivot_wider(id_cols = Accession, names_from = Biological.name , values_from = Abundances) 

```

```{r}

# Modify the order of the columns of TMTs to fit to the order of the design matrix
TMTs <- TMTs %>% dplyr::select(Accession, sort(Labels$Biological.name))

```

# Limma preparation

```{r}

# Keep the Accession order in a variable
Protein.IDs <- TMTs$Accession

# create a group variable to construct the design matrix
# The order of the condition groups should be the same as the order of the TMTs columns
group <-  sort(Labels$Day)
group <- factor(group)

# Construct the design matrix
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)

 # Modify the design into dataframe
design <- cbind(design)

# Let's have a look to the design 
design_to_view = design
rownames(design_to_view) = colnames(TMTs %>% dplyr::select(-Accession))
design_to_view

```

```{r}
# Construct Matrix of custom contrasts
contrast <- makeContrasts(
  D1vsD0 = D1 - D0, 
  D3vsD0 = D3 - D0,
  D3vsD1 = D3 - D1,
  levels = design
)

# Let's have a look
contrast

```




```{r}

# do the linear model fitting
# ---------------------------

fit <- lmFit(TMTs %>% dplyr::select(-Accession), design) # Define the parameters
fit1 <- contrasts.fit(fit, contrast) # To compare the 2 conditions
fit2 <- eBayes(fit1) # To modify the stats

```


# Summary of the Limma significant protein numbers

```{r}

# Let's have a look to the results
datatable(summary(decideTests(fit2)))

```

# D1 vs D0

```{r}
# num = numerator of the ratio
num = "D1"

# den = denominator of the ratio
den = "D0"

# Extract Limma DE results
DE.D1vsD0 <- topTable(fit2, coef = paste0(num, "vs", den), number= nrow(TMTs), sort.by = "none", adjust.method="BH")

```

```{r, fig.width = 7, fig.height=5}

# Add some information
DE.D1vsD0_Results <- get_results(limma_Results = DE.D1vsD0, d = TMTs, Protein.IDs, 
                                  Gene.names = Proteins %>% dplyr::select(Accession, Protein.Group.Accessions,Gene.names, Protein.Desc))


```

# D3 vs D0

```{r}
# num = numerator of the ratio
num = "D3"

# den = denominator of the ratio
den = "D0"

# Extract Limma DE results
DE.D3vsD0 <- topTable(fit2, coef = paste0(num, "vs", den), number= nrow(TMTs), sort.by = "none", adjust.method="BH")

```

```{r, fig.width = 7, fig.height=5}

# Add some information
DE.D3vsD0_Results <- get_results(limma_Results = DE.D3vsD0, d = TMTs, Protein.IDs, 
                                  Gene.names = Proteins %>% dplyr::select(Accession, Protein.Group.Accessions,Gene.names, Protein.Desc))


```

# D3 vs D1

```{r}

# num = numerator of the ratio
num = "D3"

# den = denominator of the ratio
den = "D1"

# Extract Limma DE results
DE.D3vsD1 <- topTable(fit2, coef = paste0(num, "vs", den), number= nrow(TMTs), sort.by = "none", adjust.method="BH")

```

```{r, fig.width = 7, fig.height=5}


# Add some information
DE.D3vsD1_Results <- get_results(limma_Results = DE.D3vsD1, d = TMTs, Protein.IDs, 
                                  Gene.names = Proteins %>% dplyr::select(Accession, Protein.Group.Accessions,Gene.names, Protein.Desc))

```

```{r}
# Write the results into an excel file
# ------------------------------------

# Get general imformation about proteins
All <- Proteins %>% dplyr::select(Accession, Master,Description, Number.of.Peptides, 
                                           Number.of.Unique.Peptides, Number.of.PSMs,
                                            Biological.Process, Cellular.Component, Molecular.Function, Gene.Symbol,
                                            WikiPathways,KEGG.Pathways, Reactome.Pathways)

# Add TMT abundances
All <- All %>% dplyr::left_join(TMTs, by = "Accession")

# Define the different comparisons done with limma
comparison = c("DE.D1vsD0_Results", "DE.D3vsD0_Results", "DE.D3vsD1_Results" )

# Remove the "DE," from the comparison variable to define titles of these comparisons
comparison.title <- gsub("DE.", "", comparison)

# Get the results of the all the different Limma analyses
results_list <- map(1:3, function(x){
  # get results
  df <- get(comparison[x]) %>% dplyr::select(Accession, logFC, p_value, fdr, candidate)
  
  # Add the name of the comparison to the different data (logFC, pvalue, fdr and status)
  Accession <- df$Accession
  df <- df %>% dplyr::select(-Accession)
  colnames(df) <- paste0(comparison.title[x], "_", colnames(df))
  df$Accession <- Accession
  df
})

# Add all the previous information to the global dataframe
All <- All %>% left_join(unique(results_list %>% purrr::reduce(full_join, by = "Accession")), by = "Accession") 

# remove the string "_Results" from different column names
colnames(All) <- sub("_Results", "", colnames(All))

```


```{r}
# Replace NA by "" for a better readability in excel
All <- All %>% replace_na(list(""))

# Write the data in a txt file 
write.table(All, "QuantifiedProteins.txt", sep="\t", row.names = FALSE)

```

# Two-sample Hotelling's T-Square test statistic

<div id="info">
<div id= "info-icon">`r fa("info-circle", fill = "steelblue", height="18px")`</div>
<div id= "info-text">

The multivariate extension of the common two-group Student's t-test
In a t-test, differences in the mean response between two populations are studied.  
T2 is used when the number of response variables are two or more.
The null hypothesis is that the group means for all response variables are equal.

</div>
</div>


```{r}
# Modify the order of the columns of TMTs to fit to the order of the design matrix
TMTs <- TMTs %>% dplyr::select(Accession, "D0_1", "D1_1", "D3_1", "D0_2", "D1_2", "D3_2", "D0_3", "D1_3", "D3_3")

# Put the Accessions into the row names
TMTs <- TMTs %>% column_to_rownames("Accession")

```

```{r}

# Define the different parameters for the mb.long test
replicates <- c("Rep1", "Rep2", "Rep3")
times <- 3
replicate.nb <- 3
assay <- rep(replicates, each=times)
time.grp <- rep(c(1:times), replicate.nb)
size <- rep(replicate.nb, nrow(TMTs))

# mb.long test 
MB.1D <- mb.long(TMTs, method = "1D", times = times, reps=size, rep.grp=assay, time.grp=time.grp, prior.df=NULL)

```

```{r}

print(paste("percent of top ranked proteins after T2 :", MB.1D$percent[2], "%")) # 60 %

# Get the number of top ranked proteins (by the timecourse package directly)
N <- as.integer(MB.1D$percent[2] /100 * nrow(MB.1D$M)) # 2751 = 60 % of the 4585 proteins 
# N

# And get the corresponding T value
Tval <- sort(MB.1D$HotellingT2, decreasing = TRUE)[N] # T > 2.07
# Tval

# Get the top ranked data by the timecourse package    
top.ranked.data <- as.data.frame(MB.1D$M[MB.1D$pos.HotellingT2[1:N],], stringsAsFactors=FALSE)
top.ranked.data <- rownames_to_column(top.ranked.data, var="Accession")

paste("Number of proteins - top ranking after T2 :", nrow(top.ranked.data), "\n")   ##  2751
    
# Add genereal information about the proteins and the time course ranking.
top.ranked.data_4writing <- top.ranked.data %>% mutate(rank = seq(1, nrow(top.ranked.data), 1)) %>%
                          left_join(Proteins %>% dplyr::select(Accession, Description, 
                                                               Biological.Process, Cellular.Component, Molecular.Function , 
                                                                Pfam.IDs, GO.Accessions , Entrez.Gene.ID , Ensembl.Gene.ID , 
                                                                Gene.Symbol , KEGG.Pathways, Reactome.Pathway.Accessions ,
                                                                Reactome.Pathway.Accessions.All , Reactome.Pathways,
                                                               Number.of.Protein.Pathway.Groups), by = "Accession")

# Write the time course results into a txt file
write.table(top.ranked.data_4writing, "TimeCourse_top_ranked_data.txt", row.names = FALSE, sep="\t")

```

## Normalisation according to D0

```{r}
### T2 Hotelling test after a normalization following D0
D0Norm_TMTs <- TMTs %>% mutate(D0vsD0_1_ratio = D0_1 - D0_1, 
                               D1vsD0_1_ratio = D1_1 - D0_1, 
                               D3vsD0_1_ratio = D3_1 - D0_1, 
                               D0vsD0_2_ratio = D0_2 - D0_2, 
                               D1vsD0_2_ratio = D1_2 - D0_2, 
                               D3vsD0_2_ratio = D3_2 - D0_2, 
                               D0vsD0_3_ratio = D0_3 - D0_3, 
                               D1vsD0_3_ratio = D1_3 - D0_3, 
                               D3vsD0_3_ratio = D3_3 - D0_3)

# Keep only the normalized data according to D0 and remove the string "_ratio"
D0Norm_TMTs <- D0Norm_TMTs  %>% dplyr::select(matches("ratio"))
colnames(D0Norm_TMTs) <- sub("_ratio", "", colnames(D0Norm_TMTs))


```

```{r}

# Define the different parameters for the mb.long test
replicates <- c("Rep1", "Rep2", "Rep3")
times <- 3
replicate.nb <- 3
assay <- rep(replicates, each=times)
time.grp <- rep(c(1:times), replicate.nb)
size <- rep(replicate.nb, nrow(TMTs))

# mb.long test 
Norm_MB.1D <- mb.long(D0Norm_TMTs, method = "1D", times = times, reps=size, rep.grp=assay, time.grp=time.grp, prior.df=NULL)

```

```{r}

print(paste("percent of top ranked proteins after T2 :", Norm_MB.1D$percent[2], "%"))

# Get the number of top ranked proteins (by the timecourse package directly)
N <- as.integer(Norm_MB.1D$percent[2] /100 * nrow(Norm_MB.1D$M)) # 2751 = 60 % of the 4585 proteins  
# N

# And get the corresponding T value
Tval <- sort(Norm_MB.1D$HotellingT2, decreasing = TRUE)[N] # T > 2.0777
# Tval

# Get the top ranked data by the timecourse package    
top.ranked.data <- as.data.frame(Norm_MB.1D$M[Norm_MB.1D$pos.HotellingT2[1:N],], stringsAsFactors=FALSE)
top.ranked.data <- rownames_to_column(top.ranked.data, var="Accession")

paste("Number of proteins - top ranking after T2 :", nrow(top.ranked.data), "\n")   ##  2180    

# Add general information about the proteins and the time course ranking.
top.ranked.data_4writing <- top.ranked.data %>% mutate(rank = seq(1, nrow(top.ranked.data), 1)) %>% 
                                left_join(Proteins %>% dplyr::select(Accession, Description,
                                        Biological.Process, Cellular.Component, Molecular.Function, 
                                        Pfam.IDs, GO.Accessions, Entrez.Gene.ID, Ensembl.Gene.ID,
                                        Gene.Symbol, KEGG.Pathways, Reactome.Pathway.Accessions,
                                        Reactome.Pathway.Accessions.All , Reactome.Pathways ,Number.of.Protein.Pathway.Groups), 
                                        by = "Accession")

# Write the time course results into a txt file
write.table(top.ranked.data_4writing, "TimeCourse_top_ranked_data_NormtoD0.txt", row.names = FALSE, sep="\t")
    
```


```{r}
# Get the different information about the R session
sessionInfo()
```
