The goal of mspms is provide a R Package that can be used to easily and reproducibly analyze data resulting from the Multiplex Substrate Profiling by Mass Spectrometry (MSP-MS) method.
Additionally, we provide a graphical user interface powered by shiny apps that allows for a user to utilize the method without requiring any R coding knowledge.
Multiplex Substrate Profiling by Mass Spectrometry (MSP-MS) is a powerful method used to determine the substrate specificity of proteases. This method is of interest for groups interested in the study of proteases and their role as regulators of biological pathways, whether applied to the study of disease states, the development of diagnostic and prognostic tests, generation of tool compounds, or rational design of protease targeting therapeutics. Analysis of the MS based data produced by MSP-MS is a multi-step process involving detection and quantification of peptides, imputation, normalization, cleavage sequence identification, statistics, and data visualizations. This process can be challenging, especially for scientists with limited programming experience. To overcome these issues, we developed the mspms R package to facilitate the analysis of MSP-MS data utilizing good software/data analysis practices.
mspms differs from existing proteomics packages that are available in the Bioconductor project in that it is specifically designed to analyze MSP-MS data. This involves unique preprocessing steps and data visualizations tailored to this specific assay.
In order to do so, mspms uses several excellent packages from the Bioconductor project to provide a framework to easily and reproducibly analyze MSP-MS data. These include:
-
QFeatures: management and processing of peptide quantitative features and sample metadata.
-
MsCoreUtils: log2 transformation, imputation (QRILC), and normalization (center.median).
Some excellent non-Bioconductor packages are also utilized.
-
ggplot2: data visualizations
-
heatmaply: interactive heatmaps
-
rstatix: conduct statistics
-
ggseqlogo: iceLogo motif visualization.
-
downloadthis: gives users the ability to download data from standard .html report generated by mspms::generate_report().
Lastly, mspms takes advantage of iceLogos (Colaert, N. et al. Nature Methods 6, 786-787 (2009)) to viusalize overrepresented amino acid motifs relative to a background set by implementing components of the Java software in R.
You can install the development version from github.
devtools::install_github("baynec2/mspms")To generate a general report using your own data, run the following code. It requires data that has been prepared for mspms data analysis by a converter function. For more information see subsequent sections.
mspms::generate_report(
prepared_data = mspms::peaks_prepared_data,
outdir = "../Desktop/mspms_report"
)The above command will generate a .html file containing a generic mspms analysis.
There is much more that can be done using the mspms package- see the following sections for more information.
Functions in this package are logically divided into 4 broad types of functions.
-
Pre-processing data. These functions are focused on getting all of the data necessary for a mspms analysis together and in a consistent format.
-
Data processing. These functions allow the user to process the proteomics data (imputation, normalization, etc).
-
Statistics. These methods allow the user to perform basic statistics on the normalized/processed data.
-
Data visualization. These functions allow the user to visualize the data.
Note that only a subset of functions are exported to the user. Helper functions can be found in the preprocessing_helper_functions.R, processing_helper_functions.R,icelogo_helper_functions.R, statistics_helper_functions.R, and iceLogo_helper_functions.R files.
These functions were not intended to be directly used by the user and are not exported in an effort to make the package API cleaner and more intuitive.
Pre-processing data here refers to standardizing the data exported from different software solutions for analyzing mass spectrometry-based proteomics and adding useful information specific to the MSP-MS assay.
Specifically, the data from the file exported from the upstream software is parsed to contain the peptide (detected peptide sequence, cleavages are represented as “_”), library_id (name of the peptide in the peptide library each detected peptide maps to), and intensities for each sample.
Information about each peptide relative to the peptide it was derived from is then added. This includes the cleavage_seq (the motif the peptide a user specified n residues up and down from the cleavage event), cleavage_pos (position of the peptide library the cleavage product was cleaved at ), peptide_type (whether the detected peptide is a cleavage product of the peptide library, or a full length member of the peptide library).
All of the information is then loaded as a QFeatures object along with colData (metadata about each sample) where it is stored as a SummarizedExperiment object named “peptides”.
As briefly described above, the colData contains metadata describing the samples you are trying to analyze.
It is very important that this contains at least the following columns: quantCols”,“group”,“condition”,and “time”.
The quantCols field contains the name of each sample. These names must match whatever the samples are named in the files exported from the upstream proteomics software.
The group field contains text specifying what group each sample belongs to. All replicate samples should have the same group name
The condition field contains the experimental conditions that vary in your MSP-MS experiment (this depends on the experimental question, could be type of protease, the type of inhibitor, etc)
The time field contains the duration of incubation of your protease of interest with the peptide library for each sample.
An example of a valid colData file is shown below.
colData <- readr::read_csv(system.file("extdata/colData.csv",
package = "mspms"
))
#> Rows: 12 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): quantCols, group, condition
#> dbl (1): time
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(colData)
#> # A tibble: 6 × 4
#> quantCols group condition time
#> <chr> <chr> <chr> <dbl>
#> 1 CatA_0000_1 CatA_0000 CatA 0
#> 2 CatA_0000_2 CatA_0000 CatA 0
#> 3 CatA_0000_3 CatA_0000 CatA 0
#> 4 CatA_0000_4 CatA_0000 CatA 0
#> 5 CatA_0030_1 CatA_0030 CatA 30
#> 6 CatA_0030_2 CatA_0030 CatA 30Analysis using a .csv exported from PEAKS software is supported.
an exported .csv file can be generated from PEAKS as follows:
- Create a New Project in PEAKS.
- Select data to add to PEAKS project.
- Add data, rename samples, set correct parameters. • Highlight all
the added data.
• Select Create new sample for each file.
• Select appropriate instrument and fragmentation method.
• Set Enzyme to “None”.
• Rename samples according to enzymes, time points, and replicates.
• Can also ensure appropriate data is selected for these.
• Select “Data Refinement” to continue. - Data Refinement.
• On the top right, select MSP-MS as the predefined parameters.
• Parameters shown are what should be used for MSP-MS data analysis.
• Select “Identification” to continue. - Peptide Identification.
• On the top right, select MSP-MS as the predefined parameters.
• Parameters shown are what should be used for MSP-MS data analysis.
• Identification -> Select Database -> TPD_237.
• Identification -> no PTMs.
• Can remove the PTMs by highlighting them and selecting Remove - Label Free Quantification.
• Make sure Label Free is selected.
• Group samples -> add quadruplicates of samples to new group.
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(mspms)
#>
#> Attaching package: 'mspms'
#> The following object is masked _by_ '.GlobalEnv':
#>
#> colData
# File path of peaks lfq file
lfq_filepath <- system.file("extdata/peaks_protein-peptides-lfq.csv",
package = "mspms"
)
colData_filepath <- system.file("extdata/colData.csv", package = "mspms")
# Prepare the data
peaks_prepared_data <- mspms::prepare_peaks(lfq_filepath,
colData_filepath,
quality_threshold = 0.3,
n_residues = 4
)
#> Rows: 2444 Columns: 32
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (6): Protein Accession, Peptide, Used, Candidate, Sample Profile (Ratio...
#> dbl (25): Protein Group, Protein ID, Quality, Significance, Avg. ppm, Avg. A...
#> lgl (1): PTM
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> 841 peptides were removed because they had a quality score < 0.3 (34%)
#>
#> Rows: 12 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): quantCols, group, condition
#> dbl (1): time
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Checking arguments.
#>
#> Loading data as a 'SummarizedExperiment' object.
#>
#> Formatting sample annotations (colData).
#>
#> Formatting data as a 'QFeatures' object.
peaks_prepared_data
#> An instance of class QFeatures containing 1 assays:
#> [1] peptides: SummarizedExperiment with 2071 rows and 12 columnsData exported from Fragpipe as a .tsv file is supported.
To generate this file, run a LFQ Fragpipe workflow using your desired parameters (we recommend using match between runs).
You will find a file named “combined_peptide.tsv” in the output folder from Fragpipe.
This can be loaded into mspms as follows:
combined_peptide_filepath <- system.file("extdata/fragpipe_combined_peptide.tsv",
package = "mspms"
)
colData_filepath <- system.file("extdata/colData.csv", package = "mspms")
fragpipe_prepared_data <- prepare_fragpipe(
combined_peptide_filepath,
colData_filepath
)
#> Rows: 1847 Columns: 63
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr (19): Peptide Sequence, Prev AA, Next AA, Protein, Protein ID, Entry Nam...
#> dbl (40): Start, End, Peptide Length, CatA_0000_1 Spectral Count, CatA_0000_...
#> num (1): Charges
#> lgl (3): Gene, Mapped Genes, Mapped Proteins
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Rows: 12 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): quantCols, group, condition
#> dbl (1): time
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Checking arguments.
#>
#> Loading data as a 'SummarizedExperiment' object.
#>
#> Formatting sample annotations (colData).
#>
#> Formatting data as a 'QFeatures' object.
fragpipe_prepared_data
#> An instance of class QFeatures containing 1 assays:
#> [1] peptides: SummarizedExperiment with 1694 rows and 12 columnsAnalysis of a PeptideGroups.txt file exported from proteome discoverer is supported.
To generate this .txt file: 1. Ensure that non-normalized/ non-scaled
abundances are included in your PeptideGroup table.
2. File > export > To Text (tab delimited) > Items to be exported -
check the Peptide Groups box > Export.
Note that the fileIDs of each sample in proteome discoverer (By default F1 - Fn) must match the quantCols in colData.
peptide_groups_filepath <- system.file(
"extdata/proteome_discoverer_PeptideGroups.txt",
package = "mspms"
)
colData_filepath <- system.file("extdata/proteome_discover_colData.csv",
package = "mspms"
)
prepared_pd_data <- prepare_pd(peptide_groups_filepath, colData_filepath)
#> Rows: 1207 Columns: 68
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr (17): Confidence, Annotated Sequence, Master Protein Accessions, Positio...
#> dbl (48): Peptide Groups Peptide Group ID, Qvality PEP, Qvality q-value, # P...
#> lgl (3): Checked, Modifications, Modifications in Master Proteins
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Rows: 12 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (3): quantCols, group, condition
#> dbl (1): time
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Checking arguments.
#>
#> Loading data as a 'SummarizedExperiment' object.
#>
#> Formatting sample annotations (colData).
#>
#> Formatting data as a 'QFeatures' object.
prepared_pd_data
#> An instance of class QFeatures containing 1 assays:
#> [1] peptides: SummarizedExperiment with 1161 rows and 12 columnsData processing includes imputation and normalization of the proteomics data. This is all done in a QFeatures object using MSCoreUtils methods. Data from each step is stored as a SummarizedExperiment within the QFeatures object with the indicated name.
- Data is first log2 transformed (“peptides_log”).
- Data is normalized using the center.median method (“peptides_log_norm”).
- Missing values are then imputed using the QRILC method (“peptides_log_impute_norm”).
- Data is then reverse log2 transformed (“peptides_norm”).
set.seed(2)
processed_qf <- process_qf(peaks_prepared_data)
#> Loading required namespace: imputeLCMD
#> Imputing along margin 2 (samples/columns).
#> Your row data contain missing values. Please read the relevant
#> section(s) in the aggregateFeatures manual page regarding the effects
#> of missing values on data aggregation.
processed_qf
#> An instance of class QFeatures containing 5 assays:
#> [1] peptides: SummarizedExperiment with 2071 rows and 12 columns
#> [2] peptides_log: SummarizedExperiment with 2071 rows and 12 columns
#> [3] peptides_log_norm: SummarizedExperiment with 2071 rows and 12 columns
#> [4] peptides_log_impute_norm: SummarizedExperiment with 2071 rows and 12 columns
#> [5] peptides_norm: SummarizedExperiment with 2071 rows and 12 columnsMSP-MS analysis has traditionally used T-tests to what peptides are significantly different relative to T0. The authors of the assay suggest using a T-test derived FDR adjusted p value threshold <= 0.05 and log2 fold change threshold >= 3.
We provide a user-friendly implementation of this statistical test as follows:
log2fc_t_test_data <- mspms::log2fc_t_test(processed_qf = processed_qf)
#> Joining with `by = join_by(quantCols)`
#> Joining with `by = join_by(quantCols)`Please note that statistical analyses are one of the strengths of the Bioconductor project and several packages provide alternative statistical approaches that are more flexible than t-tests and could be used with the QFeatures objects produced by mspms.
This includes packages using empirical Bayesian approaches (see limma) and mixed modelling (see msqrob2), that are suited for proteomics data.
Users are encouraged to use these excellent packages.
We provide a number of data visualizations as part of the mspms package.
A critical initial step in MSP-MS data analysis is a rigorous quality control assessment. Given the reliance on a specific peptide library, calculating the percentage of library peptides detected in each sample provides a straightforward yet powerful indicator of data integrity. Anomalously low detection rates signal potential instrument or experimental problems.
We can see the percentage of the peptide library that were undetected among each group of replicates as follows (for full length peptides a threshold of ~< 10 is considered good while for cleavage products a threshold ~< 5 is good):
plot_qc_check(processed_qf,
full_length_threshold = 10,
cleavage_product_threshold = 5
)
#> Joining with `by = join_by(quantCols)`
We can also see in what percentage of samples each peptide in the peptide library were undetected in.
plot_nd_peptides(processed_qf)
#> Joining with `by = join_by(quantCols)`
We can then convert our QFeatures object into a “long” tibble using the mspms_tidy function.
This makes the data easy to manipulate prior to ggplot2 or plotly based plotting functions.
mspms_tidy_data <- mspms_tidy(processed_qf)
#> Joining with `by = join_by(quantCols)`We can inspect our data using an interactive heatmap as follows:
plot_heatmap(mspms_tidy_data)note that the interactive plot is not returned in this vignette due to size constraints.
We can also inspect our data using a PCA as follows.
plot_pca(mspms_tidy_data, value_colname = "peptides_norm")
We can generate volcano plots showing the significant peptides from each condition as follows:
plot_volcano(log2fc_t_test_data)
We can also create plots showing the number of times we observe a cleavage at each position of the peptide_library among a set of peptides of interest.
We can do this as follows:
sig_cleavage_data <- log2fc_t_test_data %>%
dplyr::filter(p.adj <= 0.05, log2fc > 3)
p1 <- mspms::plot_cleavages_per_pos(sig_cleavage_data)
p1
Icelogos are a method for visualizing conserved patterns in protein sequences through probability theory.
As part of the mspms R package, we have implemented the logic to generate iceLogos in R.
To use this, we first have to define a vector of cleavage sequences that we are interested in. We will do this for the significant CatA cleavages below.
catA_sig_cleavages <- log2fc_t_test_data %>%
dplyr::filter(p.adj <= 0.05, log2fc > 3) %>%
dplyr::filter(condition == "CatA") %>%
dplyr::pull(cleavage_seq) %>%
unique()We also need to know all of the possible cleavages we could see among our peptide library. We can generate that as below:
all_possible_8mers_from_228_library <- calculate_all_cleavages(
mspms::peptide_library$library_real_sequence,
n_AA_after_cleavage = 4
)We can then generate an iceLogo relative to the background of all possible cleavages in our peptide library like so:
plot_icelogo(catA_sig_cleavages,
background_universe = all_possible_8mers_from_228_library
)
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
We also provide a function that conveniently plots all icelogos from all conditions from the experiment all together.
sig_cleavage_data <- log2fc_t_test_data %>%
dplyr::filter(p.adj <= 0.05, log2fc > 3)
plot_all_icelogos(sig_cleavage_data)
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
We can also plot the intensities of peptides of interest over time.
First we need to define which peptides we are interested in. Here we will use the t-test results, and look at the peptide with the max log2fc.
max_log2fc_pep <- log2fc_t_test_data %>%
dplyr::filter(p.adj <= 0.05, log2fc > 3) %>%
dplyr::filter(log2fc == max(log2fc)) %>%
pull(peptide)We can then filter our tidy mspms data to only include this peptide and plot
p1 <- mspms_tidy_data %>%
dplyr::filter(peptide == max_log2fc_pep) %>%
plot_time_course()
p1
sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 15.1
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/Los_Angeles
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mspms_0.99.5 dplyr_1.1.4
#>
#> loaded via a namespace (and not attached):
#> [1] gridExtra_2.3 sandwich_3.1-1
#> [3] rlang_1.1.4 magrittr_2.0.3
#> [5] clue_0.3-66 matrixStats_1.4.1
#> [7] compiler_4.4.1 vctrs_0.6.5
#> [9] reshape2_1.4.4 stringr_1.5.1
#> [11] ProtGenerics_1.36.0 pkgconfig_2.0.3
#> [13] crayon_1.5.3 fastmap_1.2.0
#> [15] backports_1.5.0 XVector_0.44.0
#> [17] labeling_0.4.3 utf8_1.2.4
#> [19] rmarkdown_2.29 tzdb_0.4.0
#> [21] UCSC.utils_1.0.0 purrr_1.0.2
#> [23] bit_4.5.0 xfun_0.49
#> [25] MultiAssayExperiment_1.30.3 zlibbioc_1.50.0
#> [27] ggseqlogo_0.2 GenomeInfoDb_1.40.1
#> [29] jsonlite_1.8.9 gmm_1.8
#> [31] DelayedArray_0.30.1 broom_1.0.7
#> [33] parallel_4.4.1 cluster_2.1.6
#> [35] R6_2.5.1 stringi_1.8.4
#> [37] car_3.1-3 GenomicRanges_1.56.2
#> [39] Rcpp_1.0.13-1 SummarizedExperiment_1.34.0
#> [41] knitr_1.49 zoo_1.8-12
#> [43] readr_2.1.5 IRanges_2.38.1
#> [45] Matrix_1.7-1 igraph_2.1.1
#> [47] tidyselect_1.2.1 rstudioapi_0.17.1
#> [49] abind_1.4-8 yaml_2.3.10
#> [51] lattice_0.22-6 tibble_3.2.1
#> [53] plyr_1.8.9 Biobase_2.64.0
#> [55] withr_3.0.2 evaluate_1.0.1
#> [57] tmvtnorm_1.6 norm_1.0-11.1
#> [59] pillar_1.9.0 BiocManager_1.30.25
#> [61] ggpubr_0.6.0 MatrixGenerics_1.16.0
#> [63] carData_3.0-5 stats4_4.4.1
#> [65] generics_0.1.3 vroom_1.6.5
#> [67] S4Vectors_0.42.1 hms_1.1.3
#> [69] ggplot2_3.5.1 munsell_0.5.1
#> [71] scales_1.3.0 BiocStyle_2.32.1
#> [73] glue_1.8.0 lazyeval_0.2.2
#> [75] tools_4.4.1 QFeatures_1.14.2
#> [77] ggsignif_0.6.4 imputeLCMD_2.1
#> [79] mvtnorm_1.3-2 cowplot_1.1.3
#> [81] grid_4.4.1 impute_1.78.0
#> [83] tidyr_1.3.1 MsCoreUtils_1.16.1
#> [85] colorspace_2.1-1 GenomeInfoDbData_1.2.12
#> [87] Formula_1.2-5 cli_3.6.3
#> [89] fansi_1.0.6 S4Arrays_1.4.1
#> [91] AnnotationFilter_1.28.0 pcaMethods_1.96.0
#> [93] gtable_0.3.6 rstatix_0.7.2
#> [95] digest_0.6.37 BiocGenerics_0.50.0
#> [97] SparseArray_1.4.8 farver_2.1.2
#> [99] htmltools_0.5.8.1 lifecycle_1.0.4
#> [101] httr_1.4.7 bit64_4.5.2
#> [103] MASS_7.3-61