In MSstatsTMTPTM: Post Translational Modification (PTM) Significance Analysis in shotgun mass spectrometry-based proteomic experiments with tandem mass tag (TMT) labeling

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=8, 
  fig.height=8
)

library(MSstatsTMTPTM)
library(MSstatsTMT)
library(MSstats)
library(dplyr)

This Vignette provides an example workflow for how to use the package MSstatsTMTPTM. It also provides examples and an analysis of how adjusting for global protein levels allows for better interpretations of PTM modeling results.

Installation

To install this package, start R (version "4.0") and enter:

``` {r, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")

BiocManager::install("MSstatsTMTPTM")

## 1. Workflow

### 1.1 Raw Data Format

The first step is to load in the raw dataset for both the PTM and Protein 
datasets. This can be the output of the MSstatsTMT converter functions: 
`PDtoMSstatsTMTFormat`, `SpectroMinetoMSstatsTMTFormat`, and 
`OpenMStoMSstatsTMTFormat`. Both the PTM and protein datasets must include the 
following columns: `ProteinName`, `PeptideSequence`, `Charge`, `PSM`, `Mixture`,
`TechRepMixture`, `Run`, `Channel`, `Condition`, `BioReplicate`, and 
`Intensity`.

#### 1.1.1 Raw PTM Data

```r
# read in raw data files
# raw.ptm <- read.csv(file="raw.ptm.csv", header=TRUE)
head(raw.ptm)

It is important to note the ProteinName column in the PTM dataset represents the modification sites. The location of the modification must be added into the ProteinName. For example the first row shows Protein_12_S703 for ProteinName, with Y474 being the modificaiton site. This can be done as shown above, or by adding the PeptideSequence into the ProteinName, ex. Protein_12_Peptide_491 for the first row.

1.1.1 Raw Protein Data

head(raw.protein)
# raw.protein <- read.csv(file="raw.protein.csv", header=TRUE)

The raw Protein dataset looks similar to the PTM dataset, however the ProteinName column does not contain a modification site.

1.2 proteinSummarization

After loading in the input data, the next step is to use the proteinSummarization function from MSstatsTMT. This provides the summarized dataset needed to model the protein/PTM abundance. The summarization for PTM and Protein datasets should be done separately. The function will summarize the Protein dataset up to the protein level and will summarize the PTM dataset up to the PTM level. The different summarizations are caused by adding the PTM site into the ProteinName field. For details about normalization and imputation options in proteinSummarization please review the package documentation here: MSstatsTMT Package.

``` {r, results='hide', message=FALSE, warning=FALSE}

Use proteinSummarization function from MSstatsTMT to summarize raw data

quant.msstats.ptm <- proteinSummarization(raw.ptm, method = "msstats", global_norm = TRUE, reference_norm = FALSE, MBimpute = TRUE)

quant.msstats.protein <- proteinSummarization(raw.protein, method = "msstats", global_norm = TRUE, reference_norm = FALSE, MBimpute = TRUE)

``` {r}
head(quant.msstats.ptm)
head(quant.msstats.protein)

1.3 groupComparisonTMTPTM

After the two datasets are summarized, both the summarized PTM and protein datasets should be used in the modeling function groupComparisonTMTPTM. First a full pairwise comparison is made between all conditions in the experiment.

``` {r, results='hide', message=FALSE, warning=FALSE}

test for all the possible pairs of conditions

model.results.pairwise <- groupComparisonTMTPTM(data.ptm=quant.msstats.ptm, data.protein=quant.msstats.protein)

Optionally, a specific contrast matrix can be defined. Below is an example of a 
contrast matrix and how it is passed into the `groupComparisonTMTPTM` function.

``` {r,  message=FALSE, warning=FALSE}
# Specify comparisons
comparison<-matrix(c(1,0,0,-1,0,0,
                     0,1,0,0,-1,0,
                     0,0,-1,0,0,-1,
                     1,0,-1,0,0,0,
                     0,1,-1,0,0,0,
                     0,0,0,1,0,-1,
                     0,0,0,0,1,-1,
                     .25,.25,-.5,.25,.25,-.5,
                     1/3,1/3,1/3,-1/3,-1/3,-1/3),nrow=9, ncol=6, byrow=TRUE)

# Set the names of each row
row.names(comparison)<-c('1-4', '2-5', '3-6', '1-3', 
                         '2-3', '4-6', '5-6', 'Partial', 'Third')
# Set the column names
colnames(comparison)<- c('Condition_1', 'Condition_2', 'Condition_3', 
                         'Condition_4', 'Condition_5', 'Condition_6')

comparison

``` {r, results='hide', message=FALSE, warning=FALSE}

test for specified condition comparisons only

model.results.contrast <- groupComparisonTMTPTM(data.ptm=quant.msstats.ptm, data.protein=quant.msstats.protein, contrast.matrix = comparison)

``` {r}
names(model.results.contrast)
ptm_model <- model.results.contrast[[1]]
protein_model <- model.results.contrast[[2]]
adjusted_model <- model.results.contrast[[3]]

head(adjusted_model)

The modeling function will return a list consisting of three dataframes.One each for the PTM-level, Protein-level, and adjusted PTM-level group comparison result.

1.4 Example Volcano Plot

