Introduction

While 1, 2, and 3D pathway analyses are useful for data generated from experiments with different treatment/conditions, analysis designed for time-course data may be better suited to analysis experiments that profile multiple time points.

Here, we will apply ClueR which is an R package specifically designed for time-course proteomic and phosphoproteomic data analysis Yang et al. 2015.

Loading packages and data

We will load the PhosR package with few other packages we will use for this tutorial.

suppressPackageStartupMessages({
  library(parallel)
  library(ggplot2)
  library(ClueR)
  library(reactome.db)
  library(org.Mm.eg.db)
  library(annotate)
  library(PhosR)
})

We will load a dataset integrated from two time-course datasets of early and intermediate insulin signalling in mouse liver upon insulin stimulation to demonstrate the time-course phosphoproteomic data analyses.

data("phospho.liver.Ins.TC.ratio.RUV.pe")
ppe <- phospho.liver.Ins.TC.ratio.RUV.pe
ppe
#> class: PhosphoExperiment 
#> dim: 800 90 
#> metadata(0):
#> assays(1): Quantification
#> rownames(800): LARP7;256; SRSF10;131; ... SIK3;493; GSK3A;21;
#> rowData names(0):
#> colnames(90): Intensity.Liver_Ins_0s_Bio7 Intensity.Liver_Ins_0s_Bio8
#>   ... Intensity.Liver_Ins_10m_Bio5 Intensity.Liver_Ins_10m_Bio6
#> colData names(0):

Gene-centric analyses of the liver phosphoproteome data

Let us start with gene-centric analysis. Such analysis can be directly applied to proteomics data. It can also be applied to phosphoproteomic data by using the phosCollapse function to summarise phosphosite information to proteins.

# take grouping information
grps <- sapply(strsplit(colnames(ppe), "_"), 
            function(x)x[3])

# select differentially phosphorylated sites
sites.p <- matANOVA(ppe@assays@data$Quantification, 
                    grps)
ppm <- meanAbundance(ppe@assays@data$Quantification, grps)
sel <- which((sites.p < 0.05) & (rowSums(abs(ppm) > 1) != 0))
ppm_filtered <- ppm[sel,]

# summarise phosphosites information into gene level
ppm_gene <- phosCollapse(ppm_filtered, 
            gsub(";.+", "", rownames(ppm_filtered)), 
                stat = apply(abs(ppm_filtered), 1, max), by = "max")

# perform ClueR to identify optimal number of clusters
pathways = as.list(reactomePATHID2EXTID)

pathways = pathways[which(grepl("R-MMU", names(pathways), ignore.case = TRUE))]

pathways = lapply(pathways, function(path) {
    gene_name = unname(getSYMBOL(path, data = "org.Mm.eg"))
    toupper(unique(gene_name))
})

RNGkind("L'Ecuyer-CMRG")
set.seed(123)
c1 <- runClue(ppm_gene, annotation=pathways, 
            kRange = seq(2,10), rep = 5, effectiveSize = c(5, 100), 
            pvalueCutoff = 0.05, alpha = 0.5)

# Visualise the evaluation results
data <- data.frame(Success=as.numeric(c1$evlMat), Freq=rep(seq(2,10), each=5))
myplot <- ggplot(data, aes(x=Freq, y=Success)) + 
    geom_boxplot(aes(x = factor(Freq), fill="gray")) +
    stat_smooth(method="loess", colour="red", size=3, span = 0.5) +
    xlab("# of cluster") + 
    ylab("Enrichment score") + 
    theme_classic()
myplot


set.seed(123)
best <- clustOptimal(c1, rep=5, mfrow=c(2, 3), visualize = TRUE)

Site-centric analyses of the liver phosphoproteome data

Phosphosite-centric analyses will perform using kinase-substrate annotation information from PhosphoSitePlus.

data("PhosphoSitePlus")
RNGkind("L'Ecuyer-CMRG")
set.seed(1)
PhosphoSite.mouse2 = mapply(function(kinase) {
  gsub("(.*)(;[A-Z])([0-9]+;)", "\\1;\\3", kinase)
}, PhosphoSite.mouse)

# perform ClueR to identify optimal number of clusters
c3 <- runClue(ppm_filtered, annotation=PhosphoSite.mouse2, kRange = 2:10, rep = 5, effectiveSize = c(5, 100), pvalueCutoff = 0.05, alpha = 0.5)

# Visualise the evaluation results
data <- data.frame(Success=as.numeric(c3$evlMat), Freq=rep(2:10, each=5))
myplot <- ggplot(data, aes(x=Freq, y=Success)) + geom_boxplot(aes(x = factor(Freq), fill="gray"))+
  stat_smooth(method="loess", colour="red", size=3, span = 0.5) + xlab("# of cluster")+ ylab("Enrichment score")+theme_classic()
myplot


set.seed(1)
best <- clustOptimal(c3, rep=10, mfrow=c(2, 3), visualize = TRUE)


