1 Load PRONE Package

library(PRONE)

2 Load Data (TMT)

Here, we are directly working with the SummarizedExperiment data. For more information on how to create the SummarizedExperiment from a proteomics data set, please refer to the “Get Started” vignette.

The example TMT data set originates from (Biadglegne et al. 2022).

data("tuberculosis_TMT_se")
se <- tuberculosis_TMT_se

3 Normalize Data

In order to compare the performance of different normalization methods on their ability to detect differentially expressed proteins, we performed some normalization with batch effect correction technqiues here. For more details about how to normalize data and evaluate the normalization approaches quantitatively and qualitatively using PRONE, please refer to the “Normalization” vignette.

se_norm <- normalize_se(se, c("IRS_on_RobNorm", "IRS_on_Median", 
                              "IRS_on_LoessF", "IRS_on_Quantile"), 
                        combination_pattern = "_on_")
#> IRS normalization performed on RobNorm-normalized data completed.
#> IRS normalization performed on Median-normalized data completed.
#> LoessF normalization not yet performed. Single LoessF normalization performed now.
#> LoessF completed.
#> IRS normalization performed on LoessF-normalized data completed.
#> Quantile normalization not yet performed. Single Quantile normalization performed now.
#> Quantile completed.
#> IRS normalization performed on Quantile-normalized data completed.

4 Differential Expression Analysis

After having performed normalization and evaluated the different normalization methods via qualitative and quantitative analysis, differential expression analysis can be used to further analyze the differences of the normalization methods and assess the impact of normalization on downstream analyses.

However before, you need to remove the reference samples in case of a TMT experiment. This can be easily done with the function remove_reference_samples().

se_norm <- remove_reference_samples(se_norm)
#> 2 reference samples removed from the SummarizedExperiment object.

4.1 Run DE Analysis

4.1.1 Specific Comparisons

First, you need to specify the comparisons you want to perform in DE analysis. For this, the function specify_comparisons() was developed which helps to build the right comparison strings.

However, you can also just simply create a vector of comparisons to ensure the correct order and handle this vector over to the DE analysis method.

comparisons <- specify_comparisons(se_norm, condition = "Group", 
                                   sep = NULL, control = NULL)

comparisons <- c("PTB-HC", "TBL-HC", "TBL-PTB", "Rx-PTB")

The function specify_comparisons() becomes handy when having a lot of sample groups and foremost samples of multiple conditions (not the case for this dataset).

