Integration of pan-cancer transcriptomics with RPPA proteomics reveals mechanisms of epithelial-mesenchymal transition

Abstract

Integrating data from multiple regulatory layers across cancer types could elucidate additional mechanisms of oncogenesis. Using antibody-based protein profiling of 736 cancer cell lines, along with matching transcriptomic data, we show that pan-cancer bimodality in the amounts of mRNA, protein, and protein phosphorylation reveals mechanisms related to the epithelial-mesenchymal transition (EMT). Based on the bimodal expression of E-cadherin, we define an EMT signature consisting of 239 genes, many of which were not previously associated with EMT. By querying gene expression signatures collected from cancer cell lines after small-molecule perturbations, we identify enrichment for histone deacetylase (HDAC) inhibitors as inducers of EMT, and kinase inhibitors as mesenchymal-to-epithelial transition (MET) promoters. Causal modeling of protein-based signaling identifies putative drivers of EMT. In conclusion, integrative analysis of pan-cancer proteomic and transcriptomic data reveals key regulatory mechanisms of oncogenic transformation.

Fig 1. Pan-cancer cell line data from CCLE transcriptomic and reverse phase protein arrays (RPPA) cluster by tissue of origin and E-cadherin expression but not by prior metastasis classification.

(A) Overlap of available RPPA and CCLE data with regard to cancer cell lines (left), measured transcripts and proteins (middle), and proteins measured for both basal expression and phosphorylation levels (right). The colored areas indicate data used to calculate and compare Euclidean distances between cell lines. (B) t-SNE plots of overlapping cancer cell lines based on protein, transcript, and equally weighted combined data. Each point represents a cell line and is colored by the tissue of origin (top), E-cadherin expression (middle), or tumor classification (bottom). NS: not specified. (C) Comparing pairwise distances between all cell lines using a linear model at the mRNA or protein levels. The red points show the top-100 highest residuals of cell line pairs, and the blue points the top-100 lowest residuals. (D) Dendrograms of breast cancer cell lines mapped for transcriptomic and RPPA data. The leaves of the trees were arranged to minimize the number of crossing lines between leaves of the two trees. L1-5 represents clusters found within the luminal subtype of breast cancer cell lines.

Fig 2. Bimodal protein expression and phosphorylation detected across cancer types associate with known oncogenic processes including EMT.

(A) Two-component Gaussian mixture model fit to E-cadherin protein expression. The lines indicate the probability density contribution from the low (-) and high (+) expression components. The histogram represents the RPPA measurements for the cell lines. (B) By comparing a two- versus one-component fit using the Bayesian Information Criterion (BIC), 260 out of 450 RPPA measurements supported bimodal expression. (C) Heat map of the posterior probabilities of each cell line belonging to the low (-, blue) or high (+, red) mixture component for the top-20 most bimodal proteins. The posterior probabilities can be thought of as soft assignments for the cell lines to low or high expression. Shannon entropy of the tissues assigned to low and high expression quantify the tissue diversity giving rise to the bimodal fits. (D) Overview of classification approach of proteins in terms of bimodality, tissue diversity (Shannon entropy), and frequency of cell lines assigned to the fitted distributions. (E) Significant GO terms for common bimodal proteins that were not found to be significant for non-bimodal proteins (p < 0.05, Benjamini-Hochberg).

Fig 3. Bimodal coupling between regulatory layers.

(A) Comparison of mRNA and protein bimodality colored by the coupling (Spearman’s correlation) between posterior probabilities of two-component Gaussian mixture models. The red percentage indicates the fraction of compared genes with coupled (r_b > 0.5) bimodalities at the transcript- and protein level. (B) Bimodal coupling of phosphosites and protein expression. (C) High confidence assignments (p < 0.1) to low or high expression for selection of bimodally coupled mRNA-protein pairs. (D) Scatter plot of E-cadherin mRNA and protein expression, indicating in red the 30 cell lines assigned to high mRNA but low protein expression (+–). r is Pearson’s correlation and r_b the bimodal coupling coefficient. (E) Tissue of origin and CDH1 (E-cadherin) mutational status of 30 cell lines with high E-cahderin mRNA but low protein expression. Of these cell lines, 3 out of 4 cell lines genotyped for CDH1 in COSMIC, all had mutations in the coding sequence. fs: frameshift, *: missense.

Fig 4. Proteins and phosphosites with coupled bimodality form network communities associated with EMT and intermediate transitions.

(A) Pan-cancer protein communities detected by Spearman’s correlation of the posterior probabilities of cell lines having low or high expression (|r_b| > 0.3). Only RPPA measurements associated with bimodal fits with high tissue diversity were included. Network communities were detected by calculating the leading non-negative eigenvector according to Newman’s method. Only edges within identified communities are shown, colored by the magnitude of the bimodal coupling coefficients. The size and color of the nodes represent the fitted mixing parameters from the Gaussian mixture models, quantifying whether underlying switches are common or rare in cancer cell lines. Each community was manually named according to plausible biological mechanisms by conducting a literature search for their protein members. Asterisks (*) indicate proteins with reported mechanisms linked to EMT. (B) Proposed interpretation of a two-step transition from the endothelial–E, to the mesenchymal–M states through two identified modules: EMT1 and EMT3. (C) Supporting protein expression data showing that Claudin7 and E-cadherin are correlated. Cell lines are colored by the tissue of origin (see Fig 3 for tissue labels).

Fig 5. Pan-cancer bimodal coupling between E-cadherin protein expression and genome-wide transcripts defines an EMT signature, predicting EMT- and MET-inducing small-molecules.

(A) Top-25 transcripts in CCLE with the strongest positive and top 24 negative bimodal coupling coefficients (r_b) to E-cadherin protein expression. Red squares indicate previous EMT signature genes in non-small cell lung carcinoma published by Byers et al. [51]. To define an EMT signature, we considered transcripts with |r_b| > 0.5, resulting in 215 epithelial and 24 mesenchymal markers. (B) Distribution of bimodal coupling coefficients, showing that E-cadherin coupling coefficients are shifted towards negative values compared to all measured proteins. (C) Overlap of EMT signature with previously published transcriptomic EMT signatures. The ‘mesenchymal’ bar plot is for the inversely correlated (coupled) genes and the ‘epithelial’ for the positively correlated genes. (D-E) Gene set enrichment analysis of epithelial part of the EMT signature. The TF enrichment analysis used ChIP-seq data to predict TFs involved in the regulation of the epithelial genes. The pie charts indicate the fraction of the signature genes associated with significantly enriched terms or TFs. (F) Small-molecule perturbations predicted to induce EMT and MET based on L1000 cell line data and the L1000CDS² method. The top-50 signatures are shown with results from multiple cell lines or concentrations aggregated by boxplots. PK: protein kinase.

Fig 6. Bayesian networks of proteins and phosphosites inferred from pan-cancer cell lines identify drivers of EMT and correlate to tumor networks.

All Bayesian network structures were inferred by a Fast Greedy Search algorithm. (A) Network centrality statistics of the directed causal graph over all measured proteins pertaining to the influence of proteins on cancer signaling. (B) In and out degree distributions of proteins and phosphosites from the inferred network. (C) Causal neighborhood (1st neighbors) of EMT markers E-cadherin, Rab25, and Claudin7. Tissue-specific correlations in support for each edge are shown as bars along the edges. The layout was determined using the hierarchical Sugiyama algorithm with all edges oriented downwards. (D) Network comparisons between cell line and tumor data. (E) Distribution of average connectivity in bootstrapped Bayesian networks. (F) Distribution of networks of bimodal coupling coefficients.