A pan-cancer proteomic perspective on The Cancer Genome Atlas.
Journal: 2015/October - Nature Communications
ISSN: 2041-1723
Abstract:
Protein levels and function are poorly predicted by genomic and transcriptomic analysis of patient tumours. Therefore, direct study of the functional proteome has the potential to provide a wealth of information that complements and extends genomic, epigenomic and transcriptomic analysis in The Cancer Genome Atlas (TCGA) projects. Here we use reverse-phase protein arrays to analyse 3,467 patient samples from 11 TCGA 'Pan-Cancer' diseases, using 181 high-quality antibodies that target 128 total proteins and 53 post-translationally modified proteins. The resultant proteomic data are integrated with genomic and transcriptomic analyses of the same samples to identify commonalities, differences, emergent pathways and network biology within and across tumour lineages. In addition, tissue-specific signals are reduced computationally to enhance biomarker and target discovery spanning multiple tumour lineages. This integrative analysis, with an emphasis on pathways and potentially actionable proteins, provides a framework for determining the prognostic, predictive and therapeutic relevance of the functional proteome.
Relations:
Content
Citations
(122)
References
(57)
Diseases
(1)
Chemicals
(2)
Organisms
(1)
Processes
(4)
Affiliates
(4)
Similar articles
Articles by the same authors
Discussion board
Nature communications. Dec/31/2013; 5: 3887-3887
Published online May/28/2014

A pan-cancer proteomic perspective on The Cancer Genome Atlas

+20 authors

Abstract

Introduction

