Pan-Cancer Molecular Classes Transcending Tumor Lineage Across 32 Cancer Types, Multiple Data Platforms, and over 10,000 Cases

This article is featured in Highlights of This Issue, p. 2027

Cancer is not a single disease, and at the molecular level there is widespread heterogeneity that may be observed from patient to patient. Nevertheless, unsupervised classification of tumors on the basis of molecular profiling data can reveal major subtypes existing within a given cancer type according to tissue of origin (1, 2). Such molecular-based subtypes can reflect altered pathways within different cancer subsets, which could have important implications for applying existing therapies or for developing new therapeutic approaches (3). The Cancer Genome Atlas (TCGA), a large-scale effort to comprehensively characterize over 10,000 human cancers at the molecular level, provides a common platform for the study of diverse cancer types, with multiple levels of data including mRNA, miRNA, protein, DNA methylation, copy number, and mutation (2). For most cancer types represented in TCGA, an individual study of the molecular landscape of that cancer type was carried out (2).

With data generation completed, there is opportunity for systematic analyses of the entire TCGA pan-cancer cohort (4), including defining molecular subtypes and associated pathways relevant to multiple cancer types. One challenge, in the identification of cancer subsets transcending the tissue of origin, is that widespread molecular patterns are associated with tumor lineage and histology (5–7). Although TCGA datasets have been harmonized to allow for cross-cancer type comparisons—which is useful for addressing a host of questions—an alternative approach to molecular classification on the basis of these data would be to first computationally subtract the molecular differences between cancer types (8, 9). This alternative approach would have the effect of consolidating the individual subtypes that might be discoverable in individual cancer types into super-types or pan-cancer “classes” that transcend tissue or histology distinctions.

The aim of this study was to define pan-cancer molecular-based subtypes or “classes,” which would transcend tumor lineage across the over 10,000 human cancers and 32 cancer types profiled by TCGA, using data from multiple molecular profiling platforms to characterize these classes in terms of associated pathways.

