Aggregate_PSMs_to_Sites.Rmd

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

Data input:
1) modified PSM-table (output of "IM.Rmd" script)
2) MaxQuant site table (e.g. "Acetyl (K)sites.txt")

Data output:
1) modified site table (aggregated from PSMs)


This script performs aggregation of PSMs to higher-level site features as listed in MaxQuant site tables 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., message=FALSE, warning=FALSE, echo=FALSE}

## Load packages
library(tidyverse)
library(rlist)


## 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 PSM table containing: A) reporter ion intensity columns (impurity-corrected and normalized); B) columns providing feature-wise EIL-values (estimated interference level) and PPF-values (precursor purity fraction); C) a unique ID by which other tables (i.e. Site table) refer to it.
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 site table refers to the PSM table.
id_columnname_PSM_table = "id"

## Specify the column denoting the integer number of modifications of interest for the respective precursor peptide of the PSM.
number_modifications_columnname = "Acetyl..K." 

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

## Specify if remaining NAs in PSM's reporter ion intensities should be imputed with the minimum MS2 intensity. Note that this can overestimate the real intensity! The recommended setting is FALSE (even though it leaves NAs in the data).
substitute_NAs_with_min_intensity = FALSE

## Specify whether to split site entries (i.e. rows in the site table) by number of modifications. If TRUE, this will create additional rows per unique site id if PSMs for that site list different unique numbers of modifications (e.g. 1-time acetylated, 2-time acetylated, etc.). The recommended setting is TRUE.
split_sites_by_number_of_mods = TRUE




## 2) Parameters specific to the site table:

## Specify filepath tosite table (e.g. MaxQuant's "Acetyl (K)Sites.txt" or Phospho (STY)Sites.txt") that lists features of higher aggregation level than PSMs. 
filepath_feature_table =  "./Acetyl (K)Sites.txt"   

## Specify the PSM ID column name in the site table (e.g. MS MS IDs" column in MaxQuant's site table). 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" 

## Specify minimum site score threshold. Sites with lower score will be filtered out.
score_threshold_feature_table = 0

## Specify minimum top 3 average intensity filter threshold.
top3_threshold_feature_table = 0


