Skip to content

Commit

Permalink
draft for phospho notebooks
Browse files Browse the repository at this point in the history
  • Loading branch information
@TomSmithCGAT
TomSmithCGAT committed Jul 2, 2024
1 parent 596e62a commit b3e9820
Show file tree
Hide file tree
Showing 8 changed files with 499 additions and 174 deletions.
114 changes: 80 additions & 34 deletions Markdowns/TMT_phospho.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,15 +10,16 @@ output:
bibliography: bib.json
---

```{r, message=FALSE, warning=FALSE}
```{r, message=FALSE, warning=FALSE, include=FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```

### Preamble
### Load dependencies

Load the required libraries.

```{r, message=FALSE}
library(camprotR)
Expand All @@ -29,13 +30,47 @@ library(dplyr)
library(here)
```

### Preamble

The study of phosphorylation by MS-proteomics typically involves enrichment of phosphopeptides using titanium dioxide, since phospho-peptides are relatively low abundance and missed without enrichment.

Since phospho-peptides inform us about the phosphorylation status of a single amino acid, or multiple amino acides within a peptide, we don't want to summarise the phosphoproteomics data to the protein level. Rather, we want to perform a phosphosite-orientated analysis. This is significantly different to standard quantitative proteomics and requires an alternative data processing workflow, which is presented here.

Interpretation of the changes in phosphorylation are aided by quantification of the total protein abundance, since we are usually interested in the change in the proportion of the protein which is phosphorylated. There are many ways to design a phosphoproteomics experiment. The ideal experimental design is for the samples to be labelled with TMT and for the phospho-enrichment to be performed on the pooled TMT-labelled samples, leaving some of the TMT-lablled pool for the total protein abundance quantification. This approach avoids the technical variance that arises with separate phospho-enrichment for each sample and also ensures the total and phospho samples originate from the same material. The limitation of this approach is that the phospho-enrichment input is limited to the amount of material which can be TMT labelled.


### Input data

We start by reading in a file containing PSM-level output from Proteome Discoverer (PD). This data comes from a published benchmark experiment where yeast peptides were spiked into human peptides at 3 known amounts to provide ground truth fold changes (see below). For more details, see [@http://zotero.org/users/5634351/items/LG3W8G4T]
We start by reading in files containing PSM-level output from Proteome Discoverer (PD).

The total data comes from a published benchmark experiment where yeast peptides were spiked into human peptides at 3 known amounts to provide ground truth fold changes (see below). For more details, see [@http://zotero.org/users/5634351/items/LG3W8G4T]

<img src="https://github.com/CambridgeCentreForProteomics/notch/blob/master/images/toc.png?raw=true" width="70%"/>

The data we will use is available through the `Proteomics.analysis.data` package.
Alongside the published total data shown above, we also performed a yeast/human mix at different ratios, from which phosphopeptides were enriched, to simulate an experiment with changes in phoshorylation. The spike in ratios for the total and phosphopeptide TMT plexes are shown below.


|Tag | Total | Phopsho | Ratio |
|------|-------|---------|-------|
| 126 |1x |1x |1 |
| 127N |1x |1x |1 |
| 127C |1x |1x |1 |
| 128N |1x |1x |1 |
| 128C |2x |6x |6 |
| 129N |2x |6x |6 |
| 129C |2x |6x |6 |
| 130N |6x |6x |1 |
| 130C |6x |6x |1 |
| 131 |6x |6x |1 |

From the above, we can see that if we compare tags 126-128N (total 1x; phospho 1x) to tags 128C-129C (total 2x; phospho 6x), we will have a 2-fold increase total protein, but an 6-fold increase in phosphorylation. Whereas, comparing tags 126-128N to 130N-131, we have an identical 6-fold increase in total protein and phosphoprotein.

Note, this experimental design is not quite the same as a phospho-enrichment from the same TMT plex used for the total proteomics, since we needed to spike in at different ratios for the two TMT plexes. Nonetheless, given each sample was simply a defined a spike in between two sample, we expect this data to closely simulate an experimental design with phospho-enrichment from the total TMT plex.


## Processing and QC of Phosphoproteomics data

The data from the phospho-enrichment of the phospho TMT plex is unpublished, but both the total and phospho data is is available through the `Proteomics.analysis.data` package.

```{r}
psm_total_data <- read.delim(
Expand All @@ -47,7 +82,9 @@ psm_phospho_data <- read.delim(
package = "Proteomics.analysis.data"))
```

The first step is to remove contaminant proteins. These were defined using the cRAP database. Below, we parse the cRAP fasta to extract the IDs for the cRAP proteins, in both 'cRAP' format and Uniprot IDs for these proteins.
### Parse PeptideGroups.txt files
To simply the process of reading in the data and performing initial filtering, we will use `camprotR::parse_features`. This function will read in the data and remove contaminant proteins and features without quantification data. Contaminant proteins were defined using the [cRAP database](https://www.thegpm.org/crap/) and provided to PD. We need to obtain their accessions and provide these to `camprotR::parse_features`. Below, we parse the cRAP FASTA to extract the IDs for the cRAP proteins, in both 'cRAP' format and Uniprot IDs for these proteins.


```{r}
crap_fasta_inf <- system.file(
Expand Down Expand Up @@ -116,13 +153,18 @@ We can also add the sequence around the phosphosite, which can be useful to e.g
psm_phospho_data_parsed <- add_site_sequence(psm_phospho_data_parsed, protein_fasta_inf)
```

Finally, we need to add the positions of the peptide for the total peptides, so we can intersect the positions of the phosphosite and total peptides for the statistical testing. We can use `add_peptide_positions` for this, though note that this code chunk takes a long time to run (the function needs to be optimised).
Finally, we need to add the positions of the peptide for the total peptides, so we can intersect the positions of the phosphosite and total peptides for the statistical testing. We can use `add_peptide_positions` for this, though note that this code chunk will take around 1 minute time to run - the function could likely be further optimised.

```{r}
psm_total_data_flt <- add_peptide_positions(psm_total_data_flt, proteome_fasta=protein_fasta_inf)
```

We now store the filtered PSM data for total and phospho in MSnSets, the standard data object for proteomics in R. See the `vignette("msnset", package="camprotR")` for more details.
### Creating an `MSnSet`

We now store the filtered PSM data for total and phospho peptides in MSnSets, the standard data object for proteomics in R. See the `vignette("msnset", package="camprotR")` for more details.For this, we need to generate a matrix of PSM-level abundances, where rows are features (PSMs), and a separate data.frame with information about the features, e.g master protein assignment. See the `vignette("msnset", package="camprotR")` for more details about MSnSets.

We need to create two separate MSnSets for the total and phospho peptides. In both cases, we will manually define the experimental conditions.

```{r}
# Abundance columns for TMT PD-output start with Abundance
Expand All @@ -135,9 +177,7 @@ psm.total.f <- psm_total_data_flt[, setdiff(colnames(psm_total_data_flt), abunda
# update the column names to remove the 'Abundance.` prefix
colnames(psm.total.e) <- gsub('Abundance.', '', colnames(psm.total.e))
# we don't have 'phenotype' data to add so we just define the
# 'expression' data and 'feature' data
# Manually define the 'phenotype' data (experimental details)
psm.total.p <- data.frame(spike=rep(factor(c('x1', 'x2', 'x6')), c(4,3,3)),
condition=rep(1:3, c(4,3,3)),
row.names=colnames(psm.total.e))
Expand All @@ -157,9 +197,7 @@ psm.phospho.f <- psm_phospho_data_parsed[, setdiff(colnames(psm_phospho_data_par
# update the column names to remove the 'Abundance.` prefix
colnames(psm.phospho.e) <- gsub('Abundance.', '', colnames(psm.phospho.e))
# we don't have 'phenotype' data to add so we just define the
# 'expression' data and 'feature' data
# Manually define the 'phenotype' data (experimental details)
psm.phospho.p <- data.frame(spike=rep(factor(c('x1', 'x6')), c(4,6)),
condition=rep(1:3, c(4,3,3)),
row.names=colnames(psm.phospho.e))
Expand All @@ -168,10 +206,14 @@ psm.phospho <- MSnbase::MSnSet(exprs = psm.phospho.e, fData = psm.phospho.f, pDa
```

To simply the downstream plotting and filtering, we create a list containg the two MSnSets.

```{r}
psm_res <- list('Total'=psm.total, 'Phospho'=psm.phospho)
```


Below, we assess the quantification distributions. There are no clear outlier samples or any other concerns.
```{r}
for(x in names(psm_res)){
Expand All @@ -190,46 +232,45 @@ for(x in names(psm_res)){
}
```

Note that where the signal:noise > 10, there are far fewer missing values.
```{r}
for(x in psm_res){
plotNA(x)
}
```
### Removing low quality PSMs
We want to remove low Signal:Noise (S:N) PSMs, since the quantification values will be less accurate and there will be more missing values. We can inspect the relationship between S:N and missing values using the `plot_missing_SN` function.

We use Note that where the signal:noise > 10, there are far fewer missing values.

```{r, fig.height=5, fig.width=7, fig.fullwidth=TRUE, fig.cap="Missing values per PSM, in relation to the signal:noise ratio"}
for(x in names(psm_res)){
p <- psm_res[[x]] %>% plot_missing_SN(bins=50) +
p <- psm_res[[x]] %>% plot_missing_SN(bins=40) +
ggtitle(x)
print(p)
}
```

We can also look into this relationship at the tag level using `plot_missing_SN_per_sample`. In this case, there is no tag which appears to have a high proportion of missing values when signal:noise > 5. If there were, this may warrant further exploration, e.g was one of the sample preparations inadequate such that fewer peptides were labeled?
We can also look into this relationship at the tag level using `plot_missing_SN_per_sample`. In this case, there is no tag which appears to have a high proportion of missing values when signal:noise > 10. If there were, this may warrant further exploration, e.g was one of the sample preparations inadequate such that fewer peptides were labeled?

```{r}
for(x in names(psm_res)){
p <- psm_res[[x]] %>% plot_missing_SN_per_sample(bins=50) +
p <- psm_res[[x]] %>% plot_missing_SN_per_sample(bins=40) +
ggtitle(x)
print(p)
}
```

Based on the above, we will filter the PSMs to only retain those with S:N > 10 using `filter_TMT_PSMs`. Using the same function, we will also remove PSMs with interference/co-isolation >10%.
Based on the above, we will filter the PSMs to only retain those with S:N > 10 using `filter_TMT_PSMs`. Using the same function, we will also remove PSMs with interference/co-isolation >10%, since the mixed species design means our abundances are particularly vulnerable to issues with co-isolation.

```{r}
psm_filt_sn_int <- psm_res %>% lapply(function(x){
filter_TMT_PSMs(x, inter_thresh = 10, sn_thresh = 5)})
```
Since we have merged together material from two species, we will perform a couple more

Finally, we will also filter our PSMs using the Delta.Score, e.g the difference between the score for the top two hits. Where this is

Below, we will also filter our PSMs using the Delta.Score, e.g the difference between the score for the top two hits. Where this is the difference between the top score in the MS2 spectrum matching and the second top score. For a mixed species sample, it's likely some peptides will be difficult to assign to the protein from the correct species, so we are being extra cautious here. For a typical phospho TMT proteomics experiment, it would likely be reasonable to use a more relaxed filter above for the interference (e.g 50%) and to not filter by Delta Score at all.

```{r}
psm_filt_sn_int_delta <- psm_filt_sn_int %>% lapply(function(x){
Expand All @@ -242,8 +283,8 @@ psm_filt_sn_int_delta <- psm_filt_sn_int %>% lapply(function(x){
```


For a typical phospho TMT proteomics experiment, it would likely be reasonable to use a more relaxed filter for the interference (e.g 50%) and to not filter by Delta Score at all.

Finally, we summarise our phosphopeptide quantification values, using the combination of the protein ID and the PTM position(s) as the key. This means that if two peptides cover the same phospho site(s) due to missed cleavage, they will be summarised into a single phospho_site quantification.
```{r}
phospho_sites <- psm_filt_sn_int_delta$Phospho %>%
MSnbase::combineFeatures(
Expand All @@ -256,14 +297,15 @@ print(nrow(phospho_sites))
```

Note that our filtering has removed a lot of PSMs, especially for the phosphosites. We've gone from:
- 5582 PSMs in the input file
- 5335 PSMs after removing contaminants and PSMs without a unique master protein
- 4301 PSMs after removing those without a phosphorylation site
- 2852 PSMs after removing those with a phosphoRS score < 75
- 1524 PSMs after removing those with interference > 10 or Signal:Noise < 10
- 1008 PSMs after removing those with Delta Score < 0.2.
- 666 phosphosites after combining PSMs for the same phosphosite
- `r nrow(psm_phospho_data)` PSMs in the input file
- `r nrow(psm_total_data_flt)` PSMs after removing contaminants and PSMs without a unique master protein
- `r nrow(psm_phospho_data_parsed)` PSMs after removing those without a phosphorylation site, a phosphoRS score < 75, interference > 10 %, Signal:Noise < 10 or Delta Score < 0.2.
- `r nrow(phospho_sites)` after combining PSMs for the same phosphosite

As stated previously, the stringent interference threshold and the use of a delta score threshold are because of the mixed species design. In a more typical experiment, a 50% interference threshold may be suitable and a Delta score threshold would likely not be required.


We also need to summarise our total peptide quantification data. As we shall see in the statistical testing notebook, we shall perform further summarisation for the total peptides, but this is better performed later, for reasons which will become clear.
```{r}
total <- psm_filt_sn_int_delta$Total %>%
MSnbase::combineFeatures(
Expand All @@ -272,11 +314,15 @@ total <- psm_filt_sn_int_delta$Total %>%
```


We save the object to disk so we can read it back into memory when we need it
Finally, we save the object to disk so we can read it back into memory when we need it

```{r, eval=FALSE}
saveRDS(psm_filt_sn_int, here('results/benchmark_tmt_phospho_psm_filt.rds'))
saveRDS(total, here('results/total.rds'))
saveRDS(phospho_sites, here('results/phospho_sites.rds'))
```


For an example how to perform statistical testing for changes in phosphorylation, see
[Intersecting phosphosites and total peptides and statistical testing](https://mrctoxbioinformatics.github.io/Proteomics_data_analysis/Markdowns/TMT_phospho_stats.html)
Loading

0 comments on commit b3e9820

Please sign in to comment.