The models from the groupComparisonTMTPTM function can be used in the model visualization function, groupComparisonPlots, from the base MSstats. Below is a Volcano Plot for the Adjusted PTM model. Note: the input for groupComparisonPlots should be one data.frame from output of groupComparisonTMTPTM.

groupComparisonPlots(data = adjusted_model,
                     type = 'VolcanoPlot',
                     ProteinName = FALSE,
                     which.Comparison = '1-4',
                     address = FALSE)

2. Analysis

2.1 How to adjust PTMs for changes in global protein levels?

In order to adjust the PTM abundance for global protein abundance there are two main steps. First the PTM and global protein datasets must be modeled and then the resulting model parameters are combined. The combination is done using the formulas below:

Log2FC: $\log_2FC_{PTM} - \log_2FC_{Protein}$

SE: $\sqrt{SE_{PTM}^2 + SE_{Protein}^2}$

DF: $(SE_{PTM}^2 + SE_{Protein}^2)^2 \biggm/ (\frac{SE_{PTM}^2}{DF_{PTM}} + \frac{SE_{Protein}^2}{DF_{Protein}})$

Please see the package MSstatsPTM for further explanation of adjustment strategy and formulas.

2.2 Example PTM

dataProcessPlotsTMTPTM(data.ptm = raw.ptm,
                       data.protein = raw.protein,
                       data.ptm.summarization = quant.msstats.ptm,
                       data.protein.summarization = quant.msstats.protein,
                       type = 'ProfilePlot',
                       which.Protein = 'Protein_2391_Y40',
                       address = FALSE)

In the plots above the data points for the PTM Protein_2391_Y40 and Protein Protein_2391 are shown for all conditions. To take a look at a specific comparison, the input data can be filtered.

dataProcessPlotsTMTPTM(data.ptm = raw.ptm %>% filter(
  Condition %in% c('Condition_1', 'Condition_4')),
                       data.protein = raw.protein %>% filter(
                         Condition %in% c('Condition_1', 'Condition_4')),
                       data.ptm.summarization = quant.msstats.ptm %>% filter(
                         Condition %in% c('Condition_1', 'Condition_4')),
                       data.protein.summarization = 
    quant.msstats.protein %>% filter(Condition %in% c(
      'Condition_1', 'Condition_4')),
                       type = 'ProfilePlot',
                       which.Protein = 'Protein_2391_Y40',
                       originalPlot = FALSE,
                       address = FALSE)

model_df <- rbind(adjusted_model %>% filter(
  Protein == 'Protein_2391_Y40' & Label == '1-4') %>% select(-Tvalue),
                  ptm_model %>% filter(
                    Protein == 'Protein_2391_Y40' & Label == '1-4'
                    ) %>% select(-issue),
                  protein_model %>% filter(
                    Protein == 'Protein_2391' & Label == '1-4'
                    ) %>% select(-issue))
model_df <- data.frame(model_df)
rownames(model_df) <- c('Adjusted PTM', 'PTM', 'Protein')
model_df

The example above shows the added insight of adjusting for Protein level. Originally the PTM model shows a small positive log2FC of 0.144. However, the global protein abundance changed strongly in the negative direction, -0.508. When the PTM is adjusted for the global the log2FC is increased to 0.652, which is a much larger absolute change than without adjustment. The change in abundance of this PTM may have been marked as insignificant without adjusting for protein levels. Note the change in SE and DF once the PTM model is adjusted.

Additionally, there are other situations where the PTM log2FC will be reduced after adjustment.

dataProcessPlotsTMTPTM(data.ptm = raw.ptm %>% filter(
  Condition %in% c('Condition_2', 'Condition_5')),
                       data.protein = raw.protein %>% filter(
                         Condition %in% c('Condition_2', 'Condition_5')),
                       data.ptm.summarization = quant.msstats.ptm %>% 
    filter(Condition %in% c('Condition_2', 'Condition_5')),
                       data.protein.summarization = quant.msstats.protein %>% 
    filter(Condition %in% c('Condition_2', 'Condition_5')),
                       type = 'ProfilePlot',
                       which.Protein = 'Protein_1076_Y67',
                       originalPlot = FALSE,
                       address = FALSE)

model_df <- rbind(adjusted_model %>% filter(
  Protein == 'Protein_1076_Y67' & Label == '2-5') %>% select(-Tvalue),
                  ptm_model %>% filter(
                    Protein == 'Protein_1076_Y67' & Label == '2-5'
                    ) %>% select(-issue),
                  protein_model %>% filter(
                    Protein == 'Protein_1076' & Label == '2-5'
                    ) %>% select(-issue))

model_df <- data.frame(model_df)
rownames(model_df) <- c('Adjusted PTM', 'PTM', 'Protein')
model_df

Note - In this example Protein_1076 was only available in one run of the global protein experiment.

In this example the PTM model originally showed a large log2FC of 0.242, however most of this log2FC was due to the global protein abundance change which was 0.356. Once adjusted, the PTM shows a much smaller log2FC of -0.114. Without adjustment a potentially incorrect inference could have been drawn for this PTM.

From these two examples the added insight into PTM abundance change from protein adjustment can be observed. If the goal of the experiment is to find the true abundance change due to PTMs than removing the effect of the protein is very important.

Session information

sessionInfo()

MSstatsTMTPTM documentation built on Feb. 18, 2021, 2 a.m.