Aggregate_PSMs_to_Proteins.Rmd

---
title: "Aggregate PSMs to MaxQuant proteinGroups"
author: "Moritz Madern"
date: "2023-02-13"
output: html_document
---

Data input:
1) modified PSM-table (i.e. the output of "IM.Rmd" script termed "PSM.txt")
2) MaxQuant proteins table "proteinGroups.txt" 

Data output:
1) modified MaxQuant proteinGroups.txt


This script performs aggregation of PSMs to higher-level protein features as listed in MaxQuant "proteinGroups.txt" by a) summation of normalized reporter intensities (both interference-corrected and non interference-corrected), b) calculating weighted averages to obtain aggregated EIL values, c) and more (e.g. filtering, aggregating PPF, visualization). 
Note: This script currently supports MaxQuant search engine output only.


```{r Load required packages etc.}

## Load packages
library(tidyverse)


## Create Results folder
if (!file.exists("Results")){
  dir.create("Results")
}


## Extract current working directory
wd <- getwd()

```


```{r Specify required parameters}

## 1) Parameters specific to the PSM table:

## Specify file path to modified PSM table (output of Interference Modeling script "IM.Rmd").
filepath_PSM_table = "Results/modified_PSM.txt"

## Specify non interference-corrected reporter ion column name pattern in PSM table (expected to be isotopic impurity-corrected and normalized).
reporter_pattern_PSM = "reporters_[0-9]*(|N|C)_norm$"

## Specify interference-corrected reporter ion column name pattern in PSM table (expected to be isotopic impurity-corrected and normalized).
reporter_pattern_PSM_interference_corrected = "reporters_[0-9]*(|N|C)_norm__interference_corrected$"

## Specify EIL column name in PSM table.
eil_columnname = "EIL"

## Specify PPF column name in PSM table.
ppf_columnname = "PPF"

## Specify minimum PPF filter threshold. The script filters out PSMs with lower PPF values.
ppf_threshold = 0

## Specify the name of the ID column by which the protein table refers to the PSM table.
id_columnname_PSM_table = "id"

## Specify the name of the column containing the minimum recorded intensity in a PSM's MS2 scan.
min_ms2_int_columnname = "minIntensity_MS2"

## Specify if remaining NAs in PSM's reporter ion intensities should be substituted with the minimum MS2 intensity
substitute_NAs_with_min_intensity = FALSE




## 2) Parameters specific to the site table and subsequent filtering of sites:

## Specify filepath to the protein table (e.g. MaxQuant's "proteinGroups.txt") that lists features of a higher aggregation level than PSMs.
filepath_feature_table =  "proteinGroups.txt"

## Specify the PSM ID column name in the protein table (e.g. "MS MS IDs" column in MaxQuant's "proteinGroups.txt"). This column points to the IDs of every PSM that counts to the aggregated feature (separated by ";"). Based on this information, aggregation is performed.
msms_id_columnname_feature_table = "MS.MS.IDs" 

```


```{r Read in PSM table and derive relevant objects}

## Read in PSM table
df_PSM <- read.delim(filepath_PSM_table, header=TRUE, sep="\t")


## Print the head of PSM reporter intensities (to check if specified pattern is correct!)
df_PSM[,grepl(names(df_PSM), pattern=reporter_pattern_PSM)] %>% head()


## Rename EIL column of PSM table to "EIL"
names(df_PSM)[names(df_PSM) == eil_columnname] <- "EIL"


## Rename PPF column of PSM table to "PPF"
names(df_PSM)[names(df_PSM) == ppf_columnname] <- "PPF"


## Rename ID column name of PSM table to "id"
names(df_PSM)[names(df_PSM) == id_columnname_PSM_table] <- "id"


## Perform additional filtering of PSMs based on PPF
writeLines(paste0("PSMs before filtering based on PPF: ", nrow(df_PSM)))
df_PSM <- df_PSM[df_PSM$PPF >= ppf_threshold,]
writeLines(paste0("\nPSMs after filtering based on PPF: ", nrow(df_PSM)))


## Print the PPF threshold used for filtering
writeLines(paste0("\nThe chosen PPF threshold used for filtering: ", ppf_threshold))

```


