Imputation of phosphoproteomics data with PhosR
2 December 2022
imputation.RmdIntroduction
PhosR is a package for the all-rounded analysis of phosphoproteomic data from processing to downstream analysis. This vignette will provide a step-by-step workflow of how PhosR can be used to process and analyse a panel of phosphoproteomic datasets. As one of the first steps of data processing, we will begin by performing filtering and imputation of the data with PhosR.
Loading packages and data
First, we will load the PhosR package. If you already haven’t done so, please install PhosR as instructed in the main page.
suppressPackageStartupMessages({
library(PhosR)
})Setting up the data
We assume that you will have the raw data processed using platforms frequently used for mass-spectrometry based proteomics such as MaxQuant. For demonstration purposes, we will take a parts of phosphoproteomic data generated by Humphrey et al. with accession number PXD001792. The dataset contains the phosphoproteomic quantifications of two mouse liver cell lines (Hepa1.6 and FL38B) that were treated with either PBS (mock) or insulin.
Let us load the PhosphoExperiment (ppe) object.
data("phospho.cells.Ins.pe")
ppe <- phospho.cells.Ins.pe
class(ppe)
#> [1] "PhosphoExperiment"
#> attr(,"package")
#> [1] "PhosR"A quick glance of the object.
ppe
#> class: PhosphoExperiment
#> dim: 5000 24
#> metadata(0):
#> assays(1): Quantification
#> rownames(5000): Q7TPV4;MYBBP1A;S1321;PQSALPKKRARLSLVSRSPSLLQSGVKKRRV
#> Q3UR85;MYRF;S304;PARAPSPPWPPQGPLSPGTGSLPLSIARAQT ...
#> P28659-4;NA;S18;AFKLDFLPEMMVDHCSLNSSPVSKKMNGTLD
#> E9Q8I9;FRY;S1380;HNIELVDSRLLLPGSSPSSPEDEVKDREGEV
#> rowData names(0):
#> colnames(24): Intensity.FL83B_Control_1 Intensity.FL83B_Control_2 ...
#> Intensity.Hepa1.6_Ins_5 Intensity.Hepa1.6_Ins_6
#> colData names(0):We will take the grouping information from colnames of our matrix.
For each cell line, there are two conditions (Control vs Insulin-stimulated) and 6 replicates for each condition.
# FL38B
gsub("Intensity.", "", grps)[1:12]
#> [1] "FL83B_Control" "FL83B_Control" "FL83B_Control" "FL83B_Control"
#> [5] "FL83B_Control" "FL83B_Control" "FL83B_Ins" "FL83B_Ins"
#> [9] "FL83B_Ins" "FL83B_Ins" "FL83B_Ins" "FL83B_Ins"
# Hepa1
gsub("Intensity.", "", grps)[13:24]
#> [1] "Hepa1.6_Control" "Hepa1.6_Control" "Hepa1.6_Control" "Hepa1.6_Control"
#> [5] "Hepa1.6_Control" "Hepa1.6_Control" "Hepa1.6_Ins" "Hepa1.6_Ins"
#> [9] "Hepa1.6_Ins" "Hepa1.6_Ins" "Hepa1.6_Ins" "Hepa1.6_Ins"Note that there are in total 24 samples and 5,000 phosphosites profiled.
dim(ppe)
#> [1] 5000 24Filtering of phosphosites
Next, we will perform some filtering of phosphosites so that only phosphosites with quantification for at least 50% of the replicates in at least one of the conditions are retained. For this filtering step, we use the selectGrps function. The filtering leaves us with 1,772 phosphosites.
ppe_filtered <- selectGrps(ppe, grps, 0.5, n=1)
dim(ppe_filtered)
#> [1] 1772 24selectGrps gives you the option to relax the threshold for filtering. The filtering threshold can therefore be optimized for each dataset.
# In cases where you have fewer replicates ( e.g.,triplicates), you may want to select phosphosites quantified in 70% of replicates.
ppe_filtered_v1 <- selectGrps(ppe, grps, 0.7, n=1)
dim(ppe_filtered_v1)
#> [1] 1330 24Imputation of phosphosites
We can proceed to imputation now that we have filtered for suboptimal phosphosites. To take advantage of data structure and experimental design, PhosR provides users with a lot of flexibility for imputation. There are three functions for imputation: scImpute,tInmpute, and ptImpute. Here, we will demonstrate the use of scImpute and ptImpute.
Site- and condition-specific imputation
The scImpute function is used for site- and condition-specific imputation. A pre-defined thereshold is used to select phosphosites to impute. Phosphosites with missing values equal to or greater than a predefined value will be imputed by sampling from the empirical normal distribution constructed from the quantification values of phosphosites from the same condition.
In the above example, only phosphosites that are quantified in more than 50% of samples from the same condition will be imputed.
Paired tail-based imputation
We then perform paired tail-based imputation on the dataset imputed with scImpute. Paired tail-based imputation performs imputation of phosphosites that have missing values in all replicates in one condition (e.g. in basal) but not in another condition (e.g., in stimulation). This method of imputation ensures that we do not accidentally filter phosphosites that seemingly have low detection rate.
As for scImpute, we can set a predefined threshold to in another condition (e.g. stimulation), the tail-based imputation is applied to impute for the missing values in the first condition.
set.seed(123)
ppe_imputed <- ppe_imputed_tmp
ppe_imputed[,seq(6)] <- ptImpute(ppe_imputed[,seq(7,12)],
ppe_imputed[,seq(6)],
percent1 = 0.6, percent2 = 0, paired = FALSE)
ppe_imputed[,seq(13,18)] <- ptImpute(ppe_imputed[,seq(19,24)],
ppe_imputed[,seq(13,18)],
percent1 = 0.6, percent2 = 0, paired = FALSE)Lastly, we perform normalisation of the filtered and imputed phosphoproteomic data.
ppe_imputed_scaled <- medianScaling(ppe_imputed, scale = FALSE, assay = "imputed")Quantification plots
A useful function in PhosR is to visualize the percentage of quantified sites before and after filtering and imputation. The main inputs of plotQC are the quantification matrix, sample labels (equating the column names of the matrix), an integer indicating the panel to plot, and lastly, a color vector. To visualize the percentage of quantified sites, use the plotQC function and set panel = quantify to visualise bar plots of samples.
p1 = plotQC(ppe_filtered@assays@data$Quantification, labels=colnames(ppe_filtered),
panel = "quantify", grps = grps)
p2 = plotQC(ppe_imputed_scaled@assays@data$scaled,
labels=colnames(ppe_imputed_scaled), panel = "quantify", grps = grps)
ggpubr::ggarrange(p1, p2, nrow = 1)
By setting panel = dendrogram, we can visualise the results of unsupervised hierarchical clustering of samples as a dendrogram. The dendrogram demonstrates that imputation has improved the clustering of the samples so that replicates from the same conditions cluster together.
p1 = plotQC(ppe_filtered@assays@data$Quantification,
labels=colnames(ppe_filtered), panel = "dendrogram",
grps = grps)
p2 = plotQC(ppe_imputed_scaled@assays@data$scaled,
labels=colnames(ppe_imputed_scaled),
panel = "dendrogram", grps = grps)
ggpubr::ggarrange(p1, p2, nrow = 1)
We can now move onto the next step in the PhosR workflow: integration of datasets and batch correction.
SessionInfo
sessionInfo()
#> R version 4.2.2 (2022-10-31)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur ... 10.16
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] PhosR_1.9.1 BiocStyle_2.26.0
#>
#> loaded via a namespace (and not attached):
#> [1] ggtext_0.1.2 bitops_1.0-7
#> [3] matrixStats_0.63.0 fs_1.5.2
#> [5] RColorBrewer_1.1-3 rprojroot_2.0.3
#> [7] GenomeInfoDb_1.34.3 tools_4.2.2
#> [9] backports_1.4.1 bslib_0.4.1
#> [11] utf8_1.2.2 R6_2.5.1
#> [13] BiocGenerics_0.44.0 colorspace_2.0-3
#> [15] withr_2.5.0 tidyselect_1.2.0
#> [17] gridExtra_2.3 GGally_2.1.2
#> [19] preprocessCore_1.60.0 compiler_4.2.2
#> [21] textshaping_0.3.6 cli_3.4.1
#> [23] Biobase_2.58.0 xml2_1.3.3
#> [25] network_1.18.0 desc_1.4.2
#> [27] DelayedArray_0.24.0 ggdendro_0.1.23
#> [29] labeling_0.4.2 bookdown_0.30
#> [31] sass_0.4.4 scales_1.2.1
#> [33] proxy_0.4-27 pkgdown_2.0.6
#> [35] commonmark_1.8.1 systemfonts_1.0.4
#> [37] stringr_1.4.1 digest_0.6.30
#> [39] rmarkdown_2.18 XVector_0.38.0
#> [41] pkgconfig_2.0.3 htmltools_0.5.3
#> [43] MatrixGenerics_1.10.0 highr_0.9
#> [45] ruv_0.9.7.1 limma_3.54.0
#> [47] fastmap_1.1.0 rlang_1.0.6
#> [49] GlobalOptions_0.1.2 farver_2.1.1
#> [51] shape_1.4.6 jquerylib_0.1.4
#> [53] generics_0.1.3 jsonlite_1.8.3
#> [55] statnet.common_4.7.0 dendextend_1.16.0
#> [57] dplyr_1.0.10 car_3.1-1
#> [59] RCurl_1.98-1.9 magrittr_2.0.3
#> [61] GenomeInfoDbData_1.2.9 Matrix_1.5-1
#> [63] Rcpp_1.0.9 munsell_0.5.0
#> [65] S4Vectors_0.36.0 fansi_1.0.3
#> [67] abind_1.4-5 viridis_0.6.2
#> [69] lifecycle_1.0.3 stringi_1.7.8
#> [71] yaml_2.3.6 carData_3.0-5
#> [73] MASS_7.3-58.1 SummarizedExperiment_1.28.0
#> [75] zlibbioc_1.44.0 plyr_1.8.8
#> [77] grid_4.2.2 lattice_0.20-45
#> [79] cowplot_1.1.1 gridtext_0.1.5
#> [81] circlize_0.4.15 knitr_1.41
#> [83] pillar_1.8.1 igraph_1.3.5
#> [85] ggpubr_0.5.0 GenomicRanges_1.50.1
#> [87] markdown_1.4 ggsignif_0.6.4
#> [89] reshape2_1.4.4 stats4_4.2.2
#> [91] glue_1.6.2 evaluate_0.18
#> [93] pcaMethods_1.90.0 BiocManager_1.30.19
#> [95] vctrs_0.5.1 gtable_0.3.1
#> [97] purrr_0.3.5 tidyr_1.2.1
#> [99] reshape_0.8.9 cachem_1.0.6
#> [101] ggplot2_3.4.0 xfun_0.35
#> [103] broom_1.0.1 e1071_1.7-12
#> [105] coda_0.19-4 rstatix_0.7.1
#> [107] ragg_1.2.4 class_7.3-20
#> [109] viridisLite_0.4.1 tibble_3.1.8
#> [111] pheatmap_1.0.12 memoise_2.0.1
#> [113] IRanges_2.32.0