# Finding enriched pathways from each cluster
best$enrichList
#> $`cluster 1`
#>      kinase   pvalue                 size
#> [1,] "PRKACA" "0.000184676866298047" "5" 
#>      substrates                                                
#> [1,] "NR1H3;196;|MARCKS;163;|PRKACA;339;|ITPR1;1755;|SIK3;493;"
#> 
#> $`cluster 3`
#>      kinase         pvalue                 size
#> [1,] "Humphrey.Akt" "0.000162969329853963" "5" 
#> [2,] "Yang.Akt"     "0.000165386907010959" "6" 
#>      substrates                                                         
#> [1,] "TSC2;939;|PFKFB2;486;|FOXO3;252;|FOXO1;316;|GSK3A;21;"            
#> [2,] "AKT1S1;247;|TSC2;939;|PFKFB2;486;|FOXO3;252;|FOXO1;316;|GSK3A;21;"

SessionInfo

sessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Catalina 10.15.7
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#> 
#> Random number generation:
#>  RNG:     L'Ecuyer-CMRG 
#>  Normal:  Inversion 
#>  Sample:  Rejection 
#>  
#> 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] stats4    parallel  stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] PhosR_1.5.1          annotate_1.70.0      XML_3.99-0.8        
#>  [4] org.Mm.eg.db_3.13.0  reactome.db_1.76.0   AnnotationDbi_1.54.1
#>  [7] IRanges_2.26.0       S4Vectors_0.30.2     Biobase_2.52.0      
#> [10] BiocGenerics_0.38.0  ClueR_1.4            e1071_1.7-9         
#> [13] ggplot2_3.3.5       
#> 
#> loaded via a namespace (and not attached):
#>   [1] colorspace_2.0-2            ggsignif_0.6.3             
#>   [3] ellipsis_0.3.2              class_7.3-19               
#>   [5] rprojroot_2.0.2             circlize_0.4.13            
#>   [7] XVector_0.32.0              GenomicRanges_1.44.0       
#>   [9] GlobalOptions_0.1.2         ggdendro_0.1.22            
#>  [11] fs_1.5.0                    proxy_0.4-26               
#>  [13] farver_2.1.0                ggpubr_0.4.0               
#>  [15] bit64_4.0.5                 fansi_0.5.0                
#>  [17] splines_4.1.2               cachem_1.0.6               
#>  [19] knitr_1.36                  broom_0.7.10               
#>  [21] png_0.1-7                   pheatmap_1.0.12            
#>  [23] BiocManager_1.30.16         compiler_4.1.2             
#>  [25] httr_1.4.2                  backports_1.3.0            
#>  [27] Matrix_1.3-4                fastmap_1.1.0              
#>  [29] limma_3.48.3                htmltools_0.5.2            
#>  [31] tools_4.1.2                 igraph_1.2.8               
#>  [33] coda_0.19-4                 gtable_0.3.0               
#>  [35] glue_1.5.0                  GenomeInfoDbData_1.2.6     
#>  [37] reshape2_1.4.4              dplyr_1.0.7                
#>  [39] Rcpp_1.0.7                  carData_3.0-4              
#>  [41] statnet.common_4.5.0        pkgdown_1.6.1              
#>  [43] vctrs_0.3.8                 Biostrings_2.60.2          
#>  [45] nlme_3.1-153                preprocessCore_1.54.0      
#>  [47] xfun_0.28                   stringr_1.4.0              
#>  [49] network_1.17.1              lifecycle_1.0.1            
#>  [51] rstatix_0.7.0               dendextend_1.15.2          
#>  [53] zlibbioc_1.38.0             MASS_7.3-54                
#>  [55] scales_1.1.1                BiocStyle_2.20.2           
#>  [57] ragg_1.1.3                  pcaMethods_1.84.0          
#>  [59] MatrixGenerics_1.4.3        SummarizedExperiment_1.22.0
#>  [61] RColorBrewer_1.1-2          yaml_2.2.1                 
#>  [63] memoise_2.0.0               gridExtra_2.3              
#>  [65] reshape_0.8.8               stringi_1.7.5              
#>  [67] RSQLite_2.2.8               highr_0.9                  
#>  [69] desc_1.4.0                  shape_1.4.6                
#>  [71] GenomeInfoDb_1.28.4         rlang_0.4.12               
#>  [73] pkgconfig_2.0.3             systemfonts_1.0.2          
#>  [75] matrixStats_0.61.0          bitops_1.0-7               
#>  [77] evaluate_0.14               lattice_0.20-45            
#>  [79] purrr_0.3.4                 ruv_0.9.7.1                
#>  [81] labeling_0.4.2              bit_4.0.4                  
#>  [83] tidyselect_1.1.1            GGally_2.1.2               
#>  [85] plyr_1.8.6                  magrittr_2.0.1             
#>  [87] R6_2.5.1                    generics_0.1.1             
#>  [89] DelayedArray_0.18.0         DBI_1.1.1                  
#>  [91] mgcv_1.8-38                 pillar_1.6.4               
#>  [93] withr_2.4.2                 KEGGREST_1.32.0            
#>  [95] abind_1.4-5                 RCurl_1.98-1.5             
#>  [97] tibble_3.1.6                crayon_1.4.2               
#>  [99] car_3.0-12                  utf8_1.2.2                 
#> [101] rmarkdown_2.11              viridis_0.6.2              
#> [103] grid_4.1.2                  blob_1.2.2                 
#> [105] digest_0.6.28               xtable_1.8-4               
#> [107] tidyr_1.1.4                 textshaping_0.3.6          
#> [109] munsell_0.5.0               viridisLite_0.4.0