```{r Read in protein table and derive relevant objects. Also Filter out reverse and CONs, + based on score. Extract gene names from fasta if empty}

## Read in feature table
df_feature <- read.delim(filepath_feature_table, header=TRUE, sep="\t")


## Filter feature table (assuming it is a MaxQuant output): remove CONs and reverse
writeLines("Filtering rows based on various criteria. Number of rows...")
writeLines(paste0("Before filtering: ", nrow(df_feature)))
df_feature <- df_feature[df_feature$Potential.contaminant == "",]
df_feature <- df_feature[df_feature$Reverse=="",]
writeLines(paste0("\nAfter filtering contaminants and reverse hits: ", nrow(df_feature)))


## Filter for at least 2 razor + unique peptides
bool_keep <- df_feature$Razor...unique.peptides >= 2
df_feature <- df_feature[bool_keep,]
writeLines(paste0("\nAfter filtering for at least 2 razor+unique peptides: ", nrow(df_feature)))


## Filter out only identified by site
bool_keep <- df_feature$Only.identified.by.site != "+"
df_feature <- df_feature[bool_keep,]
writeLines(paste0("\nAfter filtering out proteins only identified by site: ", nrow(df_feature)))


## Check gene name column
writeLines("\nExtracting Gene names from fasta header:")
writeLines(paste0("Number of Gene name entries that are empty: ", sum(df_feature$Gene.names == "")))


## Extract Fasta Headers, and within them, gene names.
gene_names <- character(nrow(df_feature))        #initialize final vector
fasta_headers <- df_feature$Fasta.headers
fasta_headers <- ifelse(substring(fasta_headers, first=1, last=1)==";",
                          no=fasta_headers, 
                          yes=substring(fasta_headers, first=regexpr(fasta_headers,pattern="[A-Za-z0-9]"), last=10000))
fasta_headers_split <- strsplit(fasta_headers, split=";")
for (i in 1:length(fasta_headers_split)){
  # in case protein i had a valid entry already, take the original entry
  if ("Gene.names" %in% names(df_feature) && !df_feature$`Gene.names`[i]== ""){
    gene_names[i] <- df_feature$Gene.names[i]
    next()
  }
  fasta_headers_split_i <- fasta_headers_split[[i]]
  if (length(fasta_headers_split_i) > 0 && grepl(fasta_headers_split_i,pattern="GN=")){
      gene_names_split_i <- substring(fasta_headers_split_i ,first= regexpr(fasta_headers_split_i,pattern="GN=") + 3, last=nchar(fasta_headers_split_i))
  } else{
      gene_names_split_i <- substring(fasta_headers_split_i ,first = 1, last=nchar(fasta_headers_split_i))
  }
  gene_names_split_i
  gene_names_split_i <- substring(gene_names_split_i, first=1, last=regexpr(gene_names_split_i, pattern=" |$")-1)
  gene_names_split_i
  if(length(gene_names_split_i)>1){
    gene_names_i <- paste0(gene_names_split_i, collapse=";")
  } else{
    gene_names_i <- gene_names_split_i
  }
  gene_names[i] <- ifelse(length(gene_names_i)>0, yes=gene_names_i, no="")
}

## Replace old vector with new gene name vector
df_feature$Gene.names <- gene_names
writeLines(paste0("Number of Gene name entries that are empty after matching with fasta header: ", sum(df_feature$Gene.names == "")))

```