Results are based upon data generated by TCGA Research Network (http://cancergenome.nih.gov/). Molecular data were aggregated from public repositories. Tumors spanned 32 different TCGA projects, each project representing a specific cancer type, listed as follows: LAML, acute myeloid leukemia; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; LGG, lower grade glioma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CHOL, cholangiocarcinoma; CRC, colorectal adenocarcinoma (combining COAD and READ projects); ESCA, esophageal carcinoma; GBM, glioblastoma multiforme; HNSC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; DLBC, lymphoid neoplasm diffuse large B-cell lymphoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumors; THYM, thymoma; THCA, thyroid carcinoma; UCS, uterine carcinosarcoma; UCEC, uterine corpus endometrial carcinoma; UVM, uveal melanoma. Cancer molecular profiling data were generated through informed consent as part of previously published studies and analyzed in accordance with each original study's data use guidelines and restrictions.

For the mRNA platform, log2-transformed expression values within each cancer type (as defined by TCGA project) were normalized to standard deviations from the median. By k-means clustering method (using the “kmeans” R function, with 10 clusters and 1000 maximum iterations and 25 random sets), cancer cases were subtyped, based on the top 2,000 features with the most variable expression on average by cancer type (using standard deviation of the log2-transformed expression values as a measure of variability, excluding genes on X or Y chromosomes).

We defined the top differential genes associated with each subtype, or pan-cancer “class.” Taking the top 2,000 mRNA features, we first computed the two-sided t test for each gene and each class, comparing expression levels of each class with that of the rest of the tumors. We then selected the top 100 genes with the lowest P-value for each subtype; however, for the c2 class only three genes were associated, and for the c7 class only 51 genes were associated, resulting in 854 top class-specific genes in all. In a similar manner, top features associated with pan-cancer class for RPPA, miRNA, and DNA methylation platforms were identified.

The Fantom datasets of gene expression by cell type (10) were analyzed using a previously utilized approach (7). Briefly, the top 2,000 most variable mRNAs (used above for the clustering analyses) were examined in Fantom. Logged expression values for each gene in the fantom dataset were centered on the median of sample profiles. For each fantom differential expression profile, the inter-profile correlation (Spearman's) was taken with that of each TCGA pan-cancer differential expression profile (with genes normalized within each TCGA project to standard deviations from the median).

We examined an external gene expression profiling dataset of multiple cancer types from the Expression Project for Oncology (expO; GSE2109), classifying each external tumor profile by pan-cancer class as defined by TCGA data. Within each cancer type, log2-transformed genes in the expO dataset were normalized to standard deviations from the median. As a classifier, the top set of 854 mRNAs distinguishing between the TCGA pan-cancer classes was used. For each pan-cancer class, the average value for each gene was computed, based on the centered TCGA expression data matrix. The Pearson's correlation between each expO profile and each TCGA pan-cancer class averaged profile was computed. Each expO case was assigned to TCGA pan-cancer class, based on which class profile showed the highest correlation with the given expO profile.

All P values were two-sided unless otherwise specified. Scoring of expression profiles for pathway-associated gene signatures was carried out as previously described (9, 11). Additional methods details are provided in supplemental. Supplementary Fig. S1 provides a schematic of the various analyses performed and their relation to supplementary data files.

Across cancer types, tumor lineage and squamous histology were found to represent major determinates of unsupervised molecular classification. In total, our study involved 11,232 human cancer cases representing 32 different major types, for which TCGA generated data on one or more of the following molecular characterization platforms (Supplementary Data S1): whole exome sequencing (WES, n = 10,224 cases), somatic DNA copy by SNP array (n = 10,845), RNA-seq (n = 10,224), miRNA-seq (n = 10,128), DNA methylation (n = 10,959), and reverse phase protein array (RPPA, n = 7,663). Using established analytical approaches (5–7), 10,662 TCGA cases (with data available for at least three platforms) were subtyped according to each of the above data platforms excluding WES, with the various subtype calls for each sample then being consolidated to define multiplatform-based molecular subtypes.

With hierarchical clustering of the platform-level subtype assignments (Supplementary Figs. S2A and S3), the cases segregated largely according to cancer type, with higher level branches of the clustering tree representing major tissue-based categories including breast, liver, kidney, brain, gastrointestinal, squamous (head and neck, lung squamous, esophageal, bladder, cervical), lung adenocarcinoma, immune-related, uterus, and skin. In this pan-cancer molecular subtyping, basal-like BRCA (breast cancer) was molecularly distinct from other BRCA. K-means clustering was applied to the platform-level subtype assignment matrix to formally define a discrete number of subtypes (Supplementary Figs. S2B, 4A, and 4B, and Supplementary Data S2). With a 25-subtype solution, for example (Supplementary Fig. S2B), most subtypes were specific to either a single cancer type or multiple types related by common tissue of origin.

Although the above pan-cancer molecular subtypes were found to be highly concordant with their tissue-of-origin counterparts, we went on to pursue an alternate analytical strategy to define subtypes that would not be driven by tumor lineage-specific markers or patterns. For expression (mRNA, miRNA, protein) and DNA methylation datasets, values were first normalized or centered within each cancer type, thereby computationally removing any tissue or histology dominant effects. For defining pan-cancer subtypes using the above datasets, we had initially considered multiplatform-based solutions (e.g., analyzing results from the five data platforms together to define subtypes; ref. 9), but instead opted to use RNA-seq platform to define subtypes and then examine the other platforms for correlates of these RNA-seq-based subtypes. RNA-seq data would represent the full range of features for cancer-relevant pathways, whereas miRNA and protein data platforms both represent expression-based data but with a more limited number of features and represented pathways as compared to RNA-seq data. In addition, DNA copy alterations (while having a normal control) vary by cancer type (Supplementary Fig. S3 and ref. 9), and so would serve to segregate samples by cancer type.

The RNA-seq platform (with values normalized within cancer type) was used to define 10 different subtypes or “classes” of cancer (Table 1) across the 10,224 cases in TCGA cohort with available data (where with solutions of more than 10 subtypes no appreciable increase in clustering consensus was observed; Supplementary Fig. S5A), with the other profiling platforms then being used to characterize these classes in terms of associated pathways as described below. These pan-cancer molecular classes, referred to here as “c1” to “c10,” were each characterized by widespread molecular patterns (Fig. 1A; Supplementary Figs. S5B–S5F and S6). For each class, the top genes most differentially expressed in the given class versus the rest of the tumors were identified (Fig. 1A; Supplementary Data S2), where the top differential patterns—involving 854 genes in all—could be observed to span cancer type. Notably, very few genes were strongly associated with the c2 class in particular, which represented one distinguishing feature of this class as compared to the other classes, although this c2 class would represent a basis of comparison. Normalized differential expression patterns for specific genes of interest—including MYC, MKI67, CTAG1B, HIF1A, CD274, VIM, and ZEB1—could further distinguish between the classes and suggested associated pathways as further explored below.

Specific miRNAs and proteins (with values normalized within cancer type) could distinguish between the pan-cancer molecular classes (Fig. 1B; Supplementary Fig. S5E and Supplementary Data S2), with, for example, miR200 family members showing the lowest differential expression in c8 and c7 classes, and with immune cell protein markers LCK and SYK being highest in the c3 class. Differential DNA methylation patterns could also distinguish between the subtypes (Fig. 1C; Supplementary Figs. S5B and S5F, and Supplementary Data S2), with distinctive patterns associated with c5 class in particular. For a subset of genes, significant anticorrelations between methylation and expression could be identified (Supplementary Fig. S5F and Supplementary Data S2). Within the top differentially expressed genes underscoring each pan-cancer molecular class, specific gene categories were overrepresented (Supplementary Fig. S7), including immune-related genes being most highly expressed in the c3 class, neuron-related genes being highly expressed in the c4 class, and cytoskeleton- and keratin-related genes being highly expressed in the c9 class. The c4, c5, and c6 classes tended to show a higher degree of genome-wide copy alterations (Fig. 1B; Supplementary Fig. S5B). Each of the pan-cancer classes were found to span cases from multiple cancer types (Figs. 1B and C), with the notable exception of c5, which consisted entirely of BRCA cases (n = 187) and represents the basal-like breast cancer molecular subtype (ref. 12; Fig. 1B). As compared to the other pan-cancer classes, the c5 class was associated with better overall patient survival (Fig. 1D); however, compared with other breast cancers, c5 has a poorer survival (3). Overall, the pan-cancer molecular classes showed significant concordances with other molecular subtype designations, which had been previously made for a subset of TCGA cases in individual cancer type studies (refs. 6, 12–24; Fig. 1E).

Across the entire TCGA pan-cancer cohort, assessment of genes within pathways demonstrated a high number of somatic alterations (mutation, copy alteration, or epigenetic silencing) involving p53 (62.7% of 10,224 cases with exome data available), PI3K/AKT/mTOR (44.1%), receptor tyrosine kinase signaling (RTK, 43.9%), chromatin modification (40.8%), SWI/SNF complex (32.0%), Wnt/β-catenin (22.3%), MYC (10.9%), NRF2-ARE (9.3%), and Hippo signaling (4.5%; Supplementary Fig. S8A). The above pathways were found to be altered in different ways involving different genes in different cancer types (Supplementary Fig. S9A). A number of pathway-level or individual gene-level alterations surveyed were highly represented within specific cancer types or pan-cancer classes (Supplementary Figs. S8B and S9B). In particular, TP53 mutations were highly represented within both c4 and c5 tumors, and MYC amplifications were highly represented in c5 tumors (Supplementary Figs. S8B and 9B). Furthermore, c10 tumors were enriched for somatic alterations involving p53 pathway (including TP53 and ATM mutations), RTKs (MET mutations), and numerous genes involving chromatin modification and SWI/SNF complex (Supplementary Figs. S8B and S9B). Some of the somatic mutation associations observed might be attributable to cancer type- or mutation rate-specific patterns (Supplementary Figs. S10 and S11); for example, TP53 and chromatin modifier mutations being enriched within c4 tumors would reflect in part the types of cancers more highly represented within c4 (BLCA, CESC, HNSC, LUSC, etc.).

We sought to examine the effects on pathway activation—as measured by mRNA or protein signature—of somatic alterations impacting specific pathways noted above. Previously-defined gene expression signatures for p53, k-ras, MTOR, Wnt/β-catenin, MYC, NRF2, and YAP1 (6, 25–27) were applied to the mRNA or protein expression profiles of the TCGA samples, whereby each sample profile was scored for each of the above pathways. For each pathway considered, relative levels of the corresponding signature were significantly different between somatically altered versus unaltered cases for that pathway, and the differences were in the anticipated direction (e.g. p53-related alterations were associated with decreased levels of p53 transcriptional targets, and MTOR-related alterations were associated with increased MTOR proteomic signaling; Supplementary Fig. S8C). These results would demonstrate a widespread level of concordance between observed DNA-level alterations and global expression patterns. At the same time, within the somatically altered and unaltered groups for each pathway, a wide range of genes signature levels were evident (Supplementary Fig. S8C), suggesting that other factors besides somatic mutation or copy alteration (e.g., tumor microenvironmental influences) may impact pathways.

To gain insight into pathways that would distinguish between the various pan-cancer molecular classes, we applied a number of gene signatures to TCGA expression profiles with values normalized within each cancer type. A number of pathways appeared more or less active for different pan-cancer classes (Fig. 2A), as further explored below. For example, gene signature scores for epithelial–mesenchymal transition (EMT), hypoxia, NRF2/KEAP1, Wnt, and Notch all were higher in c3, c7, and c8 classes, as compared to the rest of the tumors. Focusing here on EMT, representing mesenchymal features, we observed a strong negative association between the expression of miR-200 family members and their promoter methylation levels (Fig. 2B and 2C; Supplementary Fig. S12A), which demonstrates epigenetic regulation of these miRNAs across a large subset of human cancers. The miR200 family suppresses ZEB1, a key transcriptional regulator of EMT. Across the entire TCGA cohort, we examined differential expression values for a set of genes and proteins representing canonical mesenchymal or epithelial markers (11), as well as for miR200 family and ZEB1 genes (Fig. 2B). Manifestation of mesenchymal features was highest in the c7 and c8 tumors and lowest in the c6 tumors (which appeared more epithelial), and miR200 expression was lowest in c7 and c8 tumors, with associated DNA methylation being highest in c8 tumors; c3 tumors also showed mesenchymal-associated patterns but less strongly than did c7 and c8 tumors, and c9 tumors showed high expression of some mesenchymal and epithelial marker genes, as well as high miR200. The c7 and c8 pan-cancer classes each involved a wide range of cancer types (Fig. 2D).

For some tumor subsets, the observed associations with mesenchymal features could conceivably be attributable to surrounding stromal cells within the tumor sample, as well as to changes within the actual cancer cells. For example, invasive lobular carcinoma (ILC), the second most prevalent histologic subtype of invasive breast cancer, is characterized by small discohesive neoplastic cells invading the stroma in a single-file pattern (12). Immunohistochemical analysis in breast ILC has demonstrated that the neoplastic lobular cells tend to retain their epithelial identity, with mesenchymal markers being expressed by the fibroblasts in the prominent stromal component of ILC (28). However, for other cancer types, changes within cancer cells from epithelial to mesenchymal states may occur at the tumor invasive front, although not within the main tumor bulk (29). Of our 10 pan-cancer molecular classes, four were characterized by relatively lower sample purity: c3, c5 (representing basal-like breast cancer), c7, and c8 (Fig. 1B; Supplementary Figs. S12B and S12C). Interestingly, of these four classes, only the c8 class was strongly associated with the breast ILC cases in TCGA (Supplementary Fig. S12B). In contrast, renal cell carcinomas of the previously described “CC-e.3” molecular subtype manifesting EMT (30) associated with the c7 class but not with the c8 class, and the EMT-associated “SQ.1” molecular subtype of lung cancer (7) was distributed among c3, c7, and c8 classes (Supplementary Fig. S12B). Three molecular subtypes previously associated with mesenchymal features for specific cancer types, GBM:mesenchymal, OV:mesenchymal, and HNSC:mesenchymal, were each respectively associated with c3, c7, and c8 pan-cancer classes (Fig. 1E). In comparison to c7, c8 was further distinguished by high expression of fatty acid metabolism genes (Fig. 2A), and associated hypoxia-related changes in c7 and c8, for example, suggest an altered microenvironment. All of this would indicate that the molecular classes associated with lower purity would each represent distinctive biology, apart from the technical aspects involved with sample collection.

Analysis of gene expression patterns of normal tissues and cells can provide meaningful context to the widespread differential patterns observed in cancer (7, 31). We examined the top 2,000 differential mRNAs from TCGA cohort (from Supplementary Fig. S5B) in a public expression dataset from the Fantom consortium (10) of 850 profiles representing various human cell and tissue specimens. Inter-correlations between Fantom profiles and TCGA profiles revealed each pan-cancer molecular class to manifest distinctive patterns of global similarity or dissimilarity with specific categories of normal cells and tissues (Fig. 3A; Supplementary Data S5). In particular, c3, c5, and c10 classes were strongly associated with immune-related cells and tissues; the c4 class was strongly associated with cells and tissues related to the central nervous system (CNS), and c8 class was associated with CNS to a somewhat lesser degree; the c7 and c8 classes were associated with fibroblasts and other mesenchymal-related categories (reflecting the mesenchymal marker patterns observed above; Fig. 2B); and the c9 class was associated with adipose cells and heart tissues (which may be reflective of the association with fatty acid metabolism observed above; Fig. 2A). Analysis of the Fantom mouse expression dataset yielded similar associations (Supplementary Fig. S13A). In previous analyses utilizing TCGA expression profiles as normalized across all cancers, brain cancers and blood cancers associated as a group with fantom CNS and immune-related profiles, respectively (7).

Focusing here on the immune-related associations found for specific pan-cancer molecular classes, we went on to survey TCGA expression profiles for a set of genes representing potential targets for immunotherapy (6), including cancer-testis (CT) antigen genes and genes involved in immune checkpoint pathway (Fig. 3B), such as genes encoding PD1, PDL1, PDL2, and CTLA4. As a group, CT antigen genes were highest in both c4 and c5 classes, with a subset of CT antigen genes being particularly higher in c5 class (see Fig. 3B). However, immune checkpoint genes were most strongly differentially expressed within the c3 class, whereas the c8 and c10 classes also showed increased though relatively lower expression of these genes. Differential protein expression of T-cell marker LCK and B-cell marker SYK indicated the presence of T cells in c3, c8, and c10 classes, and the presence of B cells in c3 and c10 classes (Fig. 3B). Analysis of gene expression signatures from Bindea and colleagues (32) suggested that levels of infiltrating immune cell types were highest within the c3 class, followed by the c8 class, whereas the c10 class showed signatures for some but not all immune cell types (Supplementary Fig. S13B).

The above associations of c4 tumors with CNS tissues and cells suggested features of neuroendocrine tumors (NET), which arise from cells of the endocrine and nervous systems and which occur in many different tissues throughout the human body (33). Genes encoding canonical markers of NETs—including CHGA (chromogranin A), SYP (synaptophysin), NCAM1 (CD56), ENO2 (neuron-specific enolase), and CDX2—were all differentially higher in c4 tumors as compared with the other pan-cancer classes (Fig. 3C). In addition, a set of 51 genes in a panel of NET markers previously uncovered using gene expression profiling was examined (34), with the majority of these being differentially higher within the c4 class (Fig. 3C). TCGA LUAD cases previously determined to represent neuroendocrine cancer were also significantly enriched within the c4 class (Supplementary Fig. S12B, 6/14 cases, P < 1E−5, one-sided Fisher exact test). The c3, c10, and c4 pan-cancer classes each involved a wide range of cancer types (Fig. 3D), with c4 class in particular being composed of sizable percentages of head and neck, cervical, lung, bladder, ovarian, and gastrointestinal cancers.

The above pathway-associated gene signatures (Fig. 2A) indicated differential activation of specific pathways among the pan-cancer molecular classes. A number of these signatures were related to metabolism, including fatty acid metabolism, glycolysis/gluconeogenesis, pentose phosphate pathway, TCA cycle, and oxidative phosphorylation. The individual genes that comprised these signatures can be represented in a pathway diagram (Fig. 4A), whereby the c1, c6, and c8 classes were each observed to show evidence for the altered utilization of specific metabolic pathways, as compared to the rest of the cancer cases. Both c1 and c6 tumors showed increased expression of genes involved with oxidative phosphorylation, whereas c6 tumors also showed higher expression of genes involved in the TCA cycle and of genes involved in fatty acid synthesis. However, c8 tumors appeared to downregulate TCA cycle, oxidative phosphorylation, and fatty acid synthesis pathways, and to upregulate genes involving in fatty acid metabolism.

Deregulated pathways were also reflective of tumor microenvironmental effects at work in distinct subsets of cancers. Differential expression patterns involving the c3, c7, or c8 pan-cancer molecular classes were consistent with a consensus model (35) of tumor-associated macrophage (TAM) roles in the tumor microenvironment (Fig. 4B), whereby monocytes are first recruited to the tumor microenvironment (e.g., by CSF1 and CCL2). Tumor-secreted cytokines then have the potential to polarize recruited monocytes into TAMs, which play vital roles in a number of processes including immune suppression and EMT. EMT may also be induced by either hypoxia or Wnt signaling pathway (36–38), both of which appear increased in c3, c7, and c8 tumors (Figs. 2A and 4B); NRF2/KEAP1 and Notch pathways also appeared increased in these tumors. Immune suppression may involve multiple redundant immune checkpoints, which were most associated with c3, c5, c8, and c10 classes (Figs. 3B and 4C); these tumors showed higher expression of ligands such as PDL1 and PDL2, along with higher expression of the corresponding receptors associated with T cells.

The pan-cancer molecular class associations, as observed in TCGA datasets, were also examined in an external multicancer expression dataset. From the expO dataset, mRNA expression profiles of 2041 cancer cases representing 26 different types were obtained. Within each cancer type, expression values for each gene in the expO dataset were first normalized, and then each expO tumor profile was classified by pan-cancer molecular class as defined by TCGA data (Fig. 5A; Supplementary Figs. S14A and S14B), where the top set of 854 mRNAs distinguishing between our 10 classes (Fig. 1B) was used as the classifier. In the same manner as carried out above for TCGA datasets, expO expression profiles were also scored for mRNA signatures related to specific pathways or normal cell types (Fig. 5B; Supplementary Fig. S15), where similar overall trends of pathway-level differences between classes as originally observed in TCGA cohort could also be observed in the expO cohort. In particular, c5 tumors were composed almost entirely of breast cancer cases; anticipated pathway-level differences involving metabolism, YAP1, MYC, EMT, hypoxia, NRF2/KEAP1, WNT, NOTCH, and immune checkpoint were observed in expO cohort; CT antigen genes were higher as a group in c4, with the same subset also being high in c5 class as observed for TCGA; and neuroendocrine-associated global patterns and markers were associated with c4 class.

In a similar manner to that of the expO dataset, cell line profiles from the cancer cell line encyclopedia (CCLE) expression dataset (39) were each assigned to a pan-cancer class (Supplementary Fig. S16). Not all patterns of interest observable in human tumor data were as apparent in the CCLE results, which could be attributed to a number of factors including growth conditions of cell lines lacking tumor microenvironmental effects and differences in the types of cancers represented in CCLE; however, a number of associations were identifiable in CCLE, including c7 and c8 classes with EMT and hypoxia, and c4 class with neuroendocrine-related patterns (e.g., involving small cell lung cancers).

Although previous studies have greatly served to elucidate the molecular landscape of the individual cancer types represented in TCGA, our pan-cancer molecular classes provide an excellent framework for examining pathways or processes that would cut across these individual types. All but one of these molecular classes each involves a wide range of cancer types and would therefore be relevant to the study of multiple diseases. It is also remarkable that basal-like breast cancer forms its own molecular class—made up only of breast cancers—entirely distinct from all of the other cancer cases examined. Although basal-like breast cancer is already understood to represent a fundamentally different disease than other types of breast cancer, questions remain as to the origins of the observed molecular differences (40). The distinct patterns associated with basal-like breast cancer at the epigenetic and transcriptional levels, in this study, could suggest that this disease would have a different cell-of-origin from that of other breast cancers.

The landscape of pan-cancer–associated biological processes and pathways as uncovered here would include many that were well-established as having a functional role in the experimental setting, but for which the full extent of their involvement in human cancers may have been unclear. Processes and pathways such as metabolism, immune checkpoint, hypoxia, NRF2-ARE, HIPPO, Wnt, EMT, and Notch have been well characterized in the experimental setting, using models such as cell lines and mice (41). This study finds each of the above to involve large subsets of human cancers, as observed based on the coordinate expression patterns involving large numbers of genes. DNA mutation events alone (and tumor evolution by proxy) would not appear to represent the sole driver of these processes in cancer, where tumor microenvironment influences likely play a major role. For example, mutations in VHL or in RTK genes may accelerate processes of hypoxia or RTK signaling, respectively, or microenvironmental conditions such as lack of oxygen or availability of growth factors could conceivably achieve the same respective results. Individual molecular markers of our pan-cancer molecular classes that might seem most relevant from a therapeutic standpoint—including markers of processes and pathways noted above—could potentially be evaluated in the setting of patient care in future studies.

Multiple pan-cancer classes manifested patterns that we could ascribe to the noncancer cellular component of the tumor, including immune infiltrates and cancer-associated fibroblasts. Our finding of different stroma-associated tumor subsets, each with distinctive features that would distinguish it from the other pan-cancer classes, suggests various biological roles for the stromal component in human cancer. Three of our pan-cancer classes (c3, c5, and c10) showed particularly strong patterns of immune cell infiltrates, whereas one of these (c3) showing stronger patterns of immune checkpoint pathway genes and of genes involved in specific metabolic pathways. Three of our pan-cancer classes (c3, c7, and c8) showed patterns of mesenchymal cells, with two of these (c7 and c8) showing strong patterns of fibroblasts in particular, and with one of these (c8) showing patterns associated with fatty acid metabolism. The tumor microenvironment—which consists of a mixture of fibroblasts, myofibroblasts, endothelial cells, immune cells, other cells, and altered extracellular matrix—is understood to play an important role in the initiation and progression of various cancers (42, 43). For example, EMT of the cancer cells can occur in areas of fibrosis and may account for up to 40% of tumor-associated fibroblasts (43). As reflected in TCGA data, the various ways in which the tumor microenvironment can influence cancer may be involved within distinct subsets of human tumors.

Our discovery of a molecular class of cancers manifesting a differential expression profile associated with both NETs and with normal cells and tissues of the CNS would potentially have implications for a large subset of human cancers, with ∼4% of TCGA cases belonging to this “c4” class. Only through our analytical approach of subtracting out tissue-level molecular differences were we able to uncover this neuroendocrine-associated pan-cancer class (where otherwise such CNS-related patterns would be more strongly associated with TCGA brain cancers). NETs are understood to occur throughout the body, including in breast, cervical, lung, pancreas, ovarian, and gastrointestinal tissues, although it is believed that the incidence of these tumors is relatively rare, estimated at about 1% of cancers diagnosed in the United States (44, 45). If on the order of 4% of human cancers could be considered as NETs (or at least manifesting a molecular profile associated with these tumors), this would be well outside the range of previous estimates, suggesting that a large proportion of NETs are not diagnosed as such, but rather are grouped in with other cancers sharing the same tissue of origin. The tumor samples originally contributed to TCGA came from many different treatment facilities, with no uniform practices for determining neuroendocrine features being in place. Current strategies for identifying NETs would include staging at surgery, pathological grading, blood Chromogranin A (CgA) measurements, and detection of circulating tumor cells, with there being a clear need for additional and better biomarkers (34). The more accurate identification of neuroendocrine-associated cancers in particular would have important implications regarding the treatment of this disease (46).

No potential conflicts of interest were disclosed.

Conception and design: C.J. Creighton

Development of methodology: C.J. Creighton

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): F. Chen, Y. Zhang, B. Deneen, D.J. Kwiatkowski, C.J. Creighton

Writing, review, and/or revision of the manuscript: D.L. Gibbons, D.J. Kwiatkowski, M. Ittmann, C.J. Creighton

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C.J. Creighton

Study supervision: C.J. Creighton

This work was supported in part by Cancer Prevention and Research Institute of Texas (CPRIT) grants RP120713 C2 (to C.J. Creighton), RP150405 (to D.L. Gibbons), and RP120713 P2 (to D.L. Gibbons), and by the NIH grant P30CA125123 (to C. Creighton).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.