For instance, you have a data set with diabetic and healthy samples (condition) that were measured at day 3 (3d), day 7 (7d), and day 17 (d14) after operation (timepoints. In DE analysis, you want to perform comparisons between the diabetic and healthy samples at each time point but also between the different time points for a fixed condition (e.g. diabetic). Instead of writing all comparisons manually, the function helps to build the right comparison strings by only considering comparisons where at least one of the groups (timepoint or condition) remains static. For this you only need to have a column in your metadata named for instance “Condition” that combines the two groups by “”, with values diabetic_3d, diabetic_7d, diabetic_14d, etc. And in the function you specify the column name ”Condition” and the separator (sep = ””).

4.1.2 Perform DE Analysis

The function run_DE() performs the DE analysis on the selected SummarizedExperiment and comparisons. A novel feature of PRONE is to compare the DE results of different normalization methods as it has been shown that normalization has an impact on downstream analyses.

Therefore, DE analysis can be performed on multiple assays (normalization methods) at once using the already known “ain” parameter.

In the following, a more detailed explanation of the other parameters is given:

  • The condition of the SummarizedExperiment object, specified at the beginning, can be used (condition = NULL) or any other column of the meta data can specified.
  • Three methods are available for DE analysis: limma (Ritchie et al. 2015), DEqMS (Zhu et al. 2020), and ROTS (Suomi et al. 2017).
  • The logFC parameter specifies if the log2 fold change should be calculated.
  • The logFC_up and logFC_down parameters specify the log2 fold change thresholds for up- and down-regulated proteins.
  • The p_adj parameter specifies if the p-values should be adjusted.
  • The alpha parameter specifies the significance level for the p-values.

In addition, the following parameters are specific for the DE analysis methods:

  • Limma (DOI: 10.18129/B9.bioc.limma):
    • The covariate parameter can be used to include a covariate in the DE analysis in case limma is performed. This can be any column of the meta data.
    • The trend parameter specifies if the trend test should be performed in case limma is used.
    • The robust parameter specifies if the robust method should be used in case limma is used.
  • ROTS (DOI: 10.18129/B9.bioc.ROTS):
    • B is the number of bootstrapping and K is the number of top-ranked features for reproducibility optimization.
  • DEqMS (DOI: 10.18129/B9.bioc.DEqMS)
    • The DEqMS_PSMs_column need to be specified. This can be any column in the rowData(se) that represents the number of quantified peptides or PSMs. It is used in DEqMS to estimate prior variance for proteins quantified by different number of PSMs.
de_res <- run_DE(se = se_norm, 
                     comparisons = comparisons,
                     ain = NULL, 
                     condition = NULL, 
                     DE_method = "limma", 
                     logFC = TRUE, 
                     logFC_up = 1, 
                     logFC_down = -1, 
                     p_adj = TRUE, 
                     alpha = 0.05,
                     covariate = NULL, 
                     trend = TRUE,
                     robust = TRUE,
                     B = 100,
                     K = 500
                 ) 
#> Condition of SummarizedExperiment used!
#> All assays of the SummarizedExperiment will be used.
#> DE Analysis will not be performed on raw data.
#> log2: DE analysis completed.
#> RobNorm: DE analysis completed.
#> IRS_on_RobNorm: DE analysis completed.
#> Median: DE analysis completed.
#> IRS_on_Median: DE analysis completed.
#> LoessF: DE analysis completed.
#> IRS_on_LoessF: DE analysis completed.
#> Quantile: DE analysis completed.
#> IRS_on_Quantile: DE analysis completed.

If you want to apply other logFC or p-value threshold, there is no need to re-run the DE analysis again. With apply_thresholds(), you can simply change the threshold values.

For instance, if you want to not apply a logFC threshold and only consider proteins with an adjusted p-value of 0.05 as DE, just set the logFC parameter to FALSE. In this case, the proteins are not classified into up- and down-regulated but only as significant change or no change.

new_de_res <- apply_thresholds(de_res = de_res, logFC = FALSE, p_adj = TRUE, 
                               alpha = 0.05)

However, if you still want to see if a protein has a positive or negative logFC, you can set the logFC parameter to TRUE with the logFC_up and logFC_down parameters to 0.

new_de_res <- apply_thresholds(de_res = de_res, logFC = TRUE, 
                               logFC_up = 0, logFC_down = 0, 
                               p_adj = TRUE, alpha = 0.05)

4.2 Visualize DE Results

4.2.1 Barplot

To get an overview of the DE results of the different normalization methods, you can visualize the number of significant DE proteins per normalization method in a barplot using plot_overview_DE_bar(). This plot can be generated in different ways by specifying the “plot_type” parameter:

  • single: plot a single barplot for each comparison (default)
  • facet_comp: facet the barplot by comparison
  • stacked: stack the number of DE per comparison
  • facet_regulation: stack the number of DE per comparison but facet by up- and down-regulated
plot_overview_DE_bar(de_res, ain = NULL, comparisons = comparisons, 
                     plot_type = "facet_regulation")
#> All normalization methods of de_res will be visualized.
Barplot showing the number of DE proteins per normalization method colored by comparison and facetted by up- and down-regulation.

Figure 4.1: Barplot showing the number of DE proteins per normalization method colored by comparison and facetted by up- and down-regulation.

You can also just visualize two specific comparisons:

plot_overview_DE_bar(de_res, ain = NULL, comparisons = comparisons[seq_len(2)], 
                     plot_type = "facet_comp")
#> All normalization methods of de_res will be visualized.
Barplot showing the number of DE proteins per normalization method faceted by comparison.

Figure 4.2: Barplot showing the number of DE proteins per normalization method faceted by comparison.

4.2.1.1 Tile Plot

You can also get an overview of the DE results in form of a heatmap using the plot_overview_DE_tile().

plot_overview_DE_tile(de_res)
#> All comparisons of de_res will be visualized.
#> All normalization methods of de_res will be visualized.