```{r Perform aggregation of non-interference-corrected reporter ion intensities and EIL + PPF values from PSM-table to protein table based on PSM ID reference column}

## Extract msms ID column in feature table
msms_id_pointer <- df_feature[,msms_id_columnname_feature_table]


## Initiate aggregated reporter intensity matrix
m_agg_reporterInt <- matrix(numeric(1),
                            ncol=sum(grepl(names(df_PSM), pattern=reporter_pattern_PSM)),  # number of reporter ion columns in PSM-table
                            nrow=nrow(df_feature))
colnames(m_agg_reporterInt) <- paste0(grep(names(df_PSM), pattern=reporter_pattern_PSM, value = TRUE), "__aggregated")
 

## Initiate aggregated EIL vector and PPF vector
v_agg_eil <- numeric(nrow(df_feature))
v_agg_ppf <- numeric(nrow(df_feature))


## Initiate vector indicating whether PSMs for the respective site could be found in the PSM table
v_agg_control <- rep("+", times=nrow(df_feature))


## Go over each site i in df_feature and aggregate their reporter intensities + EIL values  + PPF values from the PSM-table
for (i in 1:nrow(df_feature)){
  
  # extract corresponding PSM IDs
  psm_ids_i <- msms_id_pointer[i]
  
  # split them at ";"
  psm_ids_i <- strsplit(psm_ids_i, split=";")[[1]]
  
  # check if the feature really points to PSMs in the PSM-table
  df_PSM_i <- df_PSM[df_PSM$id %in% psm_ids_i,]
  if (nrow(df_PSM_i)==0){ 
    v_agg_control[i] <- "-"
    next()
  }
  
  # extract reporter intensity columns
  m_PSM_i <- df_PSM_i[,grepl(names(df_PSM_i), pattern=reporter_pattern_PSM)] %>% as.matrix(.)
  
  # substitute 0s with NAs
  m_PSM_i[m_PSM_i==0] <- NA
  
  # calculate rowSums to infer weights in EIL aggregation + PPF aggregation
  rowsums_i <- rowSums(m_PSM_i, na.rm=TRUE)
  weights_i <- rowsums_i/sum(rowsums_i)
  
  # aggregate EIL and PPF
  v_agg_eil[i] <- sum(weights_i * df_PSM_i$EIL)
  v_agg_ppf[i] <- sum(weights_i * df_PSM_i$PPF)
  
  # substitute NAs and zeros with minimum intensity in MS2 scan if specified
  if (substitute_NAs_with_min_intensity){
    for (j in 1:nrow(m_PSM_i)){
      row_j <- m_PSM_i[j,]
      min_j <- df_PSM_i[j,min_ms2_int_columnname]
      row_j[is.na(row_j) | row_j == 0] <- min_j
      m_PSM_i[j,] <- row_j
    }
  }
  
  # calculate column sums (i.e. aggregation from PSMs to feature by summation)
  m_agg_reporterInt[i,] <- colSums(m_PSM_i, na.rm = TRUE)
}


## Plot total intensities per channel
barplot(colSums(m_agg_reporterInt), las=2, cex.names=0.8, col="grey", border="grey", main="Total intensities after aggregation")


## Merge the results with df_PSM
df_feature <- cbind(df_feature, m_agg_reporterInt)
df_feature$EIL <- v_agg_eil
df_feature$PPF <- v_agg_ppf
df_feature$agg_control <- v_agg_control

```