The Cancer Genome Atlas (TCGA) is generating comprehensive molecular profiles for each of at least 33 different human tumor types (http://cancergenome.nih.gov). The overarching goal is to elucidate the landscape of DNA and RNA aberrations within and across tumor lineages and integrate the information with clinical characteristics, including patient outcome.

Previous studies have indicated only a partial concordance between genomic copy number, RNA levels and protein levels in both patient samples and cell lines1, 2, 3 at least, in part, because protein levels and, in particular, phosphoprotein levels represent an integration of the complex genomic and transcriptomic aberrations accumulated in each tumor combined with translational and post-translational regulation that cannot be fully captured by genomic and transcriptomic analysis. Hence, functional protein analysis using reverse-phase protein arrays (RPPA), which are highly applicable to study the large numbers of TCGA samples, was added to the TCGA effort to integrate proteomic characterization of tumors with already available genomic, transcriptomic and clinical information. The Clinical Proteomic Tumor Analysis Consortium (CPTAC, http://proteomics.cancer.gov/programs/cptacnetwork) is starting to use mass spectrometry to analyze a large fraction of the human proteome for a select subset of TCGA tumors. However, a comprehensive mass spectrometry analysis across all TCGA samples is not likely to be available in the near future. Thus, while earlier TCGA analyses were primarily based on genomic and transcriptomic characteristics4, 5, 6, 7, 8, 9, 10, the current study is driven by proteomic processes within and across cancer types.

Here we report an RPPA-based proteomic analysis using 181 high-quality antibodies that target total (n=128), cleaved (n=1), acetylated (n=1) and phosphorylated forms (n=51) of proteins in 3,467 TCGA patient samples across 11 “Pan-Cancer” tumor types. The function space covered by the antibodies used in the RPPA analysis includes proliferation, DNA damage, polarity, vesicle function, EMT, invasiveness, hormone signaling, apoptosis, metabolism, immunological, and stromal function as well as transmembrane receptors, integrin, TGFβ, LKB1/AMPK, TSC/mTOR, PI3K/Akt, Ras/MAPK, Hippo, Notch, and Wnt/beta-catenin signaling. Thus, the function space encompasses major functional and signaling pathways of relevance to human cancer. The TCGA tumor types included are those with mature RPPA data: breast cancer (BRCA, n=747), colon (COAD, n=334) and rectal (READ, n=130) adenocarcinoma, renal clear cell carcinoma (KIRC, n=454), high-grade serous ovarian cystadenocarcinoma (OVCA, n=412), uterine corpus endometrial carcinoma (UCEC, n=404), lung adenocarcinoma (LUAD, n=237), head and neck squamous cell carcinoma (HNSC, n=212), lung squamous cell carcinoma (LUSC, n=195), bladder urothelial carcinoma (BLCA, n=127) and glioblastoma multiforme (GBM, n=215) 4, 5, 6, 7, 8, 9, 10. We show that the functional proteome gives important, independent insights in TCGA data that are not captured by genomics or transcriptomics. Although samples predominantly cluster by tumor lineage, we also show that part of the tissue dominant effects can be removed computationally to elucidate common processes driving cellular behavior across tumor lineages. We present proteins and pathways that correlate with outcomes within certain tumor lineages and we identify multiple protein links and proteins that are associated with pathway activation. Taken together, the data and analytical resources presented in this manuscript are aimed at facilitating future research for targeted therapies that span multiple tumors.

Results

Correlations between protein and other data types

Protein data for 3,467 samples across 11 diseases were compared to mRNA, miRNA, copy number, and mutation data for the same samples. A novel approach, called “replicates-based normalization” (RBN, Methods), mitigated batch effects facilitating creation of a single Pan-Cancer protein dataset merging samples across 6 different batches. The RBN output is equivalent to all 3,467 samples being run in a single batch. In contrast to random (trans) protein:mRNA pairs (mean Spearman’s ρ = −0.006), almost half of matched (cis) protein:mRNA pairs in the RBN set demonstrated correlation beyond that expected by chance (mean Spearman’s ρ = 0.3) in both the overall Pan-Cancer dataset (t-test P < 2.2e-16, n=206 matched protein:mRNA pairs) and within particular diseases (Fig. 1a, Supplementary Fig. 1, Supplementary Data 1,2). Approximately 44% of matched (cis) protein:mRNA pairs had a correlation >= 0.3. For micro-RNAs, as expected, (trans) protein:miRNA correlations were much weaker with a mean positive Spearman’s ρ = 0.07, and a mean negative Spearman’s ρ = −0.07 (Supplementary Data 3). On the other hand, (trans) protein:protein correlations, including phosphoproteins, were higher (mean positive Spearman’s ρ = 0.15, mean negative Spearman’s ρ = −0.13, Supplementary Data 4). Detailed protein:protein and phosphoprotein:protein correlations across the total dataset and in particular diseases are available at the TCPA portal11. The results show, not surprisingly, that matched (cis) mRNA:protein correlations were the highest on average (ρ = 0.3), followed by (trans) protein:protein correlations (ρ ≈ ±0.15), whereas (trans) protein:miRNA correlations were lowest on average (ρ = ±0.07).

A similar analysis for CNV vs. protein fold change showed a mean fold change of 1.05 for amplifications and 0.95 for deletions in cis (Supplementary Data 5,6). Mutation vs. protein (cis) analysis showed a mean fold change of 1.2 for mutations that increased expression, and 0.9 for mutations that decreased expression (Supplementary Data 7,8), showing that mutations, in general, are associated with greater average fold changes than copy number variations, perhaps due to nonsense mediated RNA degradation. Complete tables are available at: (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

HER2 analysis as an example

We then focused on HER2 as an illustrative example. A comparison of relative HER2 (ERBB2) protein levels across tumor types illustrates the potential utility of a pan-cancer proteomic analysis. While the overall HER2 protein:mRNA correlation was 0.53 (P = 5e-177), the correlation was 0.61 (P = 1e-69) in BRCA, where HER2-targeted therapy has been demonstrated to be effective (Spearman’s correlations Fig. 1, Supplementary Data 1). Importantly, phosphoHER2Y1248 protein:mRNA correlation was 0.552 (P = 3e-54) and HER2:phosphoHER2Y1248 protein:protein correlation was 0.67 (P = 4e-98) in breast cancer consistent with ability of RPPA to capture both total and phosphoprotein levels from TCGA samples (n=2,503 for overall and n=674 for BRCA correlations and P-value computations using t-distribution test and adjusted for multiple hypotheses testing using Benjamini Hochberg adjustment. n=2,479 in Fig. 1). Based on correlations with DNA, RNA and protein levels in HER2-positive breast cancers, HER2 protein levels were defined as elevated if the relative HER2 level was ≥1.46 (see Methods) (Fig. 1b-d). We also set a cutoff at the relative protein level of 1.00 (which is roughly equivalent to 3+ staining on clinical immunohistochemistry analysis of the breast cancer samples and represent the top 12% of patient samples, see Methods). Using either cutoff, 10–15% of breast cancers demonstrated elevated HER2 by DNA copy number, RNA and protein consistent with clinical data12, 13 (Fig. 1b). Based on those cutoffs, approximately 25% of serous endometrial cancers had coordinated elevation of HER2 DNA, RNA, and protein levels, an even higher frequency than breast cancer. BLCA, colorectal cancer and LUAD demonstrated a higher frequency of elevated protein levels than predicted by mRNA and DNA levels. In an independent cohort of 26 LUAD cell lines using the same cutoffs, 7 of the cell lines had high HER2 protein levels, whereas only 2 cell lines had high mRNA levels, consistent with our observation of elevated protein levels occurring at a higher frequency than elevated RNA levels (Supplementary Table 1, Supplementary Fig. 2)14.

Discordance between HER2 DNA copy number and protein levels has been observed in multiple individual tumors types previously15, 16, 17, 18, 19, 20. Besides diversity in methodology, a number of cancer specific hypotheses, including post-translational regulation of HER2 expression, cytoplasmic HER2 localization16, intratumoral heterogeneity of HER2 amplification19 or polysomy 1717, 20 have been suggested. This clearly contrasts breast cancer, where HER2 levels are usually highly correlated at the DNA, RNA and protein level21, 22, 23, 24. With the advent of TDM1 toxin conjugate therapy (trastuzumab emtansine)25, 26, the higher frequency of elevated HER2 protein levels in BLCA, LUAD, endometrial, and colorectal cancers supports the (pre)clinical exploration of TDM1, which binds HER2 to deliver a potent cell-cycle toxin (a mechanism of activity independent from trastuzumab, a drug with limited activity in endometrial cancer in previous studies27) in these tumor lineages.

Unsupervised clustering analysis

Unsupervised clustering identified eight robust clusters (Clusters A-H, Fig. 2a) when batch effects were mitigated by RBN. Not surprisingly, RBN cluster membership is defined primarily by tumor type with the exception of cluster_E and cluster_F, which include multiple diseases (Fig. 2b). Bladder cancer, however, did not generate a dominant cluster but, rather, was co-located with other tumor lineages in multiple clusters. To identify potential discriminators of clusters, we compared the ability of proteins, RNAs, miRNAs and mutations for each cluster to different samples from those in all other clusters (top 25 discriminators, Supplementary Tables 2-5, all the discriminators at http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA). Supplementary Table 2 highlights the contribution of individual proteins in driving the different clusters. Associations of specific mutations and copy number changes with the clusters were primarily based on known associations of mutations and copy number changes with tumor lineage.4, 5, 6, 7, 8, 9, 10

Cluster_E includes 70% of basal-like breast cancers, the majority of HER2 positive breast cancers (87%) and the largest group of bladder cancers (35%), including many with amplified HER2 (Fig. 2a,b). Cluster_E is defined by TP53 mutations, elevated HER2, cyclinB1 and Rab25 protein levels and low ER and PR levels (Supplementary Table 2). Cluster_F includes smoking-related, upper aerodigestive tract cancers (HNSC, LUAD, and LUSC) and subsets of other tumor types. Cluster_F contains the majority of a “squamous cancer” subset (94%), P<0.0001, Chi-squared test), recently identified through other Pan-Cancer subtype analyses (Hoadley K, personal communication). However, cluster_F also contains an equally large number of non-squamous tumors, predominantly LUAD (58% of the non-squamous tumors in cluster_F). Membership in cluster_F is associated with TP53 mutations and elevated total and phosphorylated EGFR (EGFRp1068 and EGFRp1173), phosphorylated SRC (SRCpY527) and low ER and PR levels. Although TP53 mutations are usually associated with copy number changes and a limited number of recurrent mutations in cancer genes7, cluster_F is unexpectedly enriched in recurrent cancer gene mutations (Supplementary Table 6). Within the group of current smokers in cluster_F (Supplementary Fig. 3), tumors with TP53 mutations show significantly higher rates of co-mutations in the top-25 driver mutations (Methods, P<0.0001, t-test, n=162).

Hormonally responsive ‘women’s cancers’ (luminal BRCA, OVCA, UCEC) form a major tumor super cluster. Basal-like breast cancers and HER2-positive breast cancers are distinct from luminal breast cancers, being located in cluster_E (the majority of HER2 (87%) and basal-like (70%)) and cluster_F (subset of basal-like (25%)). This is consistent with previous data suggesting that HER2 and basal-like breast cancer are distinct from luminal breast cancer5. In light of the recent identification of a “reactive” breast cancer subtype5, we split the luminal cluster into two (reactive breast cluster_A1 and non-reactive ER-positive breast cluster_A2).

For some tumor lineages, localization to different clusters reflects differences in prognosis. Breast cancers located in different clusters demonstrate distinct outcomes: tumors in cluster_E and cluster_F are associated with the worst outcome, probably due to the inclusion of HER2-positive and basal-like tumors. Reactive cluster_A1 shows a better outcome than cluster_A2 (Fig. 2c). The poor outcome associated with KIRC in cluster_F (Fig. 2d) may be due to the absence of VHL mutations (Fisher's exact test (FE), P = 0.008, n=454), which has been associated with a worse outcome in kidney cancer28. Bladder cancers in cluster_B show worse survival compared to all other BLCA, which may be due to associations with TP53 mutation (FE, P<0.001) and cMYC amplification (FE, P = 0.042) (n=127) (Fig. 2e).

We evaluated the concordance between RBN protein clusters and mRNA clusters derived from the same sample set (Supplementary Table 7). Most of the protein clusters predominantly corresponded to a single respective mRNA cluster despite the mRNA clusters being defined with a pool of about 20,000 mRNAs, whereas only 181 proteins and phosphoproteins were used to generate the protein clusters. Therefore, many of the features defining the mRNA clusters were captured by just a few proteins. This agreement between RNA and protein based clustering provides validation of the quality of the protein data, as well as the selection of protein targets in the arrays. However, clusters E and F were noticeably different from their mRNA counterparts. Unlike protein cluster_E that contains BLCA and BRCA, bladder cancer formed a separate cluster in mRNA data, distinct from HER2 and basal-like breast cancers. LUAD also formed a separate mRNA cluster, distinct from the LUSC/HNSC mRNA cluster, unlike protein cluster_F that contains LUAD as well as LUSC and HNSC.

Reduction of tissue-specific proteomic signatures

Tumor lineage represents the dominant determinant of protein clustering using the RBN approach (Fig. 2). We, therefore, investigated whether further transforming the RBN data to reduce tissue signatures by median centering within tissue types (MC, see Methods) would identify clinically or biologically relevant protein patterns that span multiple tumor lineages (Fig. 3a). Using MC, we obtained 7 clusters (I-VII) that were no longer strongly correlated with tumor lineage, as evident from the top annotation bar in Fig. 3a (Supplementary Fig. 4), and from the tissue vs. cluster cross-tabulation (Fig. 3b). This allowed exploration of molecular events that spanned multiple tissues, which was not possible with the RBN approach. Supplementary Table 8 shows a contingency table the distribution of samples across RBN vs. MC clusters, highlighting the differences between the clusters. Supplementary Tables 9-12 show the top 25 proteins, mRNAs, miRNAs, and mutations that discriminated different MC clusters (full table available at http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

Cluster_I was primarily driven by phosphoPEA15, YB1, EEF2 and ETS1 proteins (Supplementary Table 9), which were markedly elevated in a subset of colorectal tumors (18%). Cluster_I exhibited enrichment of APC and KRAS mutations, very few HER2 amplifications, but moderately high HER2 protein levels (Fig. 3a, Supplementary Tables 9,12). It also had evidence for suppressed DNA damage response, apoptosis, and mTOR and MAPK pathway levels (Fig. 4b). Cluster_II was divided into two further sub-clusters, one primarily driven by HER2 (IIa) and one by EGFR (IIb) (Supplementary Table 9). Interestingly, a subset of OVCA, UCEC, BLCA and LUAD samples that had HER2 amplification and HER2 protein levels comparable to breast HER2+ samples were located in cluster_IIa, raising intriguing opportunities for (pre)clinical investigation of HER2 targeted therapy and particularly TDM1 therapy as noted above. Cluster_IIa also had activated RTK and cell cycle pathways, but suppressed hormonal signaling pathways (Fig. 4b). Similarly, a subset of HNSC and lung samples that had EGFR levels comparable to a subset of GBM samples (28%) were located in cluster_IIb, warranting exploration of potential benefit from EGFR pathway-targeted drugs29. Tumors in cluster_IIb were enriched in EGFR mutations, contained few PTEN mutations, and had elevated RTK pathway and suppressed mTOR pathway signatures. Clusters III-VII consisted of a mixture of all tissue types. Cluster_V was the most distinctive, exhibiting a strong “reactive” signature5, with elevated MYH11, RICTOR, Caveolin1, and Collagen VI, and an activated EMT signature. Cluster_V also exhibited low cell cycle, Wnt-signaling and DNA damage response pathway signatures. Cluster_V contained the majority of the breast reactive samples along with multiple other tumors with a “reactive” signature consistent with the reactive phenotype being a pan-cancer characteristic. Cluster_III was the antithesis of “reactive” cluster_V and was primarily driven by elevated BRAF, ER-alpha and E-cadherin (Fig. 3b). In contrast to cluster_V, cluster_III had low EMT, apoptosis, and MAPK pathway signatures, but high DNA damage and hormonal pathway signatures. Patients in cluster_III may potentially benefit from (pre)clinical hormone targeting therapies. Cluster_III also had high beta-catenin levels, suggesting activation of the canonical Wnt-signaling pathway. Cluster_IV also had high beta-catenin, as well as activated AKT, MAPK and mTOR pathways, but suppressed DNA damage, apoptosis, EMT and cell cycle pathways. Cluster_IV and cluster_VII were antitheses. The high levels of phosphoAKT and phosphoMAPK in cluster_IV, suggested evaluation of (pre)clinical benefit from kinase-targeted therapies. Cluster_VI showed high EMT, cell cycle, apoptosis, mTOR and MAPK pathway signatures, also suggesting further evaluation of kinase-targeted therapies. Cluster_VI had low beta-catenin, consistent with suppressed Wnt-signaling. Cluster_VII also showed low beta-catenin, with suppressed AKT, MAPK, mTOR and RTK pathways.

Interestingly, clinical outcomes correlated with MC cluster membership, indicating the power to identify important tissue-independent processes. COAD in cluster_V had better outcome compared to COAD located in other clusters (Fig. 3g) (n=334), which may, in part, be due to depletion of mutations in TP53 (6% vs. 15%, Fisher’s Exact (FE) P = 0.05), APC (14% vs. 25%, FE P = 0.044) and KRAS (5% vs. 16%, FE P = 0.013), consistent with previous literature showing these are associated with a worse outcome30, 31, 32. The poor outcome for KIRC in cluster_VII may be partly due to enrichment of TP53 mutations (6% vs. 0.8%, FE P = 0.005, n=454) (Fig. 3c). In contrast, KIRC in cluster_IV are associated with better prognosis (Fig. 3e). For OVCA, membership in cluster_VII is associated with improved survival (Fig. 3d). LUSC in cluster_V appear to have worse prognosis, which may be related to elevated EMT pathway activity compared to LUSC in other clusters (Supplementary Fig. 5)33, 34, as well as low E-cadherin protein levels (Fig. 3f). Thus, reduction of tissue-specific signatures reveals a number of processes that transcend tissue boundaries and may represent cross-tissue biological, prognostic, and therapeutic opportunities.

Analysis of pathways and targets

To capitalize on the RPPA data, we developed a series of pathway predictors (see Methods), based on member proteins selected by literature review (Supplementary Table 13). TSC/mTOR signaling, which integrates information from the PI3K/Akt, Ras/MAPK and LKB1/AMPK pathways35, was treated as a separate pathway, as was the hormone_a (ER, pER and PR) and a series of downstream components of the hormone signaling pathway (hormone_b36, 37, 38). All proteins and genomic events with a Spearman’s ρ>0.3 or ρ<−0.3 for association with the pathway score are also presented (See methods, Fig. 4, Supplementary Figs.6-9, Supplementary Table 13) providing additional information on potential pathway membership.

In general in the RBN analysis, pathway scores were associated with tumor lineage (Fig. 4a, Supplementary Fig. 10). In Figures 4a,b, each cell in the heatmap represents the mean pathway score for that cluster or tumor lineage. Blue represents a suppressed pathway, red means an activated pathway, and white representing a score that does not differ across the set (see Methods). As expected, individual RBN clusters (Fig. 4a) show similar pathway scores to their dominant constituent tumor lineages, e.g. GBM is similar to cluster_H, KIRC is similar to cluster_G, etc. However, as clusters E and F do not consist of a single predominant lineage, their pathway score pattern is not concordant with any one tumor lineage. Similarly, the MC heatmap (Fig. 4b) shows that MC clusters, in which tissue specific effects are removed, do not reflect a single tumor type. This emergent phenotype illustrates the mitigation of tissue-specific signatures by MC, and the emergence of new, pan-cancer patterns that span multiple tumor types. In Supplementary Fig. 10, the data is transformed so that the color spectrum in the heatmaps represents absolute values of pathway scores (where only score magnitude is considered) and thus reflects ‘distance from the global pathway mean’, rather than relative protein level (see Methods). This emphasizes that both low (e.g. inhibitors) and high protein levels can be markers of pathway activity. Thus in Supplementary Fig. 10, UCEC and HNSC have a near identical hormone_a score, caused by a high (UCEC) and low (HNSC) protein score respectively. The pathway-based analyses benefit hugely from the large dataset providing sufficient power to identify associations that could otherwise not be robustly identified.

Focusing on individual pathway analysis (Fig. 4c-f, Supplementary Figs. 6-9), the high degree of correlation between pathway members, including phosphoproteins, supports the ability of RPPA to capture high-quality information including phosphoprotein levels from TCGA samples. Unexpectedly, the proteins driving the pathway signatures varied across individual tumors and tumor lineages, as did the associated proteins and genomic aberrations (Fig. 4, Supplementary Figs. 6,8). This suggests that intrinsic gene expression patterns or mutational patterns provide important contributions to convergent functional pathway output. The EMT signature, which may also represent reactive stroma, showed the greatest variation, being markedly elevated in GBM and reactive BRCA tumors (Fig. 4c,e). Significant variation in EMT was also observed within disease type and RBN clusters. For example, Cluster_F (HNSC, LUAD, LUSC) showed a separation into distinct epithelial and mesenchymal groups based on the EMT score and related protein EMT markers. RTK and downstream signaling signatures were elevated in GBM, likely due to EGFR amplification and activation of downstream signaling events (Fig. 2). Endometrial, ovarian and most breast cancers demonstrated a high hormone_a signature (Fig. 4d,f). However, an elevated hormone_b signature, indicative of functional downstream activation, was restricted to luminal, reactive, and HER2-positive breast cancers (Supplementary Fig. 11) suggesting differential “wiring” of hormonal signaling across tumor lineages. HER2-positive breast cancers, whether ER-positive or -negative, demonstrated elevated levels of GATA3, INPP4B, and AR (hormone_b signature) suggestive of active downstream hormonal signaling despite low levels of ER, pER and PR in many of the HER2-positive tumors (Fig. 2, Supplementary Fig. 11). A subset of endometrial cancers had massively elevated pAkt levels, likely due to the high frequency of coordinated genomic aberrations in the PI3K pathway, in particular, the loss/mutation of PTEN10, 39 which is consistent with responsiveness of endometrial cancers to PI3K pathway inhibitors40, 41.

We analyzed a number of potentially actionable proteins (n=25, Fig. 5a,b), selected based on a literature review (Supplementary Methods) for associations with proteomic and genomic events as well as for potential ability of proteomics to identify patients likely to benefit from targeted therapies. Luminal breast cancers (including AR-positive triple-negative breast cancers which cluster with luminal breast cancers) demonstrated selective elevation of AR, BCL2, FASN and pACC, suggesting these molecules or their associated pathways as potential therapeutic targets. The elevation of HER3 in KIRC may represent a therapeutic opportunity. SRC is activated in all but the hormone-responsive and bladder cancers, offering another potential therapeutic opportunity. EGFR activity, in general, parallels SRC activity, but in GBM is associated with NOTCH1 and HER3 activation, suggesting an interesting opportunity for exploration of combination therapy in (pre)clinical studies. PhosphoSRC, which is a downstream target of EGFR, was highly expressed in a subset of HNSC tumors, suggesting that these may be more sensitive to EGFR targeting strategies. As noted above, HER2 levels are elevated in a subset of UCEC, BLCA, BRCA and colorectal cancers and may represent responsiveness to HER2 targeted therapy. MYC, which may become targetable by emerging therapeutic approaches42, is selectively amplified and expressed in high-grade serous ovarian cancer and may represent an important target in this disease that currently lacks targeted opportunities7.

To determine whether protein levels, including phosphoproteins, can predict patient outcome, we determined correlations with overall survival (see TCPA11) for a comprehensive analysis) using Cox Proportional Hazards (CoxPH) models. In the complete Pan-Cancer dataset, 80 (including 24 phosphoproteins) of the 181 proteins demonstrated a significant (corrected for multiple comparisons) correlation with outcome. Importantly, 57 proteins, including 19 phosphoproteins, showed a multiple comparisons corrected correlation with outcome in KIRC. However, with the exception of breast cancer (13 candidates), this approach showed five or fewer proteins correlating with outcome in other tumor lineages. Why kidney cancer shows such strong correlations is not completely understood, but may reflect the maturity of the outcome data in this dataset. For some of the other diseases included in the Pan-Cancer dataset, the associated outcome data are immature, for example, the low number of events in the BRCA and endometrial cancer datasets limits the ability to detect the prognostic and predictive value of protein markers.

To extend the single protein analysis available in TCPA, we performed a formal training/test set analysis of pathways and potentially actionable proteins. As indicated in Supplementary Table 14; 17 predictors (4 pathways, 9 total proteins and 4 phosphoproteins) passed a rigorous training/test set approach and showed a robust correlation with outcome in at least one disease. As expected from the analysis of single proteins, most surviving correlations were in kidney cancer. Several pathway predictors that survived the training/test set approach demonstrated marked associations with patient outcome in the overall sample sets (Supplementary Fig. 12). PhosphoSRC (SRCpY416) and the transferrin receptor (TFRC) showed an association in three diseases suggesting particular importance for outcome. However, the effects of the TFRC on patient outcomes were different across diseases suggesting an interaction with lineage-specific events. TFRC expression was associated with a significantly worse prognosis in LUAD and KIRC. These findings have potential implications for clinical targeting using TFRC for targeted delivery of chemotherapy or other agents43. Comparing the performance of the optimized cutoff approach with medians, quartiles or tertiles, more often applied in literature, we note that up to 50% of the predictions from the optimized cutoff approach were confirmed using these alternative cutoffs. However the optimized cutoff approach, combined with a rigorous training and test set evaluation, performed better in 17 out of 21 (81%) cases (as indicated by lower P-values) compared to the use of median, tertiles or quartiles.

Network visualization

Based on the availability of protein data across a large number of samples, we used a probabilistic graphical models approach44, 45 without the inclusion of prior knowledge to create an unbiased signaling network (Fig. 6, see Methods). We used the relatively large number of samples per tumor lineage to elucidate links in specific cancers and across multiple cancers, inferring networks using tumor lineage-specific samples. Interplay between nodes was quantified using scores from the graphical model analysis (see Methods) that identify links between nodes whilst controlling for the effects of all other observed nodes. Several expected links were observed across most tumor types, including pMEK with pERK, beta-catenin with E-cadherin and pPKCdelta with pPKCalpha and pPKCbeta, supporting the ability of RPPA analysis to yield high-quality signaling information from TCGA samples. Other expected links were seen in only a subset of tumors such as pAKT with pPRAS40 and pTSC2 (tuberinPT1462), consistent with differential wiring of signaling pathways in different cancers. A number of other links such as MYH11 with Rictor, cyclinB1 with FOXM1, and pACC with FASN were not expected and warrant further exploration. The interplay between p85 and PTEN is consistent with our demonstration that p85 is a key determinant of PTEN stability39, 46. The negative link between pAKT and PTEN was expected, but the one between p85 and claudin7 in LUSC was not and may be worthy of further exploration. PI3K/AKT signaling does not link clearly to mTOR, which appears to primarily be downstream of MAPK signaling47, 48, 49. The relatively weak links in the PIK3K/AKT pathway are striking given the degree of antibody representation for this pathway in the RPPA analysis. Key nodes such as CDK1 unexpectedly linked a wide range of protein pathways. Overall, the data suggest that the EGFR receptor family, together with the linked MEK and MAPK pathways, is the dominant determinant of signaling across the cancer lineages in the Pan-Cancer analysis. Using independent datasets in breast cancer, ovarian cancer and endometrial cancer, as well as published research, many of the strongest protein links in the network could be validated (Supplementary Fig. 13 and Supplementary Table 15), supporting the notion that large RPPA-based protein datasets can be used to “learn” networks in an unbiased manner.

Discussion

Cellular biology is effectuated in considerable part by proteins, and, unfortunately neither DNA copy number nor mRNA expression are able to faithfully predict protein level and in particular the post-translational modifications of proteins that are necessary for function (Fig. 1, Supplementary Fig. 1)1, 2, 3, 50, 51. Hence, evaluation of the functional proteome offers the ability to complement genomic and transcriptomic analysis in projects like the TCGA for identification of biomarkers and elucidation of underlying biological mechanisms both within and across diseases. The availability of high-quality proteomic data across large numbers of samples makes the case more compelling. In sum, a proteomic view of TCGA data yields insights that cannot be acquired through analysis driven solely by genomics or transcriptomics. The high degree of correlation between proteins, including phosphoproteins, in signaling pathways (Figs. 2,4,6, TCPA11) supports the applicability of RPPA analysis to TCGA samples. Further, the ability to construct de novo signaling networks (Fig. 6) that capture many known relationships supports the contention that proteomic data derived from the RPPA analysis of TCGA samples can be used to inform systems-level analyses of signaling pathways and networks. Full integrative analysis of the DNA, RNA and protein relationships embodied in the several TCGA datasets will require additional analysis, but a number of interesting observations are immediately apparent.

Analysis of this large dataset demonstrates that, in general, tumor type and subtype are the dominant determinants of protein levels. This observation highlights the risk inherent in disease-specific studies that commonalities, differences, and themes that emerge across tumor types will remain undiscovered. We therefore implemented a computational approach, MC, to decrease the dominant effect of tissue-specific protein expression. This approach allowed for the discovery of processes that drive cellular behavior across tumor types and made it possible to identify tumor characteristics that warrant exploration as therapeutic opportunities. The analysis of individual therapeutically relevant proteins (e.g., HER2, Figs. 1,5) and pathways (Fig. 4) permitted classification of patient samples based on pathway activity and therapeutic tractability across different tumor types. The ability of the Pan-Cancer analysis to identify the discordance between HER2 CNV, mRNA expression and protein expression in colorectal and serous endometrial cancers (Fig. 1) argues that a broad protein-based analysis of patient samples across multiple diseases can highlight potential therapeutic opportunities not obvious from studies within single diseases or driven by RNA and DNA analysis alone.

The pathway analysis (Fig. 4, Supplementary Figs. 6-9) identifies multiple protein changes that are associated with the same functional outcome (i.e., pathway activation) in different samples and tumor types (Fig. 4). A number of proteins and genomic events correlate with pathway scores, developed using proteins defined by literature review (Fig. 4). Although some of those relationships could be identified by including members of upstream or downstream signaling or interacting pathways, many of the associations would not be predicted a priori, demonstrating that these approaches offer the potential for discovery of novel pathway connections. The ability to identify unexpected correlations was particularly clear in the network analysis (Fig. 6). For example, the strong links between MYH11 and RICTOR and between ETS1 and pPEA15 across tumor types offer opportunities for discovering new functional relationships. Some associations we reported, such as that of the mTOR pathway with MEK and MAPK, while supported by the literature47, 48, 49, 52 do not currently receive adequate consideration. Although molecular pathways often seem “set in stone”, the identification of unbiased signaling networks using large datasets can provide a powerful tool to identify tissue-specific networks, as well as to demonstrate the importance of “non-canonical” interplay, allowing for re-conceptualization of networks and the role they play in specific diseases.

A major goal of the molecular characterization of tumors, is the identification of tumor subsets and specific aberrations that can be used in the clinic as biomarkers and/or for targeted therapy (either single-agent or in combination). A bird’s eye view of the functional proteome of large sample sets encompassing multiple tumor lineages may help to suggest potential unexpected targets that are applicable to disease subsets or across diseases. The ability to identify many biomarkers associated with patient outcome (TCPA) and the ability of a set of biomarkers to pass a rigorous training/test set approach (Supplementary Table 14) suggest that additional Pan-Cancer analyses, as well as mechanistic analyses of the current proteomics study will improve our ability to understand tumorigenesis and identify new markers and targets.

Methods

Description of the protein data

Proteomic data were generated by RPPA across 3,467 patient tumors obtained from TCGA, including 747 breast (BRCA), 464 colon and rectal adenocarcinoma (COAD and READ), 454 renal clear cell carcinoma (KIRC), 412 high-grade serous ovarian cystadenocarcinoma (OVCA), 404 uterine corpus endometrial carcinoma (UCEC), 237 lung adenocarcinoma (LUAD), 212 head and neck squamous cell carcinoma (HNSC), 195 lung squamous cell carcinoma (LUSC), 127 bladder urothelial carcinoma (BLCA) and 215 glioblastoma multiforme (GBM). Those were all the samples we could obtain from TCGA and no samples were excluded. The result is, to our knowledge, the largest and most diversified database of tissue protein levels yet available, an unparalleled basis for rich functional analysis.

RPPA methodology has been described in 4, 5, 6, 7, 8, 9, 10 and is also provided in the Supplementary Methods. In total 181 high-quality antibodies targeting total (n=128), cleaved (n=1), acetylated (n=1) and phosphoproteins (n=51) were used (detailed in Supplementary Data 9). In the RPPA assay, antibodies to phosphoHER2 and phosphoEGFR have been noticed to cross-react, especially when the opposite molecule is present at very high levels. This mainly concerns EGFRpY1068 (but not EGFRpY1173), which cross-reacts with overexpressed HER2pY1248. Taking into account their favorable signal:noise ratio (10:1), useful information is contributed by both if expressed differentially, and they are thus both included. The antibodies encompass major functional and signaling pathways of relevance to human cancer. Pathways included are proliferation, DNA damage, polarity, vesicle function, EMT, invasiveness, hormone signaling, apoptosis, immunological, stromal, TGFα/β, transmembrane receptors, metabolism, LKB1/AMPK, TSC/mTOR, PI3K/Akt, Ras/MAPK, Hippo, Notch, and Wnt/beta-catenin signaling (Supplementary Fig. 14 and Fig. 6) with minimal redundant information (Supplementary Fig. 15). Supplementary Fig. 16 shows a representative image of a typical antibody slide.

The numbers of patient samples and antibodies are greater than those presented in previous TCGA marker papers4, 5, 6, 7, 8, 9, 10 based on the availability of additional samples as well as validation of additional antibodies. The detailed TCGA datasets are available online (https://tcga-data.nci.nih.gov/tcga) and combined with a number of visualization and analytic tools from TCPA (http://app1.bioinformatics.mdanderson.org/tcpa/_design/basic/index.html). High-resolution images of all heatmaps and the network are available online (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA). Some key clinical variables are shown in Supplementary Tables 16-17; extensive clinical information for all lineages is available online (https://tcga-data.nci.nih.gov/tcga) and available in the various TCGA marker papers4, 5, 6, 7, 8, 9, 10.

Protein correlations

To match the 181 antibodies available, 162 unique mRNAs were selected from downloaded RNASeqV2 data (https://tcga-data.nci.nih.gov/tcga), resulting in 184 matched and 24,282 random protein:mRNA pairs. Spearman’s rank correlations were computed on both the random and matched pairs, with associated P-values (Supplementary Data 1,2 and at http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA). The ρ values of the matched pairs were plotted in histogram form; the ρ values of the random pairs are represented as a background curve (Fig. 1a, Supplementary Fig. 1). Student’s t-test was used to compare the ρ values of all matched pairs with the ρ values of all random pairs (mean matched pairs: 0.3; mean random pairs: −0.006) and showed a significant difference (P<2.2e-16). For miRNA vs. protein, because the number of miRNAs and proteins were small, we computed all pair-wise Spearman’s rank correlations with t-test P-values (Supplementary Data 3). For the CNV vs. protein expression analysis, we divided the samples into groups of amplified vs. copy number neutral and deleted vs. copy number neutral, and computed the mean fold changes in protein expression. Similarly, to compare mutation vs. protein, we divided the samples into mutated vs. wildtype and computed the fold changes in protein expression. We then used t-tests to evaluate statistical significance of the fold changes (Supplementary Data 5-8 and at ttp://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA

Furthermore, we computed all pair-wise protein:protein correlations using the entire Pan-Cancer dataset, in total 16,290 correlations (Supplementary Data 4 and at http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA). The top-10% had a Spearman’s rank correlation coefficient magnitude of 0.3 or higher (Bonferroni-adjusted P≤3.67e-67). Consequently, we considered a correlation magnitude of 0.3 or higher (sign independent) as a reasonable cutoff threshold for the analysis presented in the pathway sections of this study (Fig. 4, Supplementary Figs. 6-10, Supplementary Table 13).

Discriminator selection

To detect the discriminating biomarkers for each cluster (obtained by hierarchical clustering using the RPPA data normalized by either RBN or MC), LIMMA53 was used for the continuous data (protein, mRNA, miRNA) by comparing samples in each cluster with samples in all the other clusters together; information gain54 was used to select the categorical discriminators (mutation). The resulting data were sorted by decreasing order of the log-odds for the former and by decreasing information gain for the latter method. The top-25 most significant discriminators are shown in Supplementary Tables 2-5 (RBN) and Supplementary Tables 9-12(MC). The complete overview of protein, mRNA, miRNA and mutation discriminators can be accessed online (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

BRCA and UCEC subdivision

For BRCA subtypes, first the reactive subtypes were classified according to the method described in the TCGA marker paper5. The other subtypes were then classified based on PAM50. For UCEC subtype classification, serous samples were first selected based on the integrative cluster (serous-like) reported in the TCGA marker paper10. Clinical histopathological subtype (https://tcga-data.nci.nih.gov/tcga) was used in any remaining cases.

HER2 cutoffs

Normal tissues for the lineages studied in this paper have been reported to have low or medium HER2 levels (http://www.proteinatlas.org/ENSG00000141736/tissue/staining+overview)55. To identify the threshold of HER2 mRNA and protein expression in breast cancer that could classify tumors as HER2-positive, we obtained PAM50 classifications for all the TCGA breast cancer samples, and divided them into two groups; HER2-positive and non-HER2-positive samples. We used the conjunctive rule algorithm in Weka software56 to determine the best HER2 total protein cutoff that separated the HER2-positive from the non- HER2-positive samples based on HER2 (ERBB2) copy number. The best protein threshold was found to be 1.46, which yielded 93% accuracy of prediction and a receiver operator characteristic (ROC) area under the curve (AUC) of 0.81. We did a similar analysis using HER2 mRNA and found a best cutoff of 14.26 (in log2 frame), which yielded an accuracy of 93% and a ROC AUC of 0.82. In addition to trastuzumab, other drugs targeting HER2 have entered clinical trials, such as TDM1, for which HER2 expression on the cell surface is sufficient to achieve preferential binding to the cell and therapeutic impact. Since data for TDM1 response is not readily available, a threshold of HER2 expression that may be sufficient to expect a response could not be calculated. We therefore compared samples in which HER2 was amplified vs. not amplified, aiming to find a threshold that might be reasonable to test. Using the dot plots, a protein threshold of ≥1.00 was adopted; roughly equivalent to 3+ on immunohistochemistry of clinical samples in breast cancer. The crosstab in Figure 1b gives the breakdown of percentage of samples above these thresholds for each tumor. If a tumor lineage had more than 5% HER2-positive samples according to any of the cutoffs, this is indicated in red.

General heatmap section

A two-way unsupervised hierarchical clustering analysis was used to discover the groups of biological objects sharing common characteristics57, 58, and a two-dimensional heat map was drawn to visualize protein expression patterns. We used Ward linkage as the agglomeration rule and 1-Pearson correlation as the dissimilarity metric. Based on protein expression patterns and guided by the clustering dendrogram, we divided the RBN data set into 8 clusters and the MC dendrogram into 7 clusters. As seen in the RBN heatmap, most clusters represented one major disease. Exceptions were clusters_E and _F. Based on the recent TCGA marker paper5, the hormone-responsive breast cancer cluster (cluster_A) in the RBN dendrogram was further divided into 2 subclusters, A1 (reactive breast cancers) and A2 (remaining luminal breast cancers). Based on marked enrichment with clinically relevant proteins, cluster_II in the MC dendrogram was further divided into two subclusters, cluster_IIa (HER2 elevated) and cluster_IIb (EGFR elevated). Hierarchical clustering analysis was performed using R, version 2.15.1 (http://www.r-project.org/). Heatmaps were generated using an NGCHM R-package59. Annotation bars were added to the heatmap that included tumor lineage, purity and ploidy; stromal and immune scores; BLCA subtype and PAM50 classification (BRCA). Significantly mutated genes (present in more than 5% of tumors in the dataset, resulting in 16 genes) are included as are the two most frequently observed amplifications. Statistical significance for the annotation bars on top of the various heatmaps was calculated by Chi-squared test (tumor lineage, mutations and amplifications), ANOVA’s F test (purity, ploidy stromal and immune score), and Fisher’s exact test (PAM50 and BLCA subtypes). Data are missing for BLCA subtype (15/127), BRCA subtype (52/747) and HER2 and MYC amplification (64/3,467).

Batch effects removal

The 3,467 RPPA Pan-Cancer samples were run in 6 batches in total, resulting in potential batch effects upon merging the sets. Batch effects in RPPA data are a known concern, even when controlling for critical materials such as the treated glass slides, antibodies, enzymes and suppliers60. A new algorithm, replicates-based normalization (RBN), was therefore developed, using replicate samples run across multiple batches to adjust the data for batch effects. The underlying hypothesis is that any observed variation between replicates in different batches is primarily due to linear batch effects plus a component due to random noise. Given a sufficiently large number of replicates, the random noise is expected to cancel out (mean=zero by definition). Remaining differences are treated as systematic batch effects. We can compute those effects for each antibody and subtract them out. In one batch, many samples with duplicates in the other 5 batches were run, and could therefore serve as anchor for all batches. The number of duplicate samples with each batch varied between 71 and 207. This batch was designated “anchor” batch and was used unchanged. We then computed the means and standard deviations of the common samples in the anchor batch and each of the other batches. The difference between the means of each antibody in the two batches and the ratio of the standard deviations provided an estimate of the systematic effects between the batches for that antibody (both location-wise and scale-wise). Each data point in the non-anchor batch was adjusted by subtracting the difference in means and multiplying by the inverse ratio of the standard deviations to cancel out those systematic differences. Whether RBN could successfully integrate batches while preserving known biological variation, was tested on TCGA breast cancer samples. As breast cancer subtypes (luminal, HER2-positive and basal-like) are well established13, we expected the subtypes from different batches to cluster together. Without RBN, the batches clustered by batch. After RBN, the batches clustered by subtypes spanning multiple batches (Fig. 2). Details of these experiments have been published previously61.

Reducing tissue differences to cluster across tumors

Using RBN, batches of RPPA data could be merged successfully. However, as protein levels of different tumors are (usually) quite distinct from each other, most samples clustered by tumor lineage (Fig. 2). Normal cells differentiate into different tissues by turning on or off different sets of genes. When cells become malignant, they retain many tissue-specific expression characteristics. We hypothesized that tissue-specific effects exist because of those expression differences and equalizing the median expression of genes across tumors might reduce those effects. A gene that is turned off in all the samples of a tumor lineage will have little variation in expression, similar to a gene that is always turned on, which will also have little variation, but an overall high level. To compare across tumor lineages, we started with the batch corrected RBN data and took sets of all samples belonging to each tumor lineage. We subtracted the median protein expression across all the samples from a single lineage (median centering, MC), making the median expression of all proteins within any given tumor equal to zero. That removed the fixed, bias component from that tissue lineage but retained the variable component found in each tumor. Since the tissue specific component had been removed, we could then compare the variable component (which was relative in scale) in each tumor sample across different tissues. That allowed for the comparison between samples with high/low expression in one tumor and samples with high/low expression in another tumor, such as HER2 or EGFR expression. Basal-like breast cancer was treated as a separate tumor lineage from the other breast cancer samples due to its expression profile being so different that it did not merge with any other tumor or even other breast samples during RBN clustering.

Tumor purity and ploidy

We obtained tumor purity and ploidy data based on the ABSOLUTE algorithm62 from TCGA Pan-Cancer working group. We calculated stromal and immune scores based on the ESTIMATE algorithm using the TCGA Pan-Cancer gene expression dataset (syn1695373, https://www.synapse.org/#!Synapse:syn169537363).

Pathway analysis

For each pathway, members, illustrated in Supplementary Table 13, were predefined based on a Pubmed literature search on review articles describing the various pathways in detail. RBN RPPA data were median-centered and normalized by standard deviation across all samples for each component to obtain the relative protein level. The pathway score is then the sum of the relative protein level of all positive regulatory components minus that of negative regulatory components in a particular pathway. We averaged antibodies targeting different phosphorylated forms of the same protein with ρ >0.85 (Pearson’s correlation). The pathway scores are visualized in the bar just above the heatmap and as a dotplot below the heatmap (median and inter-quartile range indicated, Fig. 4c-f, Supplementary Figs. 6-9). Subsequently, for each version of the pathway scores (RBN or MC derived), a Spearman’s rank correlation test was performed between each pathway score and every protein. If the ρ was >0.3 or <−0.3, the protein was included in the heatmap. Regardless of the ρ, pathway members for the given pathways were included. Annotation bars (from Fig. 2) were included if they were statistically significantly associated with the pathway (P<0.05, Kruskal-Wallis test, n=3,467), corrected for multiple testing64. Tumor lineage and cluster were included to facilitate interpretation. For each pathway heatmap (RBN or MC), the samples were first sorted by the alphabetic order of either cluster or tumor, and then by the increasing order of a pathway score.

Using the heatmap method described above, two additional summary heat maps for the pathway scores (RBN and MC) were generated (Fig. 4a,b) to provide an overall view of the relationships between tumors, unsupervised clusters and pathways. Mean pathway scores were calculated for each tumor as well as cluster variables, and the combined mean pathway scores were standardized for each pathway across all tumor and clusters. In both the individual pathway plots and the heatmap summary plots, hierarchical clustering was based on Pearson’s correlation-based distance matrices65 and Ward linkage. The dynamic heat maps were generated using the R-package NGCHM59. Each cell in the heatmap represents the mean pathway score of all the samples in that cluster or tumor lineage, with blue representing a suppressed pathway, red representing an activated pathway, and white representing neither.

Supplementary Figure 10a,b shows a similar measure of pathway activity, but on the absolute scale. The supplementary figure is derived as follows. First, the RPPA dataset (either RBN or MC normalized) is globally scaled so that the protein expression level measurements have zero mean and unit standard deviation over all samples. Next, for each cluster and tumor lineage, we calculate the mean (scaled) protein expression level for each protein. We then convert these means to their absolute value (as low or high mean protein levels could both be markers of pathway activity), obtaining an absolute mean protein level for each protein in each cluster or tumor lineage. Finally, for each pathway, we calculate the average of the absolute mean levels over the proteins that participate in the pathway. This value is designated the differential pathway activity score, as it indicates the deviation from the mean expression of a pathway in a given cluster or tumor lineage, and can thus be seen as a proxy for pathway activation/deactivation.

Actionable protein analysis

The analysis focused on the potential ability of proteomics to predict response to proteins currently of increased interest, due to proposed targetability or potentiality as a drug target in the drug development stage. The list of proteins is not exhaustive, but rather includes many different processes and pathways with varying importance in different tumor lineages included in this study. In the Supplementary Methods, registered trials targeting many of these proteins are included.

To visualize the expression pattern of these 25 proteins, heatmaps were generated59 using the RBN dataset (Fig. 5). Proteins were ordered by unsupervised hierarchical clustering and samples were ordered by cluster (disease) membership and within each, ordered by unsupervised hierarchical clustering. Ward’s method and 1-Pearson correlation were used as a dissimilarity metric and linkage.

Network analysis

Networks were estimated using statistical models known as a probabilistic graphical models (specifically Gaussian graphical models)45. These models use an undirected graph or network to describe probabilistic relationships between variables. In contrast to pair-wise correlation analysis, the networks are rooted in a global, multi-dimensional approach that identifies links between nodes whilst controlling for the effects of all other observed nodes.

Statistical inference of networks is a so-called ‘high-dimensional’ problem because network descriptions require a large number of parameters relative to available sample sizes (especially at the disease or cluster level). This motivates a need for regularization to learn sparse, parsimonious networks and thereby control over-fitting. We used l1-penalization for this purpose, specifically via an algorithm known as graphical lasso44, as implemented in the R-package huge66. A parameter λ that controls the strength of penalization was set by 10-fold cross-validation in all cases. To prevent artifacts that can arise due to duplicated nodes, related nodes that were relatively highly correlated were merged prior to network analysis. In each such case, only one of the set of correlated nodes was used for network inference and the remaining merged nodes are shown in white. Since protein levels are measured in arbitrary units (depending on affinity and avidity of specific antibodies), for each network the data were standardized before applying the graphical lasso, such that each protein had zero mean and unit variance.

Outcome analysis

A training-test approach was adopted for survival analysis. In each of the 11 tumor lineages, samples with survival data available were randomly divided into training (2/3) and test (1/3) sets with balanced events in both sets. The training set was used to obtain an optimized cutoff, which was “locked” (i.e. used without change) on the test set. Essentially, samples were sorted based on the protein expression of the interesting gene or pathway score. Each possible cutoff in the middle 60% of samples was checked using Cox’s regression model. The cutoff with lowest P-value was chosen as the optimized cutoff. In the test set, samples were divided into high and low groups according to this optimized cutoff by either percentage or absolute value. Then the hazard ratio, Wald’s test P-value and Kaplan-Meier survival curves of the two groups were examined by Cox’s regression analysis. Only the predictors that were successfully validated in the test set are shown in Supplementary Table 14 Kaplan–Meier survival curves were generated to illustrate the survival differences in the four significant pathways using the whole sample set (Supplementary Fig. 12).

Supplementary Material

Figure 1

HER2 RPPA correlations with copy number and mRNA

a Histogram of Spearman’s rank correlation (ρ values) for 206 pairs of proteins and matched mRNAs across all tumor types. The black curve represents the background of ρ values using 28,960 random protein-mRNA pairs in the same dataset.

b Crosstab identifying HER2-positive tumors by copy number, mRNA expression and protein expression across 11 tumor types. Cutoffs are defined in Methods. BRCA and UCEC are subdivided for clinical relevance regarding HER2 protein levels. Total sample numbers with analyses for all three platforms (CNV, mRNA and protein) are indicated in parentheses. Percentages ≥5% are highlighted (red).

c Relationship between HER2 copy number and HER2 protein level by RPPA across all tumor types (n=2,479). The box represents the lower quartile, median and upper quartile, whereas the whiskers represent the most extreme data point within 1.5 × interquartile range from the edge of the box. Each point represents a sample, color-coded by tumor type or subtype. As expected, ERBB2 amplified samples have much higher HER2 protein levels than non-amplified samples.

d Relationship between HER2 mRNA and protein expression across all tumor types (n=2,479). Each protein represents a sample, color-coded by tumor type or subtype. Spearman’s correlation between HER2 protein and mRNA is 0.53.

Figure 2

Unsupervised clustering and analyses based on the RBN dataset

a Heatmap depicting protein levels after unsupervised hierarchical clustering of the RBN dataset consisting of 3,467 cancer samples across 11 tumor types and 181 antibodies. Protein levels are indicated on a low-to-high scale (blue-white-red). Eight clusters are defined. Cluster_A has been subdivided into two clusters (A1 and A2), based on the differences between BRCA reactive and remaining luminal subtypes5. Annotation bars include tumor type (BRCA-basal separately indicated); purity and ploidy (ABSOLUTE algorithm); stromal and immune scores (ESTIMATE algorithm); BRCA (PAM50 classification) and BLCA subtype; 16 significantly mutated genes and two frequently observed amplifications. The statistical significance of correlations between the clusters and each variable is indicated to the left of each annotation bar (n=3,467, chi-squared, Fisher’s Exact, and ANOVA’s F test. See Methods).

b Crosstab showing the number of tumor samples in each cluster.

c-e Kaplan Meier curves showing overall survival of (c) the BRCA located in four separate clusters (A1, A2, E and F, n=740), (d) KIRC in cluster_F vs. KIRC in other clusters (n=454) and (e) BLCA in cluster_B vs. BLCA in other clusters (n=127). Follow-up was capped at 60 months due to limited number of events beyond this time. Statistical difference in outcome between groups is indicated by P-value (log-rank test). A high-resolution, interactive version of the heatmap with zooming capability, can be found at (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

Figure 3

Unsupervised clustering and analyses based on the MC dataset

a Heatmap showing protein expression after unsupervised hierarchical clustering of 3,467 cancer samples across 11 tumor types and 181 antibodies. Protein levels are indicated on a low-to-high scale (blue-white-red). Seven clusters were defined. Cluster_II has been subdivided manually into two clusters (IIa and IIb) based on significant difference in expression of the proteins of interest (HER2 and EGFR). Annotation bars include tumor lineage (BRCA-basal separately indicated), purity and ploidy (ABSOLUTE algorithm); stromal and immune scores (ESTIMATE algorithm); BRCA (PAM50 classification) and BLCA subtype; 16 significantly mutated genes and two frequently observed amplifications. Statistical significance of the correlations between the clusters and each variable is indicated left of the annotation bars (n=3,467, chi-squared, Fisher’s Exact, and ANOVA’s F test. See Methods).

b Crosstab showing the number of tumor samples in each cluster.

c-g Kaplan Meier curves showing overall survival in (c) the KIRC in cluster_VII vs. in all other clusters (n=454), (d) OVCA in cluster_VII vs. in all other clusters (n=412), (e) KIRC in cluster_IV vs. in all other clusters (n=454), (f) LUSC in cluster_V vs. in all other clusters (n=195) and (g) COAD in cluster_V vs. in all other clusters (n=334). Follow-up has been capped at 60 months months, due to limited number of events beyond this time. Statistical difference in outcome between groups is indicated by P-value (log-rank test). A high-resolution, interactive version of the heatmap with zooming capability, can be found at (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

Figure 4

Pathway analyses

Pathway analyses of the dataset by RBN clusters, MC clusters and tumor type. For pathway predictor members see Supplementary Table 13.

a-b Heatmaps depicting mean pathway scores after unsupervised hierarchical clustering on tumor lineages and protein clusters based on the (a) RBN and (b) MC datasets. The heatmaps were clustered on both axes. As expected, RBN clusters show a strong association with tumor lineages, with very similar patterns between them, whereas MC clusters do not associate with any particular tumor lineage.

c-f The heatmaps, supervised on the sample axis, depict the protein levels of the pathway members and of proteins with a high correlation (ρ>0.3/ ρ<−0.3, Spearman’s correlation) to the pathway predictor across RBN clusters (c-d) and tumor lineages (e-f). The EMT pathway (c and e) and the hormone_a pathway (d and f) are shown. Samples are first sorted by either cluster (c-d) or tumor lineage (e-f), then by pathway score (from low to high) within cluster or tumor lineage. Dotplots (lower panel) represent the pathway score for each sample. Each box represents the lower quartile, median and upper quartile, whereas the whiskers represent the most extreme data point within 1.5 × inter-quartile range from the edge of the box. Annotation bars (selected from Fig. 2) are included if statistically associated with the pathway score (P <0.05, Kruskal-Wallis test, n=3,467). Pathway members are marked in red on the left hand side. High-resolution images of the heatmaps can be found online (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

Figure 5

Analyses of selected potentially actionable proteins

a-b Heatmaps, supervised on the sample axis, depicting protein level of 25 proteins that are (potentially) actionable based on the RBN dataset. Proteins were ordered by unsupervised hierarchical clustering and samples were ordered by (a) cluster and (b) tumor lineage membership and within each ordered by unsupervised hierarchical clustering. Annotation bars include tumor lineage, purity and ploidy (ABSOLUTE algorithm); stromal and immune scores (ESTIMATE algorithm); BRCA (PAM50 classification) and BLCA subtype; 16 significantly mutated genes and two frequently observed amplifications. High-resolution images of the heatmaps can be found online (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

Figure 6

Unbiased data-driven signaling network

Unbiased signaling network based on a probabilistic graphical models analysis, visualizing all 11 tumor lineages individually. Interplay between nodes was quantified using scores from the graphical model analysis (see Methods), that identify links between nodes whilst controlling for the effects of all other observed nodes. The strength of links is indicated by the thickness of the line whilst the color indicates the tumor lineage in which the link was observed; only the strongest links are shown. Nodes in white are related nodes that were highly correlated and therefore merged prior to network analysis. The adjacent correlated (green) node was then used for network generation. Positive (negative) correlations are indicated with continuous (dotted) lines. A high-resolution image of the network can be found online (http://bioinformatics.mdanderson.org/main/TCGA/Pancan11/RPPA).

Footnotes

Conflict of interest

The authors have no conflicts of interest to declare.

Author contributions

R.A. and G.B.M. lead the analysis in the paper. G.B.M. supervised the project and manuscript writing. R.A., P.K.S.N, and H.M.J.W. wrote much of the manuscript, supplement, and response to reviewers. J.N.W. edited the manuscript and response to reviewers. M.S., G.B.M., P.K.S.N., Y.L., H.M.J.W., J.N.W., F.M.B. and R.A. helped in interpreting the results. F.Z., Z.J., W.L., S.L. and J.Y. performed analyses included in the paper and in the response to reviewers. K.Y. and R.G.W.V. performed purity analysis. J.L., H.L., B.M.B. and P.K.S.N. made the data and analyses available to the public at the TCPA portal and at the project home page. E.G.S. and P.T.R. performed function space analysis of the proteins. J.D.M., L.D., P.T., J.V.H., and L.A.B. provided unpublished lung cancer cell line data for validation. S.M.H., F.D., N.S. and S.M. performed network analyses. Y.L. generated the RPPA data under G.B.M.’s supervision.

Acknowledgements

We acknowledge contributions from the TCGA Research Network and its TCGA Pan-Cancer Analysis Working Group coordinated by J.M. Stuart, C. Sander and I. Shmulevich. This study was supported by the National Institutes of Health (TCGA CA143883 and CCSG grant P30 CA016672 to J.N.W. and G.B.M., NCI P50CA70907 to J.D.M. and J.V.H., NCI U54 CA112970 to S.M.); and by the Chapman Foundations and the Michael & Susan Dell Foundation (Lorraine Dell Program in Bioinformatics for Personalization of Cancer Medicine) to J.N.W. Additional support was provided to S.M. by the Cancer Systems Biology Center grant from the Netherlands Organisation for Scientific Research. L.A.B. acknowledges support from the UTMDACC Physician Scientist Award, LUNGevity Foundation, the North Carolina Chapter of National Lung Cancer Partnership, and The Sidney Kimmel Foundation for Cancer Research.

References

  • 1. MyhreSInfluence of DNA copy number and mRNA levels on the expression of breast cancer related proteinsMolecular oncology20137704718[PubMed][Google Scholar]
  • 2. ParkESIntegrative analysis of proteomic signatures, mutations, and drug responsiveness in the NCI 60 cancer cell line setMolecular cancer therapeutics20109257267[PubMed][Google Scholar]
  • 3. ShankavaramUTTranscript and protein expression profiles of the NCI-60 cancer cell panel: an integromic microarray studyMolecular cancer therapeutics20076820832[PubMed][Google Scholar]
  • 4. Cancer Genome AtlasNComprehensive molecular characterization of human colon and rectal cancerNature2012487330337[PubMed][Google Scholar]
  • 5. Cancer Genome AtlasNComprehensive molecular portraits of human breast tumoursNature20124906170[PubMed][Google Scholar]
  • 6. Cancer Genome Atlas ResearchNComprehensive genomic characterization defines human glioblastoma genes and core pathwaysNature200845510611068[PubMed][Google Scholar]
  • 7. Cancer Genome Atlas ResearchNIntegrated genomic analyses of ovarian carcinomaNature2011474609615[PubMed][Google Scholar]
  • 8. Cancer Genome Atlas ResearchNComprehensive genomic characterization of squamous cell lung cancersNature2012489519525[PubMed][Google Scholar]
  • 9. Cancer Genome Atlas ResearchNComprehensive molecular characterization of clear cell renal cell carcinomaNature20134994349[PubMed][Google Scholar]
  • 10. Cancer Genome Atlas ResearchNIntegrated genomic characterization of endometrial carcinomaNature20134976773[PubMed][Google Scholar]
  • 11. LiJTCPA: a resource for cancer functional proteomics dataNature methods20131010461047[PubMed][Google Scholar]
  • 12. PayneSJBowenRLJonesJLWellsCAPredictive markers in breast cancer--the presentHistopathology2008528290[PubMed][Google Scholar]
  • 13. SorlieTGene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implicationsProceedings of the National Academy of Sciences of the United States of America2001981086910874[PubMed][Google Scholar]
  • 14. CardnellRJProteomic markers of DNA repair and PI3K pathway activation predict response to the PARP inhibitor BMN 673 in small cell lung cancerClinical cancer research : an official journal of the American Association for Cancer Research20131963226328[PubMed][Google Scholar]
  • 15. AwayaHTakeshimaYFuronakaOKohnoNInaiKGene amplification and protein expression of EGFR and HER2 by chromogenic in situ hybridisation and immunohistochemistry in atypical adenomatous hyperplasia and adenocarcinoma of the lungJournal of clinical pathology20055810761080[PubMed][Google Scholar]
  • 16. BlokEJKuppenPJvan LeeuwenJESierCFCytoplasmic Overexpression of HER2: a Key Factor in Colorectal CancerClinical Medicine Insights Oncology201374151[PubMed][Google Scholar]
  • 17. CanerVNo strong association between HER-2/neu protein overexpression and gene amplification in high-grade invasive urothelial carcinomasPathology oncology research : POR200814261266[PubMed][Google Scholar]
  • 18. FleischmannARotzerDSeilerRStuderUEThalmannGNHer2 amplification is significantly more frequent in lymph node metastases from urothelial bladder cancer than in the primary tumoursEuropean urology201160350357[PubMed][Google Scholar]
  • 19. GrobTJHeterogeneity of ERBB2 amplification in adenocarcinoma, squamous cell carcinoma and large cell undifferentiated carcinoma of the lungModern pathology : an official journal of the United States and Canadian Academy of Pathology, Inc20122515661573[Google Scholar]
  • 20. SlomovitzBMHer-2/neu overexpression and amplification in uterine papillary serous carcinomaJournal of clinical oncology : official journal of the American Society of Clinical Oncology20042231263132[PubMed][Google Scholar]
  • 21. CuadrosMVillegasRSystematic review of HER2 breast cancer testingApplied immunohistochemistry & molecular morphology : AIMM / official publication of the Society for Applied Immunohistochemistry20091717[PubMed][Google Scholar]
  • 22. GrimmEESchmidtRASwansonPEDintzisSMAllisonKHAchieving 95% crossmethodological concordance in HER2 testing: causes and implications of discordant casesAmerican journal of clinical pathology2010134284292[PubMed][Google Scholar]
  • 23. PressMFHER-2 gene amplification, HER-2 and epidermal growth factor receptor mRNA and protein expression, and lapatinib efficacy in women with metastatic breast cancerClinical cancer research : an official journal of the American Association for Cancer Research20081478617870[PubMed][Google Scholar]
  • 24. YazijiHHER-2 testing in breast cancer using parallel tissue-based methodsJAMA : the journal of the American Medical Association200429119721977[PubMed][Google Scholar]
  • 25. BarginearMFJohnVBudmanDRTrastuzumab-DM1: a clinical update of the novel antibody-drug conjugate for HER2-overexpressing breast cancerMolecular medicine20121814731479[PubMed][Google Scholar]
  • 26. HurvitzSAPhase II randomized study of trastuzumab emtansine versus trastuzumab plus docetaxel in patients with human epidermal growth factor receptor 2-positive metastatic breast cancerJournal of clinical oncology : official journal of the American Society of Clinical Oncology20133111571163[PubMed][Google Scholar]
  • 27. FlemingGFPhase II trial of trastuzumab in women with advanced or recurrent, HER2-positive endometrial carcinoma: a Gynecologic Oncology Group studyGynecologic oncology20101161520[PubMed][Google Scholar]
  • 28. YaoMVHL tumor suppressor gene alterations associated with good prognosis in sporadic clear-cell renal carcinomaJournal of the National Cancer Institute20029415691575[PubMed][Google Scholar]
  • 29. VivancoIDifferential sensitivity of glioma- versus lung cancer-specific EGFR mutations to EGFR kinase inhibitorsCancer discovery20122458471[PubMed][Google Scholar]
  • 30. BacolodMDBaranyFMolecular profiling of colon tumors: the search for clinically relevant biomarkers of progression, prognosis, therapeutics, and predispositionAnnals of surgical oncology20111836943700[PubMed][Google Scholar]
  • 31. ImamuraYSpecific mutations in KRAS codons 12 and 13, and patient prognosis in 1075 BRAF wild-type colorectal cancersClinical cancer research : an official journal of the American Association for Cancer Research20121847534763[PubMed][Google Scholar]
  • 32. MalhotraPAlterations in K-ras, APC and p53-multiple genetic pathway in colorectal cancer among IndiansTumour biology : the journal of the International Society for Oncodevelopmental Biology and Medicine20133419011911[PubMed][Google Scholar]
  • 33. BremnesRMHigh-throughput tissue microarray analysis used to evaluate biology and prognostic significance of the E-cadherin pathway in non-small-cell lung cancerJournal of clinical oncology : official journal of the American Society of Clinical Oncology20022024172428[PubMed][Google Scholar]
  • 34. ZhangHClinical significance of E-cadherin, beta-catenin, vimentin and S100A4 expression in completely resected squamous cell lung carcinomaJournal of clinical pathology2013[Google Scholar]
  • 35. LaplanteMSabatiniDMmTOR signaling in growth control and diseaseCell2012149274293[PubMed][Google Scholar]
  • 36. HodgsonMCDecreased expression and androgen regulation of the tumor suppressor gene INPP4B in prostate cancerCancer research201171572582[PubMed][Google Scholar]
  • 37. LiuSExpression of autotaxin and lysophosphatidic acid receptors increases mammary tumorigenesis, invasion, and metastasesCancer cell200915539550[PubMed][Google Scholar]
  • 38. PratAAdamoBCheangMCAndersCKCareyLAPerouCMMolecular characterization of basal-like and non-basal-like triple-negative breast cancerThe oncologist201318123133[PubMed][Google Scholar]
  • 39. CheungLWHigh frequency of PIK3R1 and PIK3R2 mutations in endometrial cancer elucidates a novel mechanism for regulation of PTEN protein stabilityCancer discovery20111170185[PubMed][Google Scholar]
  • 40. SalvesenHBHaldorsenISTrovikJMarkers for individualised therapy in endometrial carcinomaThe lancet oncology201213e353e361[PubMed][Google Scholar]
  • 41. SlomovitzBMColemanRLThe PI3K/AKT/mTOR pathway as a therapeutic target in endometrial cancerClinical cancer research : an official journal of the American Association for Cancer Research20121858565864[PubMed][Google Scholar]
  • 42. HoriuchiDMYC pathway activation in triple-negative breast cancer is synthetic lethal with CDK inhibitionThe Journal of experimental medicine2012209679696[PubMed][Google Scholar]
  • 43. DanielsTRDelgadoTHelgueraGPenichetMLThe transferrin receptor part II: targeted delivery of therapeutic agents into cancer cellsClinical immunology2006121159176[PubMed][Google Scholar]
  • 44. FriedmanJHastieTTibshiraniRSparse inverse covariance estimation with the graphical lassoBiostatistics20089432441[PubMed][Google Scholar]
  • 45. RueHHeldLGaussian Markov Random Fields: Theory and Applications2005
  • 46. ChagparRBDirect positive regulation of PTEN by the p85 subunit of phosphatidylinositol 3-kinaseProceedings of the National Academy of Sciences of the United States of America201010754715476[PubMed][Google Scholar]
  • 47. MemmottRMDennisPAAkt-dependent and -independent mechanisms of mTOR regulation in cancerCellular signalling200921656664[PubMed][Google Scholar]
  • 48. SerraVRSK3/4 mediate resistance to PI3K pathway inhibitors in breast cancerThe Journal of clinical investigation201312325512563[PubMed][Google Scholar]
  • 49. ShawRJCantleyLCRas, PI(3)K and mTOR signalling controls tumour cell growthNature2006441424430[PubMed][Google Scholar]
  • 50. CanelMOverexpression of focal adhesion kinase in head and neck squamous cell carcinoma is independent of fak gene copy numberClinical cancer research : an official journal of the American Association for Cancer Research20061232723279[PubMed][Google Scholar]
  • 51. MyllykangasSIntegrated gene copy number and expression microarray analysis of gastric cancer highlights potential target genesInternational journal of cancer Journal international du cancer2008123817825[PubMed][Google Scholar]
  • 52. MaLChenZErdjument-BromageHTempstPPandolfiPPPhosphorylation and functional inactivation of TSC2 by Erk implications for tuberous sclerosis and cancer pathogenesisCell2005121179193[PubMed][Google Scholar]
  • 53. SmythGKLimma: linear models for microarray dataBioinformatics and Computational Biology Solutions using R and Bioconductor2005New YorkSpringer[Google Scholar]
  • 54. Robnik-SikonjaMSavickyPCORElearn: CORElearn - classification, regression, feature evaluation and ordinal evaluationR package version 0.9.412013(ed^(eds)[Google Scholar]
  • 55. UhlenMTowards a knowledge-based Human Protein AtlasNature biotechnology20102812481250[Google Scholar]
  • 56. HallMFrankEHolmesGPfahringerBReutemannPWittenIHThe WEKA Data Mining Software: An UpdateSIGKDD Explorations200911[Google Scholar]
  • 57. EisenMBSpellmanPTBrownPOBotsteinDCluster analysis and display of genomewide expression patternsProceedings of the National Academy of Sciences of the United States of America1998951486314868[PubMed][Google Scholar]
  • 58. HartiganJAClustering Algorithms1975John Wiley & Sons
  • 59. BroomBNGCHM: Utilities for creating Next Generation Clustered Heat MapsR package version 0.5.12013(ed^(eds)[Google Scholar]
  • 60. NeeleyESKornblauSMCoombesKRBaggerlyKAVariable slope normalization of reverse phase protein arraysBioinformatics20092513841389[PubMed][Google Scholar]
  • 61. HennessyBTA Technical Assessment of the Utility of Reverse Phase Protein Arrays for the Study of the Functional Proteome in Non-microdissected Human Breast CancersClinical proteomics20106129151[PubMed][Google Scholar]
  • 62. CarterSLAbsolute quantification of somatic DNA alterations in human cancerNature biotechnology201230413421[Google Scholar]
  • 63. YoshiharaKInferring tumour purity and stromal and immune cell admixture from expression dataNature communications201342612[Google Scholar]
  • 64. BenjaminiYHochbergYControlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society Series19955712[Google Scholar]
  • 65. CoombesKRClassDiscovery: Classes and methods for "class discovery" with microarrays or proteomicsR package version 2.13.42012(ed^(eds)[Google Scholar]
  • 66. ZhaoTiuHRoederKLaffertyJWassermanLThe huge Package for Highdimensional Undirected Graph Estimation in RThe Journal of Machine Learning Research2012134[Google Scholar]
Collaboration tool especially designed for Life Science professionals.Drag-and-drop any entity to your messages.