A Curated Resource for Phosphosite-specific Signature Analysis

Abstract

Signaling pathways are orchestrated by post-translational modifications (PTMs) such as phosphorylation. However, pathway analysis of PTM data sets generated by mass spectrometry (MS)-based proteomics is typically performed at a gene-centric level because of the lack of appropriately curated PTM signature databases and bioinformatic tools that leverage PTM site-specific information. Here we present the first version of PTMsigDB, a database of modification site-specific signatures of perturbations, kinase activities and signaling pathways curated from more than 2,500 publications. We adapted the widely used single sample Gene Set Enrichment Analysis approach to utilize PTMsigDB, enabling PTMSignature Enrichment Analysis (PTM-SEA) of quantitative MS data. We used a well-characterized data set of epidermal growth factor (EGF)-perturbed cancer cells to evaluate our approach and demonstrated better representation of signaling events compared with gene-centric methods. We then applied PTM-SEA to analyze the phosphoproteomes of cancer cells treated with cell-cycle inhibitors and detected mechanism-of-action specific signatures of cell cycle kinases. We also applied our methods to analyze the phosphoproteomes of PI3K-inhibited human breast cancer cells and detected signatures of compounds inhibiting PI3K as well as targets downstream of PI3K (AKT, MAPK/ERK) covering a substantial fraction of the PI3K pathway. PTMsigDB and PTM-SEA can be freely accessed at https://github.com/broadinstitute/ssGSEA2.0.

Keywords: Computational Biology; Database design; Pathway Analysis; Phosphorylation; Post-translational modifications*.

Graphical abstract

Fig. 1.

Pathway analysis of phosphoproteome data sets. A, The majority of proteins have multiple phospho-acceptor sites that can be differentially occupied and vary in abundance levels measured by MS. The four examples illustrate three proteoforms with varying number of phosphorylation sites (represented by amino acid residue and position in protein sequence) which are quantified with different fold changes (represented as numbers in parenthesis). Protein A exemplifies the presence of two protein isoforms carrying an isoform-specific phosphorylation site. B, Gene-centric pathway analysis typically involves combining fold changes of multiple sites mapping to the same gene symbol by calculating a center value of abundance (mean or median) or by choosing a single site characterized by high degree of variance (thus information) across a sample cohort. Additional information provided by multiple phosphorylation sites on a single protein as well as different sites on protein isoforms are not considered. Resulting gene-centric expression matrix can then be queried against gene-centric pathway databases such as Reactome, KEGG or MSigDB, in which each pathway is represented as a collection of gene symbols. C, Site-centric pathway analysis takes all quantified phosphorylation sites into account and requires a database of pathways annotated using individual phosphorylation sites. For that purpose, we developed PTMsigDB, annotating signatures of pathways, kinases and perturbations at the level of sites. Each site is annotated with the direction of regulation in a signature, i.e. whether its abundance is decreased or elevated, exemplified by blue and red arrows, respectively. Quantitative, site-centric phosphoproteomic data can be directly queried against PTMsigDB to identify signatures of phosphosites that correlate with annotated signature sets in PTMsigDB.

Fig. 2.

PTMsigDB signature sets and scoring scheme. Signatures sets in PTMsigDB are divided into three categories (perturbations, kinases and signatures of molecular pathways) derived from four major sources (PhosphoSitePlus, NetPath, WikiPathways, LINCS). A, Perturbation signatures from PhosphoSitePlus were assembled from a total of 2483 publications. To reduce experimental noise we require that each site be consistently reported by at least two independent studies (consensus signatures). Experimental noise can be introduced by use of different cell systems (e.g. in vivo versus in vitro), different protocols (e.g. dosage of perturbation) or different technologies (e.g. low-throughput versus high throughput). B, Bar chart depicting the distribution of unidirectional (all sites of a signature are annotated with decreased OR increased abundance levels) and bidirectional (sites of a signature include both decreased AND increased abundance levels) signature sets. C, Modified scoring scheme of the single sample Gene Set Enrichment Analysis (ssGSEA) algorithm enables scoring of directional signature sets. For each direction a running-sum statistic (y axis) is calculated by walking down the ranked list of PTM sites in the data set (x axis). An enrichment score (ES) reflected by the area under the resulting curve is calculated for each direction (ES_u for “up” sites and ES_d for “down” sites) separately and combined by subtracting ES_d from ES_u. D, Heatmap depicting the overlap between phosphorylation sites commonly affected by kinase-inhibitors. The similarity matrix is based on the number of shared phosphosites between signatures. The upper and lower triangular matrices are normalized by the total number of sites of signatures listed in rows and columns, respectively.

