Differential Expression Analysis
Lis Arend
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_se3 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.
- The DEqMS_PSMs_column need to be specified. This can be any column in the
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.
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.
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.
Figure 4.3: Heatmap showing the number of DE proteins per comparison and per normalization method.
4.2.1.2 Volcano Plots
Another option is to generate volcano plots for each comparison. The function plot_volcano_DE() generates a grid of volcano plots for each normalization techniques (facet_norm = TRUE) or for each comparison (facet_comparison = TRUE). A list of volcano plots is returned.
plot_volcano_DE(de_res, ain = NULL, comparisons = comparisons[1],
facet_norm = TRUE)
#> All normalization methods of de_res will be visualized.
#> $`PTB-HC`
Figure 4.4: Volcano plots per normalization method for a single comparison.
4.2.1.3 Heatmap of DE Results
Furthermore, you can visualize the DE results in form of a heatmap. The function plot_heatmap_DE() generates a heatmap of the DE results for a specific comparison and normalization method.
plot_heatmap_DE(se_norm, de_res, ain = c("RobNorm", "IRS_on_RobNorm"),
comparison = "PTB-HC", condition = NULL, label_by = NULL,
pvalue_column = "adj.P.Val")
#> Label of SummarizedExperiment used!
#> Condition of SummarizedExperiment used!
#> <simpleError in stats::hclust(stats::dist(as.matrix(data))): NA/NaN/Inf in foreign function call (arg 10)>
#> <simpleError in stats::hclust(stats::dist(as.matrix(data))): NA/NaN/Inf in foreign function call (arg 10)>
#> $RobNorm
Figure 4.5: Individual heatmap of the DE results for a specific comparison and a selection of normalization methods. The adjusted p-values are added as row annotation, while the condition of each sample is shown as column annotation.
#>
#> $IRS_on_RobNorm
Figure 4.6: Individual heatmap of the DE results for a specific comparison and a selection of normalization methods. The adjusted p-values are added as row annotation, while the condition of each sample is shown as column annotation.
4.2.2 Intersection Analysis of DE Results
Moreover, you can also intersect the DE results of different normalization methods to see how many DE proteins overlap. You can either plot for each requested comparison an individual upset plot (plot_type = “single”) or stack the number of overlapping DE proteins per comparison (“stacked”). Not only the upset plot(s) are returned, but also a table with the intersections is provided by the plot_upset_DE()function.
intersections <- plot_upset_DE(de_res, ain = NULL,
comparisons = comparisons[seq_len(3)],
min_degree = 6, plot_type = "stacked")
#> All normalization methods of de_res will be visualized.
# put legend on top due to very long comparisons
intersections$upset[[2]] <- intersections$upset[[2]] +
ggplot2::theme(legend.position = "top", legend.direction = "vertical")
intersections$upset
Figure 4.7: Upset plot showing the overlapping DE proteins of different normalization methods colored by comparison.
Additionally, the Jaccard similarity index can be calculated to quantify the similarity of the DE results between the different normalization methods. A individual heatmap can be generated for each comparison (“plot_type = single”), a single heatmap facetted by comparison (“facet_comp”) or a single heatmap taking all comparisons into account (“all”) can be generated.
plot_jaccard_heatmap(de_res, ain = NULL, comparisons = comparisons,
plot_type = "all")
#> All normalization methods of de_res will be visualized.
Figure 4.8: Heatmap showing the Jaccard similarity indices of the DE results between different normalization methods for all comparisons.
PRONE offers the functionality to extract a consensus set of DEPs based on a selection of normalization methods and a threshold for the number of methods that need to agree on the DE status of a protein. The function get_consensus_DE() returns a list of consensus DE proteins for either each comparison separately or for all comparisons combined.
DT::datatable(extract_consensus_DE_candidates(de_res, ain = NULL,
comparisons = comparisons,
norm_thr = 0.8,
per_comparison = TRUE),
options = list(scrollX = TRUE))5 Session Info
utils::sessionInfo()
#> R Under development (unstable) (2025-01-20 r87609)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] PRONE_1.1.2 SummarizedExperiment_1.37.0
#> [3] Biobase_2.67.0 GenomicRanges_1.59.1
#> [5] GenomeInfoDb_1.43.3 IRanges_2.41.2
#> [7] S4Vectors_0.45.2 BiocGenerics_0.53.5
#> [9] generics_0.1.3 MatrixGenerics_1.19.1
#> [11] matrixStats_1.5.0
#>
#> loaded via a namespace (and not attached):
#> [1] rlang_1.1.5 magrittr_2.0.3
#> [3] GetoptLong_1.0.5 clue_0.3-66
#> [5] compiler_4.5.0 png_0.1-8
#> [7] vctrs_0.6.5 reshape2_1.4.4
#> [9] stringr_1.5.1 ProtGenerics_1.39.2
#> [11] shape_1.4.6.1 pkgconfig_2.0.3
#> [13] crayon_1.5.3 fastmap_1.2.0
#> [15] magick_2.8.5 XVector_0.47.2
#> [17] labeling_0.4.3 rmarkdown_2.29
#> [19] UCSC.utils_1.3.1 preprocessCore_1.69.0
#> [21] purrr_1.0.2 xfun_0.50
#> [23] MultiAssayExperiment_1.33.6 cachem_1.1.0
#> [25] jsonlite_1.8.9 DelayedArray_0.33.4
#> [27] BiocParallel_1.41.0 parallel_4.5.0
#> [29] cluster_2.1.8 R6_2.5.1
#> [31] RColorBrewer_1.1-3 bslib_0.8.0
#> [33] stringi_1.8.4 ComplexUpset_1.3.3
#> [35] limma_3.63.3 jquerylib_0.1.4
#> [37] Rcpp_1.0.14 bookdown_0.42
#> [39] iterators_1.0.14 knitr_1.49
#> [41] BiocBaseUtils_1.9.0 Matrix_1.7-2
#> [43] igraph_2.1.4 tidyselect_1.2.1
#> [45] abind_1.4-8 yaml_2.3.10
#> [47] doParallel_1.0.17 ggtext_0.1.2
#> [49] codetools_0.2-20 affy_1.85.1
#> [51] lattice_0.22-6 tibble_3.2.1
#> [53] plyr_1.8.9 withr_3.0.2
#> [55] evaluate_1.0.3 xml2_1.3.6
#> [57] circlize_0.4.16 pillar_1.10.1
#> [59] affyio_1.77.2 BiocManager_1.30.25
#> [61] DT_0.33 foreach_1.5.2
#> [63] MSnbase_2.33.2 MALDIquant_1.22.3
#> [65] ncdf4_1.23 ggplot2_3.5.1
#> [67] munsell_0.5.1 scales_1.3.0
#> [69] glue_1.8.0 lazyeval_0.2.2
#> [71] tools_4.5.0 data.table_1.16.4
#> [73] mzID_1.45.0 QFeatures_1.17.1
#> [75] vsn_3.75.0 mzR_2.41.1
#> [77] XML_3.99-0.18 Cairo_1.6-2
#> [79] grid_4.5.0 impute_1.81.0
#> [81] tidyr_1.3.1 crosstalk_1.2.1
#> [83] MsCoreUtils_1.19.0 colorspace_2.1-1
#> [85] patchwork_1.3.0 GenomeInfoDbData_1.2.13
#> [87] PSMatch_1.11.0 cli_3.6.3
#> [89] S4Arrays_1.7.1 ComplexHeatmap_2.23.0
#> [91] dplyr_1.1.4 AnnotationFilter_1.31.0
#> [93] pcaMethods_1.99.0 gtable_0.3.6
#> [95] sass_0.4.9 digest_0.6.37
#> [97] SparseArray_1.7.4 htmlwidgets_1.6.4
#> [99] rjson_0.2.23 farver_2.1.2
#> [101] htmltools_0.5.8.1 lifecycle_1.0.4
#> [103] httr_1.4.7 GlobalOptions_0.1.2
#> [105] statmod_1.5.0 gridtext_0.1.5
#> [107] MASS_7.3-64References
Biadglegne, Fantahun, Johannes R. Schmidt, Kathrin M. Engel, Jörg Lehmann, Robert T. Lehmann, Anja Reinert, Brigitte König, Jürgen Schiller, Stefan Kalkhof, and Ulrich Sack. 2022. “Mycobacterium Tuberculosis Affects Protein and Lipid Content of Circulating Exosomes in Infected Patients Depending on Tuberculosis Disease State.” Biomedicines 10 (4): 783. https://doi.org/10.3390/biomedicines10040783.
Ritchie, Matthew E., Belinda Phipson, Di Wu, Yifang Hu, Charity W. Law, Wei Shi, and Gordon K. Smyth. 2015. “Limma Powers Differential Expression Analyses for Rna-Sequencing and Microarray Studies.” Nucleic Acids Research 43 (7): e47–e47. https://doi.org/10.1093/nar/gkv007.
Suomi, Tomi, Fatemeh Seyednasrollah, Maria K. Jaakkola, Thomas Faux, and Laura L. Elo. 2017. “ROTS: An R Package for Reproducibility-Optimized Statistical Testing.” Edited by Timothée Poisot. PLOS Computational Biology 13 (5): e1005562. https://doi.org/10.1371/journal.pcbi.1005562.
Zhu, Yafeng, Lukas M. Orre, Yan Zhou Tran, Georgios Mermelekas, Henrik J. Johansson, Alina Malyutina, Simon Anders, and Janne Lehtiö. 2020. “Deqms: A Method for Accurate Variance Estimation in Differential Protein Expression Analysis.” Molecular & Cellular Proteomics 19 (6): 1047–57. https://doi.org/10.1074/mcp.TIR119.001646.