```{r Perform aggregation of interference-corrected reporter ion intensities and EIL values from PSM table to protein table based on PSM ID reference column}

## Extract MSMS ID column in feature table
msms_id_pointer <- df_feature[,msms_id_columnname_feature_table]


## Initiate aggregated reporter intensity matrix
m_agg_reporterInt <- matrix(numeric(1),
                            ncol=sum(grepl(names(df_PSM), pattern=reporter_pattern_PSM_interference_corrected)),  # number of reporter ion columns in PSM-table
                            nrow=nrow(df_feature))
colnames(m_agg_reporterInt) <- paste0(grep(names(df_PSM), pattern=reporter_pattern_PSM_interference_corrected, value = TRUE), "__aggregated")
 

## Initiate aggregated EIL vector and PPF vector
v_agg_eil <- numeric(nrow(df_feature))
v_agg_ppf <- numeric(nrow(df_feature))


## Initiate vector indicating whether PSMs for the respective site could be found in the PSM table
v_agg_control <- rep("+", times=nrow(df_feature))


## Go over each site i in df_feature and aggregate their reporter intensities + EIL values  + PPF values from the PSM-table
for (i in 1:nrow(df_feature)){
  
  # extract corresponding PSM IDs
  psm_ids_i <- msms_id_pointer[i]
  
  # split them at ";"
  psm_ids_i <- strsplit(psm_ids_i, split=";")[[1]]
  
  # check if the feature really points to PSMs in the PSM-table
  df_PSM_i <- df_PSM[df_PSM$id %in% psm_ids_i,]
  if (nrow(df_PSM_i)==0){ 
    v_agg_control[i] <- "-"
    next()
  }
  
  # extract reporter intensity columns
  m_PSM_i <- df_PSM_i[,grepl(names(df_PSM_i), pattern=reporter_pattern_PSM_interference_corrected)] %>% as.matrix(.)
  
  # substitute 0s with NAs
  m_PSM_i[m_PSM_i==0] <- NA
  
  # calculate rowSums to infer weights in EIL aggregation + PPF aggregation
  rowsums_i <- rowSums(m_PSM_i, na.rm=TRUE)
  weights_i <- rowsums_i/sum(rowsums_i)
  
  # aggregate EIL and PPF
  v_agg_eil[i] <- sum(weights_i * df_PSM_i$EIL)
  v_agg_ppf[i] <- sum(weights_i * df_PSM_i$PPF)
  
  # substitute NAs and zeros with minimum intensity in MS2 scan if specified
  if (substitute_NAs_with_min_intensity){
    for (j in 1:nrow(m_PSM_i)){
      row_j <- m_PSM_i[j,]
      min_j <- df_PSM_i[j,min_ms2_int_columnname]
      row_j[is.na(row_j) | row_j == 0] <- min_j
      m_PSM_i[j,] <- row_j
    }
  }
  
  # calculate column sums (i.e. aggregation from PSMs to feature by summation)
  m_agg_reporterInt[i,] <- colSums(m_PSM_i, na.rm = TRUE)
}


## Plot total intensities per channel
barplot(colSums(m_agg_reporterInt), las=2, cex.names=0.8, las=2, col="grey", border="grey", main="Total intensities after aggregation")

## Merge the results with df_PSM
df_feature <- cbind(df_feature, m_agg_reporterInt)
df_feature$EIL <- v_agg_eil
df_feature$PPF <- v_agg_ppf
df_feature$agg_control <- v_agg_control

```


```{r Filter out proteins arising from previously filtered PSMs}

## Note: This step filters out proteins that wholly derive from PSMs which have previously been filtered out in the PSM table (due to, for example: 0 isobaric labels (i.e. unlabeled); second peptides; insufficient PPF, insufficient valid values per Group)


## Get an overview of how much is filtered out
barplot(table(df_feature$agg_control), border="grey", col="grey", ylab = "frequency", xlab="passed filter", main="Overview of filtering")


## Perform filtering
writeLines("Number of proteins...")
writeLines(paste0("Before integrating PSM-filters: ", nrow(df_feature)))
df_feature <- df_feature[df_feature$agg_control == "+",]
writeLines(paste0("After integrating PSM-filters: ", nrow(df_feature)))

```


```{r Export results}

## If no specifc output name was given, name resulting table based on input file's name
ind_substring <- rev(unlist(gregexpr(filepath_feature_table, pattern="[\\]|[/]")))[1]
outputname_path <- paste0("Results/aggregated_",substring(filepath_feature_table, first = ind_substring + 1, last=nchar(filepath_feature_table))) 


## Generate output table
write.table(df_feature, file = outputname_path, sep = "\t", col.names = TRUE, row.names=FALSE, quote=FALSE)


## Report about the export
writeLines("The result table can be found in:")
print(outputname_path)

```