Fig. 3.

Phosphoproteome signatures of EGF and nocodazole treatment in HeLa cells. We applied PTM-SEA to a phosphoproteome data set from the literature. Sharma et al. (18) studied phosphoproteome dynamics of EGF treatment in HeLa S3 cells. To mitotically arrest the cells, the antineoplastic agent nocodazole was used. Because PTMsigDB contained signatures for both perturbations, this data set was used as a benchmark to evaluate PTM-SEA. A, Heatmap depicting normalized enrichment scores (NES) of signatures (rows) in PTMsigDB that were consistently enriched or depleted at FDR < 0.01 in replicate measurements of EGF/nocodazole treatments (marked by asterisks). Hierarchical clustering of enrichment scores separated samples (columns) by experimental condition (EGF, nocodazole and control). Enrichment scores of nocodazole and EGF signature showed highest magnitude in the respective experiments and are highlighted in red. B, Silhouette analysis comparing hierarchical clustering of site-centric (left panel) and gene-centric (right panel) signature enrichment scores into three clusters. The bar chart depicts silhouette scores (x axis) of each sample (y axis) colored according to the assigned cluster. Average silhouette scores for each cluster are depicted at the right side of the bar charts. C, Signature enrichment scores of the nocodazole perturbation signature calculated using site-centric (left panel) and gene-centric (right panel) approach. Numbers indicate the median NES and FDR in replicate measurements of respective perturbations (x axis).

Fig. 4.

MoA-specific kinase signatures of cell cycle inhibitors. Heatmaps depicting perturbation and kinase signature enrichment scores of three human cancer cell lines (MCF7, HL60, PC3) treated with four different cell cycle inhibitors measured in duplicates. Asterisks indicate signature scores at FDR < 0.05. Cartoons above the heatmap graphically illustrate the increasing levels (blue to red) of CDK1/2 substrate phosphorylations. The density plot illustrates cyclin A/B/E concentrations across cell cycle stages that correlate with CDK1/2 activity levels. Black T-arrows indicate stage of cell cycle the corresponding inhibitor is active. Signatures of CDK1 and CDK2 activity strongly correlate with the mode of action of each compound.

Fig. 5.

Signatures of PI3K inhibition in breast cancer cells. T47D cells were treated with the PI3Ka inhibitor BYL719 and DMSO as control for 6h and 24h and deep phosphoproteomes were acquired resulting in ∼24K phosphosites localized to a specific Ser/Thr/Tyr residue. A, Volcano plot depicting enrichment of phosphoproteome signatures in the 6h time point. The x axis represents the normalized enrichment score (NES) between DMSO (left side) and drug treatment (right side). The size of the dots scale with the relative number of scored phosphorylation sites in a signature. The gray area contains signatures that did not significantly change upon drug treatment (permutation-based FDR ≥ 5%). PTM signature sets of inhibitors are annotated with an “i” after the drug target. B, Volcano plot depicting phospho signatures after treatment for 24 h. C, Schematic representation of the PI3K-AKT-mTOR and Ras-Raf-MEK-ERK pathways. Circles indicate kinases, boxes indicate specific kinase inhibitors. Highlighted in colors are significant phospho signatures shown in the volcano plots in A and B. The pathway representation is based on Fig. 1 in (44).