```


```{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 minimum MS2 inteisity of PSM table to "minIntensity_MS2"
names(df_PSM)[names(df_PSM) == minIntensity_MS2_columnname] <- "minIntensity_MS2"


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


## Rename number of modifications colname
names(df_PSM)[names(df_PSM) == number_modifications_columnname] <- "number_modifications"


## Perform filtering of PSMs that do not contain the specified modification
writeLines(paste0("PSMs before filtering out PSMs without modification of interest: ", nrow(df_PSM)))
df_PSM <- df_PSM[df_PSM$number_modifications >0,]
writeLines(paste0("PSMs after filtering out PSMs without modification of interest: ", nrow(df_PSM)))


## Perform additional filtering of PSMs based on PPF
writeLines(paste0("\nPSMs before filtering based on PPF: ", nrow(df_PSM)))
df_PSM <- df_PSM[df_PSM$PPF >= ppf_threshold,]
writeLines(paste0("PSMs 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 site table and derive relevant objects. Also Filter out reverse and CONs, + based on score. Extract gene names from fasta if empty etc.}

## 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("After filtering contaminants and reverse hits: ", nrow(df_feature)))


## Filter feature table based on score
df_feature <- df_feature[df_feature$Score >= score_threshold_feature_table, ]
writeLines(paste0("After further filtering based on minimum score threshold: ", 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 == "")))

```


Note that MaxQuant's site table lists the intensities of sites (i.e. rows) in three distinct columns per channel (ending with "___1", "___2" and "___3", respecively) which correspond to different numbers of modifications on measured peptides. This script will also treat different numbers of modifications per site separately, but instead will add extra rows in the feature table for each unique recorded number of modification per site (see next code block). Consequently, some sites (i.e. rows) will be duplicated to create multiple rows, each resulting row listing its own set of PSM IDs that match the site but which all possess distinct numbers of modifications. 


```{r Prepare aggregation and update PSM-ID column (MS.MS.IDs) in case of splitting sites by number of modifications}

## Prepare output list
list_extendedFeatureTable <- list()


## Go over each site (row) and add extra rows in case of multiple modifications (info comes from PSM table). One extra row for each unique number of modifications per site.
for (i in 1:nrow(df_feature)){
  
  # extract feature i in site table
  df_feature_i <- df_feature[i,]
  
  # extract corresponding PSM ids of site i
  psm_ids_i <- df_feature_i[,msms_id_columnname_feature_table]
  
  # split them at ";"
  psm_ids_i <- strsplit(psm_ids_i, split=";")[[1]]
  
  # check if the feature really points to any PSMs in the PSM-table (might have been filtered out already). If not, skip to next site
  df_PSM_i <- df_PSM[df_PSM$id %in% psm_ids_i,]
  if (nrow(df_PSM_i)==0){ 
    next()
  }
  
  # infer the number of unique modification counts across all PSMs pertaining to that site
  unique_mod_number_i <- df_PSM_i$number_modifications %>% unique() %>% sort()
  
  # check whether splitting by number of modifications is desired (via parameter split_sites_by_number_of_mods). 
  if (!split_sites_by_number_of_mods){
    df_feature_i_new <- df_feature_i
    df_feature_i_new$number_modifications <- paste0(unique_mod_number_i, collapse=";")
    list_extendedFeatureTable <- list.append(list_extendedFeatureTable, df_feature_i_new)
    next()
  }
  
  # go over the unique modification numbers and create new site table rows for each one. Save in list_extendedFeatureTable
  for (j in unique_mod_number_i){
    
    # extract PSM ids of PSMs with number of modifications j of site i
    msms_ids_ij <- df_PSM_i[df_PSM_i$number_modifications  == j, id_columnname_PSM_table] %>% paste0(., collapse=";")
    
    # generate new feature table entry with modified values
    df_feature_i_new <- df_feature_i
    df_feature_i_new[,msms_id_columnname_feature_table] <- msms_ids_ij
    df_feature_i_new$number_modifications <- j
    
    # store new feature table entry
    list_extendedFeatureTable <- list.append(list_extendedFeatureTable, df_feature_i_new)
  }
}


## Unlist via rbind and store as new df_feature
df_feature <- do.call(what=rbind, args=list_extendedFeatureTable)


```


```{r Perform aggregation of non-interference-corrected reporter ion intensities and EIL + PPF values from PSM table to site 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))


## 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){ 
    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

```


```{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))


## 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){ 
    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_feature
df_feature <- cbind(df_feature, m_agg_reporterInt)
df_feature$EIL <- v_agg_eil
df_feature$PPF <- v_agg_ppf

```

 
```{r Filter out sites that do not pass a minimum average top3 intensity threshold}

## Extract aggregated intensities
m_aggregated <- df_feature[,grepl(names(df_feature), pattern="norm__aggregated")] %>% as.matrix()


## Calculate average top 3 intensity per feature
mean_top3 <- apply(m_aggregated, MARGIN = 1, FUN=function(x){
                                  x_sort <- sort(x, decreasing = TRUE)
                                  if(length(x_sort)>=3){
                                    res <- mean(x_sort[1:3])
                                  } else{
                                    res <- mean(x_sort[1:length(x_sort)])
                                  }
                                  return(res)
                                  })


## Filter out sites below average top3 threshold
hist(log2(mean_top3), breaks=100, border="grey", col="grey", main="distribution of average top3 log2 intensities")
abline(v=log2(top3_threshold_feature_table), col="red", lty="dashed", lwd=2)
bool_keep <- mean_top3 >= top3_threshold_feature_table
writeLines(paste0("Number of sites after applying average top3 filter threshold: ", nrow(df_feature)))
df_feature <- df_feature[bool_keep,]
writeLines(paste0("\nNumber of sites after applying average top3 filter threshold: ", nrow(df_feature)))
writeLines(paste0("Chosen average top3 intensity threshold for filtering: ", top3_threshold_feature_table))

```


```{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)

```